98 Chapter 2. Solution of Linear Algebraic Equations x[i]=sum/p[i]; 22 A typical use of choldc and cholsl is in the inversion of covariance matrices describing the fit of data to a model;see,e.g.,15.6.In this,and many other applications,one often needs L-1.The lower triangle of this matrix can be efficiently found from the output of choldc: for(1=1;1<n;1+)f a[1][i]=1.0/p[i]; for(j=1+1;j<n;j++)( sum=0.0; for (k=i;k<j;k++)sum -a[j][k]*a[k][i]; a[i][i]=sum/p[j]; 18881892 NUMERICAL CITED REFERENCES AND FURTHER READING: Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.ll of Handbook for Automatic Com- putation (New York:Springer-Verlag),Chapter 1/1. Gill,P.E.,Murray,W.,and Wright,M.H.1991,Numerical Linear Algebra and Optimization,vol.1 (Redwood City,CA:Addison-Wesley),84.9.2. p2今 Press. Dahlquist,G.,and Bjorck,A.1974,Numerica/Methods (Englewood Cliffs,NJ:Prentice-Hall). 65.3.5. Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),84.2. IENTIFIC 2.10 QR Decomposition There is another matrix factorization that is sometimes very useful,the so-called OR (ISBN decomposition, Numer A=0R (2.10.1) Here R is upper triangular,while Q is orthogonal,that is, 431086 Q.Q=1 (2.10.2) where QT is the transpose matrix of Q.Although the decomposition exists for a general rectangular matrix,we shall restrict our treatment to the case when all the matrices are square, with dimensions N x N. Like the other matrix factorizations we have met (LU,SVD,Cholesky),OR decompo- sition can be used to solve systems of linear equations.To solve A·x=b (2.10.3) first form Q.b and then solve R.x=QT.b (2.10.4) by backsubstitution.Since OR decomposition involves about twice as many operations as LU decomposition,it is not used for typical systems of linear equations.However,we will meet special cases where OR is the method of choice
98 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). x[i]=sum/p[i]; } } A typical use of choldc and cholsl is in the inversion of covariance matrices describing the fit of data to a model; see, e.g., §15.6. In this, and many other applications, one often needs L−1. The lower triangle of this matrix can be efficiently found from the output of choldc: for (i=1;i<=n;i++) { a[i][i]=1.0/p[i]; for (j=i+1;j<=n;j++) { sum=0.0; for (k=i;k<j;k++) sum -= a[j][k]*a[k][i]; a[j][i]=sum/p[j]; } } CITED REFERENCES AND FURTHER READING: Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag), Chapter I/1. Gill, P.E., Murray, W., and Wright, M.H. 1991, Numerical Linear Algebra and Optimization, vol. 1 (Redwood City, CA: Addison-Wesley), §4.9.2. Dahlquist, G., and Bjorck, A. 1974, Numerical Methods (Englewood Cliffs, NJ: Prentice-Hall), §5.3.5. Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), §4.2. 2.10 QR Decomposition There is another matrix factorization that is sometimes very useful, the so-called QR decomposition, A = Q · R (2.10.1) Here R is upper triangular, while Q is orthogonal, that is, QT · Q = 1 (2.10.2) where QT is the transpose matrix of Q. Although the decomposition exists for a general rectangular matrix, we shall restrict our treatment to the case when all the matrices are square, with dimensions N × N. Like the other matrix factorizations we have met (LU, SVD, Cholesky), QR decomposition can be used to solve systems of linear equations. To solve A · x = b (2.10.3) first form QT · b and then solve R · x = QT · b (2.10.4) by backsubstitution. Since QR decomposition involves about twice as many operations as LU decomposition, it is not used for typical systems of linear equations. However, we will meet special cases where QR is the method of choice
2.10 QR Decomposition 99 The standard algorithm for the OR decomposition involves successive Householder transformations (to be discussed later in $11.2).We write a Householder matrix in the form 1-uu/c where c=uu.An appropriate Householder matrix applied to a given matrix can zero all elements in a column of the matrix situated below a chosen element.Thus we arrange for the first Householder matrix Q,to zero all elements in the first column of A below the first element.Similarly Q2 zeroes all elements in the second column below the second element,and so on up to Q-1.Thus R=Qn-1.Q1·A (2.10.5) Since the Householder matrices are orthogonal, Q=(Qn-1.Q1)-1=Q1…Qm-1 (2.10.6) In most applications we don't need to form Q explicitly;we instead store it in the factored form (2.10.6).Pivoting is not usually necessary unless the matrix A is very close to singular. A general OR algorithm for rectangular matrices including pivoting is given in[1].For square -00 matrices,an implementation is the following: from NUMERICAL RECIPESI 18881892 #include #include "nrutil.h" (Nort 令 void qrdcmp(float **a,int n,float *c,float *d,int *sing) Constructs the QR decomposition of a[1..n][1..n].The upper triangular matrix R is re- turned in the upper triangle of a,except for the diagonal elements of R which are returned in America server computer, University Press. one paper THE d[1..n].The orthogonal matrix Q is represented as a product of n-1 Householder matrices ART Q1...Qn-1,where Qj =1-uj uj/cj.The ith component of uj is zero for i=1,....j-1 while the nonzero components are returned in a[i][j]for i j,...,n.sing returns as for their Programs true (1)if singularity is encountered during the decomposition,but the decomposition is still completed in this case;otherwise it returns false (0). f int i,j,k; float scale,sigma,sum,tau; to dir *sing=0 for (k=1;k<n;k++){ scale=0.0; rectcustsen OF SCIENTIFIC COMPUTING(ISBN for (i=k;i<=n;i++)scale=FMAX(scale,fabs(a[i][k])) 1f(scale==0.0)[ Singular case. *s1ng=1; c[k]=d[k]=0.0; 10621 else Form Q and Qk·A for (i=k;i<=n;i++)a[i][k]/scale; @cambridge.org Numerical Recipes books or 1988-1992 by Numerical Recipes 43108 for (sum=0.0,i=k;i<=n;i++)sum +SQR(a[i][k]); sigma=SIGN(sqrt(sum),a[k][k]); a[k][k]+sigma; c[k]=sigma*a[k][k]; d[k]=-scale*sigma; for (j=k+1;j<=n;j++){ for(sum=0.0,i=k;i<=n;i++)sum+=a[1][k]*a[i][]; (outside North Amer Software. tau=sum/c[k]; for (i=k;i<=n;i++)a[i][j]-tauwa[i][k]; d[n]=a[n][n]; if (d[n]==0.0)*sing=1; The next routine,qrsolv,is used to solve linear systems.In many applications only the part (2.10.4)of the algorithm is needed,so we separate it off into its own routine rsolv
2.10 QR Decomposition 99 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The standard algorithm for the QR decomposition involves successive Householder transformations (to be discussed later in §11.2). We write a Householder matrix in the form 1 − u ⊗ u/c where c = 1 2 u · u. An appropriate Householder matrix applied to a given matrix can zero all elements in a column of the matrix situated below a chosen element. Thus we arrange for the first Householder matrix Q1 to zero all elements in the first column of A below the first element. Similarly Q2 zeroes all elements in the second column below the second element, and so on up to Qn−1. Thus R = Qn−1 ··· Q1 · A (2.10.5) Since the Householder matrices are orthogonal, Q = (Qn−1 ··· Q1) −1 = Q1 ··· Qn−1 (2.10.6) In most applications we don’t need to form Q explicitly; we instead store it in the factored form (2.10.6). Pivoting is not usually necessary unless the matrix A is very close to singular. A general QR algorithm for rectangular matrices including pivoting is given in [1]. For square matrices, an implementation is the following: #include #include "nrutil.h" void qrdcmp(float **a, int n, float *c, float *d, int *sing) Constructs the QR decomposition of a[1..n][1..n]. The upper triangular matrix R is returned in the upper triangle of a, except for the diagonal elements of R which are returned in d[1..n]. The orthogonal matrix Q is represented as a product of n − 1 Householder matrices Q1 ... Qn−1, where Qj = 1 − uj ⊗ uj/cj . The ith component of uj is zero for i = 1,...,j − 1 while the nonzero components are returned in a[i][j] for i = j, . . . , n. sing returns as true (1) if singularity is encountered during the decomposition, but the decomposition is still completed in this case; otherwise it returns false (0). { int i,j,k; float scale,sigma,sum,tau; *sing=0; for (k=1;k<n;k++) { scale=0.0; for (i=k;i<=n;i++) scale=FMAX(scale,fabs(a[i][k])); if (scale == 0.0) { Singular case. *sing=1; c[k]=d[k]=0.0; } else { Form Qk and Qk · A. for (i=k;i<=n;i++) a[i][k] /= scale; for (sum=0.0,i=k;i<=n;i++) sum += SQR(a[i][k]); sigma=SIGN(sqrt(sum),a[k][k]); a[k][k] += sigma; c[k]=sigma*a[k][k]; d[k] = -scale*sigma; for (j=k+1;j<=n;j++) { for (sum=0.0,i=k;i<=n;i++) sum += a[i][k]*a[i][j]; tau=sum/c[k]; for (i=k;i<=n;i++) a[i][j] -= tau*a[i][k]; } } } d[n]=a[n][n]; if (d[n] == 0.0) *sing=1; } The next routine, qrsolv, is used to solve linear systems. In many applications only the part (2.10.4) of the algorithm is needed, so we separate it off into its own routine rsolv
100 Chapter 2.Solution of Linear Algebraic Equations void qrsolv(float *a,int n,float c[,float d[],float b[]) Solves the set of n linear equations A.x =b.a[1..n][1..n],c[1..n],and d[1..n]are input as the output of the routine qrdcmp and are not modified.b[1..n]is input as the right-hand side vector,and is overwritten with the solution vector on output. void rsolv(float *a,int n,float d],float b[]); int i,j; float sum,tau; for (j=1;j=1;1--){ for (sum=0.0,j=1+1;j<=n;j++)sum +=a[i][j]*b[j]; b[i]=(b[i]-sum)/d[i]; 2 IENTIFIC See [2]for details on how to use OR decomposition for constructing orthogonal bases, and for solving least-squares problems.(We prefer to use SVD,$2.6,for these purposes, because of its greater diagnostic capability in pathological cases. COMPUTING (ISBN 188812920 Updating a QR decomposition Some numerical algorithms involve solving a succession of linear systems each of which 10.621 differs only slightly from its predecessor.Instead of doing O(N)operations each time Recipes Numerica to solve the equations from scratch,one can often update a matrix factorization in O(N2) Recipes 43106 operations and use the new factorization to solve the next set of linear equations.The LU decomposition is complicated to update because of pivoting.However,OR turns out to be quite simple for a very common kind of update, A+A+s☒t (2.10.7) Software. (compare equation 2.7.1).In practice it is more convenient to work with the equivalent form A=QR→A'=Q'.R'=Q(R+u⑧v) (2.10.8) One can go back and forth between equations (2.10.7)and (2.10.8)using the fact that Q is orthogonal,giving t=v and either s=Q.u or u=QT.s (2.10.9) The algorithm [2]has two phases.In the first we apply N-1 Jacobi rotations(811.1)to reduce R+u v to upper Hessenberg form.Another N-1 Jacobi rotations transform this upper Hessenberg matrix to the new upper triangular matrix R'.The matrix Q'is simply the product of Q with the 2(N-1)Jacobi rotations.In applications we usually want Q,and the algorithm can easily be rearranged to work with this matrix instead of with Q
100 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). void qrsolv(float **a, int n, float c[], float d[], float b[]) Solves the set of n linear equations A · x = b. a[1..n][1..n], c[1..n], and d[1..n] are input as the output of the routine qrdcmp and are not modified. b[1..n] is input as the right-hand side vector, and is overwritten with the solution vector on output. { void rsolv(float **a, int n, float d[], float b[]); int i,j; float sum,tau; for (j=1;j=1;i--) { for (sum=0.0,j=i+1;j<=n;j++) sum += a[i][j]*b[j]; b[i]=(b[i]-sum)/d[i]; } } See [2] for details on how to use QR decomposition for constructing orthogonal bases, and for solving least-squares problems. (We prefer to use SVD, §2.6, for these purposes, because of its greater diagnostic capability in pathological cases.) Updating a QR decomposition Some numerical algorithms involve solving a succession of linear systems each of which differs only slightly from its predecessor. Instead of doing O(N3) operations each time to solve the equations from scratch, one can often update a matrix factorization in O(N2) operations and use the new factorization to solve the next set of linear equations. The LU decomposition is complicated to update because of pivoting. However, QR turns out to be quite simple for a very common kind of update, A → A + s ⊗ t (2.10.7) (compare equation 2.7.1). In practice it is more convenient to work with the equivalent form A = Q · R → A = Q · R = Q · (R + u ⊗ v) (2.10.8) One can go back and forth between equations (2.10.7) and (2.10.8) using the fact that Q is orthogonal, giving t = v and either s = Q · u or u = QT · s (2.10.9) The algorithm [2] has two phases. In the first we apply N − 1 Jacobi rotations (§11.1) to reduce R + u ⊗ v to upper Hessenberg form. Another N − 1 Jacobi rotations transform this upper Hessenberg matrix to the new upper triangular matrix R . The matrix Q is simply the product of Q with the 2(N − 1) Jacobi rotations. In applications we usually want QT , and the algorithm can easily be rearranged to work with this matrix instead of with Q.
2.10 QR Decomposition 101 #include #include "nrutil.h" void qrupdt(float **r,float **qt,int n,float u[],float v[]) Given the QR decomposition of some n x n matrix,calculates the QR decomposition of the matrix Q.(R+uv).The quantities are dimensioned as r [1..n][1..n],gt [1..n][1..n], u[1..n],and v[1..n].Note that Q is input and returned in qt. { void rotate(float **r,float **qt,int n,int i,float a,float b); int i,j,k; for (k=n;k>=1;k--){ Find largest k such that u[k]0. if (u[k])break; if(k=1;1--){ Transform R+u v to upper Hessenberg. granted for 19881992 rotate(r,qt,n,1,u[1],-u[1+1]); if(u[i]=0.0)u[i]=fabs(u[i+1]); 100 (including this one) else if (fabs(u[i])>fabs(u[i+1])) u[i]=fabs(u[1])*sqrt(1.0+s0R(u[1+1]/u[1]): 872 else u[i]=fabs(u[i+1])*sqrt(1.0+SQR(u[i]/u[i+1])); to any /Cambridge from NUMERICAL RECIPES IN for(j=1:j #include "nrutil.h" 是 Programs void rotate(float **r,float **qt,int n,int i,float a,float b) copy for their Given matrices r[1..n][1..n]and qt [1..n][1..n],carry out a Jacobi rotation on rows i and i+1 of each matrix.a and b are the parameters of the rotation:cos=a/va2+b2 sin b/va2+62. f int j; float c,fact,s,w,yi rectcustsen OF SCIENTIFIC COMPUTING(ISBN 1f(a==0.0)[ Avoid unnecessary overflow or underflow. c=0.0; s=(b>=0.071.0:-1.0); else if (fabs(a)>fabs(b)){ @cambridge.org 1988-1992 by Numerical Recipes 10-621 fact-b/ai c=SIGN(1.0/sqrt(1.0+(fact*fact)),a); 43108 s=fact*c; }else fact=a/b; s=SIGN(1.0/sgrt(1.0+(fact*fact)),b); c=fact*s; Software. for (j=i;j<=n;j++){ Premultiply r by Jacobi rotation. (outside North Amer y=r[1][j]; w=r[1+1][j]; r[i][j]=c*y-s*w; r[i+1][j]=s*y+c*w; for (j=1;j<=n;j++){ Premultiply gt by Jacobi rotation. y=qt[i][j]; w=gt[i+1][j]; qt[1][j]=c*y-s*和； qt[1+1][j]=s*y+c*w;
2.10 QR Decomposition 101 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include #include "nrutil.h" void qrupdt(float **r, float **qt, int n, float u[], float v[]) Given the QR decomposition of some n × n matrix, calculates the QR decomposition of the matrix Q· (R+u⊗v). The quantities are dimensioned as r[1..n][1..n], qt[1..n][1..n], u[1..n], and v[1..n]. Note that QT is input and returned in qt. { void rotate(float **r, float **qt, int n, int i, float a, float b); int i,j,k; for (k=n;k>=1;k--) { Find largest k such that u[k] = 0. if (u[k]) break; } if (k =1;i--) { Transform R + u ⊗ v to upper Hessenberg. rotate(r,qt,n,i,u[i],-u[i+1]); if (u[i] == 0.0) u[i]=fabs(u[i+1]); else if (fabs(u[i]) > fabs(u[i+1])) u[i]=fabs(u[i])*sqrt(1.0+SQR(u[i+1]/u[i])); else u[i]=fabs(u[i+1])*sqrt(1.0+SQR(u[i]/u[i+1])); } for (j=1;j #include "nrutil.h" void rotate(float **r, float **qt, int n, int i, float a, float b) Given matrices r[1..n][1..n] and qt[1..n][1..n], carry out a Jacobi rotation on rows i and i + 1 of each matrix. a and b are the parameters of the rotation: cos θ = a/√a2 + b2, sin θ = b/√a2 + b2. { int j; float c,fact,s,w,y; if (a == 0.0) { Avoid unnecessary overflow or underflow. c=0.0; s=(b >= 0.0 ? 1.0 : -1.0); } else if (fabs(a) > fabs(b)) { fact=b/a; c=SIGN(1.0/sqrt(1.0+(fact*fact)),a); s=fact*c; } else { fact=a/b; s=SIGN(1.0/sqrt(1.0+(fact*fact)),b); c=fact*s; } for (j=i;j<=n;j++) { Premultiply r by Jacobi rotation. y=r[i][j]; w=r[i+1][j]; r[i][j]=c*y-s*w; r[i+1][j]=s*y+c*w; } for (j=1;j<=n;j++) { Premultiply qt by Jacobi rotation. y=qt[i][j]; w=qt[i+1][j]; qt[i][j]=c*y-s*w; qt[i+1][j]=s*y+c*w; } }
102 Chapter 2.Solution of Linear Algebraic Equations We will make use of QR decomposition,and its updating,in 89.7. CITED REFERENCES AND FURTHER READING: Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il of Handbook for Automatic Com- putation (New York:Springer-Verlag),Chapter 1/8.[1] Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),885.2,5.3,12.6.[2] 2.11 Is Matrix Inversion an N3 Process? We close this chapter with a little entertainment,a bit of algorithmic prestidig- itation which probes more deeply into the subject of matrix inversion.We start with a seemingly simple question: How many individual multiplications does it take to perform the matrix mul- tiplication of two 2 x 2 matrices. University Press. THE a12 (2.11.1) 121 a22 Programs Eight,right?Here they are written explicitly: 爱 OF SCIENTIFIC( 61 c11=a11×b11+a12×b21 c12=a11×b12+a12×b22 (2.11.2) c21=a21×b11+a22×b21 C22=a21×b12+a22×b22 Do you think that one can write formulas for the c's that involve only seven 是留会 10-621 4310-5 multiplications?(Try it yourself,before reading on.) Numerical Recipes Such a set of formulas was,in fact,discovered by Strassen [11.The formulas are: (outside Q1≡(a11+a22）×(b11+b22) North Software. Q2≡(a21+a22)×b11 Q3≡a11×(b12-b22) Q4三a22×(-b11+b21) (2.11.3) Q5三(a11+a12）×b22 Q6≡(-a11+a21)×(b11+b2) Qr≡(a12-a22)×(b21+b22)
102 Chapter 2. Solution of Linear Algebraic Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). We will make use of QR decomposition, and its updating, in §9.7. CITED REFERENCES AND FURTHER READING: Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag), Chapter I/8. [1] Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), §§5.2, 5.3, 12.6. [2] 2.11 Is Matrix Inversion an N3 Process? We close this chapter with a little entertainment, a bit of algorithmic prestidigitation which probes more deeply into the subject of matrix inversion. We start with a seemingly simple question: How many individual multiplications does it take to perform the matrix multiplication of two 2 × 2 matrices, a11 a12 a21 a22 · b11 b12 b21 b22 = c11 c12 c21 c22 (2.11.1) Eight, right? Here they are written explicitly: c11 = a11 × b11 + a12 × b21 c12 = a11 × b12 + a12 × b22 c21 = a21 × b11 + a22 × b21 c22 = a21 × b12 + a22 × b22 (2.11.2) Do you think that one can write formulas for the c’s that involve only seven multiplications? (Try it yourself, before reading on.) Such a set of formulas was, in fact, discovered by Strassen [1]. The formulas are: Q1 ≡ (a11 + a22) × (b11 + b22) Q2 ≡ (a21 + a22) × b11 Q3 ≡ a11 × (b12 − b22) Q4 ≡ a22 × (−b11 + b21) Q5 ≡ (a11 + a12) × b22 Q6 ≡ (−a11 + a21) × (b11 + b12) Q7 ≡ (a12 − a22) × (b21 + b22) (2.11.3)