Systems of linear equations - reminder Solving linear systems Special cases Examples and applications BÍ7740: Scientific computing Systems of linear equations Vlad Popovici popovici@iba.muni. cz Institute of Biostatistics and Analyses Masaryk University, Brno IBA BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Additional references: » Golub, Van Loan, Matrix Computations, Johns Hopkins Univ. Press, 3rd Ed. 1996 IBA Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Outline O Systems of linear equations - reminder • Norms • Linear systems • Conditioning » Accuracy Q Solving linear systems • Diagonal systems • Triangular systems • Gaussian elimination Q Special cases • Symmetric positive definite systems Q Examples and applications j^' Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Norms Linear systems Conditioning Accuracy Outline O Systems of linear equations - reminder 9 Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems 9 Diagonal systems • Triangular systems 9 Gaussian elimination O Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Norms Linear systems Conditioning Accuracy Outline O Systems of linear equations - reminder • Norms 9 Linear systems • Conditioning o Accuracy Q Solving linear systems 9 Diagonal systems • Triangular systems 9 Gaussian elimination O Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Vectors and norms Let x be a vector, x as [xi,..., xn]7. The p-norm is defined Zw W=1 Special cases: op = 1: (Manhattan or city-block norm) IWh =Z,|x;| » p = 2: (Euclidean norm) ||x||2 = J^-, xf o p -> oo: (oo-norm) ||x||oo = max,|x,| The unit circles. IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Vector norms - properties Vx,y e Rn and for any norm, » ||x|| > 0 with ||x|| = 0ox = 0 • Harxll = \a\ • ||x||, Va • ||x + y|| < ||x|| + ||y|| (triangle inequality); also ll|x||-||y|||<||x-y|| • ||x|h > ||x||2 > ||x|L • ||x|h < Vn||x||2, IMI2 < Vn||x|U -» norms differ by at most a constant, hence they are equivalent Matlab: norm (x, p) IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Norms Linear systems Conditioning Accuracy Matrix norms Let a-n 312 a-in 3n1 3n2 3nn be a square matrix. » defined based on a vector norm max x#0 l|Ax|| a the maximum "stretching" applied to a vector by the matrix A 9 ||A|h = maxy£/Li |a,)| (maximum absolute column sum) • ||A||oo = max/2Li la')'l (maximum absolute row sum) • IIAH2 =? (we'll see it later) Vlad Bi7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Matrix norms - properties Let A and B be two square matrices 9 ||A|| > 0 if A * 0 9 \\aA\\ = \a\ ■ ||A||, for any scalar a • ||A + B||<||A|| + ||B|| • ||A-B||<||A||-||B|| • ||Ax|| < ||A|| • ||x|| for any vector x Matlab: norm (A, p) IBA Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Norms Linear systems Conditioning Accuracy Outline O Systems of linear equations - reminder • Linear systems 9 Conditioning o Accuracy Q Solving linear systems 9 Diagonal systems • Triangular systems 9 Gaussian elimination O Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Linear systems In general, a system of linear equations has the form: anX! + a12x2 + --- + ainxn = bi a2i + a22x2 + • • • + a2nxn = b2 a^X! + am2x2 h-----h amnxn = bm or, in matrix format, Ax = b where A is an m x n matrix (say, A g Atm,n(^)), b and x are vectors with m and n elements, respectively. In other words: can the vector b be expressed as a linear combination of columns of matrix A? Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy In Matlab o general operator "matrix division" \ » this is a wrapper for various algorithms - some we will discuss Matlab: x = a \ b Vlad BÍ7740: Scientific computing IBA Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Square matrices case (m = n) A e Mn,n{^) is singular if it has any of the following equivalent properties: » A has no inverse (A-1 does not exist) • det(A) = 0 » rank(A) < n (rank: maximum number of rows or columns that are linearly independent) » Az = 0 for some vector z + 0 Otherwise, the matrix is nonsingular. If A is nonsingular, there is a unique solution; otherwise, depending on b, there might be zero or infinitely many solutions. Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Geometrical interpretation (2D): 9 a linear equation defines a line • if A is nonsingular, the two lines intersect » if A is singular, the two lines may be parallel (no solution) or identical (infinitely many solutions) If A is singular and b e span(A) the system is consistent and has infinitely many solutions. (span(A) is the vector space generated by the columns of A.) IBA Bi7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Example - nonsingular matrix , then A is nonsingular and there is a unique solution, x = 1 Naive Matlab solution: » A=[l 2; 3 4]; b=[-l;-l]; » x = inv(A)*b; » x X - 1.0000 -1.0000 o let A = 1 2 and b = 3 4 IBA Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Norms Linear systems Conditioning Accuracy Example - singular matrix 1 2 -1 and b = 2 4 -1 o let A is a unique solution, x Naive Matlab solution: then A is nonsingular and there » A=[l 2; 2 4]; b=[-l;-l]; » x - inv(A)*b Warning: Matrix is singular to working precision. -Inf -Inf IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Outline O Systems of linear equations - reminder • Linear systems • Conditioning 9 Accuracy Q Solving linear systems 9 Diagonal systems • Triangular systems 9 Gaussian elimination O Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Singularity, norm and conditioning • condition number of a nonsingular square matrix is cond(A) = ||A||-||A-1|| » convention: cond(A) = oo for singular A • ratio between maximum streching and maximum shrinking of a nonzero vector 9 large cond(A) indicates a matrix close to singularity 9 small det(A) does not imply large cond(A) *— Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Condition number - properties » cond(A) > 1 » cond(/) = 1 (/is the identity matrix - Matlab: eye(n)) » cond(aA) = cond(A), for any A and scalar a o for a diagonal matrix D = diag(dj), d, * 0 we have cond(D)= =gf • condition number is used for assessing the accuracy of the solutions to linear systems IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Condition number: • exact computation requires matrix inverse: <» ||A|| is easy to compute • computing at low cost ||A~1|| is difficult -> expensive (even more than finding the solutions to the problem) and prone to numerical instability • in practice: estimated as a byproduct of the solution process One approach: find lower bounds on ||A"1|| and, thus, on cond(A). If Ax = y it follows that ^ ||r|| = 0 » practically, small residual is not necessarily equivalent to small error 9 since IIAxll ^ .... ||r|| —77- < cond(A) small relative residual implies small relative error, only if A is well-conditioned Vlad Bi7740: Scientific computing Systems of linear equations - reminder Norms Solving linear systems Linear systems Special cases Conditioning Examples and applications Accuracy Residuals - backward error analysis 9 let D be the "delta" matrix, such that x is the exact solution of (A + D)x = b, then llrll < IID|| l|A|| • ||x|| " 11 A|| 9 large relative residual implies large backward error and indicates an unstable algorithm • stable algorithms yield small relative residuals, regardless conditioning of nonsingular A IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Outline Diagonal systems Triangular systems Gaussian elimination 9 Systems of linear equations - reminder • Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems • Triangular systems 9 Gaussian elimination O Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications General strategy Diagonal systems Triangular systems Gaussian elimination o transform the system (mainly A) such that the solution is easier to compute (but unchanged) • if M is a nonsingular matrix the systems Ax = b MAx = Mb HOMEWORK: prove it! and have the same solution. <» trivial transformations: » permutation of rows in the system: use a permutation matrix (has exactly one 1 in each row and column, rest is 0). <» diagonal scaling: may improve the accuracy IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ..... _ . ,. Gaussian elimination Examples and applications A few relevant functions in Matlab Please, use help for details! • linsoive: solves linear systems Ax = B via various methods. You can specify the properties of A. • \ operator is a wrapper for various methods • lu computes LU factorization • triu returns upper triangular part of a matrix • trii returns lower triangular part of a matrix a diag returns the diagonal of a matrix o cond, condest used for estimating the condition number Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Outline Diagonal systems Triangular systems Gaussian elimination Q Systems of linear equations - reminder • Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems • Diagonal systems • Triangular systems 9 Gaussian elimination O Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Diagonal systems Triangular systems Gaussian elimination Diagonal systems The simplest linear system is an 0 0 a.22 0 0 ... with obvious solution x = [b;/a,;];. *1 Di = b2 Xn. bn 1 function x = diagsolve(A, b) 2 % Solve A x = b for a diagonal matrix A. 3 d = diag(A); 4 if any(d == 0), error('A is singular!'), end 5 x = b . / d; 6 return Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Outline Diagonal systems Triangular systems Gaussian elimination 9 Systems of linear equations - reminder • Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems • Triangular systems 9 Gaussian elimination O Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Diagonal systems Triangular systems Gaussian elimination aiixi + ai2x2 + ai3x3 = fc>i 822X2 + 823X3 = £>2 a33x3 = t>3 which is equivalent to aiixi =öi 2X2 3X3 a22X2 = t»2 -a23x3 a33x3 = b3 IBA Vlad BJ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Triangular systems Diagonal systems Triangular systems Gaussian elimination • A is lower triangular if a,] = 0 for / < j or upper triangular if a,y = 0 for / > j a solution is obtained by back-substitution: for an ai2 ai3 • • ai„ 0 a22 a23 • • a2n A = 0 0 333 • • a3n 0 0 0 . ann Xn — bnl 3r, b<- Yj ai'X' i=/+1 /a,;,for /' = n - 1,n - 2,..., 1 Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Back-substitution algorithm Diagonal systems Triangular systems Gaussian elimination (not vectorized!) Algorithm: Back-substitution algorithm for j = n to 1 do if ass = 0 then stop; end if Xj«- bj/af for / = 1 to y - 1 do b j <- bi - auxr, end for end for IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Exercise Diagonal systems Triangular systems Gaussian elimination a derive the forward substitution method for lower triangular matrices 9 implement in Matlab the functions fwsoive and bksoive for forward and backward substitution IBA Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Diagonal systems Triangular systems Gaussian elimination Elementary elimination matrices Goal Find tranformations of nonsingular matrices that would lead to triangular systems. J Example: let z = [z-\,z2]T with z-\ * 0, then 1 0 Z1 Z1 -Z2/Z! 1 Z2 0 use linear combinations or rows IBA Bi7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . 3 _ . , Triangular systems Special cases _ . ..... _ . ,. Gaussian elimination Examples and applications In general, 1 .. 0 0 . . 0 z{ 0 .. 1 0 . . 0 Zk Zk 0 .. • -mjf+1 1 . . 0 0 o .. • -m„ 0 . . 1 - Zn .0. where m,- = Zj/Zk, for / = k + 1,..., n. 9 pivot: Zk 9 Gaussian transformation or elementary elimination transformation: Mk Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . , _ . , Triangular systems Special cases _ . ..... _ . ,. Gaussian elimination Examples and applications Properties of the Gaussian transformation 9 Mk is nonsingular (it is lower triangular, full rank matrix) 9 Mjf = I - me[, where m = [0,..., 0, rrik+:,..., mn]7 and is the k-th column of the identity matrix 9 M^1 = I + me7: just the sign is changed for the inverse. Denote Lk = Mk 9 if M7 = I — tej,j > k, then MkMj = I - me7 + te7, so the result is sort of "union" of the two matrices. Note that the order of multiplication is important. • a similar result holds for the inverses ^ IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Outline Diagonal systems Triangular systems Gaussian elimination 9 Systems of linear equations - reminder • Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems • Triangular systems • Gaussian elimination O Special cases » Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Diagonal systems Triangular systems Gaussian elimination Gaussian elimination a transform the system Ax = b into a triangular system: o choose Mi with an as pivot to eliminate the 1st column below a-n. The new system is M-|Ax = M-|b. The solution stays the same. • next choose M2 with a22 as pivot to eliminate the 2nd colum below a22. The new system is M2N\^ Ax = IV^Mtb. The solution stays the same. «... until we get a triangular system a solve the system Mn_i ... Mi Ax = Mn_i ... Mi b by back-substitution IBA Vlad BJ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications LU factorization Diagonal systems Triangular systems Gaussian elimination o let M = Mn-i ... Mi and L = M"1 • L = (Mn_1...M1)"1 =M71...M-11 =L1...Ln_1 which is unit lower triangular. • by design, U = MA is upper triangular o then A = M"1U = LU with L lower triangular and U upper triangular 9 Gaussian elimination is a factorization of a matrix as a product of two triangular matrices: LU factorization 9 LU factorization is unique up to a scaling factor of diagonal scaling of factors Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ Gaussian elimination Examples and applications if A is factorized into LU, the system becomes LUx = b and is solved by forward-substitution (reverse order of backward s.) in lower triangular system Ly = b followed by back-substitution in Ux = y Gaussian elimination and LU factorization express the same solution process Matlab example: 1 » A = [0 1 1; 2 -1 -1; 1 1 -1] ; b 2 >> [L, U] = lu (A) ; 3 >>y=L\b; x=U\y; [2 0 1] IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ Gaussian elimination Examples and applications Note: det(A) = det(L) det(U) if at any stage, the leading entry on the diagonal is zero -> cannot choose the pivot -> interchange the row with some row below with a non-zero pivot if there is no way to choose a proper pivot, the matrix U will be singular but the factorization can be performed! the back-substitution will fail however. IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Experiment Diagonal systems Triangular systems Gaussian elimination (from C. Van Loan, "Introduction to scientific computing") Consider the system 6 1 X1 "1 +6 1 1 X2 2 with the solution [1 1]T. Write a Matlab code to solve it using LU factorization, for e= 10"2,10"4,..., 10"18. Discuss the results! Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ ^ .. Gaussian elimination Examples and applications Another application of LU decomposition Consider you have to compute the scalar a = zTA"1b e R, with z,belN and A e Rnx" nonsingular. But x = A"1b is the solution of the linear system Ax = b. So, you should use LU decomposition, compute x and then a = zTx. In Matlab: 1 [L,U] = lu(A); 2 y = L\ b; x = U\ y; 3 alpha = z' * x; IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Improving stability Diagonal systems Triangular systems Gaussian elimination 9 chose the pivot to minimize error propagation • choose the entry of largest magnitude on or below the diagonal as pivot • this is called partial pivoting 9 each Mk is preceded by a permutation matrix to interchange rows • still MA = U, but M = Mn-!Pn_i . ..M^ • L = M"1 is triangular, but not necessarily lower triangular 9 in general (Pn_1...P1)A = PA = LU 9 check again previous Matlab example • try [L, U, P] = lu (A) ; fljj? Vlad Bi7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ Gaussian elimination Examples and applications • if the pivot is sought as the largest entry in the entire unreduced submatrix, then you have complete pivoting » requires permutations or rows AND columns » there are 2 permutations matrices, P, Q, such that PAQ = LU » better numerical stability, but much more expensive in computation • in general, only partial pivoting is used with Gaussian elimination • in Matlab is implemented only for sparse matrices. See help iu Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Pivoting is not required if: Diagonal systems Triangular systems Gaussian elimination o the matrix is diagonally dominant: n 2 |a,y| < \ajj\, j = 1,...,n 9 the matrix is symmetric positive definite: A = AT and xTAx > 0, Vx * 0 Examples of symmetric positive (semi-)definite matrices from practice? Vlad Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Diagonal systems Triangular systems Gaussian elimination Residuals 9 r = b - Ax where x was obtained by Gaussian elimination • it can be shown that where E is the backward error in data matrix: (A + E)x = b and p = max(uy)/ max(a,y) is the growth factor 9 without pivoting, p is unbounded so the algorithm is unstable 9 with partial pivoting, p < 2"~1 • in practice, p * 1, so jJJfL <> nemac„ IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Residuals, cont'd Diagonal systems Triangular systems Gaussian elimination • Gaussian elimination with partial pivoting yields small relative residuals, regardless of the conditioning 9 however, computed solution is close to real solution only if the system is well-conditioned a yet a smaller growth factor can be obtained with complete pivoting, but the extra cost may not be worth IBA Bi7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ Gaussian elimination Examples and applications Example: in a 3-digit decimal arithmetic, solve 0.641 0.242 x^ 0.883 0.321 0.121 x2 0.442 9 the exact solution is [1 1]T 9 the Gaussian elimination leads to x = [0.782 1.58]T * the exact residual is r = [-0.000622 - 0.000202]7 -> as small as can be expected with 3 digits precision 9 the error is large: ||x - x|| = 0.6196 which is « 62% relative error! 9 this is because of ill-conditioning, cond(A) > 4000 Vlad Bi7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ Gaussian elimination Examples and applications What happend? The Gaussian elimination led to 0.641 0.242 *1 0.883 0 0.000242 X2 -0.000383 so x2 was the result of the division of quantities below emach, yielding an arbitrary result. The x-\ is computed to satisfy the 1st eq., resulting in small residual but large error. IBA Bi7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ ^ .. Gaussian elimination Examples and applications Implementation and complexity The general form of the Gaussian elimination is for / do for j do for k do a,y <- a,y - (aik/akk)akj end for end for end for o order of the loops is not important (for the final result) • ...but, depending on the memory storage, they have different performance Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ ^ .. Gaussian elimination Examples and applications Implementation and complexity (cont'd) • there are about n3/3 floating-point operations -> the complexity is 0(n3) o the forward-/back-substitutions require about n2 multiplications and n2 additions (for a single b) • if you try to invert A, x = A-1 b, you need n3 operations -» 3x more the LU factorization • inversion is less precise: difference between 3"1 x 18 and 18/3 in fixed-precision arithmetic • matrix inversion is convenient in formulas, but in practice you do factorizations! 9 Ex: A"1B should use LU factorization of A and then forward-and back-substitutions with columns of B mL IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Gauss-Jordan elimination Diagonal systems Triangular systems Gaussian elimination o idea: for each element of the diagonal, eliminate all the elements below AND above in the column using combinations of rows o the elimination matrix has the form 1 . . 0 0 . . 0 0 . . 1 -mit-1 0 . . 0 0 . . 0 1 0 . . 0 0 . . 0 -mfc+i 1 . . 0 0 . . 0 -m„ 0 . . 1 where m,- = a,/a/< for / = 1,..., n <» do the same to the right hand side term, too Vlad Bi7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ ^ .. Gaussian elimination Examples and applications Gauss-Jordan elimination, cont'd » the result is a diagonal matrix on Ihs a the solution is obained by dividing the entries on the transformed rhs by the terms of the diagonal • it requires n3/2 multiplications and the same number of additions —> 50% more expensive than LU decomposition o despite being more expensive, it is sometimes preferred to LU decomposition for parallel implementations • if the rhs is initialized with an identity matrix, after G-J elimination the rhs becomes A-1 Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ ^ .. Gaussian elimination Examples and applications Solving series of similar problems 9 idea: try to reuse as much as possible from previous computations • if only rhs changes, LU decomposition does not have to be recomputed « if A suffers only rank one changes, one can still use pre-computed A-1 (Sherman-Morrison formula): (A-uv7)"1 = A"1 +A"1u(1 -vTA-1u)-1vTA-1 a this has a complexity of 0(n2) compared to 0(n3) that is needed by a new inversion Vlad BÍ7740: Scientific computing Systems of linear equations - reminder _ . . Diagonal systems Solving linear systems _. . ' , _ . , Triangular systems Special cases _ . ,. . 4. _ Gaussian elimination Examples and applications For a modified equation, (A - uvT)x = b the solution is x = A"1b +A"1u(1 -vTA-1u)-1vTA-1b and is solved by the following procedure 9 solve Az = u, so z = A-1 u 9 solve Ay = b, so y = A"1 b • compute x = y + ((vTy)/(1 - vTz))z If A is already factored, this approach has a complexity 0(n2) Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Comments on scaling Diagonal systems Triangular systems Gaussian elimination » theoretically, multiplying the terms on diagonal of A and corresponding entries of b would not change the solution • in practice, it affects conditioning, choice of pivot and, by consequence, accuracy • Example: 1 0 X1 1 0 e X2. e is ill-conditioned for small e, since cond(A) = 1 /e. It becomes well-conditioned if the second equation is multiplied by 1 /e. Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Iterative refinements Diagonal systems Triangular systems Gaussian elimination » let x0 be the approximate solution to Ax = b and r0 = b - Ax0 be the corresponding residual 9 let then z0 be the solution to Az = r0 9 an improved approximate solution is then Xi = x0 + z0 HOMEWORK: prove that Axt = b 9 repeat until convergence 9 the process needs higher precision for computing a useful residual 9 not often used, but sometimes useful Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems _ _ . , Symmetric positive definite systems Special cases Examples and applications Outline 9 Systems of linear equations - reminder • Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems 9 Diagonal systems • Triangular systems 9 Gaussian elimination Q Special cases 9 Symmetric positive definite systems ) Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems _ _ . , Symmetric positive definite systems Special cases Examples and applications Special forms of linear systems For some special cases of A storage and computation time can be saved. For example, if A is • symmetric: A = AT, a,y = ay; for all i,j 9 positive definite: zTAz > 0, Vz + 0 • band diagonal: a,y = 0 if |/ - j\ > j3, where 13 is the bandwidth • sparse: most of the elements of A are zero IBA Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems _ _ . , Symmetric positive definite systems Special cases Examples and applications Outline 9 Systems of linear equations - reminder • Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems 9 Diagonal systems • Triangular systems 9 Gaussian elimination Q Special cases • Symmetric positive definite systems ) Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems _ _ . , Symmetric positive definite systems Special cases Examples and applications Symmetric positive definite systems • Cholesky decomposition: A = LLT where L is lower triangular. • A admits a Cholesky decomposition if and only if it is symmetric positive definite » if the decomposition exists, it is unique IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems _ _ . , Symmetric positive definite systems Special cases Examples and applications Cholesky decomposition algorithm with overwriting of A Algorithm: Cholesky decomposition algorithm for j = 1 to n do for k = 1 to j'•- 1 do for / = j to n do *~ &ij ~ &ik3jk, end for end for aii *~ yßs' for k = j'•+ 1 to n do &kj ^~ &kjl 3//! end for end for IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems _ _ . , Symmetric positive definite systems Special cases Examples and applications Cholesky decomposition - properties o does not need pivoting to maintain stability a only n3/6 multiplications and n3/6 additions are required o for the algorithm presented, only the lower triangle of A is modified, and can be restored, if needed, from the upper triangle » requires about half the computations and half of the memory compared with LU factorization o there are variations of Cholesky decomposition for banded matrices, for positive semi-definite matrices (semi-Cholesky decomposition) and for symmetric indefinite matrices Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems _ _ . , Symmetric positive definite systems Special cases Examples and applications Suggestions of methods to use If A is a real dense square matrix... » ...use LU decomposition with partial pivoting: A = PLU » ...and is a band matrix, use LU decomposition with pivoting and row interchanges 9 ...and is tridiagonal, use Gaussian elimination 9 ...and is symmetric positive definite, use Cholesky decomposition • ...and is symmetric tridiagonal, use special Cholesky with pivoting, A = LDLT 9 ...and is symmetric indefinite, use special Cholesky In Matlab, check the functions: choi, ichoi, ldi, lu, iiu. IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Outline 9 Systems of linear equations - reminder • Norms • Linear systems • Conditioning o Accuracy Q Solving linear systems 9 Diagonal systems • Triangular systems 9 Gaussian elimination Q Special cases 9 Symmetric positive definite systems Q Examples and applications IBA Vlad BÍ7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Polynomial interpolation o a function p(x) interpolates a set of points {(*;> yi)\i = 0,..., N} if it satisfies y, = p(x,) for all / = 0,..., N. » this leads to a system of N + 1 equations. If p(x) is a polynomial of degree M, p(x) = a^xM H-----h aix + a0, the system is of the form a0 + aix0 + --- + a/Mx^ = yo a0 + aixw h-----h aMx^ = yN where the unknowns are ao,..., au-9 if M = N -> Vandermonde matrix • in Matlab check the functions poiyf it and poiyvai <» write the Matlab function to solve the interpolation problem for M = N. Do NOT use Matlab's own functions for interpolation! IBA Vlad Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications 1D Poisson problem A two-point boundary problem, -u"(x) = y(x), xe[0,1], u(0) = u(1) = 0, where y is a given continuous function on [0,1]. If y cannot be integrated exactly, approximate solutions are sought. Using finite differences, u'(x) = lim u(x+|)-u(x-|) h->0 h „, , .. u(x + h) - 2u(x) + u(x - h) u (x) = lim —---Vj!---- Divide the interval [0,1] in m + 1 equal subintervals of length h = 1 /(m + 1) and let x,- = ih be the limits of these subtintervals, / = 0,..., m + 1. Vlad Bi7740: Scientific computing Systems of linear equations - reminder Solving linear systems Special cases Examples and applications Denote y(x,) = y(ih) = y, and u(x,) = u(ih) = u,-. Then, the problem becomes Uf+1 - 2Uj + L//-1 h2 yi, i = 1,..., m, u0 = um+i = 0. This can be written as a linear system: Tu -1 2 -1 U1 y^ u2 y2 = h2 ym-i . um where the matrix T is a Toeplitz matrix. The system can be solved using the Levinson algorithm - see levinson function in Matlab. Vlad Bi7740: Scientific computing