Course compendium, Applied matrix analysis Christopher Engstr¨ om Karl Lundeng˚ ard Sergei Silvestrov November 14, 2014 Contents 1 Review of basic linear algebra 1.1 Notation and elementary matrix operations 1.2 Linear equation systems . . . . . . . . . . . 1.2.1 Determinants . . . . . . . . . . . . . 1.3 Eigenvalues and eigenvectors . . . . . . . . 1.4 Exercises . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 . 6 . 8 . 10 . 13 . 14 2 Matrix factorization and canonical forms 2.1 Important types of matrices and some useful properties . 2.1.1 Triangular and Hessenberg matrices . . . . . . . . 2.1.2 Hermitian matrices . . . . . . . . . . . . . . . . . . 2.1.3 Unitary matrices . . . . . . . . . . . . . . . . . . . 2.1.4 Positive definite matrices . . . . . . . . . . . . . . 2.2 Factorization . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Spectral factorization . . . . . . . . . . . . . . . . 2.2.2 Rank factorization . . . . . . . . . . . . . . . . . . 2.2.3 LU and Cholesky factorization . . . . . . . . . . . 2.2.4 QR factorization . . . . . . . . . . . . . . . . . . . 2.2.5 Canonical forms . . . . . . . . . . . . . . . . . . . 2.2.6 Reduced row echelon form . . . . . . . . . . . . . . 2.2.7 Jordan normal form . . . . . . . . . . . . . . . . . 2.2.8 Singular value decomposition (SVD) . . . . . . . . 2.3 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Portfolio evaluation using Monte Carlo simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 17 17 18 18 19 19 20 20 21 22 22 23 23 24 25 25 3 Non-negative matrices and graphs 3.1 Introduction to graphs . . . . . . . . . 3.2 Connectivity and irreducible matrices 3.3 The Laplacian matrix . . . . . . . . . 3.4 Exercises . . . . . . . . . . . . . . . . 3.5 Applications . . . . . . . . . . . . . . . 3.5.1 Shortest path . . . . . . . . . . 3.5.2 Resistance distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 27 30 33 35 35 35 36 . . . . 39 39 42 43 43 . . . . . . . . . . . . . . . . . . . . . 4 Perron-Frobenius theory 4.1 Perron-Frobenius for non negative matrices 4.2 Exercises . . . . . . . . . . . . . . . . . . . 4.3 Applications . . . . . . . . . . . . . . . . . . 4.3.1 Leontief input-output model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Stochastic matrices and Perron-Frobenius 46 5.1 Stochastic matrices and Markov chains . . . . . . . . . . . . . . . . . . . . 46 1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 Irreducible, primitive stochastic matrices . . . . Irreducible, imprimitive stochastic matrices . . Reducible, stochastic matrices . . . . . . . . . . Hitting times and hitting probabilities . . . . . A short look at continuous time Markov chains Exercises . . . . . . . . . . . . . . . . . . . . . Applications . . . . . . . . . . . . . . . . . . . . 5.8.1 Return to the Leontief model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 50 51 53 55 57 58 58 6 Linear spaces and projections 6.1 Linear spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Inner product spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Projections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.1 Fourier approximation . . . . . . . . . . . . . . . . . . . . . . . 6.4.2 Finding a QR decomposition using the Gram-Schmidt process 6.4.3 Principal component analysis (PCA) and dimension reduction . . . . . . . . . . . . . . . 60 60 66 70 74 74 75 77 7 Linear transformations 7.1 Linear transformations in geometry . . . . . . . . . . . . . . . . . . . . . 7.2 surjection,injection and combined transformations . . . . . . . . . . . . 7.3 Application: Transformations in computer graphics and homogeneous coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.4 Kernel and image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.5 Isomorphisms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Least Square Method (LSM) 8.1 Finding the coefficients . . . . . . . . . . 8.2 Generalized LSM . . . . . . . . . . . . . 8.3 Weighted Least Squares Method (WLS) 8.4 Pseudoinverses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 . 82 . 84 . 85 . 88 . 91 . . . . 93 95 96 97 98 9 Matrix functions 100 9.1 Matrix power function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 9.1.1 Calculating An . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 9.2 Matrix root function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 9.2.1 Root of a diagonal matrix . . . . . . . . . . . . . . . . . . . . . . . 103 9.2.2 Finding the square root of a matrix using the Jordan canonical form104 9.3 Matrix polynomials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 9.3.1 The Cayley-Hamilton theorem . . . . . . . . . . . . . . . . . . . . 106 9.3.2 Calculating matrix polynomials . . . . . . . . . . . . . . . . . . . . 109 9.3.3 Companion matrices . . . . . . . . . . . . . . . . . . . . . . . . . . 110 9.4 Exponential matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 9.4.1 Solving matrix differential equations . . . . . . . . . . . . . . . . . 114 2 9.5 General matrix functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 10 Matrix equations 116 10.1 Tensor products . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 10.2 Solving matrix equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 11 Numerical methods for computing eigenvalues and 11.1 The Power method . . . . . . . . . . . . . . . . 11.1.1 Inverse Power method . . . . . . . . . . 11.2 The QR-method . . . . . . . . . . . . . . . . . 11.3 Applications . . . . . . . . . . . . . . . . . . . . 11.3.1 PageRank . . . . . . . . . . . . . . . . . Index eigenvectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 . 121 . 124 . 126 . 129 . 129 133 3 4 Notation v vector (usually represented by a row or column matrix) |v| length (Euclidean norm) of a vector kvk norm of a vector M matrix M > transpose M −1 inverse M ∗ pseudoinverse M complex conjugate M H Hermitian transpose M k the k:th power of matrix M, see section 9.1 mi,j The element on the i:th row and j:th column of M (k) mi,j The element on the i:th row and j:th column of Mk Mi. i:th row of M M.j j:th column of M ⊗ Kroenecker (tensor) product a scalar |a| absolute value of a scalar h·, ·i Z + inner product (scalar product) set of integers Z set of positive integers R set of real numbers C set of complex numbers Mm×n (K) set of matrices with dimension m × n and elements from the set K 5 1 Review of basic linear algebra The reader should already be familiar with the concepts and content in this chapter, instead it serves more as a reminder and overview of what we already know. We will start by giving some notation used throughout the book as well as reminding ourselves of the most basic matrix operations. We will then take a look at linear systems, eigenvalues and eigenvectors and related concepts as seen in elementary courses in linear algebra. While the methods presented here works, they are often very slow or unstable when used on a computer. Much of our later work will relate to how to find better ways to find these or similar quantities or give other examples of where and why they are useful. 1.1 Notation and elementary matrix operations We will here go through the most elementary matrix operations as well as give the notation which will be used throughout the book. We remind ourselves that a matrix is a rectangular array of numbers, symbols or expressions. Some example of matrices can be seen below: 1 0 A= 0 4 1 −2 C = −2 3 3 1 B= 0 1−i 1 0 0 D= 0 3 0 We note that matrices need not be square as in C and D, we also note that although matrices can be complex valued, in most of our examples or practical applications we need only to work with real valued matrices. • We denote every element of A as ai,j , where i is it’s row number and j is it’s column number. For example looking at C we have c3,2 = 1. • The size of a matrix is the number of rows and columns (or the index of the bottom right element): We say that C above is a 3 × 2 matrix. • A matrix with only one row or column is called a row vector or column vector respectively. We usually denote vectors with an arrow v and we assume all vectors are column vectors unless stated otherwise. • The diagonal of a matrix is the elements diagonally from the top left corner towards the bottom left (a1,1 , a2,2 , . . . an,n ). A matrix where all elements except those on the diagonal is zero is called a diagonal matrix . For example A is a diagonal matrix, but so is D although we will usually only consider square diagonal matrices. • The trace of a matrix is the sum of the elements on the diagonal. 6 Matrix addition Add every element in the first matrix with correspondign element in the second matrix. Both matrices need to be of the same size for addition to be defined. 1 0 1 −2 1+1 0−2 2 −2 + = = 0 2 2 −1 0+2 2−1 2 1 Although addition between a scalar and a matrix is undefined, sometimes authors write it as such in which case they usually means to add the scalar to every element of the matrix as if the scalar was a matrix of the same size with every element equal to the constant. The general definition of matrix addition is: Let A = B + C where all three matrices have the number of rows and columns. Then ai,j = bi,j + ci,j for every element in A. Matrix multiplication Given the two matrices A and B we get the product AB as the matrix whose elements ei,j are found by multiplying the elements in row i of A with the elements of column j of B and adding the results. 11 10 7 14 3 2 0 1 3 2 0 1 1 2 3 4 1 2 3 4 5 9 16 24 0 1 2 1 0 1 2 1 = 1 4 8 9 0 1 5 10 0 0 1 3 0 0 1 3 For the colored element we have: 1 · 3 + 2 · 1 + 3 · 0 + 4 · 0 = 5 We note that we need the number of columns in A to be the same as the number of rows in B but A can have any number of rows and B can have any number of columns. Also we have: • Generally AB 6= BA, why? • The size of AB is the number of rows of A times the number of columns of B. • The Identity matrix I is the matrix with ones on the diagonal and zero elsewere. For the Identity matrix we have: AI = IA = A. Multiplying a scalar with a matrix is done by multiplying the scalar with every element of the matrix. 1 0 3·1 3·0 3 0 3 = = 0 2 3·0 3·2 0 6 The general definition of matrix multiplication is: Let A = BC then j X ai,j = bi,k ck,j k=1 for each element in A. 7 Matrix transpose The transpose of a matrix is obtained by flipping the rows and columns such that the elements of the first row is now the elements of the first column. For example given matrix A below: 1 −2 A = −1 0 3 1 We get A transpose, denoted A> as: A> = 1 −1 3 −2 0 1 The general definition of the transpose is: Let B = A> then bi,j = aj,i for all elements in B. For complex valued matrices we usually instead talk about the Hermitian transpose where we first take the transpose and then take the complex conjugate (change sign of imaginary part). We usually denote this by AH (sometimes this is denoted by A∗ but here we will reserve that notation for pseudoinverses, see section 8.4). The general definition of the Hermitian transpose (also known as Hermitian conjugate or conjugate transpose) is: Let B = A> then bi,j = a∗j,i , where ∗ denotes the complex conjugate, for all elements in B. A real matrix where A = A> is called a symmetric matrix and a complex matrix where AH = A we call a hermitian matrix . 1.2 Linear equation systems Very often in applications we end up with a system of linear equations which we want to solve or find an as good solution as possible. Either repeatedly for the same or similar systems or only once. Matrices presents presents powerful tools which we can use to solve these problems, in fact even if we have a nonlinear system, we very often approximate it with a linear system and try to solve that instead. Linear systems will be treated repeatedly throughout the book, but here we take a short look at the problem itself, when we can solve it, and how from what we should be familiar with from courses in linear algebra. Given a linear equation system a1,1 x1 + a1,2 x2 + . . . a2,1 x1 + a2,2 x2 + . . . .. .. .. . . . am,1 x1 + am,2 x2 + . . . 8 + + a1,n xn a2,n xn .. . = = b1 , b2 , .. . + am,n xn = bm , we can represent it using a matrix and column vectors b1 x1 a1,1 a1,2 . . . a1,n a2,1 a2,2 . . . a2,n x2 b2 .. .. .. .. = .. . . . . . . . . bm xm am,1 am,2 . . . am,n Or in a more compact form as Ax = b. Solving this system of equations is the same as finding an inverse matrix A−1 such that A−1 A = I and multiply both expressions from the left with this matrix. The reader should be familiar with using Gaussian elimination to write the system in a easy to solve form, or going a step further and finding the inverse A−1 using Gauss-Jordan elimination. We remember for Gaussian elimination using elementary row operations (adding and multiplying rows) trying to find pivot elements in order to change our matrix to diagonal/row echelon form . F F F F F F F 0 0 0 F F F F F 0 0 0 0 F F F F 0 0 0 0 0 F 0 F F 0 0 0 0 0 0 F 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 After we have done this we can easily solve the system as long as it is solvable. Not all systems can be solved exactly, for example the system below is obviously unsolvable: 1 2 x1 1 = 1 2 x2 2 For square matrices we can go one step further by using the Gauss-Jordan elimination reducing the elements above the diagonal to zero as well in the same manner. Similarly some systems have more than one solution. For example will any point on the line y = b−1 (c − ax) solve the linear equation system 4 −1 x1 2 = −8 2 x2 −4 We are often interested in if there is a single unique solution to a linear equation system. Theorem 1.1. The following statements are equivalent for A ∈ Mn×n (K): 1. AX = B has a unique solution for all B ∈ Mn×m where X ∈ Mn×m . 2. AX = 0 only has the trivial solution (X = 0). 3. A has linearly independent row/column vectors. 9 4. A is invertible/non-degenerate/non-singular. 5. det(A) 6= 0. 6. A has maximum rank (rank(A) = n). 7. A’s row/column space has dimension n. 8. A’s null space has dimension 0. 9. A’s row/column vectors span K Most of these statements should be familiar already, some will be looked at more in depth later, especially when we look at linear spaces in chapter 6. 1.2.1 Determinants One of the ways to guarantee that a linear system Ax = b has a unique solution is that the determinant is non-zero det(A) 6= 0. We remember the determinant as the product of the eigenvalues. For a small matrix we can easily calculate the determinant by hand, we had the formulas for 2 × 2 and 3 × 3 matrices as: a b c d = ad − bc a b c d e f = aei − af h − bdi + bf g + cdh − ceg g h i And for larger matrices we could write it as a sum of determinants of smaller matrices: a b c d e f = a e f − b d f + c d e h i g i g h g h i And in general we get the determinant function det : Mn×n (K) 7→ K det(A) = det(A.1 , A.2 , . . . , A.k , . . . , A.n ) Some properties of the determinant: • det A> = det A • det AB = det A det B but det(A + B) 6= det(A) + det(B). • det cA = cn det A • det is multilinear 10 • det is alternating • det(I) = 1 With multilinear we mean that it’s linear in each argument det(A.1 , A.2 , . . . , a · A.k + b · B, . . . , A.n ) = = a det(A.1 , A.2 , . . . , A.k , . . . , A.n ) + b det(A.1 , A.2 , . . . B, . . . A.n ) For example if A is a 3 × 3 matrix and B, C are 3 × 1 matrices. Let M = [A.1 A.2 + B A.3 + C] then det(M) = det(A.1 , A.2 + B, A.3 + C) = = det(A.1 , A.2 , A.3 + C) + det(A.1 , B, A.3 + C) = = det(A) + det(A.1 , B, A.3 ) + det(A.1 , A.2 , C) + det(A.1 , B, C) Since the determinant is alternating, switching two columns (or rows) of A switches the sign of det(A) det(B) = det(A.1 , A.2 , . . . , A.k+1 , A.k , . . . , A.n ) = − det(A) But more important is what happens when we have two linear dependent vectors: det(A.1 , A.2 , . . . , B, B, . . . , A.n ) = 0 Then the multilinear and alternating property gives: linear dependence between columns (or rows) in A ⇔ det(A) = 0 Apart from using column or row expansion to calculate the determinant there is an alternative formula for the determinant X aσ1 ,1 aσ2 ,2 . . . aσn ,n sign(σ) (1) det(A) = σ∈Sn σ is a permutation of columns (rows), Sn is the set of all possible permutations, aσi ,i is the element that ai,i is switched with. sign is a function that is equal to 1 if an even number of rows have been switched and -1 if an odd number of rows have been switched. For example for a 2 × 2 matrix A a1,1 a1,2 A= a2,1 a2,2 We get two possible permutations: σ a do not switch places of columns, aσ1a ,1 = a1,1 , aσ2a ,2 = a2,2 . σ b switch places of columns, aσb ,1 = a1,2 , aσb ,2 = a2,1 . 1 2 And we get determinant det(A) = aσ1a ,1 aσ2a ,2 sign(σ a ) + aσb ,1 aσb ,2 sign(σ b ) = a1,1 a2,2 − 1 2 a2,1 a1,2 11 The basic idea for the formula is as follows. We start by rewriting columns as sums of columns of the identity matrix n X A.k = aσk ,k I.σk σk =1 Using the multilinearity of the determinant function we get det( n X aσ1 ,1 I.σ1 , A.2 , . . . A.n ) = σ1 =1 = n X aσ1 ,1 det(I.σ1 , A.2 , . . . A.n ) = σ1 =1 = n X aσ1 ,1 σ1 =1 = n X aσ1 ,1 σ1 =1 n X aσ2 ,2 det(I.σ1 , I.σ2 , . . . A.n ) = σ2 =1 n X n X aσ2 ,2 . . . σ2 =1 aσn ,n det(I.σ1 , I.σ2 , . . . I.σn ) σn =1 Terms where σi = σk for some i 6= k, 1 ≤ i, k ≤ n will be equal to zero. det(A) = X aσ1 ,1 aσ2 ,2 . . . aσn ,n det(I.σ1 , I.σ2 , . . . I.σn ) σ∈Sn The (sign) function captures the alternating behaviour. Q 1≤i<j≤n sign(σ) = Q (σj − σi ) 1≤i<j≤n (j − i) = ±1 While the function given here i correct, it’s extremely slow for large matrices, for an n × n matrix there is n! permutations. Instead we mostly we calculate it by column (row) expansion as we are used to for 3 × 3 matrices: column i row i z }| { z }| { n n X X det(A) = ak,i Aki = ai,k Aik k=1 (2) k=1 ˜ ki ) where A ˜ ki is equal to A with the k:th row and where Aki is equal to (−1)k+i det(A i:th column removed. Aki is called the ki-cofactor of A. 12 1.3 Eigenvalues and eigenvectors We are often interested in finding eigenvalues or eigenvectors of a matrix. Definition 1.1 If Av = λv where A ∈ Mn×n (K), λ ∈ K, v ∈ Mn×1 (K) then λ is a eigenvalue of A and v is a (right) eigenvector if v 6= 0. We call an eigenvalue λ with corresponding eigenvector v an eigenpair (λ, v). A matrix can have (and usually do) have multiple eigenvalues and eigenvectors. A single eigenvalue can also be associated with multiple eigenvectors. Often we are interested in the ”dominant” eigenvalue or the spectral radius of a matrix A. Definition 1.2 The set of eigenvalues for a matrix A is called the spectrum of A and is denoted by Sp(A). The absolute value of the largest eigenvalue is called the spectral radius ρ(A) = max |λ| λ∈Sp(A) The eigenvalues whith the largest absolute value (if there is only one) is often called the dominant eigenvalue. This and it’s corresponding eigenvector is often of particular interest as we will see for example while working with stochastic matrices in chapter 5. We can find the eigenvalues by solving det(A − λI) = 0 which is an nth degree polynomial in λ. pA (λ) = det(A − λI) is called the characteristic polynomial and pA (λ) = 0 is called the characteristic equation. Some properties of eigenvalues are: • The trace of A is the sum of the eigenvalues of A: tr(A) = n X i=1 • det(A) = n Y ai,i = n X λi i=1 λi i=1 • If A ∈ Mn×n (K) have n different eigenvalues then A is diagonalizable. • If A have no eigenvalue λ = 0 then the determinant is non-zero and the linear system Ax = b have a single unique solution. Often when interested in finding the eigenvalues or eigenvectors we first find a ”similar” matrix with the same unknown eigenvalues and then solve the problem for this new matrix instead. 13 Definition 1.3 Given two n × n matrices A, B we say that they are similar if there is an invertible n × n matrix P such that B = P−1 AP Similar matrices share not only eigenvalues but also many other properties because of that. Theorem 1.2. Let matrices A, B be similar, then A, B have the same • Eigenvalues (but generally not the same eigenvectors) • Determinant, det(A) = det(B). • Rank • Trace We will look much more at eigenvalues and eigenvectors for special types of matrices in later chapters, as well as how to compute them efficiently in chapter 11. 1.4 Exercises a) Show that this linear equation system has one single −x + 2y − z = −x + y = 2x − y + 2z = b) Show that this linear equation system has infinitely x + y − z = 2x + y + 2z = 3x + 2y − z = solution. 2 3 1 many solutions. 1 2 3 c) Show that this linear equation system has no solution. x − 2y + z = 3 x + 3y − 2z = 2 2x + y − z = 1 Calculate these determinants 2 3 1 5 4 , c = 2 −1 5 a = 5 , b = 3 2 5 0 −1 14 1.4.1 Let 3 0 A= 0 0 3 1 0 0 4 2 3 0 4 3 2 4 ,B = 2 3 3 2 3 1 1 1 4 2 3 3 4 3 2 4 ,C = 2 3 3 3 3 1 1 2 4 2 3 1 4 2 , det(C) = −26 3 3 The following determinants should all be very simple to calculate. Calculate them. a) det(A) b) det(B) c) det(C> ) 1.4.2 If the linear equation systems in 1.4 were written on the form Ax = y, for which systems would A−1 exist? 1.4.3 4 1 Find the eigenvalues and eigenvectors of M = 3 2 1.4.4 Given 1 0 0 1 1 2 1 2 1 A = −1 1 0 , B = 0 3 2 , P = 0 0 1 0 1 0 −1 0 1 0 2 3 Show that A and B are ( similar). 15 2 Matrix factorization and canonical forms Matrix factorization, also known as matrix decomposition, are methods in which we write one matrix as a product of several matrices. The aim is typically to change the problem into another easier problem through the choice of these factor matrices. Matrix factorization is very often used in order to speed up and/or stabilize calculations in numerical algorithms such as when finding eigenvalues of a matrix or solving a linear system. We will give some examples of applications but primary we will introduce different factorizations which we will then use in later chapters. One example we should already be familiar with is the diagonalization of a diagonalizable matrix A by elementary row (or column) operations. Since elementary row consist of multiplying rows with non-zero scalars and adding them together the application of a series of elementary row operations on a linear equation system, Ax = y, can be written as multiplication by a matrix, S, from the left. Elementary row operations:Ax = y → SAx = Sy The matrix S is always square and since it is easy to revert any elementary row operations it will also be invertible. Thus it is always allowed to write ˆ=y ˆ, x ˆ = Sx, y ˆ = Sy SASS−1 x = Sy ⇔ SAS−1 x Suppose that the matrix A was diagonalized in this manner. SAS−1 = D, D diagonal matrix This relation makes the following definition of a diagonalizable matrix quite natural Definition 2.1 A matrix A is diagonalizable if there exists an invertible matrix S such that A = S−1 DS Where D is a diagonal matrix. This if course an interesting property since a matrix equation Ax = y has a unique solution if A is diagonalizable. But there are many other important factorizations such as: • Spectral decomposition QΛQ−1 • LU-factorization • QR-factorization • Rank factorization CF • Jordan canonical form S−1 JS 16 • Singular value factorization UΣVH • Cholesky factorization GGH Before we look at how to find these factorizations and why they are useful we will start by looking at a couple of important types of matrices with some useful properties. 2.1 Important types of matrices and some useful properties We have already looked extensively at non-negative matrices in the earlier chapters, however there are many other types of matrices which are useful when trying to find a suitable matrix factorization. 2.1.1 Triangular and Hessenberg matrices A triangular matrix as we should already know is a matrix with all zeros above or below the diagonal. A matrix where all elements above the diagonal are zero we call lower triangular and a matrix where all elements below the diagonal are zero such as the one below, we call upper triangular. F F ... F 0 F . . . F A=. .. . . .. .. . . . 0 0 ... F We note that a triangular matrix does not need to be square although that is when they are the most useful. For a non-square matrix the diagonal is always the one starting in the upper left corner. To understand why triangular matrices are desirable we will consider to important properties. First of all we can very easily solve a linear system Ax = b if A is triangular since the last row is a equation of the form axn = b, Solving this we can substitute this answer into the row above that and get a new equation in just one variable. A second important property of triangular matrices us that the eigenvalues of a triangular matrix is the same as the elements on the diagonal. You can easily check this by looking at what happens with the characteristic polynomial det(A − Iλ) when A is diagonal. This gives rise to one important question: Is it possible to rewrite a matrix as a product of a triangular matrix and another matrix which doesn’t change the eigenvalues when multiplied with this matrix and if possible, can it be done fast? Related to triangular matrices are the Hessenberg matrices. These are matrices that are ’almost’ triangular, instead of demanding that all elements below the diagonal is zeros as for an upper triangular matrix, a matrix that is an upper Hessenberg matrix if 17 we also allow the elements immediately below the diagonal to be non zero. This gives an upper Hessenberg matrix A the form: F F F ··· F F F F F F · · · F F F 0 F F · · · F F F A = 0 0 F · · · F F F .. .. .. . . .. .. .. . . . . . . . 0 0 0 · · · F F F 0 0 0 ··· 0 F F While not as immediately useful as a triangular matrix, Hessenberg matrices are useful just from the fact that they are nearly triangular. The diagonal doesn’t hold the eigenvalues as in a triangular matrix, but they often give a good approximation or initial guess for iterative methods. Multiplication of two (both upper or lower) Hessenberg matrices also result in a new Hessenberg matrix. We will see examples of how and why Hessenberg matrices are useful later when we look at the QR-method in section 11.2. 2.1.2 Hermitian matrices We start by defining the Hermitian conjugate of a matrix: Definition 2.2 The Hermitian conjugate of a matrix A is denoted AH and is defined by (AH )i,j = (A)j,i . In other words we transpose it and change the sign of the imaginary parts. So if A is real valued, the hermitian conjugate is the same as the transpose. Definition 2.3 A matrix is said to be Hermitian (or self-adjoint) if AH = A Notice the similarities with a symmetric matrix A> = A, a Hermitian matrix can be seen as the complex extension of real symmetric matrices. One property of Hermitian (and real symmetric) are that they are normal and thus diagonalizable. 2.1.3 Unitary matrices Matrix factorizations using unitary matrices are often used over other factorization techniques because of their stability properties for numerical methods. Definition 2.4 A matrix, A, is said to be unitary if AH = A−1 . We will look more at unitary matrices when we look at linear spaces and projections in chapter 6.2. Unitary matrices have a large number of ”good” properties. Theorem 2.1. Let U be a unitary matrix, then 18 a) U is always invertible. b) U−1 is also unitary. c) | det(U)| = 1 d) (UV)H = (UV)−1 if V is also unitary. e) For any λ that is an eigenvalue of U, λ = eiω , 0 ≤ ω ≤ 2π. f ) Let v be a vector, then |Uv| = |v| (for any vector norm). H g) The rows/columns of U are orthonormal, that is Ui. Uj.H = 0, i 6= j, Uk. Uk. = 1. One of the important things we can find from this is that if we multiply a matrix A with a unitary matrix A, then UA will have the same eigenvalues as A. We will meet unitary matrices again in later chapters when we will need them and easier can illustrate why they are useful. One example of a unitary matrix is the C matrix below which rotates a vector by the angle θ around the x-axis 1 0 0 C = 0 cos(θ) − sin(θ) 0 sin(θ) cos(θ) 2.1.4 Positive definite matrices Positive definite matrices often come up in practical applications such as covariance matrices in statistics. Positive definite matrices have many properties which makes them easier to work with than general matrices. Definition 2.5 We consider a square symmetric real valued n × n matrix A, then: • A is positive definite if x> Ax is positive for all non-zero vectors x. • A is positive semidefinite if x> Ax is non-negative for all non-zero vectors x. We define negative definite and negative semidefinite in the same way except demanding that the resulting vector is negative or non-positive instead. Similarly we can define it in the same way by replacing the transpose with the hermitian for matrices with complex valued elements. 2.2 Factorization As was mentioned in the introduction factorization (also called decomposition) is taking a matrix and rewriting it as a product of several matrices. This can make it easier to solve certain problems, similarly to how factorization of numbers can made division easier. 19 A well known and previously visited example of factorization is rewriting a diagonalizable matrix as a product of an invertible matrix and a diagonal matrix A = S−1 DS, D is diagonal 2.2.1 Spectral factorization Spectral factorization is a special version of diagonal factorization. It is sometimes referred to as eigendecomposition. Theorem 2.2. Let A be an square (n × n) matrix with linearly independent rows. Then A = QΛQ−1 where Q is a square n × n matrix whose columns are the eigenvectors of A and Λ is a diagonal matrix with the corresponding eigenvalues. The spectral factorization can be thought of as a ’standardized’ version of diagonal factorization. Also relates to projection, see section ??. 2.2.2 Rank factorization Theorem 2.3. Let A be an m × n matrix with rank(A) = r (A has r independent rows/columns). Then there exist a factorization such that A = CF where C ∈ Mm×r and F ∈ Mr×n We take a short look at how to find such a factorization, We start by rewriting the matrix on reduced row echelon form 1 F 0 0 0 F 0 0 0 0 1 0 0 F 0 0 0 0 0 1 0 F 0 0 0 0 0 0 1 0 F 0 B= 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 We then create C by removing all columns in A that corresponds to a non-pivot column in B. In this example we get: C = A.2 A.4 A.5 A.6 A.8 We then create F by removing all zero rows in B. In this example we get: F = B1. B2. B3. B4. B5. This gives the Rank factorization of A. 20 2.2.3 LU and Cholesky factorization The LU factorization is a factorization of the form: A = LU Where L is lower triangular and U is upper triangular. If we can find such a factorization we can for example easily find the solution to a linear system Ax = L(Ux) = b by first solving Ly = b and then solve Ux = y. Both these systems are easy to solve since L and U are both triangular. However not every matrix A have a LU factorization, not even every square invertible matrix. Instead we often talk about LUP factorization instead: Theorem 2.4. If A is a square matrix, then there exist a factorization such that: PA = LU Where P is a permutation matrix, U is a upper triangular matrix and L is a lower triangular matrix. P−1 LU is called the LUP factorization of A. We can still, provided that we found such a factorization for example easily solve the equation system Ax = b by instead solving PAx = Pb. The LU factorization can be seen as a variation of the Gaussian elimination which we often use when attempting to solve a linear system by hand. LU decomposition can also be used to calculate the determinant of A = P−1 LU. We get the determinant as: det(A) = det(P −1 ) det(L) det(U ) = det(P −1 ) n Y i=1 ! li,i n Y ! ui,i i=1 This since we remember that the eigenvalues of a triangular matrix is the values on the diagonal. The determinant of a permutation matrix can be found to be det(P −1 ) = (−1)S where S is the number of row exchanges. If A is Hermitian and positive-definite however then we can find another related factorization, namely a Cholesky factorization: Theorem 2.5. If A is a Hermitian, positive-definite matrix, then there exist a factorization such that: A = LLH Where L is a lower triangular matrix. LLH is called the Cholesky factorization of A. The Cholesky factorization is mainly used instead of the LUP factorization when possible to solve linear systems. Some example where such matrices naturally occur is. 21 • Any covariance matrix is positive-definite and any covariance matrix based on measured data is going to be symmetric and real-valued. This means it’s hermitian and we can use the Cholesky factorization. • When simulating multiple correlated variables using Monte Carlo simulation the correlation matrix can be decomposed using the Cholesky factorization. L can then be used to simulate the dependent samples. 2.2.4 QR factorization The QR-factorization is one of the most useful factorizations. While we will only show what it is and give one example of what it can be used for, section 11 is dedicated to how we can find the QR-factorization efficiently as well as it’s use in the QR-method, one of the most widely used method to find the eigenvalues of a general square matrix. Theorem 2.6. Every n × m matrix A have a matrix decomposition A = QR where • R is a n × m upper triangular matrix. • Q is a n × n unitary matrix. If A is a square real matrix then Q is real and therefore an orthogonal matrix Q> = Q. We can find a QR-factorization can be computed using for example Householder transformations or the Gram-Schmidt process. We will look more into this in ??. The QR-factorization have uses in for example: • Given a QR-factorization we can solve a linear system Ax = b by solving Rx = Q−1 b = QH b. Which is can be done fast since R is a triangular matrix. • QR-factorization can also used in solving the linear least square problem. • It plays an important role in the QR-method used to calculate eigenvalues of a matrix. 2.2.5 Canonical forms Some factorizations are referred to as canonical forms. A canonical form is a standard way of describing an object. A single object can have several different canonical forms and which one we use might depend on what we want to find or which canonical form we can easily calculate. We have already seen one example of a canonical form which we used when looking at reducible stochastic matrices. Some other examples of canonical forms for matrices are: • Diagonal form (for diagonalizable matrices) 22 • Reduced row echelon form (for all matrices) • Jordan canonical form (for square matrices) • Singular value factorization form (for all matrices) 2.2.6 Reduced row echelon form Definition 2.6 A matrix is written on reduced row echelon form when they are written on echelon form and their pivot elements are all equal to one and all other elements in a pivot column are zero. B= 0 0 0 0 0 0 0 0 1 F 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 F F F F 0 0 0 0 0 0 0 0 1 0 0 0 Theorem 2.7. All matrices are similar to some reduced row echelon matrix. We recognize this as the result of applying Gauss-Jordan elimination to the original matrix. 2.2.7 Jordan normal form The Jordan normal form, while not especially useful in practical applications because of its sensitivity to rounding errors in numerical calculations is still good to be familiar with since it is often usefor for developing mathematical theory. We start by defining the building blocks which which we build the Jordan normal form. Definition 2.7 (Jordan block) A Jordan block is a square matrix of the form λ 0 Jm (λ) = ... 0 0 1 0 ... λ 1 ... .. . . . . . . . 0 ... λ 0 ... 0 0 0 .. . 1 λ We call a matrix made up of Jordan blocks a Jordan matrix: 23 Definition 2.8 (Jordan matrix) A Jordan matrix is a square matrix of the form 0 Jm1 (λ1 ) 0 J m 2 (λ2 ) J= .. .. . . 0 0 ... ... .. . 0 0 .. . . . . Jm1 (λk ) Theorem 2.8. All square matrices are similar to a Jordan matrix. The Jordan matrix is unique except for the order of the Jordan blocks. This Jordan matrix is called the Jordan normal form of the matrix. We remember that two matrices A, B are similar if there is a unitary matrix S such that A = S−1 BS. This means that if we have found a Jordan matrix J similar to A, they have the same eigenvalues. And since Jordan matrices are upper triangular, we can easily find the eigenvalues along the diagonal of J. Theorem 2.9 (Some other interesting properties of the Jordan normal form). Let A = S−1 JS a) The eigenvalues of J is the same as the diagonal elements of J. b) J has one eigenvector per Jordan block. c) The rank of J is equal to the number of Jordan blocks. d) The normal form is sensitive to perturbations. This means that a small change in the normal form can mean a large change in the A matrix and vice versa. 2.2.8 Singular value decomposition (SVD) Singular value decomposition is a canonical form used for example in statistics and information processing in dimension reducing methods such as Principal component analysis or calculating pseudoinverses (see section 8.4). Theorem 2.10. All A ∈ Mm×n can be factorized as A = UΣVH where U and V are unitary matrices and Sr 0 Σ= 0 0 where Sr is a diagonal matrix with r = rank(A). The diagonal elements of Sr are called the singular values. The singular values are uniquely determined by the matrix A (but not necessarily their order). 24 The singular value decomposition is related to the spectral decomposition. Given a singular value decomposition A = UΣVH we have: • The columns of U are the eigenvectors of AAH • The columns of V are the eigenvectors of AH A • The singular values on the diagonal of Σ are the square root of the non-zero eigenvalues of both AAH and AH A. Having the singular value decomposition it’s easy to solve a lot of problems found in many practical applications, for example if we have a singular values decomposition A = UΣVH we can: • Calculate pseudoinverse A+ = VΣ+ UH where Σ+ is easily found by taking the reciprocal (x−1 ) of the non-zero elements of the diagonal and transposing the resulting matrix. • It is also used in Principal component analysis (PCA) as we will see later. 2.3 Applications 2.3.1 Portfolio evaluation using Monte Carlo simulation This example uses a large amount of staistics and terminology not otherwise described in this book. For the interested reader we refer to books or courses in Monte Carlo methods and stochastic processes for the details. Suppose we have two stocks A, B with prices at time t as S a (t), S b (t). Initially we start with na stocks of A and nb stocks of B. This gives initial wealth: W0 = na S a (0) + nb S b (0) We let the stocks be for a time T after which we are interested in the new wealth of the portfolio or more specifically how sure we can be that there have not occurred any great losses. We assume the stock pricing can be modelled as geometric brownian motions (GBM) with parameters S a GBM(µa , σa ) and S b GBM(µb , σb ). We have estimated correlation ρ between the return on A and B. To evaluate our portfolio we are interested in the probability that the total wealth of our portfolio WT at the end of the period T have dropped by more than 10%. We can write this as: P WT ≤ 0.9 W0 25 Finding this probability analytically is rather hard, especially if we work with not 2 stock but say 10 dependent stocks. Instead we use Monte Carlo simulation to simulate the process a large amount of times. Every time checking at the end if WT ≤ 0.9W0 . The value of stock A at time T can be written: STa = S0a exp((µa − σa2 /2)T + √ T Va ) Where Va , Vb are multivariate normal distributed variables Va , Vb N (0, Σ). 2 σa σa σb ρ Σ= σa σb ρ σb2 So in order to generate one sample of STa , STa we need to simulate one sample from N (0, Σ). We can usually rather easily simulate independent normal distributed variables N (0, I). Our goal is then to use these independent samples to generate the dependent samples N (0, Σ). Normal distribution have the property that a linear combination of normal distributed variables is itself normal. In fact if we have normal distributed variables with mean zero and variance one Z N (0, 1) then C> Z is distributed as follows. C> Z N (0, C> C) In other words, if we can find a Cholesky decomposition C> C = Σ then we need only to calculate the product C> Z to get samples from the desired distribution. Σ is obviously real valued and symmetric, but it’s also at least positive semi definite. In fact the only time it isn’t strictly positive definite is when one one of the variables is set to be an exact linear combination of the others. Since we can be almost sure Σ is a real, symmetric, positive definite matrix we can use Cholesky factorization to find a decomposition C> C = Σ. We thus generate a large amount of samples C> C and check how often WT ≤ 0.9W0 to find the probability P (WT ≤ 0.9W0 ). Some things to think about: • Can you think of another situation where you would like to take samples from another multivariate normal distribution? • Can you think of why the covariance matrix Σ is always a real, symmetric, positive (semi)definite matrix? 26 3 Non-negative matrices and graphs Non-negative matrices are common in many applications such as when working with stochastic processes and Markov chains as well as when working with things we can represent with a graph, such as electrical networks, water flows, random walks, and so on. Non-negative and positive matrices (not to be confused with positive (semi) definite matrices) have many known properties, especially concerning their eigenvalues and eigenvectors which makes them easier to work with than general matrices. This and the next to chapters are very much related in that we will start by looking at some graph theory, how to represent a graph with a matrix and some things which we can calculate here. After which we will take a look at the properties of non-negative matrices themselves in chapter 4 and then last look at a special type of non-negative matrix related to random walks on graphs in chapter 5. 3.1 Introduction to graphs A graph is a representation of a set of objects (vertices) where some pair of objects have links (edges) connecting them. The vertices can represent a state in a Markov chain, a point on a grid, a crossing in a road network among other things. Likewise the edges could in the same examples represent the possible transitions in the Markov chain, the lines on the grid or the actual roads in a road network. Sometimes in litterature vertices are called nodes and edges are called links. An example of two graphs can be seen in Figure 1. B A B C D A C D Figure 1: A directed graph (left) and a undirected graph (right) We will concern ourselves with mainly two types of graphs, directed graphs as in the one to the left where every edge connecting two vertices have a direction. And undirected graphs as in the one to the right where the edges doesn’t have a direction but instead only denote that the two vertices are connected. Sometimes we will allow edges which goes from a vertice back to itself (so called loops) and sometimes we won’t, generally loops will be allows for directed graphs but often not when working with undirected 27 graphs. In both cases we can also assign weights to the edges as we will see later. We can represent graphs such as the ones in Figure 1 with a matrix by defining the adjacency matrix for a graph. Definition 3.1 The adjacency matrix A of a graph G with n vertices (|V | = n) is a square n × n matrix with elements ai,j such that: 1, if there is an edge from vertice i to vertice j ai,j = 0, otherwise If the graph is undirected we consider every edge between two vertices as linking in both directions. We note one important thing, we don’t define in which order we take the vertices. Sometimes it’s natural to order the vertices in a certain way, for example depending on vertice labels, however any order is valid as long as you know what order you choose. B A C 0 0 A= 0 1 0 1 1 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 0 0 D E Figure 2: A directed graph and corresponding adjacency matrix with vertices ordered {A, B, C, D, E} If the graph instead was undirected we would get the graph and adjacency matrix in ??. We note that we don’t need to order the vertices in that order, for example the matrix: 0 1 0 0 1 1 1 1 1 0 A= 0 1 0 1 0 0 1 1 0 0 1 0 0 0 0 Is also a adjacency matrix of the graph in Figure 3 where we instead have ordered the vertices as {C, B, A, D, E}. 28 B A C 0 1 A= 0 1 0 1 1 1 1 0 0 1 0 0 1 1 1 0 0 0 0 0 1 0 0 D E Figure 3: A undirected graph and corresponding adjacency matrix with vertices ordered {A, B, C, D, E} One common addition to graphs is that we also assign weights to the edges, usually we only allow the weights to be positive real values. The weights can for example be used to represent transition probabilities in a Markov chain, distance of a part of road in a road network, time needed to travel between two vertices or the capacity in a flow network. If we have a graph with weighted edges we often talk about the distance matrix rather than the adjacency matrix. Definition 3.2 The distance matrix D of a graph G with n vertices and weighted edges E is a square n × n matrix with elements di,j such that: Ei,j , if there is an edge from vertice i to vertice j di,j = 0, otherwise Where Ei,j is the weight of the edge from vertice i to vertice j. As with the adjacency matrix, if the graph is undirected we consider every edge between two vertices as linking in both directions. An example of a weighted directed graph and corresponding distance matrix can be seen in 3.2. As in the case with the adjacency matrix we don’t have a defined order of the vertices, so we can use any as long as we know what order we use. While there would be no problem using negative or non real weights on the edges in the matrix representation of the graph, in practice we still usually consider only the cases with only positive weights. This is because we have a lot of usefull theory concerning specifically nonnegative matrices (which we get if we allow only positive weights). In most real world examples where we are interested in representing the problem with a graph we also naturally end up with only positive weights, such as when the weights represent for example distance, time, probability, capacity or a flow. 29 2 B 1 1 2 A C 1 2 1 0 0 D= 0 1 0 1 2 1 0 0 0 0 0 0 0 1 2 0 0 0 0 0 2 0 0 D E Figure 4: A weighted directed graph and corresponding distance matrix with vertices ordered {A, B, C, D, E} 3.2 Connectivity and irreducible matrices It’s often interesting to see not only if there is an edge between two vertices but also if there is a connection between two vertices if we are allowed to traverse more than one edge. This is interesting both in the case of directed graphs as well as undirected ones: can we get from A to B by following the edges? And if so in how many steps? Or what is the shortest distance between A and B? We take a look at so called paths in a graph. Definition 3.3 Given a weighted or unweighted graph. • A path in a graph is a sequence of vertices e1 , e2, . . . , en such that for every vertice in the sequence there is an edge to the next vertice in the sequence. • If the first vertice in a path is the same as the last vertice we call it a cycle. • The length of a path is the number of edges in the path (counting multiple uses of the same edge). We take a look at the adjacency matrix although the same is true for the distance matrix as long as we only allow positive weights. (k) Theorem 3.1. If element ai,j of the adjacency matrix A to the power of k is non zero, there is a path from vertice i to vertice j of exactly length k. In other words: (k) ai,j ∈ Ak ⇒ There is a path of length k from vertice i to vertice j 30 The proof is left as an exercise. The problem of finding a path or the existence of a path between two vertices is an important problem in many applications regarding graphs. For example in finding the shortest path in for example a road network as well as finding if and how likely it’s to get from one state to another in a Markov chain. We divide the graph’s in two different groups depending on if we can find a path from all vertices to all other vertices. Definition 3.4 Given a weighted or unweighted graph. • Two vertices is said to be connected if there exists a path from either vertice to the other in the undirected graph. • A graph is said to be connected if it is connected for all pair of vertices in the undirected graph. • If it is also connected for all pairs in the directed graph we say that it’s strongly connected. We see that if a graph is connected or strongly connected in the case of directed graphs, we have a guarantee that we have a path from any vertex to another. Before continuing to we take a short look at connected parts of a graph. Definition 3.5 Given a (possible weighted) graph: • A connected component in a undirected graph is a maximal part of the graph where all vertices are connected with eachother. • A strongly connected component is a part of the directed graph which is strongly connected. Example of connected components in an undirected graph can be seen in Figure 5 below. We see that the the connected components can themselves be considered connected graphs. A example of strongly connected components in a directed graph can be seen in Figure 6. How many edges would you need to make the graph strongly connected? We are now ready to move over to the matrix side of things. We look at general square non-negative matrices, for which both the adjacency matrix as well as the distance matrix belong if we only allow positive weights . Theorem 3.2. A non-negative square n × n matrix A is said to be Irreducible if any of the following equivalent properties is true: 1. The graph corresponding to A is strongly connected. 2. For every pair i, j there exists a natural number m such that am i,j > 0. 31 Figure 5: A undirected graph with 3 connected components Figure 6: A directed graph with 3 strongly connected components 3. (I + A)n−1 is a positive matrix. 4. A cannot be conjugated into a block upper triangular matrix by a permutation matrix P. We leave the last property for later and give a short motivation of why the first 3 are equivalent. It’s clear that since the graph is strongly connected (1) there exist a path from any vertice to any other vertice. However if the length of the path between vertice i and vertice j is m we know that am i,j > 0. So if (1) is fulfilled so is (2) and if two is fulfilled we know there is a path between all vertices so the graph needs to be strongly connected as well (1). For the third we can take the binomial expansion and get: (I + A) n−1 n X k k = A n k=0 k Since is positive the matrix is positive if there is at least one 0 ≤ k ≤ n for every n 32 pair i, j such that aki,j > 0. However since the longest possible path without repeating any vertices is n − 1 (go through all vertices in the graph) we see that this is true only if (2) is true since we need not check any m larger than n − 1. The question regarding if a non-negative matrix is irreducible or not is a very important property which we will use repeatedly in the comming chapters both with regards to graphs as well as generall non-negative matrices. Having a clear understanding of what is and isn’t an irreducible matrix and how they relate to graphs and their distance or adjacency matrix will give a much easier understanding when we talk about Perron-Frobenius theory and Markov chains later. 3.3 The Laplacian matrix The Laplacian matrix is another matrix representation of graphs with no self loops (so called simple graphs). It can be used to among others to find the resistance between nodes in an electrical network, estimate the robustness of or find if there are any ”bottlenecks” in a network. It’s also usefull in finding some other properties of the graph as well as it appearing naturally in for example the stiffness matrix of simple spring networks. Before defining the Laplacian matrix we need to define another matrix, namely the degree matrix of a graph. Definition 3.6 The degree matrix of a undirected graph is a diagonal matrix where the elements on the diagonal is the number of edges connected with that vertex. In other words given a graph G with n vertices the degree matrix D is a n × n matrix with elements: deg(i), i = j di,j := 0, i 6= j Where deg(i) is the number of edges connected with vertice vi . • If there are loops present (vertices linking to itself) we count those twice for the degree. • If the graph is directed we use either the in-degree(number of edges pointing towards the vertice) or out-degree (number of edges pointing out from the vertice). Which one we use depend on application. • With a directed graph we also only count loops once rather than twice. You can see a undirected graph and corresponding degree matrix in Figure 7 and Figure 8 Using this and the adjacency matrix we can define the laplacian matrix: 33 2 1 4 3 0 D= 0 0 0 4 0 0 0 0 1 0 0 0 0 2 3 Figure 7: Undirected graph and corresponding degree matrix. 2 1 4 2 0 D= 0 0 0 3 0 0 0 0 1 0 0 0 0 0 3 Figure 8: Directed graph and corresponding degree matrix using in-degree. Definition 3.7 The Laplacian matrix L of a undirected graph with no self loops is defined as: L := D − A where D is the degree matrix and A is the adjacency matrix. This gives elements li,j ∈ L: i=j deg(i), −1, i 6= j, vi links to vj li,j := 0, otherwise Where deg(i) is the number of edges connected with vertice vi . An example of a graph and corresponding Laplacian matrix can be seen in Figure 9. Some properties of the Laplacian matrix L and it’s eigenvalues λ1 ≤ λ2 ≤ . . . ≤ λn . • L is positive semidefinite (x> Lx ≥ 0). • λ1 = 0 and the number of connected components is it’s algebraic multiplicity. 34 2 3 −1 −1 −1 −1 2 0 −1 L= −1 0 1 0 −1 −1 0 2 1 4 3 Figure 9: Simple undirected graph and corresponding Laplacian matrix • The smallest non zero eigenvalue is called the spectral gap. • The second smallest eigenvalue of L is the algebraic connectivity of the graph. We immediately see that it seems like the eigenvalues of the Laplacian gives us a lot of information about the graph, this is especially important since it’s positive semidefinite making more effective calculations of the eigenvalues possible, which might be important if we need to work with very large graphs. We will look more into how we can use the Laplacian matrix in the application part. 3.4 Exercises 3.1. 3.5 Applications 3.5.1 Shortest path A common problem when working with networks such as a road network or in routing on the internet is the shortest path problem. You could for example want to find the best order in which to visit a number of places (the traveling salesman problem) or maybe you want to find the shortest path between two places such as in the case of GPS in cars. There are a number of related problems and solutions depending on if we need the shortest distance between all pair of vertices or if it’s enough with a subset such as a group of vertices or a single pair of vertices. One common algorithm for finding the shortests path from one vertex to all others is Djikstras algorithm which most other similar methods are based upon. The method works by assigning a distance from the initial vertex to every other vertex and then traveling through the graph always choosing the next vertex among the unvisited vertices that minimizes the distance. At every we we update the distance to every vertex linked 35 to by the current one. More exact the method works through the following steps using the distance matrix of the resulting graph: 1. Assing distance zero to initial vertex and infinity to all other verticess. 2. Mark all verticess as unvisited and set the initial vertex as current vertex. 3. For the current vertex, calculate the distance to all it’s unvisited neighbors as the current distance plus the distance to the neighbors. Update these vertices distance if the newly found distance is lower. 4. Mark current vertex as visited. 5. If the smallest distance in among the unvisited vertices is infinity: stop! Otherwise set the unvisited vertex with the smallest distance as current vertex and go to step 3. If we are only interested in the shortest distance from the initial vertex to one specific vertex we can stop as soon as the distance to the target vertex is less than the distance to the next current vertex. However for the method to work we need to make a few assumptions, there are also a couple of properties of the graph itself that can say something about the expected result of the algorithm. Take a few minutes to think through the following questions and how they affect the result or what assumptions we need. • What could happen if there are negative weight’s in the graph? • What if the graph is directed? Should a transportation network be considered directed or undirected? • If we know the graph is or isn’t irreducible, can we say something about the result beforehand? 3.5.2 Resistance distance When planning a large electrical network its important to try and get as low resistance between vertices as possible. It is however not an easy task to use the common rules for combining resistors when the network gets large and we might want to know how to expand the network in order to lower the resistance between not just two vertices but lower it overall. The resistance distance between two vertices vi , vj in a graph G is defined as the effective resistance between them when each edge is replaced with a 1 Ω resistor. We define 1 Γ=L+ n Where L is the laplacian matrix, n is the number of vertices in G and 1 is a n × n matrix with elements equal to one. 36 The elements of the resistance distance matrix (Ω) is then: −1 −1 (Ω)i,j = Γ−1 i,i + Γj,j − 2Γi,j . For example if we have the graph below and corresponding Laplacian matrix in Figure 10 below: B 3 −1 0 −1 −1 −1 3 −1 0 −1 L= 0 −1 3 −1 −1 −1 0 −1 3 −1 −1 −1 −1 −1 4 A E C D Figure 10: Undirected graph and corresponding Laplacian matrix We get Γ = L + 1/5 and it’s inverse: Γ−1 0.4267 0.1600 = 0.0933 0.1600 0.1600 0.1600 0.4267 0.1600 0.0933 0.1600 0.0933 0.1600 0.4267 0.1600 0.1600 0.1600 0.0933 0.1600 0.4267 0.1600 0.1600 0.1600 0.1600 0.1600 0.3600 −1 −1 We get the resistance distance matrix by computing (Ω)i,j = Γ−1 i,i + Γj,j − 2Γi,j . for every element: 0 0.5333 (Ω) = 0.6667 0.5333 0.4667 0.5333 0 0.5333 0.6667 0.4667 0.6667 0.5333 0 0.5333 0.4667 0.5333 0.6667 0.5333 0 0.4667 0.4667 0.4667 0.4667 0.4667 0 And we have found the effective resistance between vertices as the value of the elements in the matrix. We can see that we have the largest resistance between the vertices A, C and B, D, which is quite natural since they are the ones furthest from eachother. Between opposite vertices we can also easily check the result since we can see it as 2 ressistors made up of 3 parallel connected resistors each in a series. This would give resistance: RA,D = R4 R5 R6 1 1 2 R1 R2 R3 + = + = R1 R2 + R1 R3 + R2 R3 R4 R5 + R4 R6 + R5 R6 1+1+1 1+1+1 3 37 Since all the resistances is 1 from the definition of resistance distance. Checking the result for the resistance between for example A and E is harder using the usual computation rules though. 38 4 Perron-Frobenius theory 4.1 Perron-Frobenius for non negative matrices Perron-Frobenius theory gives information about many of the most essential properties of non-negative square matrices. We will start by formulating item 4.1 which is essential to this chapter. While it is not immediately easy to see why these properties are useful, we will see throughout the chapter why they are important and how and when we can use them. Be sure to review what we mean by irreducible matrices if your not sure since it’s one of the recurring themes of this chapter. Theorem 4.1. Let A be a square non-negative irreducible matrix with spectral radius ρ(A) = r. Then we have: 1. r, called the Perron-Frobenius eigenvalue, is an eigenvalue of A. r is a positive real number. 2. The (right)eigenvector x belonging to r, Ax = rx, called the (right)Perron vector, is positive ( x > 0). 3. The left eigenvector y belonging to r called the left Perron vector, is positive (y> A = ry> , y > 0). 4. No other eigenvector except those belonging to r have only positive elements. 5. The Perron-Frobenius eigenvalue r has algebraic and geometric multiplicity 1. X X 6. The Perron-Frobenius eigenvalue r satisfies: min ai,j ≤ r ≤ max ai,j . i j i j We remember the spectral radius being the largest (in absolute value) eigenvalue. The theorem says a lot about the ”largest” eigenvalue as well as the corresponding eigenvector. This is also often called the dominant eigenvalue and its value and eigenvector are in many applications one of the most important things to know about the system as can be seen in the application part. We give a short motivation for the correctness of the theorem. Given (2) we can easily prove that the eigenvalue r is positive: If x is a positive vector and A is non-negative, then Ax > 0 since A is irreducible. Since Ax = rx if r is an eigenvalue belonging to the eigenvector x we have r > 0. If r was negative, zero or complex rx would be as well, but Ax is a positive vector. We also prove (4) that the only positive eigenvectors are the ones belonging to r: let A be a non-negative irreducible square matrix, set ρ(A) = r and suppose y > 0 is the left Perron vector y> A = ry> . Assume we have an (right) eigenpair (λ, v), v ≥ 0. This gives: ry> v = y> Av = y> λv ⇒ r = λ 39 We give no further proof, but leave some of the other statements as exercises. Initially the theorem was stated for positive matrices only and then later generalised for irreducible non-negative matrices. It’s easy to see that if A is positive it’s also non-negative and irreducible. Theorem 4.2. If A is positive, in addition to all properties for irreducible non-negative matrices, the following is true as well: • The Perron-Frobenius eigenvalue r has the strictly largest absolute value of all eigenvalues. I.e., for any eigenvalue λ 6= r we have |λ| < r. This last property ,that there is only one eigenvalue on the spectral radius, is a very useful property as well. In the section on Markov chains, we will see that it determines convergence properties. However this was for positive matrices and we see that it isn’t true for general non-negative irreducible matrices as in the example below. We look at the permutation matrix: 0 1 P = 1 0 It’s obviously non-negative and irreducible. From Perron-Frobenius we get that the spectral radius ρ(A) = r = 1 since the sum of both rows are one. However the eigenvalues can easily be found to be ±1 which are both on the unit circle. While there are irreducible matrices for which there are more than one eigenvalue on the spectral radius, there are some irreducible matrices for which there is only one. We let this last property of having or not having only one eigenvalue on the spectral radius divide the non-negative irreducible matrices in two different classes. We say that a non-negative irreducible matrix is primitive if it has only one eigenvalue with absolute value equal to the spectral radius. If a non-negative irreducible matrix is not primitive, it is imprimitive. Theorem 4.3. We let A be a square non-negative irreducible matrix, then: • A is primitive if and only if there is only one eigenvalue r = ρ(A) on the spectral radius. • A is primitive if and only if lim (A/r)k exists. k→∞ xy> > 0. k→∞ y> x where x, y are the right and left Perron vectors respectively. For the limit we have: lim (A/r)k = • A is imprimitive if and only if there are h > 1 eigenvalues on the spectral radius. h is called it’s index of imprimitivity. 40 We note the easily missed part that for a matrix to be either primitive or imprimitive it needs to first be non-negative and irreducible. We also see that that only one eigenvalue on the spectral radius is equivalent to the limit lim (A/r)k existing. k→∞ We return to the permutation matrix: 0 1 P = 1 0 We consider the limit lim (A/r)k , given r = 1 which we got earlier. That means that k→∞ if lim (A)k exists then P is primitive. However P2 = I and the sequence A, A2 , A3 , . . . k→∞ alternate between P and I. Since it alternates between two matrices it obviously doesn’t converge and is therefore imprimitive. In fact all irreducible permutation matrices are imprimitive since for all of them there is a k such that Pk = I (why?) and r = 1 since every row sum to one. There are two common ways to tell if a non-negative irreducible matrix is primitive: Theorem 4.4. A non-negative irreducible matrix A is primitive if: 1. If any diagonal element ai,i > 0 A is primitive (but A could still be primitive even if there are all zeroes on the diagonal). 2. A is primitive if and only if Am > 0 for some m > 0. 3. A is primitive if and only if A is aperiodic. Proof. • If (1) is true it’s obvious that there is a m > 0 such that Am > 0. This since if we look at the graph: irreducibility guarantees there is a path from all vertices to all others and Am > 0 means there is a path of exactly length m from all vertices to all others. However since we have one diagonal element ai,i > 0 we can for all paths between two vertices going through ai,i increase the length of the path by one simply by looping within ai,i once. • If Am > 0 and A had a non positive eigenvalue λ = ρ(A). Then λm should be an eigenvalue on the spectral radius of Am . However since Am is positive the only eigenvalue on ρ(Am ) is positive and simple, hence A is primitive. Conversely, if A is primitive then the limit } lim Am /rm exists and is positive. But for this to be m→∞ possible we must that Am is positive for large enough m. Especially the first statement gives a way to often immediately conclude that a primitive matrix is in fact primitive. When the first can not be used you are left with trying the second, however you might often need to try quite large m before you can find if a 41 matrix is primitive or not. To show that a matrix is imprimitive it’s often easiest to show that it’s periodic as we did with the permutation matrix. The second statement can also be used however since it can be shown that if A is primitive, the smallest possible m which gives a positive matrix can be bounded by: m ≤ n2 − 2n + 2. So if you try a m at least as large as this and the matrix still isn’t positive, it’s imprimitive. We give a short note on imprimitive matrices as well: Theorem 4.5. If A is a imprimitive matrix with h eigenvalues λ1 , λ2 , . . . , λh on the spectral radius ρ(A). Then the following are true: • alg mult(λn ) = 1, n = 1, 2, . . . , h. • λ1 , λ2 , . . . , λh are the h−th roots of ρ(A), such that: {λ1 , λ2 , . . . , λh } = {r, rω, rω 2 , . . . , rω h−1 }, ω = ei2π/h We will not look at the imprimitive matrices more than this, except noting that the eigenvalues on the spectral radius for imprimitive matrices are actually very special. Not only do they all have algebraic multiplicity 1, they are evenly spaced around the spectral radius. 4.2 Exercises 4.1. Consider the matrices: 1 1/2 0 1/2 0 0 0 1 0 1 0 1/2 1/2 0 0 0 0 1 A= 0 1 0 0 B = 1/3 1/3 0 1/3 C = 2 3 0 0 0 1 1 0 0 0 1 3 1 2 2 1 1 3 3 2 1 1 For each matrix consider: • Is it non-negative? • Is it irreducible? • Is it primitive? • What is the maximum and minimum possible value of the Perron-Frobenius eigenvalue? • What can you say about the other eigenvalues? 4.2. Show that for all irreducible permutation matrices there is a k such that Pk = I. Prove that all irreducible permutation matrices are imprimitive. 42 4.3. Consider the matrix: 1 1 2 A = 1 1 1 2 2 1 • Calculate the eigenvalues of this matrix. • Calculate the (right) eigenvector to the Perron-Frobenius eigenvalue. 4.3 Applications 4.3.1 Leontief input-output model Suppose we have a closed economic system made up of n large industries, each producing one commodity. Let a J − unit be what industry j produces that sells for one unit of money. We note that this could be a part of one commodity, for example a car would typically sell for many K − units. We create a vector S with elements sj containing the produced commodities and a matrix A with consumption of other commodities needed to create the commodities. We have: • Elements sj ∈ S, 0 ≤ sj = number of J − units produced by industry j in one year. • Elements ai,j ∈ A,0 ≤ ai,j = number of I − units needed to produce one J − unit. • ai,j sj = number of I − units consumed by industry J in one year. • If we sum over all j we can then find the total I − units that are availible to to public (not consumed by the industries themselves), we call this di . • di = si − n X ai,j sj . j=1 • S − AS is then a vector with the total availible to the public of all commodities (elements di ). Suppose the public demands a certain amount of every commodity, we create a demand vector d = (d1 , d2 , . . . , dn ) ≥ 0. We want to find the minimal supply vector S = [s1 , s2 , . . . , sn ]> ≥ 0 that satisfy the demand. We are also interested to know if we can guarantee that there actually is a solution to the system or do we need to make further assumptions? We get the linear system: S − AS = d ⇒ S = (I − A)−1 d If (I − A) is invertible we can solve and find a solution relatively easily. However we don’t know if there is one or if the solution could have negative values in S (which would be equivalent to producing a negative amount of a commodity). 43 • We use the fact that if the Neumann series ∞ X ∞ X An converges, then (I − A)−1 = n=0 An . n=0 • If we assume ρ(A) < 1 we get lim Ak = 0, Which is the same as saying the k→∞ Neumann series converge. • Since A is non-negative, ∞ X An is as well and we have guaranteed not only a n=0 solution, but also that S ≥ 0. • The question is what assumptions we need to make in order for ρ(A) < 1? We look at the column sums X ai,j , this can be seen as the total number of units X required to create one I − unit. It seems reasonable to assume ai,j ≤ 1 with strict i i inequality for at least one industry. Since otherwise the industries with go with a loss, and we can assume at least one goes with a profit. Then there is a matrix E ≥ 0, E 6= 0 with elements ei,j such that X ai,j + ei,j = 1, i (every column sum of A + E equal to one).If u is the unit vector with all elements equal to 1, we get: u> = u> (A + E) We assume ρ(A) ≥ 1 and use the Perron vector x of A to get: u> x = u> (A + E)x ⇔ x = Ax + Ex This can’t be true since Ex > 0 (E non-negative, non-zero) and Ax ≥ x (since ρ(A)x = Ax, ρ(A) ≥ 1). So if we assume the column sums of A to be less than or equal to one with at least one less than one, we get that ρ(A) < 1 and we are finished. ∞ X We also note that if A is irreducible, then An is even positive and not only nonn=0 negative. In other words if we increase the demand of any commodity, we will need to increase the supply of all commodities and not only the one we increased the demand for. Take a few minutes to think through the following questions regarding the new assumptions of the system. Apart from assuming the industries aren’t making a loss, and at 44 least one makes a profit, we also assumed A to be irreducible. This corresponds to every commodity being directly or indirectly needed to produce every commodity. This doesn’t seem unreasonable considering we work with large industries. 45 5 Stochastic matrices and Perron-Frobenius 5.1 Stochastic matrices and Markov chains A stochastic matrix, sometimes also called transition matrix or probability matrix, is used to describe the transitions between states in a Markov chain. Markov chains are used in applications in a large number of fields, some example are: • Markov chain Monte Carlo for generating numbers from very complex probability distributions in statistics. • Modeling of market behaviour, such as market crashes or periods of high or low volatility in economics. • Population processes in for example biology, queuing theory and finance. • Other less obvious places like in speech recognition, error correction, thermodynamics and PageRank from the previous chapter. Our aim here is to get an understanding of how and why Markov chains behave as they do using what we know of Perron-Frobenius theory. We start by defining stochastic matrices and Markov chains. Definition 5.1 A square non-negative matrix is • Row stochastic if every row sums to one. • Column stochastic if every column sums to one. • If it is both row and column stochastic we say that it is doubly stochastic. We will from now only say stochastic matrix when we in fact mean row stochastic matrix. We also notice that stochastic matrices need not be irreducible, however when they are we can immediately use Perron-Frobenius to find a lot of information regarding the matrix. A common example of doubly stochastic matrices are the permutation matrices since only one element in every row and column is non-zero and that element is one. Every stochastic matrix corresponds to one Markov chain: Definition 5.2 A Markov chain is a random process, with usually a finite number of states, where the probability of transitions between states only depend on the current state and not any previous state. A Markov chain is often said to be memoryless. • Formal expression of the fact that a discrete time Markov chain is memoryless: P (Xn+1 = x|X1 = x1 , X2 = x2 , . . . , Xn = xn ) = P (Xn+1 = x|Xn = xn ) 46 • Every stochastic matrix defines a Markov chain and every Markov chain defines a stochastic matrix. You can see a discrete time Markov chain as a sequence of states X0 , X1 , X2 , . . . happening one after another in discrete time and the state in the next step depends only on the current state. So if our Markov chain has 3 different states, then this corresponds to a 3 × 3 stochastic matrix with every element ai,j corresponding to the transition probability from state xi to state xj . I.e., ai,j = P (Xn+1 = j|Xn = i)). Another way to represent a Markov chain is with a directed graph by considering the stochastic matrix as the graph’s distance matrix, where the distance is now a probability rather than a true distance. In this case, with the edges representing transition probabilities, it’s natural to only consider graphs with positive weights. As an example of a Markov chain, let us consider this simple weather model: • We know the weather on day zero to be sunny. • If it’s sunny one day it will be sunny the next with probability 0.7, cloudy with probability 0.2 and rainy with probability 0.1. • If it’s cloudy it will be sunny with probability 0.3, cloudy with probability 0.4 and rainy with probability 0.3. • If it’s rainy it will be sunny with probability 0.2, cloudy with probability 0.3 and rainy with probability 0.5. Seeing this in text makes it quite hard to follow, especially if we would like to add more states. However since the weather tomorrow only depend on the weather today, we can describe the model with a Markov chain and corresponding stochastic matrix. We order the states in the stochastic matrix from sunny to rainy and put in the transition probabilities to get the stochastic matrix P and it’s corresponding graph below. The vector X0 represents the probability distribution of the weather on day zero. sunny 1 0.7 0.2 0.1 cloudy X0 = 0 P = 0.3 0.4 0.3 rainy 0 0.2 0.3 0.5 0.2 0.7 sunny cloudy 0.4 0.3 0.1 0.3 0.2 0.3 rainy 0.5 47 When we worked with the distance matrix and adjacency matrix of graphs we saw that we could take the powers of the matrix to find out if there is a path between vertices. In the same way we can take the powers of the transition matrix P and if a element pi,j ∈ Pk 6= 0 then there it’s possible to get from state i to state j in k steps in the corresponding Markov chain. However we can say even more in the case of stochastic matrices, namely we can tell not only if it’s possible to be in a state after k steps, we can tell how likely it is as well. We can easily see that if we multiply the initial state (or distribution) X> 0 with P from the right, we get the probability to be in all the states in the next step X> 1 . Generally since we work with probabilities the initial distribution should be non-negative as well as sum to one. 0.7 0.2 0.1 > 0.3 0.4 0.3 = 0.7 0.2 0.1 X> 1 = X0 P = 1 0 0 0.2 0.3 0.5 And we see there is only a probability 0.1 that it will rain tomorrow provided it is sunny today. Now we can of course multiply this new vector X1 in the same way to get the probabilities for the weather on the day after tomorrow. 0.7 0.2 0.1 > 2 > 0.3 0.4 0.3 = 0.57 0.25 0.18 X> 2 = X0 P = X1 P = 0.7 0.2 0.1 0.2 0.3 0.5 And we get the probability for rain two days from now as 0.18 if it’s sunny today. So to get the probability distribution π k of the weather on day k we can take the transition matrix to the power of k and multiply with the initial distribution X> 0 to the left. > k X> k = X0 P Sometimes the transposed transition matrix is used instead Xk = (P> )k X0 , both obviously work equally well. • One important question in the theory of Markov chains is whether a stationary distribution π > = π > P exist. Or, in other words, if there is a distribution for which the distribution in the next step is the same as the distribution in the current step. k > • We are also interested to know if the limit π = ( lim X> 0 P ) exists. If it exists, k→∞ it means that as long as we iterate the Markov chain for long enough, we will eventually arrive at a stationary distribution π > = π > P • The stationary distribution can also be seen as the distribution of the time the Markov chain spend in every state if it is iterated a large number of times. • Last, when the limit exists, we are also interested to know wheather it converges towards the same stationary distribution π regardless of initial state X0 . 48 In order to answer these questions we divide the stochastic matrices in 4 different classes: 1. Irreducible with lim (Pk ) existing. (P is primitive) k→∞ 2. Irreducible with lim (Pk ) not existing. (P is imprimitive) k→∞ 3. Reducible and lim (Pk ) existing. k→∞ 4. Reducible and lim (Pk ) not existing. k→∞ 5.2 Irreducible, primitive stochastic matrices If P is irreducible we can use Perron-Frobenius for irreducible non-negative matrices. From this we get that there is one eigenvector with only positive values, namely the Perron- Frobenius eigenvector with eigenvalue r. From the same theorem we also have that r = 1 since r is bounded by the minimum and maximum row sum, both of which are equal to one since P is a stochastic matrix. The stationary distribution is then the same as the left Perron vector since we have: y> = y> P, if r = 1. So if P is irreducible, a stationary distribution exist, and it is unique. We remember the requirement for an irreducible matrix to be primitive is for the limit lim (P/r)k exist. Since r = ρ(P) = 1 the limit lim (Pk ) exist if P is primitive, and thus k→∞ k→∞ k > π = ( lim X> 0 P ) must exist as well. k→∞ Using Perron-Frobenius for primitive matrices, for P we get the (right) Perron vector (x = Px) as x = e/n, where e is the vector with every element equal to one and n is the number of states . You can easily verify this since r = 1 and every row sums to one. If we let π = (π1 , π2 , . . . , πn ) be the left Perron vector (π > = π > P) we get the limit: lim Pk = k→∞ (e/n)π > = eπ > > 0 π > (e/n) Calculating the probability distribution of the limit gives: k > > > lim X> 0 (P ) = X0 eπ = π k→∞ And we see that if P is irreducible and primitive, not only does the stationary distribution π exist and is unique, we also get that the Markov chain converge towards this distribution regardless of initial distribution X0 . Additionally Perron-Frobenius gives us an easy to calculate way to find the stationary distribution since it’s equal to the left Perron vector of P. We give a short example using the earlier weather model. 49 We had the transition matrix: 0.7 0.2 0.1 P = 0.3 0.4 0.3 0.2 0.3 0.5 Since P is positive, it’s obviously irreducible and primitive as well and as such we know that there is a stationary distribution and the distribution of the Markov chain will converge towards it. For the limit we get: k > > > lim X> 0 P = X0 eπ = π k→∞ To find the stationary distribution we calculate the left Perron vector π which yields π = [0.46 0.28 0.26]> which is our stationary distribution. From this we can say that on average it will be sunny 46% of the time, cloudy 28% of the time and rainy 26% of the time. 5.3 Irreducible, imprimitive stochastic matrices Since P is irreducible we know that there at least exists a stationary distribution and that it’s unique. However since it’s imprimitive the limit lim Pk doesn’t exist and the k→∞ Markov chain therefore doesn’t necessarily converge towards the stationary distribution. However we can easily see that the (right) Perron vector is the same as in the primitive case (since every row sums to one). Similarly finding the stationary distribution can also be done in the same way for imprimitive matrices as with primitive ones as the left Perron vector. While we in the primitive case would also have been able to find the stationary distribution simply by iterating the chain a large number of times, this is not the case when working with imprimitive matrices. Even though the Markov chain doesn’t converge, we are still interested in seeing how the Markov chain behaves if we let it iterate for a long time.This is useful for examle when we want to sample from a distribution using an imprimitive Markov chain in Markov chain Monte Carlo methods. While the limit lim Pk doesn’t exist, what we are k→∞ interested isn’t necessarily the limit itself but the distribution of time the Markov chain spends in the different states if we let the Markov chain iterate a large number of times. For primitive matrices this is the same as the limiting distribution. In the imprimitive case we can calculate the Ces´aro sum instead: k 1X j P k→∞ k lim j=0 This converges and we get: k 1X j P = eπ > > 0 k→∞ k lim j=0 50 We multiply the limit with the initial vector X0 and get just like in the primitive case: X0 k 1 X j > > P = X> 0 eπ = π k→∞ k lim j=0 In other words the stationary distribution π exist and is the same regardless of the choice of initial state X0 . The Markov chain doesn’t converge to it but rather ’oscillates’ around it. More on why the Ces´ aro sum always exist for these matrices can be read about in [?] We give a short example using our previous 0 P = 1 small permutation matrix: 1 0 We already saw that it’s imprimitive and therefore the limit doesn’t exist. However limit of the Ces´ aro sum does exist and we still get the stationary distribution π > as the Perron vector. This gives the stationary distribution π > = [0.5 0.5]> . So while the Markov chain won’t converge it will on average spend half the time in state 1 and half the time in state 2. In fact we can easily see that it spends exactly every second step in state 1 and the other in state 2. 5.4 Reducible, stochastic matrices Reducible matrices are not as common to work with in practice as with irreducible matrices, often because the model itself is constructed in such a way as to ensure it produces an irreducible matrix. However sometimes it’s desirable to use a reducible chain, for example using an absorbing Markov chain (one or more states which we can’t leave after the Markov chain ends up in it) in risk analysis. Since we can’t immedietly apply Perron Frobenius to reducible Markov chains, we for example don’t know if there exists any stationary distribution or maybe there even exist many. What we can do is try to get as close as possible to the irreducible case. We recall that a reducible matrix A could be permuted into upper triangular form by a permutation matrix P: X Y > P AP = 0 Z If X or Z is reducible as well we can permutate the original matrix further untill we end up with: X1,1 X1,2 · · · X1,k 0 X2,2 · · · X2k > P AP = . .. .. .. .. . . . 0 0 0 Xk,k 51 Where Xi,i is either irreducible or a single zero element. One thing we can immediately see already is that if we leave for example block X2,2 we can never come back since none of the states in any block Xi,i have any connection to any of the previous blocks Xj,j , j < i. Last we look for blocks where only the diagonal block Xi,i contains non-zero elements and permutate those to the bottom so we end up with: > P AP = X1,1 0 .. . X1,2 X2,2 .. . ··· ··· .. . X1,j X2j .. . X1,j+1 X2,j+1 .. . X1,j+2 X2,j+2 .. . 0 0 0 0 0 ··· Xj,j 0 Xj,j+1 Xj+1,j+1 Xj,j+2 0 0 .. . 0 .. . ··· 0 .. . 0 .. . Xj+2,j+2 .. . 0 0 0 0 0 ··· ··· ··· ··· ··· ··· ··· .. . .. . ··· X1,k X2,k .. . Xj,k 0 0 .. . Xk,k This is called the canonical form of the stochastic matrix and we notice a couple of things. First we remind ourselves that this is still the same Markov chain, we have just relabeled the vertices by permutating the matrix using a permutation matrix and we remind ourselves that the blocks Xi,j are matrices and not necessarily single elements. We also note that this form is generally not unique since we could swap the position of two different blocks in the bottom right part or sometimes in the top left as well. Using this form we can say some things about the behaviour of the Markov chain which might not be as apparent otherwise. We call the states in the upper left the transient states. Whenever we leave one of the diagonal blocks Xi,i , i ≤ j we can never get back to the same or any previous block since all elements below the diagonal blocks are zeroes. Also all the blocks in the top left Xi,i at least on of the blocks Xi,i+1 , . . . , Xi,k are non zero (or they would belong to the bottom right part). This means that if we start in any state in the top left part, we will eventually end up in one of the blocks in the bottom right. We call the states in the bottom right the ergodic states or the absorbing states, since once you enter any of the blocks in the bottom right part you can never leave it. Either because it’s a single vertex only linking back to itself, or because it’s in itself an irreducible stochastic matrix. The first thing to note is that the initial state obviously influences the behaviour of reducible Markov chains, for example if we start in one of the ergodic states we obviously get a different sequence than if we start in another ergodic state. So even if there would exist a stationary distribution we can expect it to depend on the initial state. Often with reducible Markov chains we aren’t as interested in finding where the chain is after a long time (especially if there is only one absorbing state), but we are often more 52 interested in seeing how long time we expect it to take in order to reach a absorbing state or the probability to end up in a certain absorbind state when there are more than one. One example would be when modelling the failure rates of a machine, we know it will fail eventually, but for how long can we expect it to run? Assuming a reducible stochastic matrix A is written in canonical form. X Y A= 0 Z We can then find the limit lim Ak if all blocks in Z are primitive matrices or the k→∞ corresponding limit of the Ces´ aro sum if at least one of them are imprimitive. We look at the case where all blocks in Z is primitive: The limit can then be written as: 0 (I − X)−1 YE lim A = 0 E k→∞ k > eπ j+1 0 ··· .. 0 . eπ > j+2 E= . . .. .. ··· 0 ··· 0 0 0 .. . eπ > k Obviously we don’t find a unique stationary distribution (since every block in Z is absorbing), however the limit is still very useful. First of all we see that E contains the individual Perron vectors for the blocks in Z. So if we start in one of the absorbing blocks we find a stationary distribution unique to that block. In fact in that case we can obviously disregard the rest of the Markov chain and only concern ourselves with that one block if we have the transition matrix in canonical form. Secondly the elements ap,q of (I − X)−1 Y holds the probability that starting in p we leave the transient states and eventually hit state q in one of the absorbing blocks. If q is a single absorbing state ap,q is thus the probability to end up in q when starting in p. This is commonly called the hitting probability to hit q starting in p. The vector (I − X)−1 e is of particular interest as well. It holds the average number of steps needed to leave the transient state, we call this the hitting times for reaching the absorbing states. We take a short look at hitting times and hitting probabilities in the next section from a slightly different perspective to understand why this is correct. 5.5 Hitting times and hitting probabilities When we have a reducible Markov chain (with a finite number of states) we know that we will eventually reach one absorbing state. We consider a stochastic matrix P with 53 transient states D and two absorbing states A, B. We want to find the hitting probability hi that we end up in A starting in d ∈ D. We illustrate how to find such a hitting probability with a short example where we consider the Markov chain describes by the graph below. 2/3 1/3 1 1/3 2 3 4 2/3 We ask ourselves: what is the probability to end up in 1 starting in 2? We write this probability as hi = Pi (hit 1). After one step we are either in state 1 with probability 1/3 or we are in 3 with probability 2/3. In other words we could write h2 as: 1 2 h2 = P2,1 h1 + P2,3 h3 = h1 + h3 3 3 But h1 is the probability to hit 1 starting in 1, since we are already in 1 this probability is obviously 1. We then continue and find the expression for h3 which can be written: 1 2 h3 = P3,4 h4 + P3,2 h2 = h4 + h2 3 3 Since state 4 is absorbing so the probability h4 to hit 1 from 4 is zero and we get the equation system: 1 2 h2 = + h3 3 3 h3 = 2 h2 3 Solving this gives the hitting probability h2 = 3/5. While this method can be done by hand with a low number of states, it is mostly used to find a recurrence relation when we have an infinite number of states such as a random walk on the integers. If we have a finite number of states and thus can write out our transition matrix P we can find the hitting probability using the following: Theorem 5.1. Consider the Markov chain with transition matrix P and n states. Let A be a subset of states. The vector of hitting probabilities hA i = Pi (hit A) is the minimal non-negative solution to the linear system: A i∈A hi = 1,n X A pi,j hA i∈ /A j hi = j=1 With minimal solution we mean that for any other solution x to the system we have xi ≥ hi . 54 By writing P in canonical form we can verify that this corresponds to the result in the previous section. The minimality is needed in order to get hA i = 0 in other absorbing states than those in the absorbing block we seek. When we solve the system we can set hA i = 0, if i is any absorbing state not in A. We note that it does not matter whether A consists of absorbing states, since we are not interested in what happens after we have entered A. In terms of the previous example we could use the theorem to find the probability that we reach state 3 at some point before reaching an absorbing state. Next we continue with our example but instead ask ourselves: What is the expected number of steps starting in 2 until we hit either of the absorbing states? We define the hitting time: ki = Ei (time to hit {1,4} ) As with the hitting probability we look at the next step but now also add a 1 for the time it takes to reach the next step: 1 2 2 k2 = P2,1 k1 + P2,3 k3 + 1 = k1 + k3 + 1 = k3 + 1 3 3 3 Where we noted that k1 = 0 since we then have already reach one of the absorbing states. Next we continue with k3 and get: 2 2 1 k3 = P3,4 k4 + P3,2 k2 + 1 = k4 + k2 + 1 = k2 + 1 3 3 3 Solving this system gives the expected number of steps until absorption as k2 = 3. In the same way as for hitting time we get in general: Theorem 5.2. Consider the Markov chain with transition matrix P, n states and absorbing states A. The vector of hitting times kiA = Pi (hit A) is the minimal non-negative solution to the linear system: A i∈A ki = 0, n X A pi,j kjA i ∈ /A ki = 1 + j=1 Once again we can write the system in canonical form to verify that this corresponds to the result in the previous section. We do note however that care needs to be taken so that we don’t have other absorbing states outside of A. This since we would then get the relation for a single absorbing state outside of A as kj = 1 + kj . 5.6 A short look at continuous time Markov chains Here we will seek to give a short introduction of continuous time Markov chains and their relation to what we already know about discrete time chains without going to much into probability theory. This is not meant to give a full understanding of continuous time 55 chains, but rather highlight similarities and relations with discrete time chains and previous chapters. A continuous time Markov chain can be seen as a discrete time chain, but instead of jumping to a new state at fixed intervals 1, 2, 3, . . . we instead jump to a new state randomly in time. With randomly in time we mean that the jump from a state is exponential distributed with intensity λ (which might be different between vertices). Exponential distributed jumps mean two important things, one if we assume there is only one jump T in the interval T ∈ (0, t) then the jump is uniformly distributed on this interval P (T = x|x ∈ (0, t)). Second the jumps are memoryless which means that if there was no jump in (0, s) then the probability of a jump at a point in time after s is independent of s: P (T ∈ (0, t)) = P (T ∈ (s, s + t)|T ≥ s). We will only concern ourselves with chains with a finite state space. With infinite state space there are more things that needs to be taken into consideration such as chances for the chain to ”explode” or when the chance to eventually get back to where we started is less then one (even if it’s irreducible). For a continuous time chain with a finite state space we can make a graph of the chain. Now however we no longer have the probabilities to jump to another state as weights on the edges, but rather the intensity λ for the jump. We also naturally have no self-loops since this is handled by the intensity of the jumps to other states. An example of what it could look like can be seen in subsection 5.6 below. 1 1 2 2 1 1 1 3 4 2 In continuous time we talk about two different matrices. First we have the Q matrix which is the negative Laplacian matrix where we take into consideration the weight’s of the edges. This means that we get elements qi,j = λi,j , i 6= j where λi,j is the intensity of jumps from state i to state j. The diagonal elements we choose such that the sum of every row is zero. For the Markov chain described by the graph above we get Q matrix: −3 1 2 0 0 −1 0 1 Q= 1 0 −2 1 0 0 2 −2 If Q is irreducible and finite then there exist a stationary distribution π which we can find as π > Q = 0. In continuous time we can’t have any periodicities and irreducibility 56 is actually enough for a finite chain to guarantee convergence towards this distribution. While we can solve this system to find the stationary distribution there is also another more familiar matrix we can work with. Namely the transition matrix P(t) = etQ . The element pi,j (t) in P (t) describes the probability of being in state j at time t if the chain started in state i at time 0. We will look more at the exponential of a matrix and how to calculate it later in ??. For the moment it is enough to know that we get : P(t) = etQ = ∞ X (tQ)k n=0 k! P(t) is then a stochastic matrix for all t ≥ 0 and we can find the stationary distribution (using any t > 0) in the same way as we did for a discrete time Markov chain. For Q above we get using t = 1 0.1654 0.2075 0.3895 0.2376 0.0569 0.3841 0.2533 0.3057 P(1) = 0.1663 0.0982 0.4711 0.2645 0.1394 0.0569 0.4720 0.3317 And we can also find the stationary distribution by solving the linear system π > P(1) = π by for example finding the Perron vector π. It’s easy to see that π = [0.1429 0.1429 0.4286 0.2857]> is a stationary distribution and that it satisfies both π > P(t) = π and π > Q = 0. > 5.7 Exercises 5.1. Consider the stochastic matrix 0 0.5 0.5 0 1 P= 0 0.5 0 0.5 • Show that P is irreducible. • If we start in state 1, what is the probability to be in state 3 after 3 steps (P (X3 = s3 |X0 = s1 )). • What is the stationary distribution of this Markov chain? 5.2. Consider the stochastic matrix 0 0.5 0 0.5 1 0 0 0 T= 0 0 1 0 0.25 0.25 0.25 0.25 57 • Write T in canonical form. • Calculate the hitting times of this Markov chain. What is the expected number of steps needed to reach the absorbing state if we start in the first state? 5.8 Applications 5.8.1 Return to the Leontief model We return to the Leontief Input-Output model in subsubsection 4.3.1 for modeling a economic system composed of n industries.. We remember we got a linear system S − AS = d. We found that if A was irreducible (every industry directly or indirectly depend on all others) and every column sum to one or less, with at least one column with sum less than one (at least one goes with a profit and none goes with a loss), we could guarantee that we could satisfy any given demand d. We are now interested to answer some more questions about the system, such as what is the gross national product (GNP) of the system? To answer this we will look at the system as a Markov chain with transition matrix Q = A> . We note that we take the transpose of the previous system matrix such that every row now sums to one or less rather than column and elements pi,j now corresponds to the amount of J−units required by industry I to create one I−unit. To account for the demand vector, which essentially takes units of produced commodities out of the system we introduce a new state d in the transition matrix. Since any produced commodity that enters the ”public” Xcan’t come back into the system we let d be an absorbing state such that Qi,d = 1 − qi,j . j In other words we modify the transition matrix by adding a elements to the column corresponding to the demand d such that every row of A> sum to one. We get a transition matrix of the form: q1,1 . . . q1,n q1,d .. .. .. .. . . . Q= . qn,1 . . . qn,n qn,d 0 ... 0 1 Since the top left n × n part is irreducible, and at least one of qi,d , . . . qn,d > 0 since at least one goes with a profit we see that we can reach the absorbing state d from any other state. Our transition matrix therefore describes an absorbing Markov chain it’s also easy to see that any initial supply S leads to supplying some demand d. GNP is the total supply S produced to both satisfy the demand as well as the internal demand of the system. If we look at the mean number of steps required to reach the absorbing state from the transient states: T = (I − A> )−1 u. • We can see the elements ti as the expected production in units to create one product unit of industry i. 58 • Multiplying the demand vector d with T: (T> d) gives the expected total produced units. i.e. the GNP. 59 6 Linear spaces and projections In this chapter, which will be a bit more theoretical, we will take a look at linear spaces and their relation to matrices. This in order to get a better understanding of some of the underlying theory as well as some commonly used methods or concepts used later. We will start by defining linear spaces and inner product spaces and later in the last part take a closer look at projections and how they can be used to for example find a basis of a matrix. 6.1 Linear spaces We start by defining a linear space: Definition 6.1 A linear space or vector space over a field F is a set V together with two binary operations (+, ∗) defined such that: u, v ∈ V ⇒ u + v ∈ V, and u ∈ V, a ∈ F ⇒ a ∗ u ∈ V With the following being true, where we let u, v, w be vectors in V and a, b be scalars in F . For vector addition (+) we have: • Associativity u + (v + w) = (u + v) + w. • Commutativity: u + v = v + u. • Existence of identity element 0 ∈ V such that v + 0 = v for all v ∈ V . • Existence of a inverse element −v ∈ V such that v + (−v) = 0. For scalar multiplication (∗) we have: • Distributivity with respect to vector addition: a ∗ (u + v) = a ∗ u + a ∗ v. • Distributivity with respect to field addition: (a + b) ∗ v = a ∗ v + b ∗ v • Compatibility with field multiplication: a ∗ (b ∗ v) = (ab) ∗ v. • Existence of identity element 1 ∈ F such that 1 ∗ v = v for all v ∈ V . Usually we omit the sign (∗) for multiplication as we are used to. In order to not go to much into algebraic structures we consider a field a set with two operations (+, ·) following the ”usual” rules such as associativity and commutativity. To read more about fields and other algebraic structures we refer to a book on algebra such as [?]. Two examples of fields are the real and complex numbers R, C. 60 We take a closer look at the definition of a linear space and the requirements on the two operations (+, ∗) for V to be a linear space. First of all we see that + is a operation between two vectors in V such that the result is a new vector also in V . We also need it to be associative and commutative such that we can add in whatever order we want and/or swap place of the two vectors. There also exist a zero(identity) and inverse element. These are all things we are used to from the real numbers. Second we have another operation ∗ between a scalar in F and a vector in V such that the result is a vector in V . Similarly to + we have a few requirements on the operation such as distributivity with respect to both vector and field addition. This means means for example we can either add two vectors and then multiply or multiply with both vectors and then add the result. We also have Compatibility with field multiplication as well as the existence of a multiplicative identity element. Some examples of linear spaces using the common addition and multiplication operations are: • Rn . • Cn . • Matrices of size n × n with elements in F : Mn×n (F ). • {polynomials} In the future, unless otherwise noted we will assume the common addition and multiplication operation is used. We take a look at R2 with vectors u = (u1 , u2 ), v = (v1 , v2 ), w = (w1 , w2 ) ∈ V and scalars a, b ∈ F . If we add two vectors we get: (u1 , u2 ) + (v1 , v2 ) = (u1 + v1 , u2 + v2 ), since the addition of two real number (for example u1 + v1 ) is a new real number this obviously belong to R2 as well. Both assiciativity and commutativity immediately follows from addition in R as well. The identity element with respect to addition is obviously the zero vector 0 = (0, 0) since (u1 , u2 ) + (0, 0) = (u1 + 0, u2 + 0) = (u1 , u2 ). Similarly we find inverse elements simply by multiplying every element of the vector with −1. For scalar multiplication a ∗ (v1 , v2 ) = (av1 , av2 ) we can show that it fullfills the requirements in a similar way and as such is left to the reader. We now introduce something called subspaces of linear spaces, these can be seen as another linear space contained within a linear space. Definition 6.2 If V is a vector space, a subset U ⊆ V is called a subspace of V if U is itself a vector space with the same operations (+, ∗) as V . We also get a way to test if a subset U is a subspace of a linear space V by checking a just a couple of different things rather than all the requirements in the original definition for a linear space. 61 Theorem 6.1. The subset U ⊆ V is a subspace of V if and only if it satisfies the following 3 conditions. • 0 lies in U , where 0 is the zero vector (additive identity) of V . • If u, v ∈ U then u + v ∈ U . • If u ∈ U and a ∈ F then au ∈ U . For example we look at the linear space R2 and want to see if the points on the line y = 2x, x, y ∈ R is a subspace of R2 . We can write the points on the line as: U = {a(1, 2)|a ∈ R} U is obviously a set in R2 so we only need to show the three requirements above. If we choose a = 0 it’s clear that 0 ∈ V also lies in U . Adding two vectors in U gives a(1, 2) + b(1, 2) = (a + b)(1, 2) which is obviously in U since a + b ∈ R if a, b ∈ R. Last we multiply with a scalar b ∗ a(1, 2) = (ba)(1, 2) which also obviously belong to U since ba ∈ R. Definition 6.3 A vector v such that: v = c1 v1 + c2 v2 + . . . + cm vm = m X cj vj j=1 is called a linear combination of vectors {v1 , v2 , . . . , vn } ∈ V with coefficients in F if {c1 , c2 , . . . , cm } ∈ F . We illustrate it with an example: Considering the two vectors u = (1, −2, 1) and v = (1, 1, 2). Then the vector (3, −3, 4) is a linear combination of u, v since 2u + v = 2(1, −2, 1) + (1, 1, 2) = (3, −3, 4). However the vector (1, 0, 1) is not a linear combination of these vectors since if we write it as an equation system su + tv = (1, 0, 1) we get: 1=s+t 0 = −2s + t 1 = s + 2t The only way to satisfy the first and third equation is if s = 1, t = 0 however then the middle equation is not fulfilled. We often talk about not one linear combination of vectors but all linear combinations of a set of vectors: Definition 6.4 If {v1 , v2 , . . . , vm } is a set of vectors in V . Then the set of all linear combinations of these vectors is called their span denoted: span(v1 , v2 , . . . , vm ) 62 • If V = span(v1 , v2 , . . . , vm ) then we call these vectors a spanning set for V . For example given the two vectors u = (1, −2, 1) and v = (1, 1, 2) above we get their spanning set: span(u, v) = {su + tv|s, t ∈ R} = (s + t, −2s + t, s + 2t) So the spanning set is all vectors which we can write in this way. We recognize that the spanning set above actually describes a plane in R3 . In fact all spanning sets in a vector space V are subspaces of V . Theorem 6.2. For a span U = span(v1 , v2 , . . . , vm ) in a vector space V we have: • U is a subspace of V containing each of v1 , v2 , . . . , vm . • U is the smallest subspace containing v1 , v2 , . . . , vm , in that all subspaces containing v1 , v2 , . . . , vm must also contain U . We prove the first statement: Proof. We use the test described earlier to see if this subset is also a subspace. • U obviously contains 0 = 0v1 + 0v2 + . . . + 0vm . • Given u = a1 v1 + a2 v2 + . . . + am vm and v = b1 v1 + b2 v2 + . . . + bm vm and c ∈ F . • u + v = (a1 + b1 )v1 + (a2 + b2 )v2 + . . . + (am + bm )vm ∈ U • cu = (ca1 )v1 + (ca2 )v2 + . . . + (cam )vm ∈ U . • And we have proven that U is a subspace of V . We take a look at the spanning set span{(1, 1, 1), (1, 0, 1), (0, 1, 1)} and want to find if this spans the whole of R3 ? span{(1, 0, 1), (1, 1, 0), (0, 1, 1)} is obviously contained in R3 and it’s clear that R3 = span{(1, 0, 0), (0, 1, 0), (0, 0, 1)}. We can then prove that (1, 0, 0), (0, 1, 0), (0, 0, 1) is contained in our spanning set and use Theorem 6.2 to prove that it contains R3 . We do this by writing (1, 0, 0), (0, 1, 0), (0, 0, 1) as a linear combination of our vectors: (1, 0, 0) = (1, 1, 1) − (0, 1, 1) (0, 1, 0) = (1, 1, 1) − (1, 0, 1) Using the two first found vectors we find the last. (0, 0, 1) = (1, 1, 1) − (1, 0, 0) − (0, 1, 0) 63 It’s clear that there are generally more than one way to write a linear space as a spanning set. We will now move on to look at two different types of spanning sets. We compare two spanning sets in R2 : span((1, −1), (1, 1), (1, 2)) = R2 span((1, 2), (2, 3)) = R2 In the first spanning set we can represent the vector (3, 2) as a linear combination of the spanning vectors in two ways: (3, 2) = 1(1, −1) + 1(1, 1) + 1(1, 2) (3, 2) = 0(1, −1) + 4(1, 1) − 1(1, 2) However in the second spanning set we get: (3, 2) = s(1, 2) + t(2, 3) = (s + 2t, 2s + 3t) ⇒ s = −5 t=4 In fact, any vector (a, b) have only one unique representation (a, b) = s(1, 2) + t(2, 3). We introduce the notion of linearly dependent or independent sets of vectors. Definition 6.5 A set of vectors {v1 , v2 , . . . , vm } is called linearly indepependent if the following is true: a1 v1 + a2 v2 + . . . + am vm = 0 ⇒ a1 = a2 = . . . = am = 0 In other words, a set of vectors {v1 , v2 , . . . , vm } are linearly independent if the only way to represent 0 in these vectors is with all coefficient zero: 0 = 0v1 + 0v2 + . . . + 0vm One important consequence of this is that it also means that any other vector v can only be representes as a linear combination of {v1 , v2 , . . . , vm } in one way as well. Otherwise we could find two representations of 0 by subtracting v − v using two different representations of v. If a set of vectors isn’t linearly independent they are linearly dependent. Definition 6.6 A set of vectors {v1 , v2 , . . . , vm } is called linearly depependent iff some vj is a linear combination of the others. We consider three vectors {(1 + x), (2x + x2 ), (1 + x + x2 )} in (P2 ) and want to find if they are linearly independent or not by finding a linear combination of them that result in the zero vector. 0 = s(1 + x) + t(2x + x2 ) + u(1 + x + x2 ) 64 We get a set of linear equations: s s + 2t t + + + u u u =0 =0 =0 It’s clear that the only solution to the system is s = t = u = 0 and we see that the vectors are linearly independent. Considering this relation between solving an equation system and linear independence we can easily proove the following usefull theorem. Theorem 6.3. If a vector space V can be spanned by n vectors. Then for any set {v1 , v2 , . . . , vm } of linear independent vectors in V , m ≤ n. In the proof which we leave to the reader we show that the m vectors are linear depenm X dent if m > n by setting up the equation system xj vj = 0 We represent the m linear j=1 independent vectors as a linear combination of a span of n vectors. This will result in a system of n equations with m unknown variables, which has a nontrivial solution. In practise the theorem means that every set of linear dependent set of vectors that span a vector space V have the same number of vectors, we call this a basis. Definition 6.7 A set of vectors {v1 , v2 , . . . , vm } is called a basis in V if the following is true: • {v1 , v2 , . . . , vm } is linear independent. • V = span{v1 , v2 , . . . , vm }. We summarize our findings of bases in a linear space. • If {v1 , v2 , . . . , vm } and {w1 , w2 , . . . , wn } are two bases in V , then m = n. • In other words, all bases in a vector space contain the same number of vectors. • If V have a basis of n vectors. We say that V have the dimension n: dim V = n • We say that a vector space is finite dimensional if either V = 0 or dim V < ∞. In fact not only for all bases in a linear space V is it true that they all have the the same amount of vectors, any set of that same amount of linear independent vectors in V is a basis in V . Theorem 6.4. If dim V = n < ∞ and {v1 , v2 , . . . , vn } is a set of linear independent vectors, then {v1 , v2 , . . . , vn } is a basis in V . 65 Proof. It’s already linear independent so we then need to prove that {v1 , v2 , . . . , vn } spans V . If we assume {v1 , v2 , . . . , vn } does not span V , we could then choose a new vector vn+1 outside span{v1 , v2 , . . . , vn }. Then {v1 , v2 , . . . , vn , vn+1 } is linear independent. But we already know there is a basis {w1 , w2 , . . . , wn } which spans V since dim V = n, and for any set of linear independent vectors {v1 , v2 , . . . , vm }, we have m ≤ n. So we get a contradiction and {v1 , v2 , . . . , vn } must span V . One common solution to problems is to fit a new more appropriate basis to the problem and use that instead, this guarantee that any set of same amount of linear independent vectors is also a basis making it unnecessary to check if our new proposed basis is actually a basis or not. In essence it makes it easy to find a new basis if we already have one basis to start with. For example changeing to spherical coordinates, finding eigenvectors of a matrix, when finding many matrix decompositions, or making a fourier approximation of a function can all be seen as more or less as finding a new basis to the problem. Before continuing we give a short relation to matrices and how they can and are often used to represent a basis. Theorem 6.5. For a (real) n × n matrix M the following properties are equivalent: • The rows of M are linearly independent in Rn • The columns of M are linearly independent in Rn • The rows of M span Rn . • The columns of M span Rn . • M is invertible. 6.2 Inner product spaces We note that this far we have only seen vectors as a set of elements in F . However when we talk about vectors we are often talking about length of or angle between vectors. To be able to do this we need to add something to our linear space, namely we add another operation called a inner product and call it a inner product space. Definition 6.8 An inner product on a vector space V over F is a map h·, ·i that takes two vectors u, v ∈ V and maps to a scalar s ∈ F . We denote the inner product of two vectors u, v: hu, vi for which the following statements must be true: • hu, vi = hv, ui, (symmetric in Rn ). • hau, vi = ahu, vi. 66 • hu + v, wi = hu, wi + hv, wi. • hu, ui > 0 for all u 6= 0 in V . • A vector space V with an inner product h·, ·i we call a inner product space. The most commonly used inner product is the dot product sometimes called scalar product: u · v = u1 v1 + u2 v2 + . . . un vn ∈ Rn Using the definition above we can show that the dot product is in fact a inner product. • For the first property we get u · v = u1 v1 + u2 v2 + . . . un vn = v1 u1 + v2 u2 + . . . vn un . Which is obviously true since xy = yx in R. • For the second we get: au·v = au1 v1 +au2 v2 +. . . aun vn = a(v1 u1 +v2 u2 +. . . vn un ). Which also immediately follows from multiplication in R. • Next we get (u + v) · w = (u1 + v1 )w1 + (u2 + v2 )w2 + . . . (un + vn )wn = (v1 w1 + v2 w2 + . . . vn wn ) + (u1 w1 + u2 w2 + . . . un wn ) = (v · w) + (u · w) and the third statement is fulfilled as well. • Last we note that hu, ui is the sum of the square of all elements. Since the square of a real number r is positive if r 6= 0 the only way for the dot product to equal zero is if u = 0. But although the dot product is the most commonly used inner product, there are also other inner products and such as for example the following on the polynomials of order n Pn : Z hp, qi = 1 p(x)q(x) dx 0 Where p, q are polynomials in Pn . It’s easy to show that this is a inner product from the definition as well for the interested reader. √ You should be familiar with the Euclidean norm |v| = the norm of a inner product space in a similar way. v> v of a vector, we define Definition 6.9 If h·, ·i is an inner product in the vector space V , the norm (sometimes called length) kvk of a vector v in V is defined as: p kvk = hv, vi We immediately see that when using the dot product we get the Euclidian norm. The norm in a inner product space share some properties with those of a vector norm: 67 • kvk ≥ 0, kvk = 0 if and only if v = 0. • kcvk = |c| kvk. • kv + wk ≤ kvk + kwk. But we also have one additional important property of norms in a inner product space called the Schwarz inequality: • If u, v are vectors in V , then: hu, vi2 ≤ kuk2 kvk2 • With equality if and only if one of u, v is a scalar multiple of the other. Now that we have a notion of length (norm) in our inner product space we might also want the distance between vectors. Once again we start by looking at how we calculate distance using the dot product and the common geometric distance (in R3 ). Using the dot product we get the distance between two vectors d(u, v) as: d(u, v) = p (u1 − v1 )2 + (u2 − v2 )2 + (u3 − v3 )2 p = (u − v) · (u − v) Which we recognize as the Euclidean norm |u − v|. We define the distance between vectors in a general inner product space in the same way: Definition 6.10 The distance between two vectors u, v in an inner product space V is defined as: d(u, v) = ku − vk For the distance d(u, v) we have the following properties for vectors u, v, w in the inner product space V . • d(u, v) ≥ 0 with equality iff u = v. • d(u, v) = d(v, u). • d(u, w) ≤ d(u, v) + d(v, w). While we can now calculate the length of a vector or the distance between two vectors we still have no idea how to find out if for example two vectors are parallel or orthogonal. To do this we need to define the angle between vectors. Using Schwarz inequality we define the angle between two vectors u, v. Definition 6.11 If 0 ≤ θ ≤ π is the angle between two vectors we have: cos θ = 68 hu, vi kuk kvk We remember Schwarz inequality saying: hu, vi2 ≤ kuk2 kvk2 This gives: ⇔ hu, vi kuk kvk 2 ≤ 1 ⇔ −1 ≤ hu, vi ≤1 kuk kvk Which we then use to define the angle between two vectors. One very important case is when the inner product is zero (hu, vi = 0) and therefor the angle π/2 radians between them. This is not only interesting geometrically but for general inner product spaces as well. Definition 6.12 Two vectors u, v in a inner product space V are said to be orthogonal if: hu, vi = 0 More important is if we have a whole set of vectors, all which are orthogonal to eachother: Definition 6.13 A set of vectors {v1 , v2 , . . . , vm } is called an orthogonal set of vectors if: • vi 6= 0 for all vi . • hvi , vj i = 0 for all i 6= j. • If kvi k = 1 for all i as well, we call it an orthonormal set. If we have a orthogonal set {v1 , v2 , . . . , vm } and want a orthonormal set we can easily find a orthonormal set simply by normalizing every vector. v1 v2 vm , ,..., kv1 k kv2 k kvm k Given a orthogonal set we have a couple of important properties, first of all it’s also linear independent. Theorem 6.6. If {v1 , v2 , . . . , vm } is an orthogonal set of vectors then it is linear independent. We now have another way to find a basis in a inner product space V , rather than trying to find a set of linear independent vectors, we could instead look for a orthogonal set of vectors. The found orthogonal set will then automatically be linearly independent as well. We will look more into how to do this in later when we take a closer look at projections. Additionally we have one additional property concerning the norm and orthogonal vectors. 69 Theorem 6.7. If {v1 , v2 , . . . , vm } is an orthogonal set of vectors then: kv1 + v2 + . . . + vm k2 = kv1 k2 + kv2 k2 + . . . + kvm k2 We give a short motivation showing it for two vectors, the proof can then easily be extended for n vectors as well. Given two orthogonal vectors v, u we get: kv + uk2 = hv + u, v + ui = hv, vi + hu, ui + 2hv, ui Since v, u are orthogonal hv, ui = 0 and we get: kv + uk2 = hv, vi + hu, ui = kvk2 + kuk2 A matrix whose rows or columns form a basis of orthonormal vectors is called a unitary matrix. We remember that for a unitary matrix we have: Theorem 6.8. For a n × n matrix U the following conditions are equivalent. • U is invertible and U−1 = UH . • The rows of U are orthonormal with respect to the inner product. • The columns of U are orthonormal with respect to the inner product. • U is a unitary matrix. If U is a real matrix using the ordinary dot product we say that U i an orthogonal matrix if it satisfies the conditions above. Worth noting is that if U is real the hermitian of a matrix is the same as the transpose and we get U−1 = UT . One example of unitary matrices are the permutation matrices, in fact we can see the unitary matrices as a sort of generalization. In the same way as multiplying with a permutation matrix PAPT can be seen as changing from one basis to another by reordering the basis vectors (relabel a graph). Multiplication with a unitary matrix UAUH can also be seen as a change of basis. This is often used in practice by given M find a suitable matrix decomposition M = UAUH and then continue working with the presumably easier A instead. 6.3 Projections We will now look closer at why we usually want to work with a orthogonal basis and how to construct such a basis. We still haven’t given a reason why we would want to find a new orthogonal basis given that we usually start with the standard orthonormal basis represented by the identity matrix. This however wil mainly have to wait until the examples and applications. One of the large advantages of orthogonal bases is in the ease of finding a representation of any vector in this basis. We have already seen that the representation is unique since it’s basis vectors are linearly independent and we therefor find the new representation by solving a linear system. However if we have an orthogonal basis {e1 , e2 , . . . , em } of a inner product space V . Then for any vector v ∈ V : 70 • v= hv, e1 i hv, e2 i hv, em i e1 + e2 + . . . + em ke1 k2 ke2 k2 kem k2 • We call this the expansion of v as a linear combination of the basis vectors {e1 , e2 , . . . , em }. This gives us an easy way to express v given an orthogonal basis without the need to solve a linear equation system. Even better, if the basis is orthonormal it’s even easier since then all kei k2 = 1 and we simply need to take the inner product of the vector v and all the basis vectors. The next step is in constructing an orthogonal basis, we will take a look at projections and how to use these to create a basis. While we have already constructed bases many times before, for example every time you use Gauss-elimination, we are now interested in how we can ensure that we create a orthogonal basis and rather than only a linear independent one. One method to do this is the Gram-Schmidt process, however before looking at the method we must first show that such a basis actually exist. Theorem 6.9. In every inner product space V with dim V < ∞ there exist an orthogonal basis. Proof. We let V be an inner product space with dimension dim V = n and show it using induction. • If n = 1 then any basis is orthogonal. • We then assume every inner product space of dimension n have an orthogonal basis. Then if dim V = n+1 and {v1 , v2 , . . . , vn+1 } be any basis of V . If U = span(v1 , v2 , . . . , vn ) and {e1 , e2 , . . . , en } is an orthogonal basis of U . It’s then enough to find any vector v such that hv, ei i = 0 for all i. since then {e1 , e2 , . . . , en , v} will be an orthogonal basis in V . We consider: v = vn+1 − t1 e1 − t2 e2 − . . . − tn en Where we need to find all ti . Since vn+1 is not in U , v 6= 0 for all choices of ti . So if we can find ti for all i such that hv, ei i = 0 we are finished. For the orthogonality with ei we get hv, ei i = hvn+1 , ei i − t1 he1 , ei i − . . . − ti hei , ei i − . . . − tn hen , ei i = hvn+1 , ei i − 0 − . . . − ti kei k2 − . . . − 0 = hvn+1 , ei i − ti kei k2 And we get hv, ei i = 0 if ti = hvn+1 , ei i . kei k2 71 We noice that the proof not only shows that there is an orthogonal basis, but also gives us a way to find one given any other basis. It’s is recommended that you take some time to really understand the proof and why it works. We use the results here in order to create the Gram-Schmidt process used in order to find an orthogonal basis. We summarize the previous finding in Theorem 6.10. Let {e1 , e2 , . . . , en } be an orthogonal set of vectors in a inner product space V and v be a vector not in span(e1 , e2 , . . . , en ). • Then if: en+1 = v − hv, e1 i hv, e2 i hv, en i e1 − e2 − . . . − en 2 2 ke1 k ke2 k ken k2 • {e1 , e2 , . . . , en , en+1 } is an orthogonal set of vectors. Using this we construct the Gram-Schmidt process: • Let V be an inner product space, and {v1 , v2 , . . . vn } be a basis of V . • We then construct the vectors e1 , e2 , . . . en of an orthogonal basis in V succesively using: • e1 = v1 . • e2 = v2 − hv2 , e1 i e1 . ke1 k2 • e3 = v3 − hv3 , e1 i hv3 , e2 i e1 − e2 . 2 ke1 k ke2 k2 • ... • en = vn − hvn , e1 i hvn , en−1 i e1 − . . . − en−1 . 2 ke1 k ken−1 k2 We note that to get a orthonormal basis (and orthogonal matrix) we need to also normalize the resulting vectors. We can do that either at the end when we have found an orthogonal basis or during the algorithm itself. While this works in theory, when implemented on a computer there can be some stability problems. When computing the vector ek , there might be some small error in precision resulting in vk being not completely orthogonal. This is then compounded when we further orthogonalize against this vector and resulting vectors repeatedly. Luckily we can make a small adjustment of the method in order to handle these problems. Modified Gram-Schmidt works like Gram-Schmidt with a small modification. • u1 = v 1 . (1) • uk = v k − 72 hvk , u1 i u1 . ku1 k2 (1) (2) (1) • uk = uk − huk , u2 i u2 . ku2 k2 • ... (n−2) • (n−1) uk = (n−2) uk hu , un−1 i − k un−1 . kun−1 k2 (j) • In every step j we orthogonalize all vectors uk , k > j with uj . We note that while we assume that {v1 , v2 , . . . vn } it doesn’t actually need to be linearly independent. Since if we apply it to a set of linear dependent vectors then for some i, ui = 0 and we can discard them and not orthogonalize further vectors with them. In this case the number of resulting vectors will be the same as the dimension of the space spanned by the original vectors. One example of where we can use Gram-Schmidt is in finding a QR-decomposition of a matrix, we look more at that in the applications part here and the QR-method later. While we have already used projections, we haven’t actually defined what a projection is in respect to a inner product space. Definition 6.14 If v = v1 + v2 with v1 ∈ U and v2 ∈ U ⊥ , where U ⊥ is orthogonal to U . • The vector v1 is called the projection of v on U , denoted: projU (v). • The vector v2 is called the component of v orthogonal to U , given by: v2 = v − projU (v). Theorem 6.11. If V is a inner product space and U is a subspace of V . Then every vector v ∈ V can be written uniquely as: v = projU (v) + (v − projU (v)) • If e1 , e2 , . . . , em is any orthogonal basis of U , then: projU (v) = hv, e1 i hv, e2 i hv, em i e1 + e2 + . . . + em ke1 k2 ke2 k2 kem k2 • Additionally, for the subspaces U, U ⊥ we have: dim V = dim U + dimU ⊥ The projection of a vector v in V on a subspace U can be seen as the vector in U most similar to v in U . Looking at the Gram-Schmidt process earlier we also see that another way to see the algorithm is by adding new vectors ei by taking the reminder after subtracting the projection of ei on the subspace spanned by the previous vectors. ei = vi − projU (v), U = span(e1 , . . . , ei−1 ) 73 Theorem 6.12. If V is a inner product space and U is a subspace of V . Then projU (v) is the vector in U closest to v ∈ V . kv − projU (v)k ≤ kv − uk for all u ∈ U . Proof. For v − u we get: v − u = (v − projU (v)) + (projU (v) − u) Since the first lies in U ⊥ and the second lies in U they are orthogonal and we can use pythagoras: kv − uk2 = (v − projU (v))2 + (projU (v) − u)2 ≥ (v − projU (v))2 We look at an example in R3 where we want to find the point on a plane closest to a point in space. We consider the plane given by the equation x − y − 2z = 0 and find the point on the plane closest to the point v = (3, 1, −2). • For the plane we get x = y + 2z and the subspace U = {s + 2t, s, t|s, t ∈ R} = span((1, 1, 0), (2, 0, 1)) • Using Gram-Schmidt we get: {e1 = (1, 1, 0), e2 = (2, 0, 1)− h(2, 0, 1), (1, 1, 0)i (1, 1, 0). k(1, 1, 0)k2 • e2 = (2, 0, 1) − (1, 1, 0) = (1, −1, 1). • Then we find the closest point by finding the projection of v on U : projU (v) = hv, e1 i hv, e2 i e1 + e2 ke1 k2 ke2 k2 • = 2e1 + 0e2 = (2, 2, 0) Since we have found an orthogonal base for the plane we can now easily also find the point on the plane closest to any general point v = (x, y, z)in space as: projU (v) = hv, e1 i hv, e2 i x+y x−y+z e1 + e2 = (1, 1, 0) + (1, −1, 1) 2 2 ke1 k ke2 k 2 3 6.4 Applications 6.4.1 Fourier approximation We inspect the space C[−π, π] of continuous functions on the interval [−π, π] using the inner product: Z π hf, gi = f (x)g(x) dx π . 74 • Then {1, sin x, cos x, sin 2x, cos 2x, . . .} is an orthogonal set. • We can then approximate a function f by projecting it on a finite subspace of this set: Tn = span({1, sin x, cos x, sin 2x, cos 2x, . . . , sin mx, cos mx}) This gives f (x) ≈ tnf (x) = a0 +a1 cos x+b1 sin x+a2 cos 2xb2 sin 2x+. . .+am cos mx+ bm sin mx. • The coefficients ak , bk we get as: Z π hf (x), 1i 1 • a0 = = f (x) dx. k1k2 2π π Z hf (x), cos(kx)i 1 π • ak = f (x) cos(kx) dx. = k cos(kx)k2 π π Z hf (x), sin(kx)i 1 π • ak = = f (x) sin(kx) dx. k sin(kx)k2 π π • It’s obvious that as we take more coefficients we get a better approximation of f since T1 ⊆ T2 ⊆ T3 ⊆ . . .. • If we take an infinite number of coefficients we call this the Fourier series of f . 6.4.2 Finding a QR decomposition using the Gram-Schmidt process We remember that given a matrix M we can decompose it in two matrices M = QR such that Q is a unitary matrix and R is a upper triangular matrix. The QR-decomposition can then in turn be used to solve or find a least square solution to the linear system defined by M . It can also be used to more effectively find the eigenvalues and eigenvectors of M . Both of these uses of the QR decomposition is handled elsewhere in the book however one problem first need to be solved. If we want to use QR decomposition to speed up calculations, we need to be able to find the QR decomposition itself as fast (and stable) as possible. If we look at the columns of A as a (possibly linearly dependent) set of basis vectors. Then QR can be seen as a base change Q with coefficients in R. If we look at the Gram-Schmidt process used on the column vectors v of A we get (before normalization): ei = vi − hvi , e1 i hvi , ei−1 i e1 − . . . − ei−1 2 ke1 k kei−1 k2 ⇔ vi = ei + hvi , e1 i hvi , ei−1 i e1 + . . . + ei−1 2 ke1 k kei−1 k2 So we write column i using i basis vectors. If we collect the coefficient in R and the basis vectors in Q we get a QR decomposition. You can easily check this by considering 75 one column in A as we did above. In practice we usually first calculate Q and then find R using multiplication R = QH A since Q is unitary. We give an example considering the matrix A Schmidt process: 1 −1 2 A = 2 0 2 2 2 3 below using the unmodified Gram 3 1 2 This gives : 1 −4 A.1 1 1 2 , A0.2 = A.2 − QH −2 Q.1 = = A Q = .2 .1 .1 kA.1 k 3 3 2 4 2 6 1 1 0 0 −2 , A.4 = −3 A.3 = 3 3 1 0 Continuing with the next column we get: : −2 2 0 1 1 A.2 00 0 H 0 −1 , A.3 = A.3 − Q.2 A.3 Q.2 = −2 = Q.2 = kA0.2 k 3 3 2 1 4 1 A00.4 = −4 3 2 2 A00.3 1 −2 Q.3 = = kA00.3 k 3 1 We now have a full basis and can stop since Q.4 must be zero. For the QR-decomposition A = QR we then get the matrices: 1 −2 2 1 Q = (Q.1 , Q.2 , Q.3 ) = 2 −1 −2 3 2 2 1 And 9 3 12 9 1 R = QH A = 0 6 0 −3 3 0 0 3 6 Here are a couple of things to consider and think about: • Can you recognize the Gram-Schmidt process when written in this form? Where are we doing the normalization? 76 • Can you find the QR-decomposition using the modified algorithm? • Using that Q is hermitian we can write a linear system Ax = b as QRx = b ⇒ Rx = QH b which is very easy to solve since R is triangular. How does this compare to other standard methods such as Gauss elimination? 6.4.3 Principal component analysis (PCA) and dimension reduction Let us consider n sets of mesurements at time t = 1, 2, 3, . . . , m. For example the value of n stocks, consumption of n commoditys, rainfall at n locations, etc. While making our analysis of the data might be easy if n is small, as n grows the time to make the analysis typically grows faster (for example in n2 ). As n grows large enough we reach a point where we cannot do our analysis in a reasonable time anymore. However let us assume the measurements are dependent, this seems like a rather reasonable assumption. Our goal is to find a new basis in which we can represent the n sets using fewer basis vectors. If one measurement is a linear combination of the others we can use for example Gram-Schmidts to find a new basis. However in practice this is rarely the case, especially considering as in this case typical stochastic processes. Since Gram-Schmidt is guaranteed to find a basis with the minimum number of basis vectors since it gives an orthogonal basis, findign a basis with less basis vectors that captures all the data is impossible. However what if we can find a new basis which captures most of the data? For example we might have points in 3 dimensions which all lie on or close to a plane. In this case projecting the points to this plane should give a good approximation of the data, as well as reducing the dimension of the problem. One method used to reduce the dimension of a problem is principal component analysis (PCA): Our first step is to subtract the mean of every individual data set. We use this to construct a data matrix by inserting the different sets of zero mean measurements as row vectors Xn×m with mean zero of every row. Next we compute the singular value decomposition (SVD): X = WΣVT We remember that W is a unitary matrix of the eigenvectors of XXT , Σ is a diagonal matrix with the singular values on the diagonal and V is a unitary matrix which contains the eigenvectors of XT X. The singular values Σi,i are the square roots of the eigenvalues of XXT . The eigenvalues are proportional to the variance of the data captured by corresponding eigenvector. Essentially we can see the first eigenvector corresponding to the largest eigenvalue as the direction in which our original data have the largest variance. The second eigenvector corresponds to the direction orthogonal to the first which which captures the most of the remaining variance, the third is orthogonal to the first two and 77 so on. If we would project our data onto the space projected by the first L eigenvectors we would thus get the L dimensional space that captures the most variance in the data. If the last principal components are small they can be disregarded without a large loss of information. If we let Yn×m be the projection on the space described by the eigenvectors of W we get: Y = WT X = WT WΣVT = ΣVT If we only want the space described by the first L eigenvectors we get: T YL = WL X = ΣL VT T Where WL is a L × n matrix taken as the first L rows of W, and ΣL is a L × m matrix composed of the first L rows of Σ. We can now work with our L × m matrix YL rather than X with little loss of information if the remaining singular values are small. We note that it is also possible to work with the covariance matrix C = XXT since W is the eigenvectors of XXT And as such, if we have the covariance matrix XXT we can calculate it’s eigendecomposition rather than the singular value decomposition of X. However it is recommended to use SVD decomposition since constructing the covariance matrix XXT can give a loss of precision using the covariance method. While PCA gives the ”best” orthogonal basis, there are other methods as well as other variations of PCA to reduce the dimension of a dataset. Some things to think about concerning the use of PCA is: • PCA assumes tha data are scaled equally. If the data is of temperatures for example, what happens if one dataset is in Fahrenheit while all others are in Celsius? Similar problems arise whenever one dataset have much larger/smaller variance than the others. • How does outliers (datapoints far from the others often the result of error in measurements) effect the results of PCA? Can you think of any way to handle them better? • What happens if the datasets are independent or close to independent? Why can’t we use PCA now? Is there something else you could possible do in this case? • Can you think of any examples in your main field where PCA or similar methods could prove useful? 78 7 Linear transformations A linear transformations also called a linear map is a function mapping elements from one set to another (possibly same) set fulfilling certain conditions. Namely if the two sets are linear spaces then the map preserves vector addition and scalar multiplication. In this chapter we will take a look at linear transformations, what properties makes them important and some examples of where we can use them in applications. This as with the chapter before this will have a heavier focus on theory than other chapters. The definition of a linear transformation is as follows: Definition 7.1 Given two vector spaces V, W over the same field F , then a function T : V → W is a linear transformation if: • T (v1 + v2 ) = T (v1 ) + T (v2 ), v1 , v2 ∈ V . • T (rv) = rT (v), v ∈ V, r ∈ F . We see that the first condition is equivalent to saying that T needs to preserve vector addition and the second that it preserves scalar multiplication. We note that the two vector spaces V, W need not be the same vector space, or even have the same number of elements. This also means that the addition operation need not be the same since v1 + v2 is done in V and T (v1 ) + T (v2 ) is done in W . From the definition follows immediately a couple of properties, it’s recommended to stop and think about why if it isn’t immediately clear. • T (0) = 0 • T (−v) = −T (v), for all v ∈ V . • T (r1 v1 + r2 v2 + . . . rn vn ) = T (r1 v1 ) + T (r2 v2 ) + . . . + T (rn vn ) for all v ∈ V and all r ∈ F . Many important transformations are in fact linear, some important examples are: • The identity map x → x and the zero map x → 0. • Rotation, translation and projection in geometry. • Differentiation. • Change of one basis to another. • The expected value of a random variable. The importance of many linear transformations gives a motivation why we should study them further and try to find properties of them that we can use. In fact, just by knowing a transformation is linear we already know two useful properties from the 79 definition. For example it doesn’t matter if we first add two vectors and then change their basis or if we first change both their basis and then add them. Often it’s possible to check if a transformation is linear using only the definition as in the example below. We define a function T : R2 → R3 such that: x+y x T = 2x − y y y We then want to show that it’s a linear transformation. (x1 + y1 ) + (x2 + y2 ) x1 x T + 2 = (2x1 − y1 ) + (2x2 − y2 ) y1 y2 y1 + y2 (x1 + y1 ) (x2 + y2 ) x1 x = (2x1 − y1 ) + (2x2 − y2 ) = T +T 2 y1 y2 y1 y2 And the first requirement is true. For the second requirement T (rv) = rT (v) we get: (rx + ry) rx x = (r2x − ry) =T T r ry y ry (x + y) x = r (2x − y) = rT y y And the second requirement is true as well. Next we look at how to linear transformations act on basis vectors and how this can be used in order to create a linear transformation. Theorem 7.1. We let V, W be vector spaces over the field F . Given a basis e1 , e2 , . . . , en of V and any set of vectors {w1 , w2 , . . . , wn } ∈ W . There exists a unique linear transformation T : V → W such that T (ei ) = wi . • Additionally we have: T (v1 e1 + v2 e2 + . . . + vn en ) = v1 w1 + v2 w2 + . . . vn wn Proof. We start by looking at the uniqueness. If two transformations S, T exist such that T (ei ) = S(ei ) = wi . Then they must be the same since e1 , e2 , . . . , en spans V . You can easily see that they must be the same 80 transformation by writing any S(v) = T (v) and expanding v in the basis vectors. And we can conclude that the transformation is unique if it exists. Next we want to show that such a linear transformation always exists. Given v = v1 e1 + v2 e2 + . . . + vn en ∈ V , where v1 , v2 , . . . , vn are uniquely determined since e1 , e2 , . . . , en is a basis of V . We define a transformation T : V → W as: T (v) = T (v1 e1 + v2 e2 + . . . + vn en ) = v1 w1 + v2 w2 + . . . + vn wn For all v ∈ V . We clearly have that T (ei ) = wi and showing that T is linear can be done simply by trying: T (av1 + bv2 ) = aT (v1 ) + bT (v2 ) Using this we see that it’s quite easy to define a linear transformation, since we only need to specify where to take the basis vectors. Then everything else is handled by the linearity. This result also gives an easy way to find if two linear transformations are equal. Simply check if they give the same result when applied to the basis vectors. If they do, they need to be equal because of the uniqueness. We continue by looking at matrices multiplied with vectors acting as transformations. Consider the m × n matrix A and define a function TA : Rn → Rm such that TA (v) = Av for all (column vectors) v ∈ Rn . Show that TA is a linear transformation. We use the definition and immediately get: TA (v1 + v2 ) = A(v1 + v2 ) = Av1 + Av2 = TA (v1 ) + TA (v2 ) TA (rv) = A(rv) = r(Av) = rTA (v) And TA must be a linear transformation. We see that every m × n matrix A defines a linear transformation TA : Rn → Rm . This gives us a new perspective in which to view matrices, namely we can view them as linear transformations on vector spaces. But even more importantly we will also see that every linear transformation T : Rn → Rm can be written in this way using a m × n matrix. Theorem 7.2. If T : Rn → Rm is a linear transformation, then there exists a m×n matrix A such that T (v) = Av for all column vectors v ∈ Rn . Additionally if e1 , e2 , . . . , en is the standard basis in Rn , then A can be written as: A = [T (e1 ) T (e2 ) . . . T (en )] Where T (ei ) are the columns of A. The matrix A is called the standard matrix of the linear transformation T . 81 7.1 Linear transformations in geometry The matrix representation of linear transformations is especially useful in geometry and 3d graphics. Especially since most commonly used transformations such as rotation and scaling are in fact linear. If we can find the standard matrix for a certain transformation we can then easily apply that transformation on a large number of vectors with a simple matrix-vector multiplication. This is especially important in for example computer graphics where we often work with a huge number of points even for very primitive objects, which need to be moved whenever the object or camera (view point) moves. Here we will give some examples of linear transformations in geometry and show a couple of methods which we can use to find the standard matrix. Using the previous result we find the standard matrix for a couple of common transformaitons in Rn → Rm . We start in R2 → R2 by finding the standard matrix for rotation clockwise about the origin. By making the unit circle we can easily see that if we rotate the first basis vector [1 0]> by θ degrees we end up at [cos θ sin θ]> And if we rotate the second one [0 1]> we end up at [− sin θ cos θ]> Our standard matrix is thus cos θ − sin θ A= sin θ cos θ And we can use this to rotate any point θ degrees about the origin. For example remembering that the angle between two orthogonal vectors is θ = 90◦ , given a vector v = [x y]> we can find an orthogonal vector in 2d by multiplying with: A= cos 90 − sin 90 0 −1 = sin 90 cos 90 1 0 Giving Av = [−y x]> , which can easily seen to give zero when taking the dot product together with v. But what if we now instead consider this rotation in the xy−plane in R3 ? When looking at rotation around the z axis we get the matrix: cos θ − sin θ 0 A = sin θ cos θ 0 0 0 1 Since for the first two basis vectors we rotate in the plane we had in two dimensions, and the last basis-vector in the direction of the z−axis should remain unchanged. We can find the rotation matrices around the other axes in the same way. 82 Next we take a look at projection in R2 → R2 where we want to find the standard matrix for projection on the line y = ax. We let d = [1 a]> be the direction vector for the line. Then using the projection formula with the basis vectors we get: projd (1, 0) = 1 h(1, 0), (1, a)i d= (1, a) 2 k(1, a)k 1 + a2 projd (0, 1) = h(0, 1), (1, a)i a d= (1, a) 2 k(1, a)k 1 + a2 Which gives the standard matrix 1 1 a A= 1 + a2 a a2 But what if we instead want to find the standard matrix for projection on the plane described by x − z = 0 in R3 ? We introduce a slightly different method more suitable for higher dimensions. We have an orthogonal basis U = {(1, 0, 1), (0, 1, 0)} for points in the plane. Then instead of finding the projection of each basis vector on the plane individually, we use the projection formula to find the projection of a general point onto this plane. Doing this yields: projU (x, y, z) = h(x, y, z), (1, 0, 1)i (1, 0, 1) k(1, 0, 1)k2 h(x, y, z), (0, 1, 0)i + (0, 1, 0) = k(0, 1, 0)k2 We then split ut up in the part belonging to basis vectors: 1/2 A= 0 1/2 x+z x+z , y, 2 2 x, y, z to get the columns for the individual 0 1/2 1 0 0 1/2 While we probably could have found the standard matrix for this projection regardless, we can use the same methodology for more complex subspaces in higher dimensions. • First find an orthogonal basis for the subspace. • Then find the projection of the basis vectors in the original space on the subspace. We then get a matrix for this transformation such that we can calculate this projection of any points onto the subspace with a simple matrix-vector multiplication. Before continuing with some more theory we return to R2 and look at the reflection in the line y = ax • The reflection of v is equal to 2projd (v) − v (make a picture). 83 This gives: 2 1 a 1 0 − A= 0 1 1 + a2 a a2 1 1 − a2 2a = 2a a2 − 1 1 + a2 We note that the methodology presented here is possibly since the transformations in question are all linear. While most needed transformations in this case are linear, we will see later that there are at least one critical non-linear transformation we need to be able to do in many applications. Can you think about which one? 7.2 surjection,injection and combined transformations While nice that we can for example rotate a point, what if we want to combine transformations. For example maybe we want to rotate around a point which itself rotate around another point. We will also give a way in which to make certain non-linear transformations into linear ones, a method essential to many applications in computer graphics. We start by defining what we mean by a surjection or injection. Definition 7.2 A linear transformation (or more generally a function) F : V → U is a: • Surjection, or onto if every element w ∈ W can be written w = T (v at least one vector v ∈ V . • Injection, or one-to-one function if every element in U is mapped to by at most one element in V . • If it’s both surjective and injective it’s called a bijection. One important property is that for F to be invertible it needs to be an injection since otherwise F −1 (u), u ∈ U would map to at least two elements in V for some u. Theorem 7.3. Let F : V → W and G : W → U be two linear transformations. Then the combined transformation G ◦ F : V → U is linear. If F is invertible, then the inverse F −1 : W → V is linear. Proof. For the first part we get: G ◦ F (au + bv) = G(F (au + bv)) = G(F (au)) + G(F (bv)) = aG(F (u)) + bG(F (v)) For the second part we get: F (F −1 (au + bv)) = au + bv = aF (F −1 (u)) + bF (F −1 (v)) = F (aF −1 (u) + bF −1 (v)) 84 Since F is linear and injective. From this follows that: F −1 (au + bv) = aF −1 (u) + bF −1 (v) But if the composition of two linear transformations is itself a linear transformation, how do we represent the combined linear transformation G ◦ F ? If we look at the matrix representations of the transformations we get: AF = [F (e1 ) F (e2 ) . . . F (en )] AG = [G(e1 ) G(e2 ) . . . G(en )] Then: AG◦F = [G(F (e1 )) G(F (e2 )) . . . G(F (en ))] = AG AF In other words if we can combine two linear transformations simply by multiplying their respective standard matrices. We note that since matrix multiplication is not commutative, the order in which we apply the linear transformations does impact the resulting transformation. The transformation we want to apply ”first” should be to the right in the matrix multiplication. 7.3 Application: Transformations in computer graphics and homogeneous coordinates We want to make a simple computer graphics of our solar system containing the sun, planets as well as their moons. For simplicity we assume all bodies are spherical and all bodies move in circular paths around something else. To work with we have the points of a unit sphere located in origo. This is a good example of how we can combine relatively simple transformations in order to create very complex paths such as the movement of someone on the Moons surface relative the sun. We start by assuming the sun is stationary positioned in the middle of the solar system (origo) and everything else revolves around it. Then we only need to scale our unit sphere to appropriate size as well as rotate all the points on the sphere around some axis. Since both rotation and scaling are linear transformations we can apply this transformation without further problems. We let t be a point in time, Rsun be the rotation of the sun around it’s own axis and Ssun be the standard matrix for the scaling needed to get the appropriate size, we then get standard matrix for the points on the unit sphere representing the sun as. Asun = Rsun Ssun Lets say we choose to let the sun rotate around the z−axis, we then get the rotation matrix as before and we get transformation matrix in time t as: cos θs t − sin θs t 0 ssun 0 0 ssun 0 Asun = sin θs t cos θs t 0 0 0 0 1 0 0 ssun 85 ssun cos θs t −ssun sin θs t 0 0 = ssun sin θs t ssun cos θs t 0 0 ssun We are now finished with the sun, in every frame we multiply all the points on the unit sphere with corresponding matrix to get their position in our animation. However as soon as we want to animate our first planet we immediately move into problems, moving a point a fixed amount in any direction (translation) is not a linear translation. You can easily realize this if you let T be translation in R: T (x) = x + a → T (0 + 0) = a 6= T (0) + T (0) = 2a. The same is obviously true for higher dimensions. While we could start by adding a vector to every point and then rotate around origo, we need to take great care in order to get the rotation around it’s own axis right(rotate points, then add vector, then finally rotate around origo) as well as it now requiring 2 matrix multiplications and one vector addition. Additionally as soon as we want to add a moon or want to make more complex objects, the problem becomes more and more complicated. Luckily there is a way to handle these transformations, a trick you could say in which we represent our n dimensional vector in n + 1 dimensions. We represent a vector v = (v1 , v2 , . . . , vn )in n dimensions as a vector v = (v1 , v2 , . . . , vn , 1) in n + 1 dimensions. For two dimensions this can be seen as representing the points in the xy−plane by the points in the plane 1 unit above the xy−plane in three dimensions. Now we can represent translation x0 = x + tx , y 0 = y + ty , z 0 = z + tz (in three dimensions) using matrix multiplication: 0 1 0 0 tx x x y 0 0 1 0 ty y 0 = z 0 0 1 tz z 1 1 0 0 0 1 This is commonly referred to as homogeneous coordinates. We note that the last element of every vector is always 1 and the last row of the standard matrix is always the same as well. So by using one extra dimension even if the value of every vector in this extra dimension is always one, we can make translation into a linear transformation. • We note that origo is no longer represented by the zero vector (0, 0, 0) in homogeneous coordinates since it’s there represented by (0, 0, 0, 1). • It’s also easy to write any other previous already linear transformation in homogeneous coordinates, we simply add another row and column with all zeros except for the bottom right corner which we set to one. We are now ready to return to our animation of the solar system using homogeneous coordinates instead. • We start by rewriting the transformation matrix for the points on the suns surface simply by adding the additional row and column: 86 Asun ssun cos θt −ssun sin θt 0 ssun sin θt ssun cos θt 0 = 0 0 ssun 0 0 0 0 0 0 1 Next we want to find the transformation matrix for the points on the earth’s surface from another unit sphere in the origin. A good guideline is to start by applying all transformations of the object in reference to itself, and then apply the transformation with respect to the ”closest” reference point, in this case the position of the sun. • We start by applying the scaling SE and rotation around it’s own axis RE just like we did with the sun. • Then we need to transpose the planet away from the sun in the origin. We apply the translation matrix TES so that it’s appropriately far away. • We then apply another rotation on this to rotate all of these points around the sun (origin) RES representing the rotation of earth around the sun. • Applying all of the transformations in order gives transformation matrix AE = RES TES RE SE . Especially we note that the first two RES TES denoting the centre position of the earth relative to the sun and the last two REE SE is the position of the points on the earths surface relative a unit sphere positioned at earth’s core. We can use the first (already calculated) half of the transformation when finding the transformation matrix for the moon. In the same way we get the transformation SM , RM to apply the scaling and rotation around the moons own axis. We consider the position of the earth as the origin and translate and rotate around that giving: TM E , RM E . We notice that we now have the position of the moon relative the centre of earth, we apply the transformations ARES TES previously applied on the earth’s centre to get the moons coordinates relative the origin. This gives transformation matrix: AM = RES TES RM E TM E RM SM We notice that we can interpret this as: (from right to left). • RM SM -Position of points on the moons surface relative a unit sphere at the centre of the moon, times: • RM E TM E -Position of the centre of the moon relative the centre of earth, times: • RES TES -Position of the centre of earth relative the centre of the sun (origin). This methodology is used extensively in computer graphics and robotics in order to describe very complex movements using relatively simple transformations and hierarchies. 87 7.4 Kernel and image To every linear transformation T : V → W we associate two linear spaces. One of them a subspace of V and the other a subspace of W . Definition 7.3 Given the linear transformation T . • R(T ) = {w|w = T (v), v ∈ V } ⊂ W , is called the image of T . • In other words the subspace of W which is mapped to from V • N (T ) = {v|T (v) = 0} ⊂ V . is called the Kernel of T or the Nullspace of T . • In other words the subspace of V for which the linear transformation maps to the zero-vector in W If we look at the transformation in terms of the standard matrix we see that the image R(T ) is the solutions to the matrix equation Ax = b. While the Kernel N (T )l is the solutions to Ax = 0. Theorem 7.4. Given the linear transformation T : V → W , the kernel N (T ) is a subspace of V and the image R(T ) is a subspace of W . Proof. Since T (0) = 0 they both contain the zero vector. For the kernel we get: T (av1 + bv2 ) = aT (v1 ) + bT (v2 ) = a0 + b0 = 0, v1 , v2 ∈ N (T ) And av1 + bv2 lies in N (T ) since it satisfies the condition. So N (T ) is a subspace of V . For the image we use: w1 = T (v1 ), w2 = T (v2 ) v1 , v2 ∈ V This gives: aw1 + bw2 = aT (v1 ) + bT (v2 ) = T (av1 ) + T (bv2 ) = T (av1 + bv2 ) And aw1 + bw2 lies in R(T ). So R(T ) is a subspace of W . We can use this to easily check if a linear transformation is injective,surjective and bijective (and invertible). Theorem 7.5. A linear transformation T : V → W is: • surjective ⇔ R(T ) = W • injective ⇔ N (T ) = {0} • bijective (and invertible) ⇔ N (T ) = {0} and R(T ) = W . 88 If the two vector spaces have the same dimension (square standard matrix) it’s actually enough to show only that it’s injective or surjective and it will always be bijective as well. Theorem 7.6. Assume a linear transformation T : V → W , where dim V = dim W . • If N (T ) = {0} or R(T ) = W then T is bijective. In fact we can show even more, we can show something about the combined dimension of the kernel and image: Theorem 7.7. Assume we have a linear transformation T : V → W , where V is finite dimensional, then: • dim V = dim N (T ) + dim R(T ) We move on to a short example showing how we can use what we have learned here. Hermite interpolation is a method in which to find a curve through a given number of points. This is often needed in computer graphics as well as in construction for calculating strength of materials. We will look at one subproblem where we want to fit a curve through two points in the plane where the slope is specified in both points. The easiest possible interpolation function we can think of in this case is a third order polynomial: y = a + bx + cx2 + dx3 Which we want to fit to the points such that the curve goes through the points with the given slope. In order to do this we get a system of four unknowns and four equations: y(x1 ) = a + bx1 + cx21 + dx31 = y1 y 0 (x1 ) = b + 2cx1 + 3dx21 = s1 y(x2 ) = a + bx2 + cx22 + dx32 = y2 y 0 (x2 ) = b + 2cx2 + 3dx22 = s2 Writing it in matrix form we get a linear system in the four unknowns (a, b, c, d) with matrix A: 1 x1 x21 x31 0 1 2x1 3x21 A= 1 x2 x22 x32 0 1 2x2 3x22 While we know how to solve this problem, we ask ourselves, is it always possible to solve this system? Obviously we could simply solve it and see, but doing so (especially by hand) is rather tedious work even for this relatively small matrix. • You could also calculate the determinant, less work than solving the system but still alot of work. 89 • Instead we look at the transformation a + bx1 + cx21 + dx31 y(x1 ) y 0 (x1 ) b + 2cx1 + 3dx21 T : y → T (y) = y(x2 ) = a + bx2 + cx22 + dx32 y 0 (x2 ) b + 2cx2 + 3dx22 As the transformation from V = {y|y(x) ∈ P3 } to W = R4 . Using the basis {1, x, x2 , x3 } in V this is obviously a linear transformation with matrix A above. We take a look at the Kernel of A. For a polynomial in N (T ) we have y(x1 ) = y(x2 ) = 0, y 0 (x1 ) = y 0 (x2 ) = 0 So it has a root in x1 , x2 both with multiplicity two. However y is a third order polynomial and can only have 3 roots (counting multiplicity) except for the zero polynomial. So N (T ) = {0} and T is injective. Since it’s injective and dim P3 = dim R4 T is bijective and invertible so the system have a unique solution. We remember the definition of rank for a square matrix: Definition 7.4 For A ∈ Mn×n (F) the rank(A) is equivalent to: • The number of linearly independent rows or columns. • The dimension of the column space of A. But we can give even more definitions of rank, we now move on to define the rank of a matrix Mm×n (F) using the kernel and image of A. Definition 7.5 Given a linear transformation T : V → W with standard matrix A. • Rank or the column rank r(A) is the dimension of the image of A: dim R(A), with elements Ax, x ∈ F n . • Row rank is the dimension of the image of A> : dim R(A> ), with elements A> y, y ∈ F m. • Nullity or the right nullity v(A) is the dimension of the kernel of A: dim N (A), with elements Ax = 0, x ∈ F n . • Left nullity is the dimension of the kernel of A> : dim N (A> ), with elements Ay = 0, y ∈ F m . From this immediately follows: • Nullity zero: v(A) = 0 ⇔ A have a left inverse. 90 • rank m: r(A) = m ⇔ A have a right inverse. We can of course do the same for A> , however we will see that they are related in such a way that knowing only one of these four values (R(A), R(A> ), N (A), N (A> )), we can immediately find the others. Theorem 7.8. For every matrix A we have: dim R(A) = dim R(A> ) In other words, for non square as with the square matrices the column rank and row rank are the same. As with a square matrix you can show it by writing the matrix in echelon form. We end this part with a short example: Let A be a n × m matrix with rank r. We want to show that the space of all solutions to Ax = 0 of n homogeneous equations in m variables has dimension m − r. This space is the kernel of N (T ) with T : Rn → Rm such that T (v) = Av. We use dim V = dim N (T ) + dim R(T ) We already have the dimension of the image R(T ), which is the rank r. The dimension of V is m And we immediately get the dimension of the kernel as dim N (T ) = m − r 7.5 Isomorphisms Definition 7.6 A linear transformation T : V → W is called an Isomorphism if it’s both surjective and injective (bijective). • Two vector spaces V and W are called isomorphic if there exists an isomorphism T : V → W. • We denote this by writing V ∼ = W. We see that for two vector spaces that are isomorphic, we have a pairing v → T (v), v ∈ V, T (v) ∈ W . For this pairing we preserve the vector addition and scalar multiplication in respective space. Because addition and multiplication in one space is completely determined by the same operation in the other all vector space properties of one can be determined by those of the other. In essence you could say the two vector spaces are identical apart from notation. Theorem 7.9. For two finite dimensional spaces V, W and a linear transformation T : V → W the following are equivalent: • T is an isomorphism. • If {e1 , e2 , . . . , en } is a basis of V then {T (e1 ), T (e2 ), . . . , T (en )} is a basis of W . 91 • There exists a basis {e1 , e2 , . . . , en } of V such that {T (e1 ), T (e2 ), . . . , T (en )} is a basis of W . • There exist a inverse linear transformation T −1 : W → V such that T −1 T = 1V and T T −1 = 1W Here we see one of the important relations between isomorphic vector spaces, if we have a basis in one we can given a linear isomorphic transformation also find a basis in the other. If both vector spaces have a finite dimension we can say even more. Theorem 7.10. For the finite dimensional spaces V, W the following are equivalent. • V ∼ = W , (V and W isomorphic). • dim V = dim W In other words, the dimension of a vector space decides which spaces it’s isomorphic to. For example we can immediately conclude that the linear space R3 is isomorphic to the polynomials of order two: P2 (Polynomials that can be written on the form a + bx + cx2 ). Some examples of important isomorphisms are: • The Laplace transform, mapping (hard) differential equations into algebraic equations. • Graph isomorphisms T : V → W such that if there is an edge between edges u, v in V there is an edge between T (u), T (v) in W . • Others include isomorphisms on groups and rings (such as fields). This also give a little motivation of why the study of isomorphisms are useful, if we have a problematic vector space where something we want to do is very complicated, if we can find another isomorphic vector space where the corresponding operation is easier. We can formulate our problem in this new vector space instead, do our now easier calculations and then change back, effectively sidetracking the whole issue. 92 8 Least Square Method (LSM) Many of the properties and ideas we have examined apply to square matrices, a good example of this is theorem 1.1 which list nine conditions that are all equivalent to a square matrix being invertible. If a square matrix A is invertible then the linear equation system Ax = y is solvable. In this section we will discuss a scenario where we have a system of linear equations where the A-matrix has more rows than columns. This kind of system is called an overdetermined system (we have more equations than unknowns) and if A is a m × n matrix with m > n it is said to have m − n redundant equations. An overdetermined system can be solvable if the elements of vector y represents points on a straight line. Having 7 7 6 6 5 5 4 4 3 3 2 2 1 1 0 0 2 4 6 8 (a) A set of points on a single line. 10 0 0 2 4 6 8 10 (b) A set of points near but not on a single line. a number of points near but not on a straight line is a very common situation. One of the simplest ways to model a relation between two quantities is to use a straight line, l(x) = a + bx, and often experiments are designed such that this kind of relation is expected or the quantities measured are chosen such that the relation should appear. In reality any system, be it the temperature of a piece of machinery, the amplitude of a signal or the price of a bond, will be subject to some disturbance, often referred to as noise, that will move any move the measured values of the ideal line. In this situation the two most common problems are: • Hypothesis testing: does this data fit our theory? • Parameter estimation: how well can we fit this model to our data? In both of these cases it is important to somehow define the ’best-fitting’ line. Usually this is done in what is called the least-square sense. Our problem can be described in 93 the following way: Ac = y y1 1 x1 y2 1 x2 a , y= . A = . . , c = . . b .. . . yn 1 xn Where x = [x1 x2 . . . xn ]> and y is the data and c is the vector of coefficients that we want to find. To define what we mean by least-square sense we first note that for each x there is a corresponding spot on the line, a + bx, and a corresponding measurement y. We can create a vector that describes these distances for each measurement called the residual vector or the error vector . e = Ax − y, is the residual vector We can then measure the normal (Euclidean) length of this vector. v u n √ uX |e| = t |ei |2 = e> e i=1 Minimizing this length gives us the best-fitting line in the least-square sense. This can also easily be expanded into other kinds of functions, not just straight lines. Suppose we want to fit a parabola (second degree polynomial instead, p(x) = a+bx+cx2 ) then the problem can be written: Ac = y 1 x1 x21 y1 a 1 x2 x2 y2 2 A = . . .. , c = b , y = .. .. .. . . c 2 yn 1 xn xn Here the residual vector can easily be defined in the same way as for the straight line. This can be extended in the same way for a polynomial of any degree smaller than n − 1 where n is the number of measurements. If a polynomial of degree n − 1 is chosen we get the following situation: 1 1 A = . .. x1 x2 .. . 1 xm 94 x21 x22 .. . x2m Ac = y . . . xm−1 c1 y1 1 m−1 . . . x2 c2 y2 , c = , y = . . .. .. .. .. . . yn cn . . . xm−1 m Here the matrix A will be square and the system will be solvable as long as all the x-values are different. This means that for every set of n measurements there is a polynomial of degree n − 1 that passes through each point in the measurement. This is often referred to as interpolation of the data set (and in this case we it is equivalent to Lagrange’s method of interpolation). The matrix A here is known as a Vandermonde matrix. It and similar matrices appear in many different applications and it is also known for having a simple formula for calculating its determinant. 8.1 Finding the coefficients How can we actually find the coefficients that minimize the length of the residual vector? First we can define the square of the length of the residual vector as a function > s(e) = e e = n X |ei |2 = (Ax − y)> (Ax − y) i=1 This kind of function is a positive second degree polynomial with no mixed terms and ∂s thus has a global minima where = 0 for all 1 ≤ i ≤ n. We can find the global ∂ei minima by looking at the derivative of the function, ei is determined by ci and ∂ei = Ai,j ∂ci thus n n i=1 i=1 X X ∂s ∂ei = 2ei = 2(Ai. c − yi )Ai,j = 0 ⇔ A> Ac = A> y ∂ci ∂cj The linear equation system A> Ac = A> y is called the normal equations. The normal equations are not overdetermined and will be solvable if the x-values are different. Finding the solution the the normal equations gives the coefficients that minimizes the error function s. The reason that A> Ac = A> y are called the normal equations is that B = A> A is a normal matrix. Definition 8.1 A normal matrix is a matrix for which AH A = AAH . An alternate definition would be A = UDUH where U is a unitary matrix and D is a diagonal matrix. Proof that the definitions are equivalent. Here we can use a type of factorization that we have not discussed previously in the course called Schur factorization: A = URUH 95 where U is unitary and R triangular. AH A = UUH UH URUH = URH RUH = URRH UH = AAH if and only if R is diagonal (check this yourself using matrix multiplication). AH A = UDH UH UDUH = UDH DUH = UDDH UH = AH A since diagonal matrices always commute. Normal matrices are always diagonalizable which means they are always invertible if D has full rank. If B = A> A and A has full rank (equivalent to rank(A) = n) then B = UDUH with rank(D) = n. Thus: if rank(A) = n then x = (A> A)−1 A> y minimizes the error function or ’solves’ Ax = y in the least square sense. Sometimes the normal equations can become easier to solve if the A-matrix is factorized. One useful example of a factorization that can simplify the calculation is the QR-factorization that you will learn how to do in section 11.2. Using QR-factorization we get A = QR with Q being an orthogonal matrix, Q> = Q−1 , and R an upper triangular matrix. The normal equations can then be rewritten in this way: (A> A)c = (QR)> QR c = = R> Q> QR c = R> R c = = A> y = (QR)> y = R> Q> y ⇔ ⇔ Rc = Q> y 8.2 Generalized LSM LSM can be used to fit polynomials, as shown previously, but also other functions that can be written on the form X f (x) = gi (x) i For example, if we want to fit a trigonometric function, f (x), instead of a polynomial to some data and f (x) looks like this f (x) = a sin(x) + b sin(2x) + c cos(x) then the fitting problem can be described in this way Ac = y sin(x1 ) sin(2x1 ) cos(x1 ) y1 a sin(x2 ) sin(2x2 ) cos(x2 ) y2 A= . .. .. , c = b , y = .. .. . . . c yn sin(xn ) sin(2xn ) cos(xn ) 96 For this kind of problems we do not have the same simple guarantees that the A matrix is invertible though. There are also many functions that can be rewritten on the form X f (x) = gi (x) i for example, suppose we have a trigonometric function with a polynomial argument f (x) = tan(a + bx + cx2 ) ⇒ arctan(f (x)) = a + bx + cx2 And even though we have a non-linear model we have turned it into the linear problem Ac = y 2 1 x1 x1 arctan y1 a 1 x2 x2 arctan y2 2 A = . . , c = b , y = .. . .. .. .. . c 2 arctan yn 1 xn xn that can be solved with the regular least square method. Note thought that when doing this we no longer fit the original function but instead the rewritten function. In some applications this is problematic and in some it is not. We will not concern ourselves further with this in this course but in a real application this sort of rewrite should not be used indiscriminately. 8.3 Weighted Least Squares Method (WLS) It was previously described that the least square method was useful ion situations where some data had been collected and we wanted to analyze the data or adapt a specific model to the data.In both these cases certain parts of the data can be more relevant than the other. The data can be a compilation of data from many different sources, it can be measured by a tool that is more efficient for certain x or there is a certain area where it is more important that the model fits the data better than everywhere else, for example if we describe a system that is intended to spend most time in a state corresponding to value x the it probably more important that our model fits well near that x than far away from that x. This can be taken into consideration by using the weighted least square method . When we use the weighted least square method we add a a weight in the error function. s(e) = (e)> W(e) = n X Wi,i |ei |2 i=1 Then we get a new set of normal equations A> WAx = A> Wy 97 How to choose W? For physical experiments W is dictated by physical properties of the measured system or physical properties of the measurement tools. In statistics the 1 and correct weights for independent variables are the reciprocal of the variance wi = σi for variables with the correlation matrix Σ the weights are decided by W> W = Σ−1 to get the best linear unbiased estimate. If the different variables in x are independent then W is diagonal but if there is some correlation between them then there will be off-diagonal elements as well. 8.4 Pseudoinverses For square matrices we could describe the solution of a linear equation system using the inverse matrix Ax = y ⇔ x = A−1 y Not all matrices have an inverse, diagonalizable matrices where one eigenvalue is zero for example, but the idea behind the inverse can be generalized to and made to extend to these matrices as well. For the regular inverse we have AA−1 = A−1 A = I In this section we are going to discuss pseudoinverses, which are denoted by the following rules A∗ AA∗ = A∗ , AA∗ A = A ∗ and obey The inverse of a matrix is unique but a matrix can have many different pseudo inverses, for example, if: 0 1 A= 0 0 then all of the following are pseudo-inverses 1 1 0 0 1 0 + + + , A3 = , A2 = A1 = 1 1 1 1 1 0 In engineering there is a very important pseudoinverse called the Moore-Penrose pseudoinverse. Definition 8.2 The Moore-Penrose pseudoinverse of a matrix A is defined as follows A∗ = (A> A)A> The Moore-Penrose pseudoinverse is interesting because it is closely connected to the least square method and writing Ax = y ⇔ x = A∗ y 98 is the same as solving Ax = y using the least square method. Another good property of the Moore-Penrose pseudoinverse is that if the matrix A is invertible then A∗ = A−1 . In practice calculating the Moore-Penrose pseudoinverse is often done using the singular value factorization (SVD) or QR-factorization. Showing how to calculate the Moore-Penrose pseudoinverse using QR-factorization is very similar to showing how it is done for the least square method. −1 A∗ = (A> A)−1 A> = (QR)> QR (QR)> = −1 −1 = R> Q> QR R> Q> = R> R R> Q> = = R−1 (R> )−1 R> Q> = R−1 Q> 99 9 Matrix functions Functions are a very useful mathematical concept in many different areas. Power function: f (x) = xn , x ∈ C, n ∈ Z Root function: f (x) = Polynomials: p(x) = √ n x, x ∈ R, n ∈ Z+ n X ak xk , a, x ∈ C, n ∈ Z+ k=0 Exponential function: f (x) = ex , x ∈ C Natural logarithm: f (x) = ln(x), x ∈ C Table 1: Some common functions on the real and complex numbers. Each of the functions in table ?? are based on some fundamental idea that makes them useful for solving problems involving numbers of different kinds. If new functions that have corresponding or similar properties to these functions but acts on matrices instead of numbers could be defined then there might be ways to use them to solve problems that involve matrices of some kind. This can actually be done to some extent and this section will explore this further. 9.1 Matrix power function The normal power function with positive integer exponents looks like this f (x) = xn = |x · x ·{z . . . · x}, x ∈ C, n ∈ Z+ n Inspired by this it is quite intuitive to define a corresponding matrix version of the power function Definition 9.1 For any square matrix A ∈ Mk×k , k ∈ Z+ the matrix power function is defined as An = AAA . . . A}, n ∈ Z+ | {z n + for any n ∈ Z . The regular power function is not only defined for positive integers, but also for 0 and negative integers. How can definition 9.1 be extended to let n be any integer, not just a positive one? First consider how that x0 = 1 for any complex number x. For matrices this corresponds to 100 Definition 9.2 For any square matrix A ∈ Mk×k , k ∈ Z+ A0 = I where I is the identity matrix. For negative exponents it is well known that 1 , n ∈ Z+ xn Here it is not as intuitive and simple to find a corresponding definition for matrices. But consider that x−n = x· 1 = xx−1 = 1 x 1 can be viewed as the inverse of x. This means that there is a x natural way to define A−n , n ∈ Z+ for any invertible matrix A. This means that Definition 9.3 For any square and invertible matrix A ∈ Mk×k , k ∈ Z+ + −1 −1 −1 −1 A−n = (A−1 )n = A | A A{z . . . A }, n ∈ Z n where A−1 is the inverse of A. Now a matrix version of the power function An with integer exponent n ≥ 0 has been defined for any square matrix A. If A is also invertible An is defined for all integer n. When discussing (real or complex) numbers, x, the function xn is not limited to integer n. In fact, n can be any (real or complex number). In the next section the definition of An is extended to fractional n and in section 9.5 we will show how to define matrix versions of more general functions, including An for any (real or complex) number n. 9.1.1 Calculating An Calculating the matrix An can be done in a very straightforward manner by multiplying A with itself n times (or the inverse of A with itself −n times if n is negative). But the calculations can sometimes be simplified further. Previously, see page 13 and definition 2.1, the concept of a diagonalizable matrix has been examined. If A is a diagonalizable matrix, A = S−1 DS, then the power function (for positive integer exponents) becomes −1 −1 n An = AA . . . A = S−1 D |SS{z−1} D |SS{z−1} . . . SS | {z } DS = S D S =I =I =I n where D of course becomes 101 n d1,1 0 . . . 0 dn . . . 2,2 Dn = . .. .. . . . . 0 0 ... 0 0 .. . dnm,m Thus if a matrix A is simple to diagonalize then An becomes much less tedious to calculate for high n. The same idea can be generalized to similar matrices. Consider two matrices, A and B, that are similar, in the mathematical sense given in definition 1.3. A = P−1 BP In this case the power function of A can be calculated as follows −1 −1 −1 −1 n An = AA . . . A = P−1 B |PP {z } B |PP {z } . . . PP | {z } BP = P B P =I =I =I In some cases, for example when B is a diagonal matrix and n is large, it can be easier to calculate An in this way. In theorem 2.8 we saw that even if a matrix was not diagonalizable it could be rewritten as a matrix which is almost diagonal Jm1 (λ1 ) 0 ... 0 0 Jm2 (λ2 ) . . . 0 A = P−1 JP = P−1 P .. .. . . . ... . . 0 0 . . . Jmk (λk ) Here Jmi , 1 ≤ i ≤ k are all Jordan blocks, that is matrices with the value λi on the diagonal and 1 in the element above the diagonal, see definition 2.7. k is the number of different eigenvalues of A. This means that for any square matrix A the power function can be written like this Jm1 (λ1 )n 0 ... 0 0 Jm2 (λ2 )n . . . 0 n −1 n −1 A =P J P=P P .. .. . . . . . ... 0 0 . . . Jmk (λk )n Since all the Jordan blocks, Jmi (λi ), are either small or contains a lot of zeros each Jmi (λi )n is simple to calculate and thus Jn will often be easier to calculate than An . One kind of matrix for which it can be easier to calculate An are nilpotent matrices. Definition 9.4 A square matrix An is said to be a nilpotent matrix if there is some positive integer, n ∈ Z+ , such that An = 0 From this it follows that Ak , k ≥ n. 102 For nilpotent matrices it obviously easy to calculate An for sufficiently high numbers. One example of a nilpotent matrix is a strictly triangular matrix (a triangular matrix with zeros on the diagonal). Remember any square matrix is similar to a matrix written on the Jordan canonical form, which is an upper triangular with the matrix eigenvalues along the diagonal. This means that if all the eigenvalues for a matrix are equal to zero then the matrix will be similar to a strictly triangular matrix which meant that it will be nilpotent. 9.2 Matrix root function For (real or complex numbers) √ we have the following relation between the power function xn and the root function n x y = xn ⇔ x = √ n 1 y = yn For any positive integer, n ∈ Z+ , it is easy to do the corresponding definition for matrices Definition 9.5 If the following relation holds for two square matrices A, B ∈ Mk×k A = BBB . . . B} | {z n then B is said to be the nth root of A. This is annotated B = √ n 1 A = An √ n Simply defining the function is not very useful though. In order to make A useful a few questions need to be answered √ n • How do you find A? √ n • For how many different B is B = A √ √ 2 In this course we will only consider the square root of a matrix A = A . 9.2.1 Root of a diagonal matrix It is not obvious how to find the square root of a matrix, it is not even clear whether there exist any roots at all for a given matrix. For this reason it is suitable to first consider a specific kind of matrix for which the roots exist and are easy to find. Diagonal matrices (see page 6) are generally simple to work with. It is also easy to find their roots. Consider a diagonal matrix D ∈ Mm×m . d1,1 0 . . . 0 0 d2,2 . . . 0 D= . .. .. .. .. . . . 0 0 . . . dm,m 103 Then we can define a matrix C in this way p C= d1,1 p0 ... 0 d1,1 0 . . . 0 0 d2,2 . . . 0 0 d2,2 . . . 0 2 ⇒ CC = C = . . .. .. .. . . .. . . . . . . . . . . . p . 0 0 . . . dm,m 0 0 ... dm,m (3) From the definition it follows that C is a root of D. It is not the only root though. Each of the diagonal elements di,i has two roots, one positive and one negative and regardless which root that is chosen in C equation (3) will still hold. This means that there are 2n different Cs that are roots of D. These are not the only possible roots though. For an example of how many roots a matrix can have, see the following theorem. 1 0 Theorem 9.1. All of the following matrices are roots to I2 = 0 1 1 ∓s ∓r 1 ±s ∓r ±1 0 1 0 , , , 0 1 0 ±1 t ∓r ±s t ∓r ∓s if r2 + s2 = t2 . Thus the simple matrix I2 has infinitely many roots. Clearly roots are more complicated for matrices than they are for numbers. Fortunately though it is not always necessary to find all the roots of a matrix, often a single root can be enough, and there are several ways to make it easier to construct. For a diagonal matrix D it is very simple to find a root, simply take the matrix that had the positive root of all the elements in the diagonal of D in its diagonal. This is often referred to as the principal root. 9.2.2 Finding the square root of a matrix using the Jordan canonical form The technique for finding a root for a diagonal matrix can be extended to diagonalizable matrices. Consider a diagonalizable matrix A √ √ A = S−1 DS = S−1 D DS = √ √ √ = S−1 DSS−1 DS = S−1 DS ⇒ √ √ ⇒ A = S−1 DS By an analogous argument you can find a similar relation between any square matrix and its Jordan canonical form √ √ A = P−1 JP Finding the square root for J is not as simple as it is for a diagonal matrix, but there are some similarities 104 p Jm1 (λ1 ) p 0 0 Jm2 (λ2 ) √ √ . .. A = P−1 JP = P−1 .. . 0 0 ... ... .. . 0 0 P . . . q ... Jmk (λk ) √ This means that if a square root of the Jordan blocks are known then A can be calculated. Let J be a Jordan block √ λ 1 0 ... 0 1 λ √0 . . . 0 0 λ 1 . . . 0 0 1 λ . . . 0 1 . . . 0 J = 0 0 λ . . . 0 = λ 0 0 = λ(I + K) .. .. .. . . . .. . . . . .. .. . . .. . . . .. . . 0 0 0 ... λ 0 0 0 ... 1 Consider the function f (x) = f (x) = √ 1 + x, this function can be written as a Taylor series √ f 00 (a) 2 f (3) (0) 3 f (4) (0) 4 f 0 (0) x+ x + x + x = 1 + x = f (0) + 1! 2! 3! 4! 1 1 5 4 1 x + ... = 1 + x − x2 + x3 − 2 8 16 128 It turns out that if you switch 1 for I and x for K the equality still holds (see page 115). K is also a strictly upper-triangular matrix and is therefore a nilpotent matrix (see definition 9.4). Thus the following formula will always converge √ 1 1 1 5 4 I + K = 1 + K − K2 + K3 − K + ... 2 8 16 128 which makes it possible to find a square root for any square matrix. 9.3 Matrix polynomials A polynomial of degree n over a set K (for example the real or complex numbers) is a function p(x) that can be written as p(x) = n X ak xk , ak ∈ K k=0 where n is a finite non-negative integer and x is a variable taken from some set such that xn and addition is defined in a sensible way (more specifically x can only take values from a ring). Examples of sets that the variable x can belong to are real numbers, complex numbers and matrices. 105 Definition 9.6 A matrix polynomial (with complex coefficients) is a function p : Mn×n 7→ Mn×n that can be written as follows p(A) = n X ak Ak , ak ∈ C, A ∈ Mn×n k=0 where n is a non-negative integer. Matrix polynomials has a number of convenient properties, one of them can be said to be inherited from the power function Theorem 9.2. For two similar matrices, A and B, it is true for any polynomial p that p(A) = p(P−1 BP) = P−1 p(B)P Proof. This theorem is easily proved by replacing A with P−1 BP in the expression describing the polynomial and using An = P−1 Bn P. 9.3.1 The Cayley-Hamilton theorem A very useful and important theorem for matrix polynomials is the Cayley-Hamilton theorem. Theorem 9.3 (Cayley-Hamilton theorem). Let pA (λ) = det(λI − A) then pA (A) = 0 The Cayley-Hamilton theorem can be useful both for simplifying the calculations of matrix polynomials and for solving regular polynomial equations. This will be demonstrated in the following sections. The rest of this section focuses on proving the CayleyHamilton theorem. Proving the Cayley-Hamilton theorem becomes easier using the following formula. Theorem 9.4 (The adjugate formula). For any square matrix A the following equality holds A adj(A) = adj(A)A = det(A)I (4) where adj(A) is a square matrix such that adj(A)ki = Aki where Aki is the ki-cofactor of A, (see equation (2)) Proof of the adjugate formula. Using equation (2) gives det(A) = n X ai,k Ai,k k=1 If we collect all the cofactors, Aki , in a matrix 106 (adj(A))ki = Aik Note the transpose. It can then easily be verified that n X n X a1,k A1k k=1 n X a2,k A1k A adj(A) = k=1 .. n . X a A n,k k=1 n X k=1 n X 1k k=1 k=1 n X det(A) n X a2,k A1k = k=1 .. n . X a A n,k a1,k A2k a2,k A2k .. . an,k A2k a1,k A2k det(A) .. . n X k=1 a1,k Ank k=1 n X ... a2,k Ank = k=1 .. .. . . n X ... an,k Ank k=1 k=1 1k ... n X an,k A2k ... n X a1,k Ank ... a2,k Ank k=1 . .. .. . ... det(A) k=1 n X (5) k=1 In order to get equation (4) we need to show that all off-diagonal element are equal to zero. To do this, consider a square matrix B. From this matrix we construct a new matrix, C, by replacing the i:th row with the j:th row (i 6= j). This means that in C row i and row j are identical cik = cj,k = bj,k , 1 ≤ k ≤ n The matrices B and C will also share any cofactors created by removing row i Cik = Bik Using (2) on C gives det(C) = n X k=1 ci,k Cik = n X cj,k Cik , i 6= j = k=1 n X bj,k Bik k=1 Since C has two identical rows the rows are not linearly independent which means det(C) = 0 Thus, for any square matrix B it is true that 107 n X bj,k Bik = 0 k=1 This means that all elements in matrix (5) except the diagonal elements are equal to zero. This gives the adjugate formula A adj(A) = det(AI Checking that the same holds for adj(A)A can be done in an analogous way. Now the adjugate formula can be used to prove the Cayley-Hamilton theorem Proof of the Cayley-Hamilton theorem. Set pA (λ) = det(λI − A) and let B = adj(λI − A) then the adjugate formula gives pA (λ)I = (λI − A)B Every element in B is a polynomial of degree n − 1 or lower, since each element in B is the determinant of a (n − 1) × (n − 1) matrix. Thus B can be rewritten as a matrix polynomial n−1 X λk Kk B= k=0 where Bk is a matrix containing the coefficients of the polynomial terms with degree k in B. Combining the to formulas gives (λI − A)B = (λI − A) n−1 X λk Bk = k=0 = λn Bn−1 + λn−1 (Bn−2 − ABn−1 ) + . . . + λ(B0 − AB1 ) + (−AB0 ) = = λn I + λn−1 an−1 I + . . . + a1 λ + a0 = pA (λ)I If this expression holds for all λ then the following three conditions must be fulfilled ak I = Bk−1 − ABk for 1 ≤ k ≤ n − 1 I = Bn−1 (7) a0 I = −AB0 (8) Next consider the polynomial pA (A) pA (A) = An + An−1 an−1 + . . . + Aa1 + Ia0 Combining this with condition (6)-(8) results in 108 (6) pA (A) = An + An−1 (Bn−2 − ABn−1 ) + . . . + A(B0 − AB1 ) + (−AB0 ) = = An (I − Bn−1 ) +An−1 (Bn−2 − Bn−2 ) + . . . + A (B0 − B0 ) = 0 | | {z } | {z } {z } =0 =0 =0 9.3.2 Calculating matrix polynomials In section 9.1.1 a number of different techniques and properties worth considering when calculating the power function was described. Naturally all of this is applicable to matrix polynomials as well. But there is other techniques that can be used to simplify the calculations of a matrix polynomial. On page 106 it was mentioned that the Cayley-Hamilton theorem could be used to calculate matrix polynomials. It is very simple to see hoe the following theorem follows from Cayley-Hamilton. Theorem 9.5. Let pA (λ) = det(λI − A) where A is a square matrix. Then the following formula holds pA (A) = 0 c · pA (A) = 0 where c ∈ C q(A) · pA (A) = 0 where q is a polynomial Many things that can be done with polynomials of numbers can also be with matrix polynomials. Factorization is one such thing. If p(x) = (x − λ1 )(x − λ2 ) . . . (x − λn ) then p(A) = (A − λ1 I)(A − λ2 I) . . . (A − λn I) This means that if a polynomial on normal numbers can be simplified by factorization, the corresponding matrix polynomial can be simplified in a corresponding way. It also means that polynomial division can be used for matrix polynomials as well as for normal polynomials. That is, any matrix polynomial, p(A), with degree n can be rewritten as p(A) = q(A)π(A) + r(A) where q(A), π(A) and r(A) are all polynomials of degree deg(q) = n − k, deg(π) = k and deg(r) ≤ n with 1 ≤ k ≤ n. If π(A) is chosen as one of the polynomials in theorem 9.5 then π(A) = 0 and p(A) = r(A) with potentially lowers the degree of p(A) making it easier to calculate. One especially interesting consequence of this is the following theorem Theorem 9.6. For two polynomials, a and b, with different degree, deg(a) ≥ deg(b), such that a(A) = 0 and b(A) = 0 then the following relation must apply a(A) = c(A)b(A) where c(A) is a polynomial with degree deg(c) = deg(a) − deg(b). Proof. Dividing a and b gives a(A) = c(A)b(A) + r(A). From b(A) = 0 it follows that a(A) = r(A) but a(A) = 0 thus r(A) = 0. 109 9.3.3 Companion matrices If λ is an eigenvalue of a square matrix A then det(λI − A) = 0. Since an n × n matrix has n eigenvalues (some of which can be the same) and since pA (λ) = det(λI − A) is a polynomial of λ with degree n this means that each eigenvalue corresponds to a root of the polynomial pA (λ). This can be expressed as a formula pA (λ) = det(λI − A) = (λ − λ1 )(λ − λ2 ) . . . (λ − λn ) where λi for 1 ≤ i ≤ n are eigenvalues of A. This means that if we have found the eigenvalues of some matrix we have also found the roots of some polynomial and vice versa. That solving polynomials can be used to find eigenvalues of matrices is something that we have seen before, the standard method of solving det(λI − A) = 0 is exactly this. But there are other ways of finding the eigenvalues of a matrix, in section 11 a number of different ways of calculating the roots numerically will be shown for example. There is also a simple way of constructing a matrix that has the roots of a given polynomial as eigenvalues. This matrix is called the companion matrix of a polynomial. The origin of this kind of matrix is linear differential equations. Consider a linear differential equation of order n y (n) + an−1 y (n−1) + . . . + a1 y (1) + a0 = f where y (k) is the k:th derivative of y and f is a function. When solving this kind of differential equation it is often helpful to find the roots of the characteristic polynomial p(x) = xn + an−1 xn−1 + . . . + a1 x + a0 To get the companion matrix we rewrite the differential equation using state variables x0 = y x1 = y (1) .. . xn−1 = y (n−1) xn = y (n) Using the state variables x1 , . . . xn the differential equation can be written as a system of first order differential equations. If we denote a vector containing all xi with x and a vector containing all x0i with x0 then the differential equation can be written as a matrix equation 110 x0 = C(p)> x ⇔ x00 = x1 x01 = x2 .. . x0n−1 = xn x0n = −a0 x0 − a1 x1 − . . . − an−1 xn + f With C(p) defined as follows. Definition 9.7 Let p(x) = a0 + a1 x + a2 x2 + . . . + an−1 xn−1 + xn 0 0 ... 0 1 0 . . . 0 C(p) = 0 1 . . . 0 .. .. . . .. . . . . then the companion matrix of p is −a0 −a1 −a2 .. . 0 0 . . . 0 −an−1 In some literature C(p)> is referred to as the companion matrix instead of C(p). Theorem 9.7. If C(p) is a companion matrix to the polynomial p then the eigenvalues of C(p) are roots of p. Proof. Let ei denote vectors in the standard basis, that is column vectors with a 1 in element i and 0 everywhere else. For these vectors multiplication with C(p) will result in the following C(p)ei = ei+1 for 0 ≤ i < n n X C(p)en = −ai ei (9) (10) i=1 From (9) it can be easily shown that ei = C(p)i e1 for 0 ≤ i < n n en = C(p) e1 (11) (12) Combining equation (9)-(12) gives C(p)n e1 = n−1 X −ai C(p)i e1 ⇔ i=0 n n−1 ⇔ C(p) + an−1 C(p) + . . . + a1 C(p) + a0 I = 0 ⇔ p(C(p)) = 0 111 By the Cayley-Hamilton theorem pC (C(p)) = 0 with pC (λ) = det(λI − C(p)). Thus we have two polynomials that are zero for the companion matrix. p(C(p)) = 0 pC (C(p)) = 0 Then we can use theorem (9.6) to see that that p(A) = c · pC (A) (since p and pC have the same degree). This means that p(λ) and pC (λ) must have the same roots. The eigenvalues of C(p) are roots of pC (λ) and thus the theorem is proved. Calculating roots of polynomials using companion matrices can be quite efficient, especially since a polynomial of high degree gives a sparse matrix , that is a matrix with only a few non-zero elements. Sparse matrices can be used quite efficiently by computers and many algorithms can be specifically modified to work much faster on sparse matrices than non-sparse matrices. The power method described in section 11.1 can be used to efficiently find numerical eigenvalues of sparse matrices. Many computer programs that finds roots of polynomials uses companion matrices, for example the roots command in MATLAB (see [5]). 9.4 Exponential matrix So far we have constructed matrix functions by taking the definition of a function for numbers and replaced numbers with matrices when it was possible. Now we are going to construct a matrix version of a function in a different way, we are going to identify the most important properties of a function and then construct a matrix function that has the same properties. The function that we are going to construct a matrix version of the exponential function. A common definition for the exponential function is x n ) n The exponential function is useful in many ways but most of its usefulness comes from the following three properties ex = lim (1 + n→∞ d at (e ) = aeat = eat a dt (eat )−1 = e−at ea · eb = ea+b = eb+a = eb · ea (13) (14) (15) Here a, b and t are all real or complex numbers. Multiplication of numbers is usually commutative (a · b = b · a) but matrix multiplication is often not commutative. In equation (13)-(15) the various places where commutativity comes in has been written out explicitly since we want to have commutativity in the same way for the matrix version. 112 The most complex function we have defined so far is the matrix polynomial. Can we construct a polynomial such that it fulfills (13)-(13)? Consider what happens to a general polynomial of At when we take the first time derivative p(At) = n X k=0 n n−1 k=1 l=0 X X d ak (At) ⇔ (p(At)) = ak kAk tk−1 = A al+1 (l + 1)Al tl dt k The last equality is achieved by simply taking l = k − 1. Compare this to Ap(At) = A n X al Al tl l=0 If the polynomial coefficients are chosen such that al+1 (l + 1) = al these expressions would be identical except for the final term. Let us choose the correct coefficients al+1 (l + 1) = al ⇔ al = 1 l! and simply avoid the problem of the last term not matching up by making the polynomial have infinite degree. Definition 9.8 The matrix polynomial function is defined as e At = ∞ k X t k=0 k! Ak , A ∈ Mn×n It is simple to create eA from this by setting t = 1. An important thing to note here is that since the sum in the definition is infinite we should check for convergence. As it turns out the exponential function will always converge. This follows from a more general theorem presented in section 9.5. Theorem 9.8. For the matrix exponential function the following is true d At (e ) = AeAt = eAt A dt (eAt )−1 = e−At A A e ·e =e A+B =e B+A (16) (17) B A = e · e if AB = BA (18) where A and B are square matrices. This definition of eAt has all the properties we want, with the exception that we need the extra condition on (18) compared to (15). There is no way around this unfortunately. That the definition gives a function that has property (16) is clear from the construction. Property (17) and (18) can be shown by direct calculations using the definition. Calculating the exponential matrix can be challenging. First of all we should note that the exponential function for numbers can be written like this 113 ex = ∞ X xk k=0 k! using a Taylor series. Using this we can see how to calculate the matrix exponential for a diagonal matrix D ∞ X dk1,1 0 ... 0 k! k=0 ∞ ed1,1 0 ... 0 k X d 2,2 0 ... 0 ed2,2 . . . 0 0 D k! e = = . . .. . k=0 . . . . . . . .. .. .. .. . . . . dm,m 0 0 ... e ∞ X dkn,n 0 0 ... k! k=0 Since the exponential matrix is essentially a polynomial it should not be surprising that the following relation holds. A = S−1 BS ⇒ eA = S−1 eB S On page 102 it was discussed that nilpotent matrices could make it easier to calculate the matrix power function. This is of course also true for the matrix exponential function, since it turn the infinites sum into a finite sum and thus turns the exponential function into a finite polynomial. There are also many ways of computing the exponential matrix function numerically. A nice collection can be found in [10]. 9.4.1 Solving matrix differential equations One of the most common applications of the exponential function for numbers is solving differential equations. df = kf (x) ⇒ f (x) = c · ekx dx The matrix exponential function can be used in the same way. Theorem 9.9. The exponential matrix exA is the solution to the initial value problem dU = AU dx U (0) = I Proof. From the construction of exA we already know that if fulfills the differential equation. Checking that the initial condition is fulfilled can be done simply by setting x = 0 in the definition. 114 9.5 General matrix functions There is a systematic way of creating a matrix version of a function. When the matrix exponential function was defined in the previous section (see page 113) we found a matrix polynomial of infinite degree that created the matrix version. One way of doing this is given below Definition 9.9 (Matrix function via Hermite interpolation) If a function f (λ) is defined on all λ that are eigenvalues of a square matrix A then the matrix function f (A) is defined as f (A) = h sX i −1 X f (k) (λi ) i=1 k=0 k! φik (A) where φik are the Hermite interpolation polynomials (defined below), λi are the different eigenvalues and si is the multiplicity of each eigenvalue. Definition 9.10 Let λi , 1 ≤ i ≤ n be complex numbers and si , 1 ≤ i ≤ n, be positive integers. Let s be the sum of all si . The Hermitian interpolation polynomials are a set of polynomials with degree lower than s that fulfills the conditions ( (l) 0 if i 6= j or k 6= l φik (λj ) = l! 1 if i = j and k = l Using this kind of definition to create a matrix function works as long as the function f is defined on the eigenvalues of the matrix. There is another definition af a matrix function that is equivalent to definition 9.9 but is only applicable to analytical functions, that is functions that can be written as a power series. Definition 9.11 (Matrix function via Cauchy Integral) Let C be a curve that encloses all eigenvalues of the square matrix A. If the function f is analytical (can be written as a power series) on C and inside C then Z 1 f (A) = f (z)(zI − A)−1 dz 2πi C The matrix RA (z) = (zI − A)−1 is called the resolvent matrix of A. Proofs that these two definitions makes sense and are equivalent can be found in [?] on page 1-8. The formula in definition 9.11 is a matrix version of Cauchy’s integral formula which can be used to derive the formula for a Taylor series. As a consequence of this is that we can construct the matrix version of any analytical function can by replacing the number variable in a Taylor series expansion (or other convergent power series) with a square matrix. For example of this see the formula for the matrix square root (page 105) and compare the definition of the matrix exponential function with the Taylor series for the exponential function for numbers. 115 10 Matrix equations A matrix equation is an equation that involves matrices and vectors instead of numbers. A simple example is Ax = y where A is a matrix, y is a known vector and x is an unknown vector. You are already familiar with some ways of solving this kind of equations, for example Gaussian elimination and matrix inverses. There are of course many other kinds of matrix equations. One example that we have seen before are the first order differential equation that was solved using the exponential matrix function, see page 114. There are many other kinds of matrix equations. A simple example would be exchanging the x and y vectors in the typical matrix equation for a pair of square matrices, X and Y. There are two variations of this kind of equation AX = Y XA = Y Another example of a matrix equation that is important in many applications (first and foremost in control theory) is Sylvester’s equation AX + XB = C There is also a version of Sylvester’s equation that is especially important known as the continous Liapunov1 equation. AX + XAH = −C There area many methods for solving specific matrix equations, but here we will take a look at a basic tool for examining certain properties of the solutions of matrix equations and turning the more exotic kind of matrix equations into the typical matrix-vector equation that we already know how to solve. For this a new tool will be introduced. 10.1 Tensor products The tensor product, usually denoted by ⊗ is an interesting tool that can be defined for many kinds of mathematical objects, not only matrices. It is sometimes called the outer product or direct product. A ⊗ B is usually pronounced ’A tensor B’. Definition 10.1 The tensor product for matrices is a bilinear function of matrices that has two properties: 1. Bilinearity: (µA + ηB) ⊗ C = µA ⊗ C + ηB ⊗ C A ⊗ (µB + ηC) = µA ⊗ B + ηA ⊗ C 2. Associativity: (A ⊗ B) ⊗ C = A ⊗ (B ⊗ C) 1 sometimes spelled Lyapunov or Ljapunow 116 It can be possible to find different function that fullfil these condition for at set of matrices but the most general (and most common) tensor product is the Kroenecker product. Definition 10.2 The Kroenecker product ⊗ between two matrices, as a1,1 B a1,2 B a2,1 B a2,2 B A⊗B= . .. .. . A ∈ Mm×n and B ∈ Mp×q , is defined . . . a1,n B . . . a2,n B .. .. . . am,1 B am,2 B . . . am,n B or Ab1,1 Ab2,1 A⊗B= . .. Ab1,2 Ab2,2 .. . ... ... .. . Ab1,n Ab2,n .. . Abm,1 Abm,2 . . . Abm,n these two definitions are equivalent but not equal, one is a permutation of the other. The Kroenecker product is quite simple and intuitive to work with, even if it creates large matrices. Theorem 10.1. For the Kroenecker product the following is true 1. (µA + ηB) ⊗ C = µA ⊗ C + ηB ⊗ C 2. A ⊗ (µB + ηC) = µA ⊗ B + ηA ⊗ C 3. (A ⊗ B) ⊗ C = A ⊗ (B ⊗ C) 4. (A ⊗ B)> = A> ⊗ B> 5. A ⊗ B = A ⊗ B (complex conjugate) 6. (A ⊗ B)H = AH ⊗ BH 7. (A ⊗ B)−1 = A−1 ⊗ B−1 , for all invertible A and B 8. det(A ⊗ B) = det(A)k det(B)n , with A ∈ Mn×n , B ∈ Mk×k 9. (A ⊗ B)(C ⊗ D) = (AC ⊗ BD), A ∈ Mm×n , B ∈ Mn×k , C ∈ Mp×q , D ∈ Mq×r 10. A ⊗ B = (A ⊗ Ik×k )(In×n × B), A ∈ Mn×n , B ∈ Mk×k The first 7 properties listed here shows that the Kroenecker product is bilinear, associative and react in a very intuitive way to transpose, complex conjugate, Hermitian conjugate and inverse operations. Number 8 shows that determinant are almost as intuitive but not quite (note that the determinant of A has been raised to the number of rows 117 in B and vice versa). Properties 9 and 10 are examples of how matrix multiplication and Kroenecker products relate to one another. Taking the Kroenecker product of two matrices gives a new matrix. Many of the properties if the old matrices relate to the properties of the old matrix in a fairly uncomplicated manner. One property that this applies to are the eigenvalues Theorem 10.2. Let {λ} be the eigenvalues of A and {µ} be the eigenvalues of B. Then the following is true: {λµ} are the eigenvalues of A ⊗ B {λ + µ} are the eigenvalues of A ⊗ B Sketch of proof. a) Use Shurs lemma A = S−1 RA S, B = T−1 RB T where RA and RB are triangular matrices with eigenvalues along the diagonal. Then A ⊗ B = (S−1 RA S) ⊗ (T−1 RB T) = = (S−1 RA ⊗ T−1 RB )(S ⊗ T) = = (S ⊗ T)−1 (RA ⊗ RB )(S ⊗ T) RA ⊗RB will be triangular and will have the elements ai,i bj,j = λi µj on the diagonal, thus λi µj is an eigenvalue of A⊗B since similar matrices have the same eigenvalues. b) Same argument as above gives A ⊗ I (B ⊗ I) have the same eigenvalues as A (B) adding the two terms together gives λ + µ is an eigenvalue. 10.2 Solving matrix equations The tensor product can be used to make it easy to solve equations like these AX = Y XA = Y Taking the tensor product between the identity matrix and some square matrix results in the following theorem Theorem 10.3. e=Y e AX = Y ⇔ (I ⊗ A)X e=Y e XA = Y ⇔ (A> ⊗ I)X where B.1 B.2 e= . . . B.n , B .. . B = B.1 B.2 B.n 118 This turns an equation consisting of n × n matrices into an equation with one n2 × n2 matrix and two n2 vectors. This equation can be solved the usual way. Solving the matrix this way is equivalent to assigning each element in the unknown matrix a variable and the solve the resulting linear equation system. Using the Kroenecker product can make subsequent manipulations of the matrices involved easier though, thanks to the properties shown in theorem 10.1. Often it is interesting to know if there exists a solution to a given equation and if that solution is unique. Using the Kroenecker product we can see for which matrices A and B that there is a unique solution to Sylvester’s equation. Theorem 10.4. Sylvester’s equation AX + XB = C has a unique solution if and only if A and −B have no common eigenvalues. Proof. Rewrite the equation using the Kroenecker product e=C e AX + XB = C ⇔ (I ⊗ A + B> ⊗ I) X | {z } K This is a normal matrix equation which is solvable if λ = 0 is not an eigenvalue of K. The eigenvalues of K are λ + µ where λ is an eigenvalue of A and µ is an eigenvalue of B. Thus if A and −B have no common eigenvalues the eigenvalues of K will 6= 0. As mentioned previously Sylvester’s equation is important in control theory especially the version of it that is called the continuous Liapunov equation. This will be illustrated using two theorems (neither theorem will be proved). A stable system is a system that will in some sense ’autocorrect’ itself. If left alone it will approach some equilibrium state and not increase indefinitely or behave chaotically. If the system is described by a first order differential matrix equation then the following condition gives stability. Definition 10.3 A matrix, A, is stable when all solutions to dx = Ax(t) dt converges to 0, x(t) → 0 when t → ∞. ⇔ All the eigenvalues of A have negative real component, Re(λ) < 0 for all λ that is an eigenvalue to A. Theorem 10.5. If A and B are stable matrices then AX + XB = −C has the unique solution Z X= ∞ etA CetB dt 0 119 Theorem 10.6. Let P be a positive definite matrix. If and only if there is a positive definite solution to AH X + CA = −P is A a stable matrix. 120 11 Numerical methods for computing eigenvalues and eigenvectors In practise we can rarely use a analytical solution to finding eigenvalues or eigenvectors of a system. Either because the system is two large making Gausselimination, by solving the characteristic polynomial or similar methods to slow or maybe we only need one or a few eigenvalues or eigenvectors or maybe a rough estimate is enough. You have probably seen by now how time consuming it is to find eigenvalues and eigenvectors analytically. Since there are a huge number of different problems where we want to find eigenvalues or eigenvectors there are also a large number of different methods with different strengths and weaknesses depending on the size and structure of the problem. We will look at two different related methods, the Power method and the QR-method and a couple of variations of both. While this chapter gives some insight into how to use or implement some methods, it also serves to help in showing that there is no generall method that somehow always works. If a general method doesn’t work, there might be another specified method for your type of problem. It’s then important to know what to look for both in our matrix as well as what we are actually trying to find. It also serves to show how we can start with a very primitive algorithm and then carefully build upon it in order to remove some or most of it’s limitations while still keeping it’s strengths. Some things to take into consideration when choosing which method to use is: • The size of the matrix, especially very large matrices which we maybe can’t even fit in memory. • Density of non-zero elements, Is the matrix a full matrix (few zeroes) or a sparse matrix (few non-zero elements). • Does the matrix have any kind of useful structure we could use, for example if it’s positive definite, symmetric, Hessian or a bandmatrix. • What do we seek? Do we need all eigenvalues and eigenvectors or just one or a few? Do we neeed exact result or is a good approximation enough? Or maybe we don’t need bth eigenvalues and eigenvectors? 11.1 The Power method The Power method is a very simple algorithm for finding the eigenvalue λ1 with the largest absolute value and it’s corresponding eigenvector. While the method have a couple of disadvantages making it unsuitable for general matrices, it proves it’s worth in a couple of specialized problems unsuitable for more general methods such as the QR-method which we will look at later. The Power method also acts as an introduction to many of the thoughts behind the QR-method and it’s variations. 121 The Power method in its basic form is a very simple and easy to implement method. We consider the matrix A for which we want to find the ”largest” eigenvalue and corresponding eigenvector. • We assume there is only one eigenvalue on the spectral radius and that A is diagonalizable. • The method works by iterating xk+1 = Axk . kAxk k Axk approaches a vector parallel to the eigenvector of the eigenkAxk k value with largest absolute value λ1 . • If λ1 6= 0 then • We get the eigenvalue λ1 from the relation Ak x ≈ λ1 Ak−1 x when k is large. It’s easy to see that if the Matrix A> is a stochastic matrix, running the Power method can be seen as iterating the Markov chain. This also motivates why we need there to be only one eigenvalue on the spectral radius since otherwise the Markov chain and therefor neither the Power method would converge in those cases. We will not look into how fast the convergence is here, but be content by noting that the relative error can be shown to be in the order of |λ2 /λ1 |k . In other words it’s not only other eigenvalues on the spectral radius that’s problematic, if there are eigenvalues close to the spectral radius they will slow down the convergence as well. We give a short example using the following matrix: −2 1 1 A = 3 −2 0 1 3 1 We start with x0 = [1 1 1]> and iterate: 0 1 −2 1 1 1 = 1 Ax0 = 3 −2 0 1 5 1 3 1 Next we calculate the norm of this vector: kAx0 k = √ 26 And last we use the norm kAx0 k in order to normalize Ax0 : x1 = √ Ax0 = [0, 1, 5]/ 26 kAx0 k Continuing for the next iteration we −2 1 Ax1 = 3 −2 1 3 122 get: 1 0 6 1 1 0 1 √ = −2 √ 26 26 1 5 8 x2 = √ Ax1 = [6, −2, 8]/ 104 kAx1 k And if we continue for a couple of more iterations we get: • x2 = [0.59, −0.17, 0.78]> . • x3 = [−0.25, 0.91, 0.33]> . • x4 = [0.41, −0.61, 0.67]> . • x7 = [−0.28, 0.85, −0.44]> . • x10 = [0.28, −0.79, 0.55]> . • For reference the true eigenvector is xtrue = [0.267, −0.802, 0.525]. We see that it does seem to converge, but note that it could take quite some time. However since every step gives an a little better approximation, if we only need a rough estimate we could stop considerably earlier. Because of the assumptions regarding the system as well as the problem with convergence if λ1 ≈ λ2 this simple version of the algorithm is rarely used unless we already have a good understanding of the system. In short a couple of things to consider when using the power method are: • If there is more than one eigenvalue on the spectral radius or if it’s multiplicity is larger than one, the method won’t converge. • Since the error is of the order |λ2 /λ1 |k , if λ2 is close to λ1 the convergence can be very slow. • We only get the (in absolute value) largest eigenvalue and corresponding eigenvector. While these are some considerable disadvantages (we can’t use it on all matrices and the convergence speed is uncertain), there are a couple of things we can do to handle these problems. However while the method does have a couple of negatives, it does have a couple of very strong points as well. • In every step we make only a single vector-matrix multiplication, making it very fast if we can ensure a fast convergence. • If A is a sparse matrix (few non-zero elements) the iterations will be fast even for extremely large matrices. • Since we do no matrix decomposition (as we do in most other common methods) we do not ”destroy” the sparcity of the matrix when calculating the matrix decomposition. 123 Especially we see the advantage of the method on sparse matrices since most other commonly used methods uses a matrix decomposition of the matrix, usually destroying the sparsity of the matrix in the process. This means that while the method might not be as reliable or as fast as other methods such as the QR-method in the general case, if we have a very large sparse matrix we might not even be able to use many other methods since they would be to slow to be usable. In some cases the matrix itself might not even fit in memory anymore after matrix decomposition further increasing the time needed. 11.1.1 Inverse Power method In order to handle some of the problems found earlier we will take a look at something called the inverse Power method. While this fixes some of the problems it doesn’t fix all of them, and there is a cost in that the method is no longer as suitible for sparse matrices as before. We consider the invertible matrix A with one eigenpair λ, x. From the definition of eigenvalues and eigenvectors we then have Ax = λx. This gives: Ax = λx ⇔ (A − µI)x = (λ − µ)x ⇔ x = (A − µI)−1 (λ − µ)x ⇔ (λ − µ)−1 x = (A − µI)−1 x And we see that ((λ − µ)−1 , x) is an eigenpair of (A − µI)−1 If we instead iterate using this new matrix (A − µI)−1 the method becomes: xk+1 = (A − µI)−1 xk k(A − µI)−1 xk k And just like the ordinary Power method this will converge towards the eigenvector of the dominant eigenvalue of (A − µI)−1 . The dominant eigenvalue of (A − µI)−1 is not necesarily the same as the dominant eigenvalue of A. Rather we get the dominant eigenvalue (λ − µI)−1 of (A − µI)−1 for the λ closest to µ. If we take a look at the error we get: (λ −1 k − µ) (second closest to µ) (λ(closest to µ) − µ)−1 (λ(closest to µ) − µ) k = (λ − µ) (second closest to µ) We can easy see that if µ is a good approximation of an eigenvalue of A then (λ(closest to µ) ≈ 0 and the error will be close to zero as well, this especially improves the convergence with multiple close eigenvalues. With a good approximation of λ the method will converge very fast, often in only one or a couple of iterations. Especially in the case where we know the eigenvalue exactly beforehand such as when working with irreducible stochastic matrices we get the error equal to zero and the method should converge in only 124 one step. We also see that if we change µ we can find the eigenvectors of not only the dominant eigenvalue, but the eigenvectors to multiple eigenvalues. However we still have problems with eigenvalues with multiplicity larger than one, or when more than one have the same absolute value. The inverse Power method also requires us to solve a linear system or calculate the inverse of the matrix (A − µI). This makes the method unsuitable for large sparse matrices where the basic Power method shows it’s greatest strength, since calculating the inverse matrix usually destroys the sparcity of the matrix. We give a short example using the inverse the dominant eigenvalue of: 3 3 1 1 A= 1 0 1 3 power method to find the eigenvector to 1 0 1 0 0 3 3 2 Since A is non-negative and irreducible (why?) we can use Perron Frobenius Theorem 4.1 to find bounds on the dominant eigenvalue. We remember this means the dominant eigenvalue is at least equal to the lowest row sum and at most equal to the largest row sum. Since at least one of the diagonal elements are non-zero we know that A is primitive as well so that we know that there is only one eigenvalue on the spectral radius. This ensures the convergence (when seeking the dominant eigenvalue). We choose µ as the largest possible value of the dominant eigenvalue: µ = max row sum = 7. This ensures that the dominant eigenvalue is the eigenvalue closest to µ, ensuring we converge towards the right eigenvector. We get: −4 3 1 0 1 −6 0 3 A − µI = A − 7I = 1 0 −6 3 1 3 0 −5 And we get its inverse as: (A − µI)−1 −0.4038 −0.1538 = −0.1538 −0.1731 −0.3173 −0.3590 −0.1923 −0.2788 −0.0673 −0.0256 −0.1923 −0.0288 −0.2308 −0.2308 −0.2308 −0.3846 Starting with x0 = [1 1 1 1]> gives: x1 = (A − µI)−1 x0 k(A − µI)−1 x0 k = [−0.5913 − 0.4463 − 0.4463 − 0.5020]> If we continue to iterate a couple of times and multiply the end result of those that are negative with −1 to end up with a positive vector we get: 125 • x1 = [0.5913 0.4463 0.4463 0.5020]> • x2 = [0.6074 0.4368 0.4368 0.4995]> • x3 = [0.6105 0.4351 0.4351 0.4986]> • xtrue = [0.6113 0.4347 0.4347 0.4983]> We see that we are very close to the true eigenvector after only three iterations, we in fact have a rather good approximation after only one iteration. Using the found eigenvector we can now also find a better approximation of the true dominant eigenvalue λ = 5.8445 as well from the relation Ax = λx. For example using x3 we get: Ax3 = [3.5721 2.5414 2.5414 2.9131]> → λ ≈ 3.5721 = 5.8511 0.6105 Which is very close to the true eigenvalue λ = 5.8445. The Inverse power method is often used when you have a good approximation of an eigenvalue and need to find an approximation of the corresponding eigenvector. This is often the case in real time systems where a particular eigenvector might needs to be computed multiple times during a very short time as the matrix changes. Overall both the Power method and the Inverse Power method are unsuitable for general matrices since we cannot guarantee the convergence without further inspecting the matrix. In order to solve the limitations of these methods we will need to calculate many eigenvalues at the same time as for example in the QR-method below. 11.2 The QR-method Often we might not know anything about the matrix beforehand and we might want to find not only one or a couple of eigenvalues/eigenvectors but all of them. We also don’t want to be limited to only finding the eigenpairs for simple eigenvalues. The QR-method solves all of these problems, it finds all eigenpairs, and it’s generally very fast unless we are working with very large matrices where we might need specified algorithms. Since the method is fast, reliable and finds all eigenpairs it is almost always the QR-method (or a variation there of) that is used whenever you calculate eigenvalues/eigenvectors using computational software. Unless you have very specific requirements or already know a more suitable method it’s usually a good idea to start by trying the QR-method and then continue looking at other methods if that is unsufficient. 126 Like the Power method, the QR-method is a iterative method, which after some time converge and gives us the eigenvalues and eigenvectors. Like the name suggest the method uses the QR-decomposition of the original matrix. We remember that we could find a QR decomposition A = QR where Q is a unitary matrix and R is upper triangular using for example the Gram-Schmidt process. The method works by iterating: • We start by assigning A0 := A • In every step we then compute the QR-decomposition Ak = Qk Rk . −1 H • We then form Ak+1 = Rk Qk = QH k Qk Rk Qk = Qk Ak Qk = Qk Ak Qk • Since Q is unitary all Ak are similar and therefor have the same eigenvalues. • When the method converge Ak will be a triangular matrix and therefor have it’s eigenvalues on the diagonal. The QR-method can be seen as a variation of the Power method, but instead of working with only one vector, we work with a complete basis of vectors, using the QR-decomposition to normalize and orthogonalize in every step. In practice a few modifications are needed in order to handle a couple of problems in it’s most basic form some of which are: • Computing the QR-decomposition of a matrix is expensive, and we need to do it in every step. • The convergence itself like with the Power method can also be very slow. • Like the Power method we have problems with complex and non-simple eigenvalues. However if we can handle these problems we have a very stable and fast algorithm for general matrices. Especially the stability of the algorithm is one of it’s huge selling point and the biggest reason why it’s used over other similar methods using for example the LUP decomposition. But how do we improve the algorithm? We start by looking at how we can speed up the computations of the QR-decomposition since even if we use Gram Schmidt or some other method it’s still a very expensive operation. We remember a matrix (A) is (upper) Hessenberg if all elements ai,j = 0, i > j + 1. In other words A is of the form: F F F ··· F F F F F F · · · F F F 0 F F · · · F F F A = 0 0 F · · · F F F .. .. .. . . . .. .. . . .. . . . . 0 0 0 · · · F F F 0 0 0 ··· 0 F F 127 If A is Hessenberg then we would speed up the algorithm considerably since: • Calculating the QR-factorization of a Hessenberg matrix is much faster than for a generall matrix. • Since multiplication of a Hessenberg matrix with a triangular matrix (both in the same direction) is itself a Hessenberg matrix. If A0 is Hessenberg every other Ak , k = 0, 1, 2, . . . is Hessenberg as well. • Since a Hessenberg matrix is nearly triangular, it’s convergence towards the triangular matrix will be improved, resulting in fewer total steps as well We don’t want to limit our algorithm to only Hessenberg matrices, however what if we can change the matrix we work with to one in Hessenberg form? This is possible since: Theorem 11.1. Every square matrix A have a matrix decomposition A = UH HU where • U is a unitary matrix. • H is a Hessenberg matrix. Since U is a unitary matrix, A and U have the same eigenvalues so in order to find the eigenvalues of A we can instead find the eigenvalues of H which is as we wanted Hessenberg. A Hessenberg decomposition can be calculated effectively using for example householder transformations. While we generally have a faster convergence using H, it can still be very slow. We remember we could solve this problem in the Power method by shifting the eigenvalues of the matrix. In the same way we can shift the eigenvalues in the QR-method using (Ak − σI). For Ak we then get: H Ak+1 = QH k (Ak − σI)Qk + σI = Qk (Ak )Qk We see that Ak , Ak+1 is still similar and therefor will have the same eigenvalues. As in the Power method we want σ to be as close as possible to one of the eigenvalues. Usually the diagonal elements of the hessenberg matrix is used since they prove to give a rather good approximation. Since the diagonal elements of Ak should give increasingly accurate approximation of the eigenvalues we can change our shift to another eigenvalue or improve the current one between iterations as well. As fast as a subdiagonal element is close to zero, we have either found an eigenvalue if it’s one of the outer elements. Or we can divide the problem in two smaller ones since we would have: A1,1 A1,2 A= 0 A2,2 128 This still leaves the problem with complex eigenvalues in real matrices. In order to solve this problem an implicit method using multiple shifts at the same time is used in practice. We leave this implicit method, also called the Francis algorithm by noting that it exist and works by much the same principle as the inverse Power algorithm. We give a summary of the modified QR-method using a single shift: 1. Compute the Hessenberg form (A = UH HU) using for example Householder transformations. 2. Assign A0 = H. Then in every step, choose a shift σ as one of the diagonal elements of Ak . 3. Compute the Qk Rk factorization of (Ak − σI). 4. Compute Ak+1 = Rk Qk + σI 5. Iterate 2 − 4 until Ak is a triangular matrix, optionally dividing the problem in two whenever a subdiagonal element becomes zero. 11.3 Applications 11.3.1 PageRank Pagerank is the method originally designed by Google to rank homepages in their search engine in order to show the most relevant pages first. In the ranking algorithm we consider homepages as nodes in a directed graph, and a link between two homepages as an edge between the corresponding vertices in the graph. Since the number of homepages is huge (and rapidly increasing) this does not only put high demand on the accuracy of the model (few looks at more than the first couple of pages of results) it also puts a huge demand on the speed of the algorithm. The principle used in the ranking is such that homepages that are linked to by many other homepages and/or homepages that are linked to by other important homepages are considered more important (higher rank) than others. An image illustrating the principle can be seen in Figure ??. PageRank uses the adjacency matrix 3.1 of the graph with a couple of modifications. First we do not allow any homepage to link to itself. Secondly, any homepages which links to no other page, usually called dangling nodes needs to be handled separately. Usually they are assumed to instead link to all pages (even itself even though we otherwise don’t allow it). This is the method we will use here, although there are other methods available as well. Before calculating PageRank we make a modification to the adjacency matrix. We weight the matrix by normalizing the sum of every non-zero row such that it sums to one. This is our matrix A. Pagerank is then the eigenvector normalized such that the sum is equal to one of the dominant eigenvalue of matrix M : M = c(A + gv> )> + (1 − c)ve> 129 image by Felipe Micaroni Lalli Where the 0 < c < 1 is a scalar, usually c = 0.85, e is the one vector and A as above is a non-negative square matrix whose rows all sum to one. v is a non-negative weigh vector with the sum of all elements equal to one. g is a vector with elements equal to zero for all elements except those corresponding to dangling nodes which are equal to one. We see that saying that we consider dangling nodes linking to all other nods is not strictly true, rather we assume they link according to a weight vector with sum one, usually however all or at least most of the elements of v are equal such that this is more or less true. To illustrate the method we give a short example where we consider the following system: 1 2 3 4 5 We assume all elements of v equal such that we give no node a higher weight. Using the matrix for the graph we get the matrix (A + gv> ) by normalizing every row to one and add the change on the third row corresponding to the dangling node we get: 0 1/3 1/3 1/3 0 0 0 0 1 0 > (A + gv ) = 1/5 1/5 1/5 1/5 1/5 0 0 0 0 1 0 0 0 1 0 130 Using this we can now 0 0 M= c/5 0 0 get our matrix M = c(A + gv> )> + (1 − c)ve> : > c/3 c/3 c/3 0 1 1 1 1 1 1 1 1 1 1 0 0 c 0 1−c 1 1 1 1 1 c/5 c/5 c/5 c/5 + 5 0 0 0 c 1 1 1 1 1 0 0 c 0 1 1 1 1 1 If we choose c = 0.85 and calculate the eigenvector to the dominant eigenvalue we end up with PageRank p = [0.0384, 0.0492, 0.0492, 0.4458, 0.4173]> . From the result we can see a couple of things, first we see that node one which have no other node linking to it have the lowest PageRank and node two and three which are only linked to by node one have equal and the second lowest PageRank. Node four have a very high PageRank since it’s linked to by many nodes and node five have high PageRank since it’s linked to by node four which have a very high PageRank. We see that the definition seems reasonable, however we still have a large problem to consider: how do we find the eigenvector to the dominant eigenvalue of such a huge matrix in a reasonable time? For reference 2008 − 2009 both Google and Bing reported indexing over 1.000.000.000.000 webpages which in turn means our matrix M would have about 1024 elements. Obviously something needs to be done. We note that the matrix M has been constructed such that it’s a non-negative matrix. And it also seems reasonable to assume it’s irreducible through our vector v. In fact unless some nodes have a zero-weight in v, M will be a positive matrix. We also have that the sum of every column of M equal to one, so we can say that M is actually a (column) stochastic matrix. This gives us a more intuitive way to describe pagerank looking at is from a Markov chain perspective. Consider a internet surfer surfing the web. The internet surfer starts surfing at a random page then with probability c it follows any of the links from the current page ending up in a new page or with probability 1 − c it starts over at a new random page. The Internet surfer then continue repeatedly either follows a link from the current page or starts over at a new random page. The pagerank of a node (page) can then be seen as the probability that the Internet surfer is at that page after a long time. While this by itself doesn’t make it easier, we note something about A. Since every homepage on average only links to a couple of other homepages, A is an extremely sparse matrix. While A have maybe 1024 elements only 1012 m, where m is the average number of links on a homepage is non-zero. If we have a sparse matrix, we could use the Power method to use the sparsity of the data rather than trying (and failing) using other methods such as the QR-method. We note that even though M is a full matrix, the multiplication xn+1 = Mxn in the Power method can be done fast by calculating Axn first. Take a couple of minutes to think through the following: 131 • How can we calculate xn+1 = Mxn by multiplying xn with cA> instead so that we can use that A is sparse? • How can the weightvector v be used in order to combat linkfarms, where people try to artificially increase ones PageRank? • The attribute ’nofollow’ can be used on links to not increase the PageRank of the target page. Can you think of any sites where this could be useful? • What influence does the value c have on the PageRank? For example what is important for a high PageRank if c is close to zero rather than 0.85? 132 Index adjacency matrix, 28 adjugate formula, 106 adjugate matrix, adj(A), 106 aperiodic, 41 associativity, 60 basis, 65 bijection, 84 bijective, 88 canonical form, 22 stochastic matrix, 52 Cayley-Hamilton theorem, 106 Ces´aro sum, 50 characteristic polynomial, 13 Cholesky factorization, 21 commutativity, 60 connected (graph), 31 connected component (graph), 31 continous time Markov chain, 56 covariance matrix, 22 cycle (graph), 30 degree matrix, 33 determinant, 10, 21 determinant function, 10 diagonal matrix, 6 diagonalizable, 13 distance, 68 distance matrix, 29 distributivity, 60 echelon form, 20 eigenvalue, 13, 39 eigenvector, 13, 39 error vector, 94 field, 60 Gauss elimination, 9 Hermitian conjugate, 18 Hermitian matrix, 18, 21 hermitian matrix, 8 hermitian transpose, 8 Hessenberg matrix, 17 hitting probability, 54 hitting time, 55 homogeneous coordinates, 86 Image, R(T), 88 imprimitive matrix, 42 injection, 84 injective, 88 inner product space, 66 interpolation, 95 irreducible matrix, 31, 39 Isomorphism, 91 Jordan block, 23 Jordan normal form, 23 Kernel, N(T), 88 Kroenecker product, 117 Laplacian matrix, 34 least-square sense, 93 linear combination, 62 linear space, 60 linear system, 8 linear transformation, 79 linearly dependent, 64 linearly independent, 64 LU-factorization, 21 LUP-factorization, 21 Markov chain, 46, 46 matrix decomposition, 16 matrix exponential function, eA , 113 matrix factorization, 16 matrix polynomial, 106 matrix power function,√An , 100 1 n matrix root function, A, A n , 103 Moore-Penrose pseudoinverse, 98 multilinear, 11 133 nilpotent matrix, 102 non-negative matrix, 39 normal equations, 95 Nullity, 90 Nullspace, N(T), 88 orthogonal, 69 orthogonal set, 69 orthonormal set, 69 path (graph), 30 permutation matrix, 21 Perron-Frobenius (theorem), 39 positive definite matrix, 19, 21 positive matrix, 40 Power method, 121 primitive matrix, 40 pseudoinverses, 98 QR-factorization, 22 QR-method, 126 rank, 90 rank factorization, 20 reduced row echelon form, 23 residual vector, 94 rotation matrix, 82 row echelon form, 9 Schwarz inequality, 68 similar matrix, 14 simple graph, 33 singular value decomposition, 24 span, 62 spanning set, 63 sparse matrix, 112 spectral factorization, 20 spectral radius, 13, 39 spectrum, 13 √ 1 square root of a matrix, A, A 2 , 103 stable, 119 standard matrix, 81 stationary distribution, 48 stochastic matrix, 46 irreducible and imprimitive, 50 134 irreducible and primitive, 49 reducible, 51 strongly connected (graph), 31 strongly connected component, 31 subspace, 61 surjection, 84 surjective, 88 SVD, 24 symmetric matrix, 8 trace, 6 transpose, 8 triangular matrix, 17, 21, 51 unitary matrix, 18, 22, 70 vector space, 60 weighted least square method, 97 References [1] K. Bryan and T. Leise. The $25, 000, 000, 000 eigenvector: the linear algebra behind google. SIAM Rev., 48(3):569–581, 2006. [2] M. Haugh. The monte carlo framework, examples from finance and generating correlated random variables. Course Notes, 2004. [3] J.R.Norris. Markov Chains. Cambridge university press, New York, 2009. [4] C. D. Meyer. Matrix Analysis and Applied Linear Algebra. SIAM, Boston, 2001. [5] C. Moler. Cleve’s corner: Chebfun, roots, October 2012. visited November 22 2012. [6] W. Nicholson. Elementary linear algebra with applications, second edition. PWS Publishing company, Boston, 1993. [7] P. G. D. J. L. Snell. Random walks and electric networks. Free Software Foundation, 2000. [8] S. Spanne. F¨ orel¨ asningar i Matristeori. KFS i Lund AB, Lund, 1994. [9] R. Tobias and L. Georg. Markovprocesser. Univ., Lund, 2000. [10] C. M. C. van Loan. Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Review, 45, 2003. 135

© Copyright 2018