51 Matrix Factorizations and Direct Solution of Linear Systems s of Linear Systems. Linear Syster Decomposition... 51.4 Symmetric Factorization 51-19 51.5 Orthogonalization and the QR Factoriation- References 51-20 The need to solve sys of linear ence engin ession "direct solution of lir ear ms”refers and dimension of the coefficient matrix.The evolution of computers has and continues to influence the development of these strategies and has also fostered particular styles of per- turbation analysis suited to illuminating their behavior.Some general themes have become dominant,as a result;others have been pushed aside.For example,Cramer's Rule may be properly thought of as a direct solution strategy for solving linear systems:however as normally manifested,it requires a much larger number of arithmetic operations than Gauss elimination and is generally much more susceptible to the deleterious effects of rounding Most current approaches for the direct solution of a linear system, Ax b,are patterne after Gauss elimination and favor an initial phase that partially decouples the system of equations:ze into tria 1 Find invertible matrices such that Se ..S2SA=U is triangular;then 2 a modifod right hand side SS b:and then 3.Determine the solution set to the triangular system Ux=y. The matrices S,S2,...S.are typically either row permutations of lower triangular matrices (Gauss transformations)or unitary matrices.In either case,inverses s are readily available.Evidently,A can be written as A=NU,where N=(5...S2S1)-1.A solution framework may be built around the availability of decompositions such as this: 1.Find a decomposition A =NU such that U is triangular and Ny =b is easily 1 =b;ther set to the triangular system Uxy. 51-1
51 Matrix Factorizations and Direct Solution of Linear Systems Christopher Beattie Virginia Polytechnic Institute and State University 51.1 Perturbations of Linear Systems ................... 51-2 51.2 Triangular Linear Systems.......................... 51-5 51.3 Gauss Elimination and LU Decomposition ....... 51-7 51.4 Symmetric Factorizations ........................... 51-13 51.5 Orthogonalization and the QR Factorization .... 51-16 References..................................................... 51-20 The need to solve systems of linear equations arises often within diverse disciplines of science, engineering, and finance. The expression “direct solution of linear systems” refers generally to computational strategies that are able to produce solutions to linear systems after a predetermined number of arithmetic operations that depends only on the structure and dimension of the coefficient matrix. The evolution of computers has and continues to influence the development of these strategies and has also fostered particular styles of perturbation analysis suited to illuminating their behavior. Some general themes have become dominant, as a result; others have been pushed aside. For example, Cramer’s Rule may be properly thought of as a direct solution strategy for solving linear systems; however as normally manifested, it requires a much larger number of arithmetic operations than Gauss elimination and is generally much more susceptible to the deleterious effects of rounding. Most current approaches for the direct solution of a linear system, Ax = b, are patterned after Gauss elimination and favor an initial phase that partially decouples the system of equations: zeros are introduced systematically into the coefficient matrix, transforming it into triangular form; the resulting triangular system is easily solved. The entire process can be viewed in this way: 1. Find invertible matrices {Si} ρ i=1 such that Sρ . . . S2S1A = U is triangular; then 2. Calculate a modified right-hand side y = Sρ . . . S2S1b; and then 3. Determine the solution set to the triangular system Ux = y. The matrices S1, S2, . . . Sρ are typically either row permutations of lower triangular matrices (Gauss transformations) or unitary matrices. In either case, inverses are readily available. Evidently, A can be written as A = NU, where N = (Sρ . . . S2S1) −1 . A solution framework may be built around the availability of decompositions such as this: 1. Find a decomposition A = NU such that U is triangular and Ny = b is easily solved; 2. Solve Ny = b; then 3. Determine the solution set to the triangular system Ux = y. 51-1
51-2 Handbook of Linear Algebra 51.1 Perturbations of Linear Systems s a sma pers degrade t the accuracy o 8 stems ofter befor within the pute nted within the internal foating errors that may be introduced in the course of computation often may be viewed in effectively as an additional contribution to this initial rer linear system for which a solution is computed will deviate slightly from the "true"linear system and it becomes of critical interest to determine whether such deviations will have a significant effect on the accuracy of the final computed result. Definitions: perturbations and A and b,respectively,the solution ix Cn satisfies the b(pre then that the pertubed system perturbed systen m(A+6A(文 x)= sistent). iated with as an approximate solution to the linear system Ax b is defined as r()=b-A%. For any文∈C",the ass ociated(norm-wise)relative backward error of the linear system Ax =b (with respect to the the p-norm,for 1sps oo)is there exist 6A.&b such that p(A,b:)=min (A+6A)=b+6b with I5Ap≤ellAll? For any文∈Cm,the associated co mponent-wise relative backward error of the linear system Ax =b is there exist 6A.6b such that (A+6A)=b+6b with 16A≤E4 w(A,b;刘=mi bl≤ebl where the absolute values and inequalities applied to vectors and matrices are interpreted component- wise:for example,BI A]means for all index pairs i,j. The (norm-wise)condition number of the linear system Ax =b(with respect to thethe p-norm.for1<p<oo）is The matrix condition number of A (with respect to the the p-norm,for 1 spsoo)is p(A)=IlAllpllA-'llp The Skeel condition number of the linear system Ax=b is comd(A,刘=AA三 xoo The Skeel matrix condition number is cond(A)=A-Al lloo
51-2 Handbook of Linear Algebra 51.1 Perturbations of Linear Systems In the computational environment afforded by current computers, the finite representation of real numbers creates a small but persistent source of errors that may on occasion severely degrade the overall accuracy of a calculation. This effect is of fundamental concern in assessing strategies for solving linear systems. Rounding errors can be introduced into the solution process for linear systems often before any calculations are performed — as soon as data are stored within the computer and represented within the internal floating point number system of the computer. Further errors that may be introduced in the course of computation often may be viewed in aggregate effectively as an additional contribution to this initial representation error. Inevitably, the linear system for which a solution is computed will deviate slightly from the “true” linear system and it becomes of critical interest to determine whether such deviations will have a significant effect on the accuracy of the final computed result. Definitions: Let A ∈ C n×n be a nonsingular matrix, b ∈ C n , and then denote by xˆ = A −1b the unique solution of the linear system Ax = b. Given data perturbations δA ∈ C n×n and δb ∈ C n to A and b, respectively, the solution perturbation, δx ∈ C n satisfies the associated perturbed linear system (A + δA)(xˆ + δx) = b + δb (presuming then that the perturbed system is consistent). For any x˜ ∈ C n , the residual vector associated with x˜ as an approximate solution to the linear system Ax = b is defined as r(x˜) = b − Ax˜. For any x˜ ∈ C n , the associated (norm-wise) relative backward error of the linear system Ax = b (with respect to the the p-norm, for 1 ≤ p ≤ ∞) is ηp(A, b; x˜) = min ε there exist δA, δb such that (A + δA)x˜ = b + δb with kδAkp ≤ εkAkp kδbkp ≤ εkbkp . For any x˜ ∈ C n , the associated component-wise relative backward error of the linear system Ax = b is ω(A, b; x˜) = min ε there exist δA, δb such that (A + δA)x˜ = b + δb with |δA| ≤ ε|A| |δb| ≤ ε|b| , where the absolute values and inequalities applied to vectors and matrices are interpreted componentwise: for example, |B| ≤ |A| means |bij| ≤ |aij| for all index pairs i, j. The (norm-wise) condition number of the linear system Ax = b (with respect to the the p-norm, for 1 ≤ p ≤ ∞) is κp(A, xˆ) = kA −1 kp kbkp kxˆkp . The matrix condition number of A (with respect to the the p-norm, for 1 ≤ p ≤ ∞) is κp(A) = kAkpkA −1 kp. The Skeel condition number of the linear system Ax = b is cond(A, xˆ) = k |A −1 | |A| |xˆ|k∞ kxˆk∞ . The Skeel matrix condition number is cond(A) = k |A −1 | |A| k∞
Matrix Factorizations and Direct Solution of Linear Systems 51-3 Facts:[Hig02].[SS90] 1.For any C",is the exact solution to any one of the following families of perturbed linear systems (A+6A)=b+6b. 2.(Rigal-Gaches Theorem)For anyC", (A,b:刘=Ap,+Ib可 with b+ob with data perturbations as in (1)andA and as a result 5AL=lbl=,(A,b:x刘. IAllP bllp 3.(Oettli-Prager Theorem)For any C", (A,b:刘=maxA闵+b If Ddiag)and Dding(n()theis an eact solution ear system (A +6A)=b+6b with 6A D1A D2 a to the 5A≤w(A,b:x)4andbl≤w(A,b:x)l4 and no smaller constant can be used in place of w(A,b;). J5A业A+6A is singular 西=min{ In particular,the perturbed coefficient matrix A+A is nonsingular if 5.1≤(A,x刘≤p(A)and1≤cond(A,x)≤cond(A)≤K(A), b p 8.If 6b =0 and A+6A is nonsingular,then 是≤陆 l18xllp 1 .fl6A≤d4bl≤dblp.ade<④then x2 2eKP(A)） 1-EKp(A)
Matrix Factorizations and Direct Solution of Linear Systems 51-3 Facts: [Hig02], [SS90] 1. For any x˜ ∈ C n, x˜ is the exact solution to any one of the following families of perturbed linear systems (A + δAθ)x˜ = b + δbθ, where θ ∈ C, δbθ = (θ − 1) r(x˜), δAθ = θ r(x˜)y˜ ∗ , and y˜ ∈ C n is any vector such that y˜ ∗x˜ = 1. In particular, for θ = 0, δA = 0 and δb = −r(x˜); for θ = 1, δA = r(x˜)y˜ ∗ and δb = 0. 2. (Rigal–Gaches Theorem) For any x˜ ∈ C n, ηp(A, b; x˜) = kr(x˜)kp kAkpkx˜kp + kbkp . If y˜ is the dual vector to x˜ with respect to the p-norm (y˜ ∗x˜ = ky˜kq kx˜kp = 1 with 1 p + 1 q = 1), then x˜ is an exact solution to the perturbed linear system (A + δAθ˜)x˜ = b + δbθ˜ with data perturbations as in (1) and ˜θ = kAkpkx˜kp kAkpkx˜kp+kbkp , and as a result kδAθ˜kp kAkp = kδbθ˜kp kbkp = ηp(A, b; x˜). 3. (Oettli–Prager Theorem) For any x˜ ∈ C n, ω(A, b; x˜) = max i |ri | (|A| |x˜| + |b|)i . If D1 = diag ri (|A| |x˜| + |b|)i and D2 = diag(sign(x˜)i), then x˜ is an exact solution to the perturbed linear system (A + δA)x˜ = b + δb with δA = D1 |A| D2 and δb = −D1 |b| |δA| ≤ ω(A, b; x˜)|A| and |δb| ≤ ω(A, b; x˜)|A| and no smaller constant can be used in place of ω(A, b; x˜). 4. The reciprocal of κp(A) is the smallest norm-wise relative distance of A to a singular matrix, i.e., 1 κp(A) = min kδAkp kAkp A + δA is singular . In particular, the perturbed coefficient matrix A + δA is nonsingular if kδAkp kAkp < 1 κp(A) . 5. 1 ≤ κp(A, xˆ) ≤ κp(A) and 1 ≤ cond(A, xˆ) ≤ cond(A) ≤ κ∞(A). 6. cond(A) = min {κ∞(DA)| D diagonal}. 7. If δA = 0, then kδxkp kxˆkp ≤ κp(A, xˆ) kδbkp kbkp . 8. If δb = 0 and A + δA is nonsingular, then kδxkp kxˆ + δxkp ≤ κp(A) kδAkp kAkp . 9. If kδAkp ≤ kAkp, kδbkp ≤ kbkp, and < 1 κp(A) , then kδxkp kxˆkp ≤ 2 κp(A) 1 − κp(A)
51-4 Handbook of Linear Algebra 10.If15A≤dAl,l5bl≤ebl,ande<eomd④then s x Examples: [1000999] 0A1= 「-9989991 1.Let A= 999998 999 1000 Then IlAll =A-=1999 so that K1(A)≈3.996×10.Consider associated with a solution&= 1997 Aperturbation of the right-hand side b.01 「-0.011 constitutes a relative change in the right. hand side of≈50o5×10-6 procperbo文+6x=- 「20.971 -18.99 constituting a relative change=19.9820=(A).The bound determined by the condition number is very nearly achieved.Note that the same perturbed solution+6x could be produced by a change in the coefficient matrix -0.01 6A==- ,学aab and subdiagonal entries equal to 1 (associated with a centered difference approximation to the second derivative).Let b be a vector with a quadratic variation in entries bk=(k-1)(100-k)/10,000, Then 2(4,)≈1，but2(A)≈4.1336×103 Since the elements of b do not have an exact binary representation,the linear system that is presented to any computational algorithm will be Ax b+6b with llob2 ellbl2 where e is the unit roundoff error.For example,if the linear system data is stored in IEE single precision format The matrix condition num (A),would yield bound of (6 x 10 1336 ng the digit furt 2 stantially 08t o th But in fact this clus of the orithm that is used true back ht be if s s substantially if the right-hand side is changed to bk=(-1)(k-1)(100-)/10.000, which only introduces a sign variation in b.In this case .(A.文)≈(A).and the com ponents of the computed solution can be expected to lose about 4 significant digits purely on
51-4 Handbook of Linear Algebra 10. If |δA| ≤ |A|, |δb| ≤ |b|, and < 1 cond(A) , then kδxk∞ kxˆk∞ ≤ 2 cond(A, xˆ) 1 − cond(A) . Examples: 1. Let A = " 1000 999 999 998# so A −1 = " −998 999 999 −1000# . Then kAk1 = kA −1 k1 = 1999 so that κ1(A) ≈ 3.996 × 106 . Consider b = " 1999 1997# associated with a solution xˆ = " 1 1 # . A perturbation of the right-hand side δb = −0.01 0.01 constitutes a relative change in the righthand side of kδbk1 kbˆk1 ≈ 5.005 × 10−6 yet it produces a perturbed solution xˆ + δx = 20.97 −18.99 constituting a relative change kδxk1 kxˆk1 = 19.98 ≤ 20 = κ1(A) kδbk1 kbk1 . The bound determined by the condition number is very nearly achieved. Note that the same perturbed solution xˆ +δx could be produced by a change in the coefficient matrix δA = ˜ry˜ ∗ = − " −0.01 0.01# h 1 39.96 − 1 39.96 i = (1/3996) " 1 −1 −1 1# constituting a relative change kδAk1 kAk1 ≈ 2.5 × 10−7 . Then (A + δA)(xˆ + δx) = b. 2. Let n = 100 and A be tridiagonal with diagonal entries equal to −2 and all superdiagonal and subdiagonal entries equal to 1 (associated with a centered difference approximation to the second derivative). Let b be a vector with a quadratic variation in entries bk = (k − 1)(100 − k)/10,000. Then κ2(A, xˆ) ≈ 1, but κ2(A) ≈ 4.1336 × 103 . Since the elements of b do not have an exact binary representation, the linear system that is presented to any computational algorithm will be Ax = b + δb with kδbk2 ≤ kbk2, where is the unit roundoff error. For example, if the linear system data is stored in IEEE single precision format, ≈ 6 × 10−8 . The matrix condition number, κ2(A), would yield a bound of (6 × 10−8 )(4.1336 × 103 ) ≈ 2.5 × 10−4 anticipating the loss of more than 4 significant digits in solution components even if all computations were done on the stored data with no further error. However, the condition number of the linear system, κ2(A, xˆ), is substantially smaller and the predicted error for the system is roughly the same as the initial representation error ≈ 6 × 10−8 , indicating that the solution will be fairly insensitive to the consequences of rounding of the right-hand side data—assuming no further errors occur. But, in fact, this conclusion remains true even if further errors occur, if whatever computational algorithm that is used produces small backward error, as might be asserted if, say, a final residual satisfies krk2 ≤ O() kbk2. This situation changes substantially if the right-hand side is changed to bk = (−1)k (k − 1)(100 − k)/10,000, which only introduces a sign variation in b. In this case, κ2(A, xˆ) ≈ κ2(A), and the components of the computed solution can be expected to lose about 4 significant digits purely on
Matrix Factorizations and Direct Solution of Linear Systems 51-5 the basis of errors that are made in the initial representation.Additional errors made in the course of the computation can hardly be expected to improve this situation. 51.2 Triangular Linear Systems Syseamsofieroquatiosorw时iahthotmnosnmennlforometntnci识 Such may be an be solved w. eff tems are e the for of lin they often constitute an interm Definitions: = r an upper triangular matrix (ti=0 for i>j)or a Facts:Hig02],[GV96] 1.[GV96,pp.88-90 Algorithm 1:Row-wise forward substitution for solving lower triangular system Input:L=[∈R"xn with(为=0fork<jb∈R Output:solution vector x ER"that satisfies Lx b x1←-b1/011 for k=2 ton 工k←(k-Lk.1k-1·工1k-1)/k.内 Algorithm 2:Column-wise back substitution for solving upper triangular system Input:U =ue nxn with =0 for k>i: for k=n down to 2 in steps of-1, 工k←b/ukk b1:k-1+b1:k-1-kU1:k-1.k x1←b/u1,1 2.Algorithm 1 involves as a core calculation dot products of portions of coefficient matrix rows with portions of the emerging solution vector.This can incur a performance penalty for large n from accumulation of dot products using a scalar recurrence.A "column-wise"reformulation may have better performance for large n.Algorithm 2 is such a "column-wise"formulation for upper triangular systems. 3.An efficient and reliable implementation for the solution of triangular systems is offered as part of the standard BLAs software library in xTRSz (see Chapter 92), where x=S,D,C,or Z according to whether data are single or double precision real,or single or double precision complex floating point numbers,respectively,and z=V or M according to whether a single system of equations is to be solved or multiple systems (sharing the same coefficient matrix)are to be solved,respectively
Matrix Factorizations and Direct Solution of Linear Systems 51-5 the basis of errors that are made in the initial representation. Additional errors made in the course of the computation can hardly be expected to improve this situation. 51.2 Triangular Linear Systems Systems of linear equations for which the unknowns may be solved for one at a time in sequence may be reordered to produce linear systems with triangular coefficient matrices. Such systems can be solved both with remarkable accuracy and remarkable efficiency. Triangular systems are the archetype for easily solvable systems of linear equations. As such, they often constitute an intermediate goal in strategies for solving linear systems. Definitions: A linear system of equations Tx = b with T ∈ C n×n (representing n equations in n unknowns) is a triangular system if T = [tij] is either an upper triangular matrix (tij = 0 for i > j) or a lower triangular matrix (tij = 0 for i < j). Facts: [Hig02], [GV96] 1. [GV96, pp. 88–90] Algorithm 1: Row-wise forward substitution for solving lower triangular system Input: L = [`ij] ∈ R n×n with `kj = 0 for k < j; b ∈ R n Output: solution vector x ∈ R n that satisfies Lx = b x1 ← b1/`1,1 for k = 2 to n xk ← (bk − Lk,1:k−1 · x1:k−1)/`k,k Algorithm 2: Column-wise back substitution for solving upper triangular system Input: U = [uij] ∈ R n×n with ukj = 0 for k > j; b ∈ R n Output: solution vector x ∈ R n that satisfies Ux = b for k = n down to 2 in steps of −1, xk ← bk/uk,k b1:k−1 ← b1:k−1 − xkU1:k−1,k x1 ← b1/u1,1 2. Algorithm 1 involves as a core calculation dot products of portions of coefficient matrix rows with portions of the emerging solution vector. This can incur a performance penalty for large n from accumulation of dot products using a scalar recurrence. A “column-wise” reformulation may have better performance for large n. Algorithm 2 is such a “column-wise” formulation for upper triangular systems. 3. An efficient and reliable implementation for the solution of triangular systems is offered as part of the standard blas software library in xTRSz (see Chapter 92), where x=S, D, C, or Z according to whether data are single or double precision real, or single or double precision complex floating point numbers, respectively, and z=V or M according to whether a single system of equations is to be solved or multiple systems (sharing the same coefficient matrix) are to be solved, respectively