23 Matrix Algorithms

1. Analysis of Algorithms Matrix algorithms Andres Mendez-Vazquez November 22, 2015 1 / 102

2. Outline 1 Introduction Basic Deﬁnitions Matrix Examples 2 Matrix Operations Introduction Matrix Multiplication The Inverse Determinants 3 Improving the Complexity of the Matrix Multiplication Back to Matrix Multiplication Strassen’s Algorithm The Algorithm How he did it? Complexity 4 Solving Systems of Linear Equations Introduction Lower Upper Decomposition Forward and Back Substitution Obtaining the Matrices Computing LU decomposition Computing LUP decomposition 5 Applications Inverting Matrices Least-squares Approximation 6 Exercises Some Exercises You Can Try!!! 2 / 102

4. Basic deﬁnitions A matrix is a rectangular array of numbers A = a11 a12 a13 a21 a22 a23 = 1 2 3 4 5 6 A transpose matrix is the matrix obtained by exchanging the rows and columns AT =    1 4 2 5 3 6    4 / 102

5. Basic deﬁnitions A matrix is a rectangular array of numbers A = a11 a12 a13 a21 a22 a23 = 1 2 3 4 5 6 A transpose matrix is the matrix obtained by exchanging the rows and columns AT =    1 4 2 5 3 6    4 / 102

7. Several cases of matrices Zero matrix        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        The diagonal matrix       a11 0 · · · 0 0 a22 · · · 0 ... ... ... ... 0 0 · · · ann       6 / 102

8. Several cases of matrices Zero matrix        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        The diagonal matrix       a11 0 · · · 0 0 a22 · · · 0 ... ... ... ... 0 0 · · · ann       6 / 102

9. Several cases of matrices Upper triangular matrix       a11 a12 · · · a1n 0 a22 · · · a2n ... ... ... ... 0 0 · · · ann       7 / 102

11. Operations on matrices They Deﬁne a Vectorial Space Matrix addition. Multiplication by scalar. The existence of zero. 9 / 102

15. Matrix Multiplication What is Matrix Multiplication? Given A, B matrices with dimensions n×n, the multiplication is deﬁned as C = AB cij = n k=1 aikbkj 11 / 102

16. Complexity and Algorithm Algorithm: Complexity Θ (n3 ) Square-Matrix-Multiply(A, B) 1 n = A.rows 2 let C be a new matrix n × n 3 for i = 1 to n 4 for j = 1 to n 5 C [i, j] = 0 6 for k = 1 to n 7 C [i, j] = C [i, j] + A [i, j] ∗ B [i, j] 8 return C 12 / 102

17. Matrix multiplication properties Properties of the Multiplication The Identity exist for a matrix A of m × n: ImA = AIn = A. The multiplication is associative: A(BC) = (AB)C. In addition, multiplication is distibutive A(B + C) = AB + AC (B + C)D = BD + CD 13 / 102

20. In addition Deﬁnition The inner product between vectors is deﬁed as xT y = n i=1 xiyi 14 / 102

22. Matrix inverses The inverse is deﬁned as the vector A−1 such that AA−1 = A−1 A = In Example 1 1 1 0 −1 = 0 1 1 −1 =⇒ 1 1 1 0 0 1 1 −1 = 1 · 0 + 1 · 1 1 · 1 − 1 · 1 1 · 0 + 1 · 0 1 · 1 + 0 · −1 = 1 0 0 1 Remark A matrix that is invertible is called non-singular. 16 / 102

25. Properties of an inverse Some properties are (BA)−1 = A−1B−1 A−1 T = AT −1 17 / 102

26. The Rank of A Rank of A A collection of vectors is x1, x2, ..., xn such that c1x1 + c2x2 + ... + cnxn = 0. The rank of a matrix is the number of linear independent rows. Theorem 1 A square matrix has full rank if and only if it is nonsingular. 18 / 102

27. The Rank of A Rank of A A collection of vectors is x1, x2, ..., xn such that c1x1 + c2x2 + ... + cnxn = 0. The rank of a matrix is the number of linear independent rows. Theorem 1 A square matrix has full rank if and only if it is nonsingular. 18 / 102

28. Other Theorems A null vector x is such that Ax = 0 Theorem 2: A matrix A has full column rank if and only if it does not have a null vector. Then, for squared matrices, we have Corollary 3: A square matrix A is singular if and only if it has a null vector. 19 / 102

29. Other Theorems A null vector x is such that Ax = 0 Theorem 2: A matrix A has full column rank if and only if it does not have a null vector. Then, for squared matrices, we have Corollary 3: A square matrix A is singular if and only if it has a null vector. 19 / 102

31. Determinants A determinant can be deﬁned recursively as follows det(A) =    a11 if n = 1 n j=1 (−1)1+ja1jdet(A[1j]) if n > 1 (1) Where (−1)i+jdet(A[ij]) is called a cofactor 21 / 102

32. Determinants A determinant can be deﬁned recursively as follows det(A) =    a11 if n = 1 n j=1 (−1)1+ja1jdet(A[1j]) if n > 1 (1) Where (−1)i+jdet(A[ij]) is called a cofactor 21 / 102

33. Theorems Theorem 4(determinant properties). The determinant of a square matrix A has the following properties: If any row or any column A is zero, then det(A) = 0. The determinant of A is multiplied by λ if the entries of any one row (or any one column) of A are all multiplied by λ. The determinant of A is unchanged if the entries in one row (respectively, column) are added to those in another row (respectively, column). The determinant of A equals the determinant of AT . The determinant of A is multiplied by −1 if any two rows (or any two columns) are exchanged. Theorem 5 An n × n matrix A is singular if and only if det(A) = 0. 22 / 102

40. Positive definite matrix Definition A positive definite matrix A is called positive definite if and only if xT Ax > 0 for all x = 0 Theorem 6 For any matrix A with full column rank, the matrix AT A is positive definite. 23 / 102

41. Positive definite matrix Definition A positive definite matrix A is called positive definite if and only if xT Ax > 0 for all x = 0 Theorem 6 For any matrix A with full column rank, the matrix AT A is positive definite. 23 / 102

43. Matrix Multiplication Problem description Given n × n matrices A,B and C: r s t u = a b c d e f g h Thus, you could compute r, s, t and u using recursion!!! r = ae + bg s = af + bh t = ce + dg u = cf + dh 25 / 102

44. Matrix Multiplication Problem description Given n × n matrices A,B and C: r s t u = a b c d e f g h Thus, you could compute r, s, t and u using recursion!!! r = ae + bg s = af + bh t = ce + dg u = cf + dh 25 / 102

45. Problem Complexity of previous approach T(n) = 8T n 2 + Θ(n2 ) Thus T(n) = Θ(n3 ) Therefore You need to use a diﬀerent type of products. 26 / 102

49. The Strassen’s Algorithm It is a divide and conquer algorithm Given A, B, C matrices with dimensions n × n, we recursively split the matrices such that we ﬁnish with 12 n 2 × n 2 sub matrices r s t u = a b c d e f g h Remember the Gauss Trick? Imagine the same for Matrix Multiplication. 28 / 102

50. The Strassen’s Algorithm It is a divide and conquer algorithm Given A, B, C matrices with dimensions n × n, we recursively split the matrices such that we ﬁnish with 12 n 2 × n 2 sub matrices r s t u = a b c d e f g h Remember the Gauss Trick? Imagine the same for Matrix Multiplication. 28 / 102

52. Algorithm Strassen’s Algorithm 1 Divide the input matrices A and B into n 2 × n 2 sub matrices. 2 Using Θ n2 scalar additions and subtractions, compute 14 matrices A1, B1, ..., A7, B7 each of which is n 2 × n 2 . 3 Recursively compute the seven matrices products Pi = AiBi for i = 1, 2, 3, ..., 7. 4 Compute the desired matrix r s t u by adding and or subtracting various combinations of the Pi matrices, using only Θ n2 scalar additions and subtractions 30 / 102

59. Strassen Observed that Trial and Error First , he generated Pi = AiBi = (αi1a + αi2b + αi3c + αi4d) · (βi1e + βi2f + βi3g + βi4h) Where αij, βij ∈ {−1, 0, 1} 32 / 102

60. Then r r = ae + bg = a b c d      +1 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0           e f g h      s s = af + bh = a b c d      +1 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 0           e f g h      33 / 102

63. Therefore t r = ce + dg = a b c d      0 0 0 0 0 0 0 0 +1 0 0 0 0 0 +1 0           e f g h      u u = cf + dh = a b c d      0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 +1           e f g h      34 / 102

64. Therefore t r = ce + dg = a b c d      0 0 0 0 0 0 0 0 +1 0 0 0 0 0 +1 0           e f g h      u u = cf + dh = a b c d      0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 +1           e f g h      34 / 102

65. Example Compute the s from P1 and P2 matrices Compute s = P1 + P2 Where P1 P1 = A1B1 = a (f − h) = af − ah = a b c d      0 +1 0 −1 0 0 0 0 0 0 0 0 0 0 0 0           e f g h      35 / 102

70. Example Compute the s from P1 and P2 matrices Where P2 P2 = A2B2 = (a + b) h = ah + bh = a b c d      0 0 0 +1 0 0 0 +1 0 0 0 0 0 0 0 0           e f g h      36 / 102

75. Complexity Because we are only computing 7 matrices T(n) = 7T n 2 + Θ n2 = Θ nlg 7 = O n2.81 . 38 / 102

76. Nevertheless We do not use Strassen’s because A constant factor hidden in the running of the algorithm is larger than the constant factor of the naive Θ n3 method. When matrices are sparse, there are faster methods. Strassen’s is not a numerically stable as the naive method. The sub matrices formed at the levels of the recursion consume space. 39 / 102

80. The Holy Grail of Matrix Multiplications O (n2 ) In a method by Virginia Vassilevska Williams (2012) Assistant Professor at Stanford The computational complexity of her method is ω < 2.3727 or O n2.3727 Better than Coppersmith and Winograd (1990) O n2.375477 Many Researchers Believe that Coppersmith, Winograd and Cohn et al. conjecture could lead to O n2 , contradicting a variant of the widely believed sun ﬂower conjecture of Erdos and Rado. 40 / 102

83. Exercises 28.1-3 28.1-5 28.1-8 28.1-9 28.2-2 28.2-5 41 / 102

85. In Many Fields From Optimization to Control We are required to solve systems of simultaneous equations. For Example For Polynomial Curve Fitting, we are given (x1, y1) , (x2, y2) , ..., (xn, yn) We want To ﬁnd a polynomial of degree n − 1 with structure p (x) = a0 + a1x + a2x2 + ... + an−1xn−1 43 / 102

88. Thus We can build a system of equations a0 + a1x1 + a2x2 1 + ... + an−1xn−1 1 = y1 a0 + a1x2 + a2x2 2 + ... + an−1xn−1 2 = y2 ... a0 + a1xn + a2x2 n + ... + an−1xn−1 n = yn We have n unknowns a0, a1, a2, ..., an−1 44 / 102

89. Thus We can build a system of equations a0 + a1x1 + a2x2 1 + ... + an−1xn−1 1 = y1 a0 + a1x2 + a2x2 2 + ... + an−1xn−1 2 = y2 ... a0 + a1xn + a2x2 n + ... + an−1xn−1 n = yn We have n unknowns a0, a1, a2, ..., an−1 44 / 102

90. Solving Systems of Linear Equations Proceed as follows We start with a set of linear equations with n unknowns: x1, x2, ..., xn    a11x1 + a12x2 + ... + a1nxn = b1 a21x1 + a22x2 + ... + a2nxn = b2 ... ... an1x1 + an2x2 + ... + annxn = bn Something Notable A set of values for x1, x2, ..., xn that satisfy all of the equations simultaneously is said to be a solution to these equations. In this section, we only treat the case in which there are exactly n equations in n unknowns. 45 / 102

94. Solving systems of linear equations continuation We can conveniently rewrite the equations as the matrix-vector equation:       a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... an1 an2 . . . ann             x1 x2 ... xn       =       b1 b2 ... bn       or, equivalently, letting A = (aij), x = (xj), and b = (bi), as Ax = b In this section, we shall be concerned predominantly with the case of which A is nonsingular, after all we want to invert A. 46 / 102

100. Overview of Lower Upper (LUP) Decomposition Intuition The idea behind LUP decomposition is to ﬁnd three n × n matrices L,U, and P such that: PA = LU where: L is a unit lower triangular matrix. U is an upper triangular matrix. P is a permutation matrix. Where We call matrices L,U, and P satisfying the above equation a LUP decomposition of the matrix A. 48 / 102

107. What is a Permutation Matrix Basically We represent the permutation P compactly by an array π[1..n]. For i = 1, 2, ..., n, the entry π[i] indicates that Piπ[i] = 1 and Pij = 0 for j = π[i]. Thus PA has aπ[i],j in row i and a column j. Pb has bπ[i] as its ith element. 49 / 102

108. What is a Permutation Matrix Basically We represent the permutation P compactly by an array π[1..n]. For i = 1, 2, ..., n, the entry π[i] indicates that Piπ[i] = 1 and Pij = 0 for j = π[i]. Thus PA has aπ[i],j in row i and a column j. Pb has bπ[i] as its ith element. 49 / 102

109. How can we use this in our advantage? Lock at this Ax = b =⇒ PAx = Pb (2) Therefore LUx = Pb (3) Now, if we make Ux = y Ly = Pb (4) 50 / 102

112. Thus We ﬁrst obtain y Then, we obtain x. 51 / 102

114. Forward and Back Substitution Forward substitution Forward substitution can solve the lower triangular system Ly = Pb in Θ(n2) time, given L, P and b. Then Since L is unit lower triangular, equation Ly = Pb can be rewritten as: y1 = bπ[1] l21y1 + y2 = bπ[2] l31y1 + l32 + y3 = bπ[3] ... ln1y1 + ln2y2 + ln3y3 + ... + yn = bπ[n] 53 / 102

115. Forward and Back Substitution Forward substitution Forward substitution can solve the lower triangular system Ly = Pb in Θ(n2) time, given L, P and b. Then Since L is unit lower triangular, equation Ly = Pb can be rewritten as: y1 = bπ[1] l21y1 + y2 = bπ[2] l31y1 + l32 + y3 = bπ[3] ... ln1y1 + ln2y2 + ln3y3 + ... + yn = bπ[n] 53 / 102

116. Forward and Back Substitution Back substitution Back substitution is similar to forward substitution. Like forward substitution, this process runs in Θ(n2) time. Since U is upper-triangular, we can rewrite the system Ux = y as u11x1 + u12x2 + ... + u1n−2xn−2 + u1n−1xn−1 + u1nxn = y1 u22x2 + ... + u2n−2xn−2 + u2n−1xn−1 + u2nxn = y2 ... un−2n−2xn−2 + un−2n−1xn−1 + un−2nxn = yn−2 un−1n−1xn−1 + un−1nxn = yn−1 unnxn = yn 54 / 102

117. Example We have Ax =    1 2 0 3 4 4 5 6 3    x =    3 7 8    = b 55 / 102

118. Example The L, U and P matrix L =    1 0 0 0.2 1 0 0.6 0.5 1    , U =    5 6 3 0 0.8 −0.6 0 0 2.5    , P =    0 0 1 1 0 0 0 1 0    56 / 102

119. Example Using forward substitution, Ly = Pb for y Ly =    1 0 0 0.2 1 0 0.6 0.5 1    y =    0 0 1 1 0 0 0 1 0       3 7 8    = Pb 57 / 102

120. Example Using forward substitution, we get y y =    8 1.4 1.5    58 / 102

121. Example Now, we use the back substitution, Ux = y for x Ux =    5 6 3 0 0.8 −0.6 0 0 2.5       x1 x2 x3    =    8 1.4 1.5    , 59 / 102

122. Example Finally, we get x =    −1.4 2.2 0.6    60 / 102

123. Forward and Back Substitution Given P, L, U, and b, the procedure LUP- SOLVE solves for x by combining forward and back substitution LUP-SOLVE(L, U, π, b) 1 n = L.rows 2 Let x be a new vector of length n 3 for i = 1 to n 4 yi = bπ[i] − i−1 j=1 lijyj 5 for i = n downto 1 6 xi = yi− n j=i+1 uijxj uii 7 return x Complexity The running time is Θ(n2). 61 / 102

128. Ok, if we have the L,U and P!!! Thus We need to ﬁnd those matrices How, we do it? We are going to use something called the Gaussian Elimination. 63 / 102

129. Ok, if we have the L,U and P!!! Thus We need to ﬁnd those matrices How, we do it? We are going to use something called the Gaussian Elimination. 63 / 102

130. For this We assume that A is a n × n Such that A is not singular We use a process known as Gaussian elimination to create LU decomposition This algorithm is recursive in nature. Properties Clearly if n = 1, we are done for L = I1 and U = A. 64 / 102

134. Computing LU decomposition For n > 1, we break A into four parts A =         a11 a12 · · · a1n a21 a22 · · · a2n ... ... ... ... an1 an2 · · · ann         =   a11 wT v A   (5) 66 / 102

135. Where We have v is a column (n − 1) −vector. wT is a row (n − 1) −vector. A is an (n − 1) × (n − 1). 67 / 102

139. Computing a LU decomposition Thus, we can do the following A = a11 wT v A = 1 0 v a11 In−1       a11 wT 0 A − vwT a11 Schur Complement       = 1 0 v a11 In−1 a11 wT 0 L U = 1 0 v a11 L a11 wT 0 U = LU 68 / 102

144. Computing a LU decomposition Pseudo-Code running in Θ (n3 ) LU-Decomposition(A) 1 n = A.rows 2 Let L and U be new n × n matrices 3 Initialize U with 0’s below the diagonal 4 Initialize L with 1’s on the diagonal and 0’s above the diagonal. 5 for k = 1 to n 6 ukk = akk 7 for i = k + 1 to n 8 lik = aik ukk lik holds vi 9 uki = aki uki holds wT i 10 for i = k + 1 to n 11 for j = k + 1 to n 12 aij = aij − likukj 13 return L and U 69 / 102

153. Example Here, we have this example 2 3 1 5 6 13 5 19 2 19 10 23 4 10 11 31 ⇒ 13 5 19 19 10 23 10 11 31 − 1 2 6 2 4 3 1 5 = 13 5 19 19 10 23 10 11 31 − 1 2 18 6 30 6 2 10 12 4 20 ⇒ 2 3 1 5 3 4 2 4 1 16 9 18 2 4 9 21 ⇒ 9 18 9 11 − 1 4 16 4 2 4 = 9 18 9 11 − 1 4 32 64 8 16 = 9 18 9 11 − 8 16 2 4 ⇒ 2 3 1 5 3 4 2 4 1 4 1 2 2 1 7 17 ⇒ 2 3 1 5 3 4 2 4 1 4 1 2 2 1 7 3 70 / 102

159. Thus We get the following decomposition      2 3 1 5 6 13 5 19 2 19 10 23 4 10 11 31      =      1 0 0 0 3 1 0 0 1 4 1 0 2 1 7 1           2 3 1 5 0 4 2 4 0 0 1 2 0 0 0 3      71 / 102

161. Observations Something Notable The elements by which we divide during LU decomposition are called pivots. They occupy the diagonal elements of the matrix U. Why the permutation P It allows us to avoid dividing by 0. 73 / 102

164. Thus, What do we want? We want P, L and U PA = LU However, we move a non-zero element, ak1 From somewhere in the ﬁrst column to the (1, 1) position of the matrix. In addition ak1 as the element in the ﬁrst column with the greatest absolute value. 74 / 102

167. Exchange Rows Thus We exchange row 1 with row k, or multiplying A by a permutation matrix Q on the left QA = ak1 wT v A With v = (a21, a31, ..., an1)T with a11 replaces ak1. wT = (ak2, ak3, ..., akn). A is a (n − 1) × (n − 1) 75 / 102

168. Exchange Rows Thus We exchange row 1 with row k, or multiplying A by a permutation matrix Q on the left QA = ak1 wT v A With v = (a21, a31, ..., an1)T with a11 replaces ak1. wT = (ak2, ak3, ..., akn). A is a (n − 1) × (n − 1) 75 / 102

169. Now, ak1 = 0 We have then QA = ak1 wT v A = 1 0 v ak1 In−1 ak1 wT 0 A − vwT ak1 76 / 102

170. Now, ak1 = 0 We have then QA = ak1 wT v A = 1 0 v ak1 In−1 ak1 wT 0 A − vwT ak1 76 / 102

171. Important Something Notable if A is nonsingular, then the Schur complement A − vwT ak1 is nonsingular, too. Now, we can ﬁnd recursively an LUP decomposition for it P A − vwT ak1 = L U Then, we deﬁne a new permutation matrix P = 1 0 0 P Q 77 / 102

174. Thus We have PA = 1 0 0 P QA = 1 0 0 P 1 0 v ak1 In−1 ak1 wT 0 A − vwT ak1 = 1 0 P v ak1 P ak1 wT 0 A − vwT ak1 = 1 0 P v ak1 In−1 ak1 wT 0 P A − vwT ak1 = 1 0 P v ak1 In−1 ak1 wT 0 L U = 1 0 P v ak1 L ak1 wT 0 U = LU 78 / 102

180. Computing a LUP decomposition Algorithm LUP-Decomposition(A) 1. n = A.rows 2. Let π [1..n] new array 3. for i = 1 to n 4. π [i] = i 5. for k = 1 to n 6. p = 0 7. for i = k to n 8. if |aik| > p 9. p = |aik| 10. k = i 11. if p == 0 12. error “Singular Matrix” 13. Exchange π [k] ←→ π [k ] 14. for i = 1 to n 15. Exchange aki ←→ ak i 16. for i = k + 1 to n 17. aik = aik akk 18. for j = k + 1 to n 19. aij = aij − aikakj 79 / 102

189. Computing a LUP decomposition Example 1 2 0 2 0.6 2 3 3 4 -2 3 5 5 4 2 4 -1 -2 3.4 -1 =⇒ 3 5 5 4 2 2 3 3 4 -2 1 2 0 2 0.6 4 -1 -2 3.4 -1 =⇒ 3 5 5 4 2 2 3 3 4 -2 1 2 0 2 0.6 4 -1 -2 3.4 -1 =⇒ 3 5 5 4 2 2 0.6 0 1.6 -3.2 1 0.4 -2 0.4 -0.2 4 -1 -1 4.2 -0.6 =⇒ 3 5 5 4 2 2 0.6 0 1.6 -3.2 1 0.4 -2 0.4 -0.2 4 -1 -1 4.2 -0.6 =⇒ 3 5 5 4 2 2 0.6 0 1.6 -3.2 1 0.4 -2 0.4 -0.2 4 -1 -1 4.2 -0.6 =⇒ 3 5 5 4 2 1 0.4 -2 0.4 -0.2 2 0.6 0 1.6 -3.2 4 -1 -1 4.2 -0.6 =⇒ 3 5 5 4 2 1 0.4 -2 0.4 -0.2 2 0.6 0 1.6 -3.2 4 -1 -1 4.2 -0.6 =⇒ 3 5 5 4 2 1 0.4 -2 0.4 -0.2 2 0.6 0 1.6 -3.2 4 -1 0.5 4 -0.5 80 / 102

198. Finally, you get The Permutation and Decomposition      0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0      P      2 0 2 0.6 3 3 4 −2 5 5 4 2 −1 −2 3.4 −1      A = ...      1 0 0 0 0.4 1 0 0 −0.2 0.5 1 0 0.6 0 0.4 1      L      5 5 4 2 0 −2 0.4 −0.2 0 0 4 −0.5 0 0 0 −3      U 81 / 102

199. Symmetric positive-definite matrices Lemma 28.9 Any symmetric positive-definite matrix is nonsingular. Lemma 28.10 If A is a symmetric positive-definite matrix, then every leading submatrix of A is symmetric and positive-definite. 82 / 102

200. Symmetric positive-definite matrices Lemma 28.9 Any symmetric positive-definite matrix is nonsingular. Lemma 28.10 If A is a symmetric positive-definite matrix, then every leading submatrix of A is symmetric and positive-definite. 82 / 102

201. Symmetric positive-definite matrices Definition: Schur complement Let A be a symmetric positive-definite matrix, and let Ak be a leading k × k submatrix of A. Partition A as: A = Ak BT B C Then, the Schur complement of A with respect to Ak is defined to be S = C − BA−1 k BT 83 / 102

205. Symmetric positive-definite matrices Lemma 28.11 (Schur complement lemma) If A is a symmetric positive-definite matrix and Ak is a leading k × k submatrix of A, then the Schur complement of A with respect to Ak is symmetric and positive-definite. Corollary 28.12 LU decomposition of a symmetric positive-definite matrix never causes a division by 0. 84 / 102

206. Symmetric positive-definite matrices Lemma 28.11 (Schur complement lemma) If A is a symmetric positive-definite matrix and Ak is a leading k × k submatrix of A, then the Schur complement of A with respect to Ak is symmetric and positive-definite. Corollary 28.12 LU decomposition of a symmetric positive-definite matrix never causes a division by 0. 84 / 102

208. Inverting matrices LUP decomposition can be used to compute a matrix inverse The computation of a matrix inverse can be speed up using techniques such as Strassen’s algorithm for matrix multiplication. 86 / 102

209. Computing a matrix inverse from a LUP decomposition Proceed as follows The equation AX = In can be viewed as a set of n distinct equations of the form Axi = ei, for i = 1, ..., n. We have a LUP decomposition of a matrix A in the form of three matrices L,U, and P such that PA = LU. Then we use the backward-forward to solve AXi = ei. 87 / 102

210. Complexity First We can compute each Xi in time Θ n2 . Thus, X can be computed in time Θ n3 . LUP decomposition is computed in time Θ n3 . Finally We can compute A−1 of a matrix A in time Θ n3 . 88 / 102

211. Complexity First We can compute each Xi in time Θ n2 . Thus, X can be computed in time Θ n3 . LUP decomposition is computed in time Θ n3 . Finally We can compute A−1 of a matrix A in time Θ n3 . 88 / 102

212. Matrix multiplication and matrix inversion Theorem 28.7 If we can invert an n × n matrix in time I(n), where I(n) = Ω(n2) and I(n) satisﬁes the regularity condition I(3n) = O(I(n)), then we can multiply two n × n matrices in time O(I(n)). 89 / 102

213. Matrix multiplication and matrix inversion Theorem 28.8 If we can multiply two n × n real matrices in time M(n), where M(n) = Ω(n2) and M(n) = O(M(n + k)) for any k in range 0 ≤ k ≤ n and M(n 2 ) ≤ cM(n) for some constant c < 1 2. Then we can compute the inverse of any real nonsingular n × n matrix in time O(M(n)). 90 / 102

215. Least-squares Approximation Fitting curves to given sets of data points is an important application of symmetric positive-deﬁnite matrices. Given (x1, y1), (x2, y2), ..., (xm, ym) where the yi are known to be subject to measurement errors. We would like to determine a function F(x) such that: yi = F(xi) + ηi for i = 1, 2, ..., m 92 / 102

216. Least-squares Approximation Continuation The form of the function F depends on the problem at hand. F(x) = n j=1 cjfj(x) A common choice is fj(x) = xj−1, which means that F(x) = c1 + c2x + c3x2 + ... + cnxn−1 is a polynomial of degree n − 1 in x. 93 / 102

217. Least-squares Approximation Continuation Let A =       f1(x1) f2(x1) . . . fn(x1) f1(x2) f2(x2) . . . fn(x2) ... ... ... ... f1(xm) f2(xm) . . . fn(xm)       denote the matrix of values of the basis functions at the given points; that is, aij = fj(xi). Let c = (ck) denote the desired size-n vector of coeﬃcients. Then, A =       f1(x1) f2(x1) . . . fn(x1) f1(x2) f2(x2) . . . fn(x2) ... ... ... ... f1(xm) f2(xm) . . . fn(xm)             c1 c2 ... cn       =       F(x1) F(x2) ... F(xm)       94 / 102

218. Least-squares Approximation Then Thus, η = Ac − y is the size of approximation errors. To minimize approximation errors, we choose to minimize the norm of the error vector , which gives us a least-squares solution. ||η||2 = ||Ac − y||2 = m i=1 n j=1 aijcj − yi 2 Thus We can minimize ||η|| by diﬀerentiating ||η|| with respect to each ck and then setting the result to 0: d||η||2 dck = m i=1 2 n j=1 aijcj − yi aik = 0 95 / 102

219. Least-squares Approximation Then Thus, η = Ac − y is the size of approximation errors. To minimize approximation errors, we choose to minimize the norm of the error vector , which gives us a least-squares solution. ||η||2 = ||Ac − y||2 = m i=1 n j=1 aijcj − yi 2 Thus We can minimize ||η|| by diﬀerentiating ||η|| with respect to each ck and then setting the result to 0: d||η||2 dck = m i=1 2 n j=1 aijcj − yi aik = 0 95 / 102

220. Least-squares Approximation We can put all derivatives The n equation for k = 1, 2, ..., n (Ac − y)T A = 0 or equivalently to AT (Ac − y) = 0 which implies AT Ac = AT y 96 / 102

221. Least-squares Approximation Continuation The AT A is symmetric: If A has full column rank, then AT A is positive- deﬁnite as well. Hence, (AT A)−1 exists, and the solution to equation AT Ac = AT y is c = ((AT A)−1AT )y = A+y where the matrix A+ = ((AT A)−1AT ) is called the pseudoinverse of the matrix A. 97 / 102

222. Least-Square Approximation Continuation As an example of producing a least-squares ﬁt, suppose that we have 5 data points (-1,2), (1,1),(2,1),(3,0),(5,3), shown as black dots in following ﬁgure 98 / 102

223. Least-squares Approximation Continuation We start with the matrix of basis-function values A =        1 x1 x2 1 1 x2 x2 2 1 x3 x2 3 1 x4 x2 4 1 x5 x2 5        =        1 −1 1 1 1 1 1 2 4 1 3 9 1 5 25        whose pseudoinverse is A+ =    0.500 0.300 0.200 0.100 −0.100 −0.388 0.093 0.190 0.193 −0.088 0.060 −0.036 −0.048 −0.036 0.060    99 / 102

224. Matrix multiplication and matrix inversion Continuation Multiplying y by A+ , we obtain the coeﬃcient vector c =    1.200 −0.757 0.214    which corresponds to the quadratic polynomial F(x) = 1.200 − 0.757x + 0.214x2 100 / 102

226. Exercises From Cormen’s book solve 34.5-1 34.5-2 34.5-3 34.5-4 34.5-5 34.5-7 34.5-8 102 / 102

23 Matrix Algorithms

Recommended

More Related Content

What's hot (20)

Viewers also liked (20)

Similar to 23 Matrix Algorithms (20)

More from Andres Mendez-Vazquez (20)

Recently uploaded (20)

23 Matrix Algorithms