Numerical Analysis I
Numerical Analysis I MA 580
Popular in Course
Popular in Mathematics (M)
This 33 page Class Notes was uploaded by Braeden Lind on Thursday October 15, 2015. The Class Notes belongs to MA 580 at North Carolina State University taught by Staff in Fall. Since its upload, it has received 8 views. For similar materials see /class/223723/ma-580-north-carolina-state-university in Mathematics (M) at North Carolina State University.
Reviews for Numerical Analysis I
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 10/15/15
Chapter 4 Iterative methods for solving linear system of equations The Gaussian elimination method for solving A7 I is quite ef cient if the size of A is small to medium in reference the available computers and dense matrices most of entries of the matrix are nonzero numbers But for several reasons sometimes an iterative method may be more ef cient I For sparse matrices the Gaussian elimination method may destroy the structure of the matrix and cause quot ll2391 s see for example 1 1 1 1 1 1 1 1 1 1 1 1 2 1 0 0 0 0 0 1 2 2 2 2 3 0 1 0 0 0 0 3 2 3 3 3 4 0 0 1 0 0 0 4 4 3 4 4 5 0 0 0 1 0 0 5 5 5 1 5 6 0 0 0 0 1 0 6 6 6 6 5 Obviously the case discussed above can be generalized to a general n by n matrix with the same structure I Large sparse matrices for which we may not be able to store all the entries of the matrix Below we show an example in two dimensions If If 56 Z Li 405 The central nite difference method with ve point stencil for Pois son equation Consider the Poisson equation quot11 quotmy 2 3 3 E Q 2 1th X C d 4 01 ua mm uoa y Dirichlet BC 402 If f E L2Z then the solution exists and it is unique Analytic solution is rarely available Now we discuss how to use the finite difference equation to solve the Poisson equation I Step 1 Generate a grid A uniform Cartesian grid can be used 3 xi a 17 139 012 4m ti m a 403 d c yjztjlty j012 n y 7 404 I We want to find an approximate solution Uij to the exact solution at all the grid points mi 93 where Mara 93 is unknown So there are m 1n 1 unknown for Dirichlet boundary condition I Step 2 Substitute the partial derivatives with a finite difference formula in terms of the function values at grid points to get quoty3719 EMU391 93 quotJn1 93 quotJ391 934 214miwy Ui1Jj1 x2 hy2 f7 jI Tij i1g44m 1 j1 n 1 where fig 2 fyi 93 The local truncation error satisfies x2 6 12 641 T w 7 7 40r U 12 6x4 12 694 a Define h max hi Ly 406 The finite difference discretization is consistent if l T 0 407 3ng H H Therefore the discretization is consistent and second order accurate If we remove the error term in the equation above and replace the exact solution 13 93 with the approximate solution Uij which is the solution of the linear system of equations Ui 1j UiLj 471 UM lt 2 2 39 7 4 zzy hyV 0ny ivy U1 f1 0 8 Numerical Analysis I 57 The finite difference scheme at a grid point soapy involves five grid points east north west south and the center The center is called the master grid point I Solve the linear system of equations to get an approximate solution at grid points how I Error analyr implementation visualization etc 406 Matrixvector form of the nite difference equations Generally if one wants to use a direct method such Gaussian elimination method or sparse matrix techniques then one needs to find out the matrix structure If one use an iterative method such Jacobi Gauss Seidel SORw methods then it may be not necessarily to have the matrix and vector form In the matrix vector form AU F the unknown is a one dimensional array For the two dimensional Poisson equations the unknowns U73 are a two dimensional array Therefore we need to order it to get a one dimensional array We also need to order the finite difference equations It is common practice that we use the same ordering for the equations and for the unknowns There are two commonly used ordering One is called the natural ordering that fits sequential computers The other one is called the red and black ordering that fits parallel computers Figure 41 The natural ordering left and the redblack ordering right The natural row ordering In the natural row ordering we order the unknowns equations rowwise therefore the kth equation corresponding to i j with the following relation k im 1j 1 139 12 m 1 339 12 n 1 409 We use the following example to verify the matrixvector form of the finite difference equations Assume that ha Ly It m n 4 so we will have nine equations and nine unknowns The coef cient matrix is 9 by 9 To write down the matrixvector form we use a onedimensional array x to express the unknown Uij 1 U11 302 U21 303 U31 U12 U22 4010 3639s U32 307 U13 U23 309 U33 If we order the equations the same way we order the unknowns then the nine equations from the standard central finite difference scheme using the five point stencil are 1 um 1110 i 4y y ar 7 h2 1 2 4 fn L2 1 1120 i 4 r t i h2 71 32333J le kg 1 quot30 1441 72 4333 6 f3 1 1402 ff lflk l fIkIf 9 7 h2 1 4 o flz L2 1 72 47 L39s fzz 1 I H42 72 363 4366 369 f32 72 1 um 1114 im 4m7m 7 h2 4 8 f1 L2 1 1124 72 7 4308 9 f23 72 1 1034 1143 7 43 L2 as 38 9 fix L2 Now we can write down the coefficient matrix easily It is black tridz39agonal and has the following form 13 I 0 1 A i 1 B 1 4011 01quot Numerical Analy where I is a 3 X 3 identity matrix 4 1 0 B 1 4 1 0 1 4 For a general 71 by 71 grid we will have B I 4 1 1 I B I 1 4 1 A 72 B I B 1 4 Note that A is a symmetric positive definite matrix and it is weakly diagonally dominant Therefore A is nonsingular and there is a unique solution The matrixvector form is useful to understand the structure of the linear system of equations and it may be necessary if a direct method such Gaussian elimination or sparse matrix techniques are used for solving the system However it is more convenient especially if an iterative method is used V sometimes to use the two parameters system i j to solve the system It is more intuitive and useful to visualize the data using two index system The eigenvalues and eigenvectors of A can be indexed by two parameters p and k corre sponding to wave numbers in the a and y directions The p kth eigenvector UNquot has 712 elements for a n by 71 matrix of the form above 1 sinp7rilt sink7rjlt ij 12 lt71 4012 quot3393 for p k 1 2 n The corresponding eigenvalues are 2 A 72 cosp7rlt 1 cosk7rlt 1 4013 The least dominant eigenvalue the smallest in the magnitude is A 27r hz 4014 The dominant eigenvalue the largest in the magnitude is 4 392 392 W W 772 4015 Therefore we have the following estimates 4 1 1 1144112 mame 2 1144 le 7 M It min A 277 4016 2 crmd2lt4gt1141121144112 W 00 Since the condition number is considered to be large we should use double precision to reduce the effect of round off errors 60 Z Li 41 Basic iterative methods for solving linear system of equa tions The idea of iterative methods is to start with an initial guess then improve the solution iteratively The first step is to rewrite the original equation f x 0 to an equivalent form x gx and then we can form an iteration xlk gxk For example to find is equivalent to solving the equation 303 5 0 This equation can be written a 5302 or a a 3 O For the second one the iteration 3 3 mltk1gtxltklef 201 30022 39 39 is called the Newton s iterative method The mathematical theory behind this is the xed point theory For a linear system of equations Ax b we hope to rewrite it an equivalent form x R c so that we can form an iteration Xkl Rle 3 given an initial guess x We want to choose such a R and c that limknooxm x5 A lb A common method is called the splitting approach in which we rewrite the matrix A A M K detM 7 0 411 Then A b can be written AI Kb b or AI Kxb or x AI IKxHil lb or x R c where R M quotK is called the iteration matrix and c M 1b is a constant vector The iterative process is then given an initial guess X we can get a sequence of x00 according to xWU Rxlk 3 412 We first discuss three basic iterative methods for solving Ax b To derive the three methods we rewrite the matrix A a 0 122 021 0 A D L U 033 031 032 0 04m Tan anz 39 anmil 0 0 012 lm 0 mu Numerical Analysis I 61 42 The Jacobi iterative method Solve for the diagonals The matrix vector form of the The Jacobi iterative method can de derived follows D L 1 b D L 1 b x DquotL 1 D lb xltk1gt DquotL UXk D 1b The component form can be written R kl k xi I Z aux m z 1 2 H m 421 PM The component form is useful for implementation while the matrixvector form is good for convergence analysis 43 The GaussSeidel iterative method Use the most up dated In the Jacobi iterative method when we com ute N we have alreadv com uted wk 2 1 Assume that an is a better approximation than why can not we use an update 303 instead of With this idea we get a new iterative method which is the when we GaussSeidel iterative method for solving Ax b The component form is 1 71 k n k r 1 r 1 r mgH bi 01 2 17331 a t 123 71 431 j1 ji1 To derive the matrixvector form of the GaussSeidel iterative method we write the component form above to a form k1 k The component form above is equivalent to 1 71 7 2 aingk0 7731ng 1 2 ijmgk i 1 2 4 H 41 432 j1 ji1 which is the component form of the following system of equations xltk1gt U b or xltk1gt D LquotUx D Lquotb Thus the iteration matrix of the GaussSeidel iterative method is D L 1U and the constant vector is c D Lquotb 62 Z Li 431 Implementation details The iterative methods are an infinity process However when we implement on a computer we have to stop it in finite time One of several of the following stopping criteria are used I 1xk XWH S tol ItHi k x l g ml I 39 HTXkUH S toL 39 k 2 kmaxt where ml and km are two given parameters 432 Pseudocode of the GaussSeidel iterative method Exk mygsnabx0tol error 1e5 xx0 k0 while error gt tol for i1 n xCi bCi for j1zn if j i xCi xCi aijxj end x1 xltiaii end error normxx0 defaut is 2norm x0x k k1 replace the old value add the counter end end while 433 The GaussSeidel iterative method for 2point boundary value prob lem Following section 13 we have Ur 2U U 7 71fyi 139 12 n l It If we use the same ordering for the equations and unknowns then the diagonals are always 2h2 the Jacobi iteration is simply 4k 4k 2 U U I U kl m y f i1fgf n1 Numerical Analy s39 I 63 No matrix is formed or only matrixvector multiplication is needed no ordering is neces sary assuming that the equations and unknowns have the same ordering For the Gauss Seidel iteration from kth to k 1th iteration we can use the following UWH U7 139 12 n l 1 Uk1 1 U1k1U1k1 2 30 i12w n l If UYFHI has not been updated then it use the value from kth iteration otherwise it uses the most updated one 434 The GaussSeidel iterative method for the nite difference method for Poisson equation For the Poisson equation um 11W f 31 y if we use the standard 5point central finite difference scheme UH 1413 4U 12 1 by 1 J 7 J g 1 1 fatmfg and the same ordering for the equations and unknowns then the Jacobi iteration is r k k k gtU gtUgt 2 J W 4 w 1 W L few i12w n l j12w n l if the solution is prescribed along the boundary Dirichlet BC Again no matrix is needed no ordering is necessary We do not need to transform the two dimensional array to a one V V one The i is rather simple 44 The successive overrelaxation SORw iterative method The Jacobi and GaussSeidel methods can be quite slow The SORw iterative method is an acceleration method by choosing appropriate parameter 1 The SORw iterative method is xltk1gt 1 UM wiggl 441 The component form is 171 k n k 1 merl 1 wyf 0 1 Sagan 2 was l7 7 442 F1 ji1 Note that it is incorrect to get GS result first then do the linear interpolation 64 Z Li The idea of the 80R method is to interpolate xk and to get a better approximation When u g 1 the new point of 1 wxk is between xk and and this it is called interpolation The iterative method is called under relaxation When u gt 1 the new point of 1 wxk is outside xk and and this it is called extrapolation The iterative method is called over relaxation Since the approach is used at every iteration it is called successive over relaxation 80R method To derive the matrixvector form of SORw method we write its component form 1 71 k n k i r 1 r r 175me 32 mm l 1 win w 2 mac 443 F1 jil This is equivalent to D wLxkl 1 MD amxm wb Thus the iteration matrix and constant vector of the SORw method are R go w D wal 1 wD wU 3301 wD wal 444 45 Convergence of basic iteration methods Xltkn RX c Using an iterative method we will get a vector sequence of xquot and we know how to tell whether it is convergent or not However for an iterative method we need to consider all possible initial guesses and constant vector 3 If the vector sequence of x00 converges to xquot then by taking limit on both sides of the iterative scheme we have x Rxquot c 451 The above equality is called the consistency condition De nition 451 The iteration methods calm Baum c convergent if for any initial guess 500 and constant vector 0 the vector sequence of new converges to the solution of the system of equations acquot likequot c Now we discuss a few su icient conditions that guarantee convergence of a basic iterative method Theorem 451 If there an associated matrix norm such that HRH lt 1 then the iteration method calm Baum c converges Numerical Analysis I 65 Proof Let elk xik xquot from the iterative method xlk Rle c and the consistency condition xquot Rxquot c we have 606 Rea 0 S Helm HRelk H S HRHHeWI S HRHHRHHe k WI S S HRHWHBWI Thus we conclude that lim Hem 0 or equivalently lim xik x kgt00 kgt00 Example Given R 09 005 08 01 Does the iterative method xlk Rle c converge We can easily get HRHl 17 which leads to no conclusion and HRHOOH 095 lt 1 which lead to the conclusion that the iterative method converges 451 Convergence speed In the theorem above the kth error depends on the initial one that we do not know The following error estimate does not need the initial error Theorem 452 If there an associated matrix norm such that HRH lt 1 we have the following error estimate for the iteration method calm Batik c HRHk llew S W use Wu 452 Proof From the iterative method xik Rxlk c we also have fxik Rxlk l c Subtracting the two we get xk1 x00 Rxk x0670 32xltk1 x0e H Rkx1 x0 Since xk1 x00 xk1 x 3 X 60 606 R new Combining the two equalities above we get I Rek Rkx1 x0 This leads to new HU Rr Rk W xlt gt 11 Finally from the Banachquots lemma we have Milk 0 He H s 1 HRH Hx1x0H 66 Z Li 452 Other sufficient conditions using the original matrix Theorem 453 If A strictly row diagonally dominant matrix then both Jacobi and GaussSeidel methods converge The GaussSeidel method converges faster in the sense that HRasHoo S HRJHOO lt 1 Proof The proof of the first part is easy For the Jacobi method we have R DquotL U thus HIM mgx lt1 j lija i aii The proof for the GaussSeidei method is not trivial and long we refer the readers to the book 3 W Demmei on page 287288 For general matrices it is unclear whether the Jacobi or GaussSeidei method converges faster even if they both converge Theorem 454 If A a symmetric positive de nite SPD matrix then the SOBw method converges for 0 lt w lt 2 Again we refer the readers to the book 3 W Demmei on page 290291 Theorem 455 If A a weakly row diagonally dominant matrix R 2 aijigiaiiit ti1t2t 39tnt jlij9 i with at least on inequality strictly and A irreducible that there no permutation PTAP A411 A412 0 A422 then both Jacobi and GaussSeidel methods converge The GaussSeidel method converges matrix such that faster in the sense that HRGsHoo S HRJHOO lt 1 Lemma 451 A irreducible if the graph of the matrix strongly connected Numerical Analysis I 67 453 A suf cient and necessary condition The spectral radius of a matrix A is defined 114 m1axiA 453 Note that from 114 3 from any associated matrix norms The proof is quite simple Let Ar be the eigenvalue of A such that oA AM and xquot 7 0 is the corresponding eigenvector then we have Axquot Aimquot Thus we have A x l g HAHHx l Since Hx 7 0 we get MA 2 HAM Theorem 456 Ari iteration method 5W Batik c converges for arbitrary arm and c if and only if oR lt 1 Proof Part A If the iterative method converge then oR lt 1 This can be done using counter proof method Assume that oR gt 1 then let Axquot Apxquot with oA AM gt 1 and Hx 7 0 If we set X xquot and c 0 then we have xik flxquot which dost not have a limit since A11 gt so The case of oR 1 is left an exercise Proof Part B If oR lt 1 then the iterative method converges for arbitrary xi and 3 They key in the proof is to find a matrix norm such that lt 1 From liner algebra Jordan s theorem we know that for any square matrix B it is similar to a Jordan canonical form that is there is a nonsingular matrix S such that 11 V 1 12 M 1 5455 1 JP Ai Note that HS IRSHQO oR lt 1 if all Jordan blocks are 1 by 1 matrix otherwise 115415551100 3 oR 1 Since oR lt 1 we can find lt 06 such that oR c lt 1 c 1 oR2 Consider a particular Jordan Block Ji and assume it is a k by k matrix Let gt494 68 Z Li Using this diagonal matrix we can get 1916 1916 D6 13175 Thus HD I66 1175551 91100 0R 6 lt 1 Finally we define a new matrix norm Rummy HD715SHRSDHOO 454 It can easily show that the definition above is indeed an associate matrix norm and HRHan oR c lt 1 we conclude that the iterative method converges for any initial guess and vector c 46 Discussion for the Poisson equations What is the best to For the finite difference methods for Poisson equations in 1D and 2D we can find the eigenvalues of the coef cient matrix which also leads to eigenvalues of the iteration matrix Jacobib Gauss Seidel SORw For 1D model problem a 31 f 31 with the Dirichlet boundary condition 140 and 141 are prescribed the coe icient matrix is The matrix is weakly row diagonally dominant not strictly and irreducible A is symmetric positive definite or A is symmetric negative definite Theorem 461 For the above a by a matrix we have MA 0RJ 5 1 T 461 where AAA i 1 2 7 are the eigenvalues of the matrix A R the iteration matrix of the Jacobi method Numerical Analysis I 69 Proof We know that By D 1L U and D 21 Let A be an eigenvalue of 13 then detgt1 JR 0 or detgt1 D 1L U 0 or detD lgtD L U 0 or detDdetgtD L U 0 Since detD 7 0 we have detgtD L U 0 or detgt 1D D L U 0 or det 1D A 0 or det1 AD A 0 or det 21 MI A 0 since D 21 Thus 21 A is an eigenvalue of A or X A AjRJ1g t12lt 7l From the relation above we have the theorem right way The following lemma gives the eigenvalues of a tridiagonal matrix Lemma 461 Let A be the following n by n matrix The eigenvalues of A are the following I k12 n n1 The eigenvector corresponding to M M cr 23cos kj7r n1 Xmsin j12 n Proof It is easy to check that AXk A1 When 139 2 3 1 we have k7 n 1 Ak 22cos k12 n Thus the spectral radius of the Jacobi iterative method for 1D Poisson problem is AMA 21 cosk7rn 1 r 1 7 r 1 MR 2 2 m cos 7k cos if 1 1 7T 2 I X N 7 1513 n1 n1 2 n1 We can see that n is getting larger the spectral radius is getting close to unit indi cating slower convergence Once we know the spectral radius we also know roughly the number of iterations needed to reach the desired accuracy For example if we wish to have roughly six significant digit the we should set 0Rk g 10 6 or k 2 6logm0R 70 Z Li 461 Finite difference method for the Poisson equation in two dimen sions In two space dimensions we have parallel results The eigenvalues for the matrix A IVAN of N by N matrix N n2 are i7r j7r A 4 2 3912H 42 Aw cos n1 cos n1 i y y n 6 The diagonals of the matrix are 4 and we have A A itj1t2Wthi Thus A xfl 21 u39 1 1 1RJ max 1 7 max 1 Mm BOW I 113312 4 113312 4 I m cos 7 cos jr cos if 1 1 7T 2 ax u u u w 7 1313ng n1 n 1 n 1 2 n1 We see that the results in 1D and 2D are pretty much the same To derive the best to for the SORw method we need to derive the eigenvalue relation between the original matrix and iteration matrix Note that R90 D wL 11 wD wU Theorem 462 The optimal to for the SOBw method for the system of equations derived from the nite di erence method for the Poisson equation 2 2 2 2 WW 2 1 1 oUt 2 1 1 C082 1 sinnLJ r1 N 4 63 Below is the sketch of the proof I Step 1 Show the eigenvalue relation between RJ and Rsogw Asm w 1 w szg its 1 w 464 I Step 2 Find the extreme value minimum of above a function of M which leads to the optimal to We refer the readers to the book 3 W Demmel on page 292293 for the proof Remark 461 Numerical Analysis I 71 I The spectral radius of Rsogw 1 t w2pRJ2 1 y 7 2 2 y 1 y y y RSOR Lu 2L RJ b RJ Lu 4 0ltLu ltu 0pt w l wtmgltwlt2 If we plot M1330 against to we can see it a quadratic curve between 0 lt w lt 2107 and it at as w getting closer to wisp which means it less sensitive in the neighborhood of dam while the second piece a linear function Thus we would rather choose to larger than smaller The optimal to only for the Poisson equation not for other elliptic problems How ever it gives a good indication of the best to the diffusion dominant in reference to the mesh size 72 Z Li Chapter 5 Computing algebraic eigenvalues and eigenvectors 5 1 Preliminary Given a square matrix A 6 RM if we can find a number A E C and x 7 0 such that A Ax then A is called an eigenvalue of A x is called a corresponding eigveiwectmn to A Note that if A Ax then Ayx Ax for any nonzero cnstant c in other words eigenvector can differ by a constant Often we prefer to use an eigenvector with unit length 1le 1 We call Ax an eigenpair if Ax A x 7 0 For an eigenpair A x we have A Ax 0 This means that AI Ax 0 has nonzero or not unique solutions This indicates that AIA is singular or detAIA 0 Thus A must be a root of the characteristic polynomial of the matrix A detAI A A LR71AR71 01A Lo There are n eigenvalues for an n by m square matrix The eigenvalues can be real complex numbers repeated roots If the matrix is real then the complex eigenvalues are in pairs that is if A a bi is one eigenvalue then A a bi is also an eigenvalue If A is a real and symmetric matrix then all the eigenvalues are real numbers Different eigenvectors corresponding to different eigenvalues are linear independent If an eigenvalue Aquot has multiplicity p which means that characteristic polynomial has the factor A A7 but no the factor A A7 1 the number is independent of eigenvectors corresponding to Aquot is less than or equal to p recall some examples in class If an n by m square matrix A has n linear independent eigenvectors then A is diagonalizable that is there is a nonsingular matrix S such that S IAS D where D diagA1 A2 An is a diagonal matrix For convenience of discussion we will use the following notations We arrange eigenval 73 74 Z Li ues of a matrix A according to M Z M Z 39M Z 2 ML Thus 0A Md A1 is called the dominant eigenvalue can be more than one while A is called the least dominant eigenvalue There are many applications of eigenvalue problems Below are a few of them Frequencies if vibration resonates etc Spectral radius and convergence of iterative method Discrete form of continuous eigenvalue problems for example the SturmLouville problem pau a qaua Aidan 0 lt a lt 1 After applying a finite difference or finite element method we would have A Ax The solutions are the basis for the Fourier series expansion I Stability analysis of dynamic systems or numerical methods 52 The Power s method The idea of the power method is starting from a nonzero vector xm 0 an approximation to an eigenvector then form an iteration xkl Axm Azx cil 4 4 4 4k1x0 Then under some conditions we can extract an eigenpair information from the sequence With the assumption that A1 satisfies Ad gt A2 which is an essential condition then we r 1 r k 1 can show that x N OX1 and 30 Wang Jr Xkl l Sketch of the proof A1 k is very large where While feasible in theory the idea is not practical in computation because x gt 0 if A lt 1 and x gt so if A gt 1 The solution is to rescale the vector sequence which leads to the following power method Numerical Analysis I Given xm 0 form the following iteration for k 1 until converges yum Axon k1 xltk1 yk l Hquot 12 W xkl 139Axkl end We can use the following stopping criteria marl ml lt tol or Hyik l yiklll lt tol or both Under some conditions the sequence of the pair maxim converge to the eigenpair corresponding to the dominant eigenvalue For simplicity of the proof we assume that A has a complete eigenvectors v1 vz vn vi 1 AW NW Theorem 521 Assume that Ad gt A2 2 A3 gt2 2 An and 500 221 moi with m 7 0 then the pair mg sum the power method converges to the eigenpair corresponding to the dominant eigenvalue Proof Note that k k Axltk71 4 3quot 1 1H 7 A lt y llyltkquot llz 7quot y Vk39l39kilAZDkiz quotYk 739k71 quot 1 4k71y1 CkAka where 39ymk71 7 1 and Up are some constants Since xik is parallel to yik and has unity length in 2norm we must have xua k A 3 Hz On the other hand we know we have k k A2 k A k A X A1 mm 12172 7 man 7 A1 A Thus we have k lim xm lim x7 ivl kgt00 kgt00 Hxikluz Hmvlllz and lim Mk lim xkrlh4xk A1 kgt00 kgt00 76 Z Li 521 The power method using the in nity norm If we use different scaling we can get different power method Uisng the infinity norm first we introduce the 30 notation for a vector x Given a vector x 30 is the first component such that HXHOO and p is the index For example if x 2 1 5 5 5 V then mp 5 with p 3 53 The inverse power method for the least dominant eigen value If A is invertible then 1Ai are the eigenvalues of A l and 1 will be the dominant eigenvalue of A71 However the following inverse power method for the least dominant eigenvalue without the need to have the intermediate steps Given xm 0 form the following iteration for k 1 until converges Ayltk1gt x00 k1 xk 1 yk Hquot 112 W xk1 1quot4xkl end With similar conditions PM lt Awd g 3 A1 the essential condition one can prove that 3320 M 3323 m with mm 1 Note that at each iteration we need to solve a linear system of equations with the same coef cient matrix This is the most expensive part of the algorithm An ef cient implementation is to get the matrix decomposition done outside the loop Let PA LU the algorithm can be written follows Numerical Analysis I 77 Given xm 0 form the following iteration for k 1 until converges mm pxm Uyk1 zk1 k1 xltk1 yk l W Hz W xk139rAxk1 end 54 Gershogorin theorem and the shifted inverse power method If we know a good approximate a to an eigenvalue A such that A 7 lt min a l in W 3K ii That is A a is the least dominant eigenvalue of A a We can use the following shifted inverse power method to find the eigenvalue A and its corresponding eigenvector Given xm 0 form the following iteration for k 1 until converges A anyww M k1 xltk1 y HyW le W xk139rAxk1 end Thus if we can find good approximations of any eigenvalue we can use the shifted inverse power method to compute it Now the question is how do we roughly located the eigenvalues of a matrix A Gershgorin theorem provides some useful hints De nition 541 Given a matrix A 6 CW the circle all points within the circle on the complex plane R A ml S 2 ijl 541 jlgtj9 i called the ith Gershgorin circle 78 Z Li Theorem 541 Gershgorin Theorem 1 Any eigenvalue have to be in one of Gershgorin circles 2 The union of k Gershgorin circles which do not intersect with other n k circles contains precisely k eigenvalues of A Proof For any eigenpair Am A Ax Consider the pth component of x such that HxHOO we have Th 2 WV393 Any3917 j Th A apply3977 2 0mm Flor From the expression above we get R M R M am 2 an s 2 am s 1 j m Flow Thus A is in the pth Gershgorin circle and the first part of the theorem is complete The proof for the second part is based continuation theory that roots of a polynomial are continuous functions of the coe icients of the polynomials The theorem is obviously true for the diagonal matrix If the radius of the Gershgorin circles increase continuously we change the zero off diagonal entries from 0 to aij the eigenvalues will move among union of the Gershgorin circles but can not cross to the disjoint ones Example Let A be the following matrix 5 1 0 A 1 2 12 0 1 8 Use the Gershgorin theorem to roughly locate the eigenvalues The three Gershgorin circles are IRl z5gl I R2 z 2gl5 01333 2 8Sl The do not interset with each other so each circle has one eigenvalue Since the matrix is a real matrix and complex eigenvalues have to be in pair we conclude that all the eigenvalues are real Thus we get Numerical Analy 79 I the dominant eigenvalue satisfies 7 3 A1 3 9 I the least dominant eigenvalue satis es 05 3 A3 3 35 I the middle eigenvalue satisfies 6 3 A1 3 4 If we wish to find the middle eigenvalue A2 we should choose a 5 Even for the dominant eigenvalue we would get faster convergence if we shift the matrix by taking 039 8 and then apply the shifted power method 55 How to nd a few or all eigenvalues If we know an eigenpair A1 x1 of a matrix A we can use the de ation method to reduce A to a oilediineiisioiial lower matrix whose eigenvalues are the same the rest eigenvalues of A The process is follows Assume that HXIHZ 1 we can form expand x to form an orthogonal basis of R x1x2 Hxn with xzyxj 373 Note that 373 0 ifi j and 35 1 Let Q x1x2 0cm then Q is an orthogonal matrix QTQ QQT I We can get QTAQ QTAx1sz Am A1 39139 2 1x1 x2 Xu 0 A Thus the eigenvalues of A are also those of A but Al is a onedimensional matrix compare with the original matrix A The deflation method is only used if we wish to find a few eigenvalues To find all eigenvalues the QR method for eigenvalues are often used The idea of the QR method is first to reduce a matrix to a simple form often upper Hessenberg matrix or tridiagonal matrix using similarity transformation S IAS so that the eigenvalues are unchanged Since the inverses of an orthogonal matrix is its transpose S is often chosen orthogonal matrix Orthogonal matrices also have better stability than other matrixes since HQxllz Hxllz and 1194112 114112 De nition 551 Giver a unit vector in Hsz 1 the Household matrix de ned as P I 2mm 551 71 It is a simple check that P 2 PT 2 P Such a matrix sometimes is also called a reflection matrix or transformation Theorem 551 If 112 Hsz then there a Household matrix P such that Pan y Proof Assume that P y then we have I QWWT x y x y 2wwa Note that 2w739x is a number thus w is parallel to x y Since w is also a unit vector in 2norm we conclude that w x yHx sz then it is simple manipulation to show that P y Example Find a Householder matrix P such that P 0 0 Note that in this example we need to find both 1 and P Since the orthogonal transfor mation does not change the 2norm we should have V3242 1392 21 5 1 3 i 139 W x yHx yllz 7 0 0 H 3 4 0 To avoid possible cancellation we should choose the opposite sign that 139 5 and W 8 0 ANAE 551 The QR decomposition of a matrix Start from AD 2 A for k 0 1 until converge Ak QkRK Ak1 13955ka end Theorem 552 If A 6 RM and Ad 2 AA 2 2 ML then 39 1i1 N Ai A that all Ak have the same eigenvalues I RH Rm Rlp lim Ak RA 22 kaoo Rm Numerical Analy s I 81 where RA a block upper triangular matrix whose diagonals R either a 1 by 1 or 2 by 2 matrix corresponding to complex eigenvalues in pair Proof of the first part 14k1 13955ka QzAkRk Shifted QR method Start from AD 2 A for k 0 1 until converge Ah UIJ Qle Ak1 Bika 0k end Stop criteria I max am lt tol 337SWISJSF2 Double shifted QR method If a is chosen a complex number then we should use the Double shifted QR method Start from AD 2 A for k 0 1 until converge Ak 0k QkRk Ak1 Bika 0k Ak1 5k Qk1Rk1 Ak2 Rhlel 51439 end QR method for finding all eigenvalues are quite expensive To reduce the computational cost Often we use the similarity transformation via Householder matrix first to reduce the original matrix to an upper Hessenberg matrix tridiagonal if the matrix is symmetric first then apply the QR method This requires rt Q steps PRAPRA PZPI AP1P2 PRPR2 where L21 t 1 0 1 aa 0 P1 7 P1 1 0 P1 I 04 0 for example The reason to choose P1 to keep the first row of A unchanged when we multiply P1 from the left is to ensure that the first columns of P1 A unchanged when we multiply P1 from the right so that those zeros will remain Z Li Chapter 6 Least squares and SVD solutions In this chapter we discuss numerical method for solving arbitrary linear system of equations 61 Least squares solutions Consider Ax b where A E Rm m 2 n and mnMA n One motivation is the curve fitting example in which we have many observed data but we find a simple polynomial to fit the data Below are some examples From the first two equations the solution should be a 0 from the third equation the solution should be a 1 from the last equation the solution should be a 2 In other words we can not find a single 3 that satisfy all the equations The system of equations is called overdetermined In general there is no classical solution to an overdetermined system of equation We need to find the bestquot solution By the best solution we mean that to minimize the error in some norm Since the residual is computable one approach is to minimize the residual in 2norm of all possible choices The solution xquot is the best solutionquot in 2norm if it satisfies Hb Ax lz ming 21wa 611 XE R If we use different norms rather than 2norm it will leads to different algorithm and different applications But 2norm is the one that is used the most and it is the simplest one 83 84 Z Li Note that an equivalent definition is the following Hb Axwg ming Axug 612 x6 R The above definition gives way to compute the best solutionquot the global minimum of a multivariable function of its component m1 m2 mn This can be seen from the following 2 39139 39139 39139 39139 39139 39139 39139 rgtx Hb Ax12 b Ax b Ax b b b Ax x A bx A Ax Note that bTAx xTATb one can easily get the gradient Mx which is woo 2A39139Ax ATb Since the columns of A are linearly independent mnkA n we have XTATAX HAXHZ gt 0 for any nonzero x and ATA is symmetric positive definite The only critical point of rx is the solution of the following normal equation ATAx ATb 613 whose unique solution is the least squares solution which minimizes the Qnorm of the residual in R space The normal equation approach not only provides a numerical method but also shows that the least squares solution is unique under the condition mnMA n that is the columns of A are lienarly independent A serious problem with the normal equation approach is the possible illconditioned system Note that if m n then cow1AA crmd2A2 A more accurate method is the QR method for the least squares solution For the following least squares problem R b1 x 0 b2 We have 2 2 2 Mb AXHz Hbl Rxllz WM and 1113 Axu 36 mar qu mug ubzu when b1 R 0 or x R lbl Particularly if R is an upper triangular matrix we can use the backward substitution to solve the tridiagonal system of equations ef ciently In the QR method for the least squares the idea is to reduce the original problem to the problem above using orthogonal particularly the Householder transformation which Numerical Analysis I 85 keeps the least squares solution unchanged From the process of the QR algorithm we can apply a sequence of Householder matrices that reduce A to an upper triangular form PRPR1 HP1A 15 1 or QA RT 0 gt then from A b we get QAX Qb or R 131 x N 0 b2 In practice we can apply the Householder matrix directly to the augmented matrix A E b in the Gaussian elimination method An example 62 Singular value decomposition SVD and SVD solution of AX b Singular value decomposition can be used to solve Ax b for arbitrary A and b v Theorem 621 Given any matrix A 6 0mm there are two orthogonal matrices U E mem UUH UHU I and V E 1mm VVH VHV I such that 71 A U 2 V where 2 621 71 TIIXR 71 72 H 0 gt 0 are called the singular values of A Note that they are positive numbers p rankA Furthermore of iA A iAA square root of nonzero eigenvalues of AHA or AA I illaxlgjgp 77 A Note that we often arrange singular values according to 71 2 72 2 2 7 gt 0 The proof process also gives a way for such a decomposition constrictive Proof Let 71 HAHz then there is an x HXIHZ 1 such that AHAxl xl Let y Ax1m we have Hm HAX1HzUl 1 Z Li Next we expand x1 to form an orthogonal basis in BTW to form an orthogonal matrix V Vx1x2wxn1x1 V1 WV WV I We also expand yl to form an orthogonal basis in Rmxm to form an orthogonal matrix U UlY1Y2W tyn1 y U1 U U U U 1 Then we have H 0 UHAV 2 Y1 Ax AV gt 71 U1 1 1 0 U1 AV This is because yf Axl xf lA AxlUl xf xla m 71 Ul Axl Ill39139in 0 yf AVl Ay1 VI xf AHAVI x VI 0 Thus from the mathematics induction principle we can continue this process to get the SVD decomposition Pseudoinverse of a matrix A From the SVD decomposition of a matrix A we can literally find the inversequot of the matrix to get its pseudoinverse A V 2 U where 2 622 l 0 nxm Particularly if m n and detA 0 then 21 A l The pseudoinverse matrix 21 of a matrix A has the following properties AAA A note that AA 7 I AAA 2 Air I 21er A A If mnkA n then 21 AHAYIAH I 2121 AA If mnkA m then 21 A AA 1 Solving AX b of arbitrary matrix A The SVD solution of A b is simply xquot Airb It has the the following properties Numerical Analysis I 87 I If m n and detA 7 0 then xquot 2 Air 2 A lb I If mnMA n then MAXquot sz Xlllnxgign MAXquot bug that is 1quot is the least squares solution I If there is more than one classical solution then xquot is the one with minimal 2nornl that is 5 xl zl v W W 223 W blz