Chapter 5. The Inverse; Numerical Methods

Vector Spaces in Physics8/6/2015Chapter 5. The Inverse; Numerical MethodsIn the Chapter 3 we discussed the solution of systems of simultaneous linear algebraicequations which could be written in the form (5-1)Ax Cusing Cramer's rule. There is another, more elegant way of solving this equation, usingthe inverse matrix. In this chapter we will define the inverse matrix and give anexpression related to Cramer's rule for calculating the elements of the inverse matrix. Wewill then discuss another approach, that of Gauss-Jordan elimination, for solvingsimultaneous linear equations and for calculating the inverse matrix. We will discuss therelative efficiencies of the two algorithms for numerical inversion of large matrices.A. The inverse of a square matrix.Definition of the inverse. The inverse of a scalar number c is another scalar, say d, suchthat the product of the two is equal to 1: c d 1. For instance, the inverse of the number5 is the number 0.2 . We have defined multiplication of one matrix by another in a wayvery analogous to multiplication of one scalar by another. We will therefore make thefollowing definition.Definition: For a given square matrix A , the matrix B is said to be the inverseof A if(5-2)BA AB I 1We then write B A .Notice that we have not guaranteed that the inverse of a given matrix exists. In fact,many matrices do not have an inverse. We shall see below that the condition for a squarematrix A to have an inverse is that its determinant not be equal to zero.Use of the inverse to solve matrix equations. Now consider the matrix equation justgiven, (5-1)Ax C 1We can solve this equation by multiplying on both sides of the equation by A : A 1 A x A 1C ; I x A 1C ;(5-3) x A 1C.Thus, knowing the inverse of the matrix A lets us immediately write down the solution x to equation (5-1).As an example, let us consider a specific example, where A is a 2x2 matrix.5-1

Vector Spaces in Physics8/6/2015 A x C; 1 2 x1 4 1 1 x 2 5 (5-4) If we knew the inverse of A , we could immediately calculate C A x . In this simplecase, we can guess the inverse matrix. We write out the condition for the inverse,A A 1 I ; 1 2 * * ? ? I ij 1 1 * * ? ? As a first guess we try to make I12 come out to zero; one possibility isA A 1 I ;(5-5) 1 2 * 2 ? 0 I ij 11* 1? Now we arrange for I21 to be zero:A A 1 I ;(5-6)(5-7) 1 2 1 2 ? 0 I ij 1 1 1 1 0 ? If we now look at the diagonal elements of I , they come out to be I11 3 and I22 -3.We can fix this up by changing the sign of the (1,2) and (2,2) elements of the inverse, andby multiplying it by 1/3. So we haveA A 1 I ; 1 2 1 1 2 1 0 I ij ; 1 1 3 1 1 0 1 1 1 2 A 1 3 1 1 Now that we have the inverse matrix, we can calculate the values x1 and x2: x A 1C1 1 2 4 3 1 1 5 1 6 3 9 (5-7)(5-8) 2 3 So, the solution to the two simultaneous linear equations is supposed to be x1 -2, x2 3.We will write out the two equations in long form and substitute in. 1 2 x1 4 ; 1 1 x 2 5 5-2

Vector Spaces in Physics8/6/2015 x1 2 x 2 4 ; x x 52 1 2 2 * 3 4 ; ( 2) 3 5 4 4 . 5 5 It checks out!The inverse matrix by the method of cofactors. Guessing the inverse has worked for a2x2 matrix - but it gets harder for larger matrices. There is a way to calculate the inverseusing cofactors, which we state here without proof:cof ( A) A 1 ij Aji(5-9)j i 1 M ji ( A) A(Here the minor Mpq(A) is the determinant of the matrix obtained by removing the p-throw and q-th column from the matrix A.)Note that you cannot calculate the inverse of a matrix using equation (5-9) if the matrixis singular (that is, if its determinant is zero). This is a general rule for square matrices:A 0 inverse does not existExample: Find the inverse of the matrix 1 2 (5-10)A 1 1 Here are the calculations of the four elements of A 1 . First calculate the determinant:1 2(5-11)A 1 ( 2) 3 1 1Then the matrix elements:5-3

Vector Spaces in Physics8/6/2015 A cof11 ( A ) A cof 21 ( A ) cof12 ( A ) cof 22 ( A ) 111 112 A 121 A 122A1 ;3AAAA1 1 1 A11 1 2 1 A12A1 2 1 A21 A2 2 1 A11 A2 ;3(5-12)1 ;31 ;3So,1 1 2 A 1 3 1 1 (5-13)1 1 2 1 2 1 3 0 I;A A 1 3 1 1 1 1 3 0 3 1 1 2 1 2 1 3 0 IA 1 A 3 1 1 1 1 3 0 3 (5-14)Check that this inverse works:Example: Calculate the inverse of the following 3x3 matrix using the method ofcofactors: 2 4 3 A 1 3 2 (5-15) 1 3 1 Solution: This is getting too long-winded. We will just do two representative elementsof A 1 .2 4 3A 1 3 2 2(3 6) 4(1 2) 3(3 3) 2;1 3 1 A 111 A 123 cof 11 ( A )Acof 32 ( A )A 1 1 1 3 23 1A 1 3 2 AB. Time required for numerical calculations.5-42 31 2 3 3 ; 2 2 1 1 ; 2 2(5-16)

Vector Spaces in Physics8/6/2015Let’s estimate the computer time required to invert a matrix by the method of cofactors.The quantity of interest is the number of floating-point operations required to carry outthe inverse. The inverse of a nxn matrix involves calculating n2 cofactors, each of themrequiring the calculation of the determinant of an (n-1)x(n-1) matrix. So we need toknow the number of operations involved in calculating a determinant. Let's start with a2x2 determinant. There are two multiplications, and an addition to add the two terms.n 2 gives 3 FLOPs. (FLOP Floating-Point Operation.) To do a 3x3 determinant, thethree elements in the top row are each multiplied by a 2x2 determinant and addedtogether: 3x(3 FLOPs) 2 FLOPs for addition; n 3 requires 3x3 2 FLOPs. Now wecan proceed more or less by induction. It is pretty clear that the determinant of a 4x4matrix requires 4 calculations of a 3x3 determinant: -- 4x3x3 FLOPs. And for a 5x5determinant, 5x4x3x3 operations. It is a pretty good approximation to say the following:No. of operations for nxn determinant n!(5-17)2This means that calculating the inverse by the cofactor method (n cofactors) requiresn2n! FLOPs.A fast PC can today do about 10 GigaFLOPs/sec. This leads to the table given belowshowing the execution time to invert matrices of increasing 324forfor inversetime (sec)for inversetime (sec)determinant (method of cofactors)(PC)(Gauss-Jordan)(PC)(n!)(n 2*n!)(4n 082562.56E-0812030000.0000003500 0.0000000572025920 0.0000025928648.64E-085040246960 0.0000246961372 1.372E-07403202580480 0.0002580482048 2.048E-0736288029393280 0.0029393282916 2.916E-0736288003628800000.0362884000 0.00000043991680048299328000.482993285324 5.324E-07479001600689762304006.897623046912 6.912E-0762270208001.05237E 12 105.23665158788 8.788E-078.7178E 101.70869E 13 1708.69450810976 1.0976E-061.3077E 122.94227E 14 29422.6732813500 0.000001352.0923E 135.35623E 15 535623.421116384 1.6384E-063.5569E 141.02794E 17 10279366.6719652 1.9652E-066.4024E 152.07437E 18 207436908.123328 2.3328E-061.2165E 174.39139E 19439138812527436 2.7436E-062.4329E 189.73161E 20 9731608032732000 0.00000325.1091E 192.25311E 22 2.25311E 1237044 3.7044E-061.124E 215.44016E 23 5.44016E 1342592 4.2592E-062.5852E 221.36757E 25 1.36757E 1548668 4.8668E-066.2045E 233.57378E 26 3.57378E 1655296 5.5296E-06Table 5-1. Floating-point operations required for calculation of n x n determinants andinverses of n x n matrices, and computer time required for the matrix inversion. Results5-5

Vector Spaces in Physics8/6/2015are given for two different numerical methods. (As a useful conversion number, thenumber of seconds in a year is about 3.14 x 107.)It can be seen from the table that the inversion of a 24x24 matrix could take a time on afast computer about equal to the age of the Universe. This suggests that a moreeconomical algorithm is desirable for inverting large matrices!Teasing Mathematica: Try this calculation of a determinant.n 500m Table[Random[],{n},{n}];Det[m]Does this suggest that the algorithm used for Table 5-1 is not the fastest known?C. The Gauss-Jordan method for solving simultaneous linear equations.There is a method for solving simultaneous linear equations that avoids the determinantsrequired in Cramer's method, and which takes many fewer operations for large matrices.We will illustrate this method for two simultaneous linear equations, and then for three.Consider the 2x2 matrix equation solved above, A x C;(5-4) 1 2 x1 4 1 1 x 2 5 This corresponds to the two linear equationsx1 2 x2 4(5-18) x1 x2 5A standard approach to such equations would be to add or subtract a multiple of oneequation from another to eliminate one variable from one of the equations. If we add thefirst equation to the second, we get x1 2 x 2 4 addt eq. (1) to eq.(2) x1 2 x 2 4 (5-19) x x 50 3x 9122 Now we eliminate x2 from the top equation, by subtracting 2/3 x the bottom equation: x1 2 x 2 4 add eq. (1) to eq.(2) x1 2 x2 4 x1 x 2 5 0 3x 2 9 (5-20)x 0 2 eq. (2) from eq.(1) subtract (2/3)x 1 0 3x 2 9 And finally, multiply the second equation by 1/3: x1 2 x 2 4 add eq. (1) to eq.(2) x1 2 x2 4 x1 x 2 5 0 3x 2 9 (5-21) x1 0 2 multiply eq. (2) by 1/3 x1 0 2 subtract (2/3)x eq. (2) from eq.(1) 0 3x 2 9 0 x2 3 So we have found that x1 -2 and x2 3, as determined earlier in the chapter using theinverse.5-6

Vector Spaces in Physics8/6/2015Note that the same operations could have been carried out using just the coefficients ofthe equations, and omitting x1 and x2, as follows. The assembly of the coefficients of x1and x2 and the constants on the right of the equation is refered to as the augmentedmatrix. 1 2 4 add eq. (1) to eq.(2) 1 2 4 1 1 5 039 (5-22) 1 0 2 multiply eq. (2) by 1/3 1 0 2 subtract (2.3)x eq. (2) from eq.(1) 039013 The results for x1 and x2 appear in the column to the right.Example: Use the Gauss-Jordan method to solve the system of linear equationsrepresented by A x C; 2 4 3 x1 1 1 3 2 x 2 1 1 3 1 x 4 3 (5-23)Solution: We set up the augmented matrix, and then set about making the matrix part ofit diagonal, with ones on the diagonal. This is done in the following systematic fashion.First use the first equation to eliminate A21 and A22. Next use the second equation toeliminate A32. 2 4 3 1 2 431 subtract 1/2 x (1) from (2) and from (3) 1 3 2 1 0 1 1 / 2 1 / 2 1 3 1 4 0 1 1/ 2 7 / 2 (5-24) 2 4 31 (3) subtract (2) from 0 1 1 / 2 1 / 2 0 0 13 Next we work upwards, using equation (3) to eliminate A23 and A13. After that, equation(2) is used to eliminate A12. At this point the matrix is diagonal. the final step is tomultiply equations (1) and (3) by a constant which makes the diagonal elements of Abecome unity.5-7

Vector Spaces in Physics8/6/2015 2 4 31 add 1./2 * (3) to (2), 2 4 0 10 add 3 * (3) to (1) 011/21/2 0102 0 0 1 3 0 0 1 3 2 0 0 2 multiply (1) by 1.2 1 0 0subtract 4 x (2) from (1)-1 0 1 0 2 and (3) by 0 1 0 0 0 1 3 0 0 1 1 (5-25) 2 3 The solution for the unknown x's is thus x1 1, x2 2, x3 -3.SUMMARY: Work your way through the matrix, zeroing the off-diagonal elements, INTHE ORDER SHOWN BELOW, zeroing ONE, then TWO, then THREE, etc. If you tryto invent your own scheme of adding and subtracting rows, it may or may not work.SIXFIVE . ONE FOUR TWO THREE D. The Gauss-Jordan method for inverting a matrix.There is a very similar procedure which leads directly to calculating the inverse of asquare matrix. Suppose that B is the inverse of A . ThenAB I ; A11 A12 A13 B11 B12 B13 1 0 0 A21 A22 A23 B21 B22 B23 0 1 0 A 31 A32 A33 B31 B32 B33 0 0 1 This can be thought of us three sets of three simultaneous linear equations: A11 A12 A13 B11 1 A21 A22 A23 B21 0 , A 31 A32 A33 B31 0 A11 A21 A 31A12A22A32A13 B12 0 A23 B22 1 ,A33 B32 0 (5-26)(5-27) A11 A12 A13 B13 0 A21 A22 A23 B23 0 , A 31 A32 A33 B33 1 These three sets of equations can be solved simultaneously, using a larger augmentedequation, as follows:5-8