Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 29
Текст из файла (страница 29)
The upper triangular matrix R is returned in the upper triangle of a, except for the diagonal elements of R which are returned ind[1..n]. The orthogonal matrix Q is represented as a product of n − 1 Householder matricesQ1 . . . Qn−1 , where Qj = 1 − uj ⊗ uj /cj . The ith component of uj is zero for i = 1, . . . , j − 1while the nonzero components are returned in a[i][j] for i = j, . .
. , n. sing returns astrue (1) if singularity is encountered during the decomposition, but the decomposition is stillcompleted 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 thepart (2.10.4) of the algorithm is needed, so we separate it off into its own routine rsolv.100Chapter 2.Solution of Linear Algebraic Equationsvoid 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] areinput as the output of the routine qrdcmp and are not modified. b[1..n] is input as theright-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<n;j++) {Form QT · b.for (sum=0.0,i=j;i<=n;i++) sum += a[i][j]*b[i];tau=sum/c[j];for (i=j;i<=n;i++) b[i] -= tau*a[i][j];}rsolv(a,n,d,b);Solve R · x = QT · b.}void rsolv(float **a, int n, float d[], float b[])Solves the set of n linear equations R · x = b, where R is an upper triangular matrix stored ina and d.
a[1..n][1..n] and d[1..n] are input as the output of the routine qrdcmp andare not modified. b[1..n] is input as the right-hand side vector, and is overwritten with thesolution vector on output.{int i,j;float sum;b[n] /= d[n];for (i=n-1;i>=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 decompositionSome numerical algorithms involve solving a succession of linear systems each of whichdiffers only slightly from its predecessor. Instead of doing O(N 3 ) operations each timeto solve the equations from scratch, one can often update a matrix factorization in O(N 2 )operations and use the new factorization to solve the next set of linear equations. The LUdecomposition is complicated to update because of pivoting. However, QR turns out to bequite 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 formA = 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 Qis orthogonal, givingt = v and either s = Q · u oru = QT · s(2.10.9)The algorithm [2] has two phases. In the first we apply N − 1 Jacobi rotations (§11.1) toreduce R + u ⊗ v to upper Hessenberg form. Another N − 1 Jacobi rotations transform thisupper Hessenberg matrix to the new upper triangular matrix R . The matrix Q is simply theproduct of Q with the 2(N − 1) Jacobi rotations. In applications we usually want QT , andthe algorithm can easily be rearranged to work with this matrix instead of with Q.2.10 QR Decomposition101#include <math.h>#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 thematrix 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) k=1;for (i=k-1;i>=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<=n;j++) r[1][j] += u[1]*v[j];for (i=1;i<k;i++)Transform upper Hessenberg matrix to upper trirotate(r,qt,n,i,r[i][i],-r[i+1][i]);angular.}#include <math.h>#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 rowsi 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;}}102Chapter 2.Solution of Linear Algebraic EquationsWe 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 HopkinsUniversity Press), §§5.2, 5.3, 12.6. [2]2.11 Is Matrix Inversion an N 3 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 startwith a seemingly simple question:How many individual multiplications does it take to perform the matrixmultiplication of two 2 × 2 matrices,a11a21a12a22 b11·b21b12b22=c11c21c12c22(2.11.1)Eight, right? Here they are written explicitly:c11 = a11 × b11 + a12 × b21c12 = a11 × b12 + a12 × b22c21 = a21 × b11 + a22 × b21(2.11.2)c22 = a21 × b12 + a22 × b22Do you think that one can write formulas for the c’s that involve only sevenmultiplications? (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 ) × b11Q3 ≡ a11 × (b12 − b22 )Q4 ≡ a22 × (−b11 + b21 )Q5 ≡ (a11 + a12 ) × b22Q6 ≡ (−a11 + a21 ) × (b11 + b12 )Q7 ≡ (a12 − a22 ) × (b21 + b22 )(2.11.3)2.11 Is Matrix Inversion an N 3 Process?103in terms of whichc11 = Q1 + Q4 − Q5 + Q7c21 = Q2 + Q4c12 = Q3 + Q5(2.11.4)c22 = Q1 + Q3 − Q2 + Q6What’s the use of this? There is one fewer multiplication than in equation(2.11.2), but many more additions and subtractions.
It is not clear that anythinghas been gained. But notice that in (2.11.3) the a’s and b’s are never commuted.Therefore (2.11.3) and (2.11.4) are valid when the a’s and b’s are themselvesmatrices. The problem of multiplying two very large matrices (of order N = 2m forsome integer m) can now be broken down recursively by partitioning the matricesinto quarters, sixteenths, etc. And note the key point: The savings is not just a factor“7/8”; it is that factor at each hierarchical level of the recursion. In total it reducesthe process of matrix multiplication to order N log2 7 instead of N 3 .What about all the extra additions in (2.11.3)–(2.11.4)? Don’t they outweighthe advantage of the fewer multiplications? For large N , it turns out that there aresix times as many additions as multiplications implied by (2.11.3)–(2.11.4). But,if N is very large, this constant factor is no match for the change in the exponentfrom N 3 to N log2 7 .With this “fast” matrix multiplication, Strassen also obtained a surprising resultfor matrix inversion [1].
Suppose that the matricesa11 a12c11 c12and(2.11.5)a21 a22c21 c22are inverses of each other. Then the c’s can be obtained from the a’s by the followingoperations (compare equations 2.7.22 and 2.7.25):R1 = Inverse(a11 )R2 = a21 × R1R3 = R1 × a12R4 = a21 × R3R5 = R4 − a22R6 = Inverse(R5 )c12 = R3 × R6c21 = R6 × R2R7 = R3 × c21c11 = R1 − R7c22 = −R6(2.11.6)104Chapter 2.Solution of Linear Algebraic EquationsIn (2.11.6) the “inverse” operator occurs just twice. It is to be interpreted as thereciprocal if the a’s and c’s are scalars, but as matrix inversion if the a’s and c’s arethemselves submatrices. Imagine doing the inversion of a very large matrix, of orderN = 2m , recursively by partitions in half. At each step, halving the order doublesthe number of inverse operations. But this means that there are only N divisions inall! So divisions don’t dominate in the recursive use of (2.11.6). Equation (2.11.6)is dominated, in fact, by its 6 multiplications.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
