Press, Teukolsly, Vetterling, Flannery - Numerical Recipes in C (523184), страница 23

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 23)

Consult [2,3] for references on this. The NAGlibrary [4] has an analyze/factorize/operate capability. A substantial collection ofroutines for sparse matrix calculation is also available from IMSL [5] as the YaleSparse Matrix Package [6].You should be aware that the special order of interchanges and eliminations,732.7 Sparse Linear Systemsprescribed by a sparse matrix method so as to minimize fill-ins and arithmeticoperations, generally acts to decrease the method’s numerical stability as comparedto, e.g., regular LU decomposition with pivoting. Scaling your problem so as tomake its nonzero matrix elements have comparable magnitudes (if you can do it)will sometimes ameliorate this problem.In the remainder of this section, we present some concepts which are applicableto some general classes of sparse matrices, and which do not necessarily depend ondetails of the pattern of sparsity.Sherman-Morrison FormulaSuppose that you have already obtained, by herculean effort, the inverse matrixA−1 of a square matrix A.

Now you want to make a “small” change in A, forexample change one element aij , or a few elements, or one row, or one column.Is there any way of calculating the corresponding change in A−1 without repeatingyour difficult labors? Yes, if your change is of the formA → (A + u ⊗ v)(2.7.1)for some vectors u and v. If u is a unit vector ei , then (2.7.1) adds the componentsof v to the ith row. (Recall that u ⊗ v is a matrix whose i, jth element is the productof the ith component of u and the jth component of v.) If v is a unit vector ej , then(2.7.1) adds the components of u to the jth column.

If both u and v are proportionalto unit vectors ei and ej respectively, then a term is added only to the element aij .The Sherman-Morrison formula gives the inverse (A + u ⊗ v)−1 , and is derivedbriefly as follows:(A + u ⊗ v)−1 = (1 + A−1 · u ⊗ v)−1 · A−1= (1 − A−1 · u ⊗ v + A−1 · u ⊗ v · A−1 · u ⊗ v − .

. .) · A−1= A−1 − A−1 · u ⊗ v · A−1 (1 − λ + λ2 − . . .)= A−1 −(A−1 · u) ⊗ (v · A−1 )1+λ(2.7.2)whereλ ≡ v · A−1 · u(2.7.3)The second line of (2.7.2) is a formal power series expansion. In the third line, theassociativity of outer and inner products is used to factor out the scalars λ.The use of (2.7.2) is this: Given A−1 and the vectors u and v, we need onlyperform two matrix multiplications and a vector dot product,z ≡ A−1 · uw ≡ (A−1 )T · vλ=v·z(2.7.4)to get the desired change in the inverseA−1→A−1 −z⊗w1+λ(2.7.5)74Chapter 2.Solution of Linear Algebraic EquationsThe whole procedure requires only 3N 2 multiplies and a like number of adds (aneven smaller number if u or v is a unit vector).The Sherman-Morrison formula can be directly applied to a class of sparseproblems. If you already have a fast way of calculating the inverse of A (e.g., atridiagonal matrix, or some other standard sparse form), then (2.7.4)–(2.7.5) allowyou to build up to your related but more complicated form, adding for example arow or column at a time.

Notice that you can apply the Sherman-Morrison formulamore than once successively, using at each stage the most recent update of A−1(equation 2.7.5). Of course, if you have to modify every row, then you are back toan N 3 method. The constant in front of the N 3 is only a few times worse than thebetter direct methods, but you have deprived yourself of the stabilizing advantagesof pivoting — so be careful.For some other sparse problems, the Sherman-Morrison formula cannot bedirectly applied for the simple reason that storage of the whole inverse matrix A−1is not feasible. If you want to add only a single correction of the form u ⊗ v,and solve the linear system(A + u ⊗ v) · x = b(2.7.6)then you proceed as follows.

Using the fast method that is presumed available forthe matrix A, solve the two auxiliary problemsA·y=bA·z=ufor the vectors y and z. In terms of these,v·yx=y−z1 + (v · z)(2.7.7)(2.7.8)as we see by multiplying (2.7.2) on the right by b.Cyclic Tridiagonal SystemsSo-called cyclic tridiagonal systems occur quite frequently, and are a goodexample of how to use the Sherman-Morrison formula in the manner just described.The equations have the form  x1r1b 1 c1 0 · · ·β a 2 b 2 c2 · · ·  x2   r2   ··· ·  · · ·  =  · · ·  (2.7.9)  xN−1rN−1· · · aN−1 bN−1 cN−1α···0aNbNxNrNThis is a tridiagonal system, except for the matrix elements α and β in the corners.Forms like this are typically generated by finite-differencing differential equationswith periodic boundary conditions (§19.4).We use the Sherman-Morrison formula, treating the system as tridiagonal plusa correction.

In the notation of equation (2.7.6), define vectors u and v to be 1γ 0 0 . .. .(2.7.10)v=u= . .  0 0β/γα2.7 Sparse Linear Systems75Here γ is arbitrary for the moment. Then the matrix A is the tridiagonal part of thematrix in (2.7.9), with two terms modified:b1 = b1 − γ,bN = bN − αβ/γ(2.7.11)We now solve equations (2.7.7) with the standard tridiagonal algorithm, and thenget the solution from equation (2.7.8).The routine cyclic below implements this algorithm. We choose the arbitraryparameter γ = −b1 to avoid loss of precision by subtraction in the first of equations(2.7.11).

In the unlikely event that this causes loss of precision in the second ofthese equations, you can make a different choice.#include "nrutil.h"void cyclic(float a[], float b[], float c[], float alpha, float beta,float r[], float x[], unsigned long n)Solves for a vector x[1..n] the “cyclic” set of linear equations given by equation (2.7.9). a,b, c, and r are input vectors, all dimensioned as [1..n], while alpha and beta are the cornerentries in the matrix. The input is not modified.{void tridag(float a[], float b[], float c[], float r[], float u[],unsigned long n);unsigned long i;float fact,gamma,*bb,*u,*z;if (n <= 2) nrerror("n too small in cyclic");bb=vector(1,n);u=vector(1,n);z=vector(1,n);gamma = -b[1];Avoid subtraction error in forming bb[1].bb[1]=b[1]-gamma;Set up the diagonal of the modified tridibb[n]=b[n]-alpha*beta/gamma;agonal system.for (i=2;i<n;i++) bb[i]=b[i];tridag(a,bb,c,r,x,n);Solve A · x = r.u[1]=gamma;Set up the vector u.u[n]=alpha;for (i=2;i<n;i++) u[i]=0.0;tridag(a,bb,c,u,z,n);Solve A · z = u.fact=(x[1]+beta*x[n]/gamma)/Form v · x/(1 + v · z).(1.0+z[1]+beta*z[n]/gamma);for (i=1;i<=n;i++) x[i] -= fact*z[i];Now get the solution vector x.free_vector(z,1,n);free_vector(u,1,n);free_vector(bb,1,n);}Woodbury FormulaIf you want to add more than a single correction term, then you cannot use (2.7.8)repeatedly, since without storing a new A−1 you will not be able to solve the auxiliaryproblems (2.7.7) efficiently after the first step.

Instead, you need the Woodbury formula,which is the block-matrix version of the Sherman-Morrison formula,(A + U · VT )−1= A−1 − A−1 · U · (1 + VT · A−1 · U)−1 · VT · A−1(2.7.12)76Chapter 2.Solution of Linear Algebraic EquationsHere A is, as usual, an N × N matrix, while U and V are N × P matrices with P < Nand usually P N . The inner piece of the correction term may become clearer if writtenas the tableau,U −1   · 1 + VT · A−1 · U ·  VT(2.7.13)where you can see that the matrix whose inverse is needed is only P × P rather than N × N .The relation between the Woodbury formula and successive applications of the ShermanMorrison formula is now clarified by noting that, if U is the matrix formed by columns out of theP vectors u1 , .

. . , uP , and V is the matrix formed by columns out of the P vectors v1 , . . . , vP ,        U ≡ u 1  · · · u P       V ≡ v1  · · · vP   then two ways of expressing the same correction to A arePA+uk ⊗ vk = (A + U · VT )(2.7.14)(2.7.15)k=1(Note that the subscripts on u and v do not denote components, but rather distinguish thedifferent column vectors.)Equation (2.7.15) reveals that, if you have A−1 in storage, then you can either make theP corrections in one fell swoop by using (2.7.12), inverting a P × P matrix, or else makethem by applying (2.7.5) P successive times.If you don’t have storage for A−1 , then you must use (2.7.12) in the following way:To solve the linear equationPA+(2.7.16)uk ⊗ vk · x = bk=1first solve the P auxiliary problemsA · z 1 = u1A · z 2 = u2···(2.7.17)A · z P = uPand construct the matrix Z by columns from the z’s obtained,      Z ≡ z1  · · · zP   (2.7.18)Next, do the P × P matrix inversionH ≡ (1 + VT · Z)−1(2.7.19)2.7 Sparse Linear Systems77Finally, solve the one further auxiliary problemA·y =b(2.7.20)In terms of these quantities, the solution is given byx = y − Z · H · (VT · y)(2.7.21)Inversion by PartitioningOnce in a while, you will encounter a matrix (not even necessarily sparse)that can be inverted efficiently by partitioning.

Suppose that the N × N matrixA is partitioned intoA=PRQS(2.7.22)where P and S are square matrices of size p × p and s × s respectively (p + s = N ).The matrices Q and R are not necessarily square, and have sizes p × s and s × p,respectively.If the inverse of A is partitioned in the same manner,−1A=PQRS!(2.7.23)then P, Q, R, S, which have the same sizes as P, Q, R, S, respectively, can befound by either the formulasP = (P − Q · S−1 · R)−1Q = −(P − Q · S−1 · R)−1 · (Q · S−1 )R = −(S−1 · R) · (P − Q · S−1 · R)−1(2.7.24)S = S−1 + (S−1 · R) · (P − Q · S−1 · R)−1 · (Q · S−1 )or else by the equivalent formulasP = P−1 + (P−1 · Q) · (S − R · P−1 · Q)−1 · (R · P−1 )Q = −(P−1 · Q) · (S − R · P−1 · Q)−1R = −(S − R · P−1 · Q)−1 · (R · P−1 )(2.7.25)S = (S − R · P−1 · Q)−1The parentheses in equations (2.7.24) and (2.7.25) highlight repeated factors thatyou may wish to compute only once. (Of course, by associativity, you can insteaddo the matrix multiplications in any order you like.) The choice between usingequation (2.7.24) and (2.7.25) depends on whether you want P or S to have the78Chapter 2.Solution of Linear Algebraic Equationssimpler formula; or on whether the repeated expression (S − R · P−1 · Q)−1 is easierto calculate than the expression (P − Q · S−1 · R)−1 ; or on the relative sizes of Pand S; or on whether P−1 or S−1 is already known.Another sometimes useful formula is for the determinant of the partitionedmatrix,det A = det P det(S − R · P−1 · Q) = det S det(P − Q · S−1 · R)(2.7.26)Indexed Storage of Sparse MatricesWe have already seen (§2.4) that tri- or band-diagonal matrices can be stored in a compactformat that allocates storage only to elements which can be nonzero, plus perhaps a few wastedlocations to make the bookkeeping easier.

What about more general sparse matrices? When asparse matrix of dimension N × N contains only a few times N nonzero elements (a typicalcase), it is surely inefficient — and often physically impossible — to allocate storage for allN 2 elements. Even if one did allocate such storage, it would be inefficient or prohibitive inmachine time to loop over all of it in search of nonzero elements.Obviously some kind of indexed storage scheme is required, one that stores only nonzeromatrix elements, along with sufficient auxiliary information to determine where an elementlogically belongs and how the various elements can be looped over in common matrixoperations.

Характеристики

Список файлов книги