Nash - Compact Numerical Methods for Computers (523163), страница 53
Текст из файла (страница 53)
Bradbury and Fletcher (1966) address thisdifficulty by setting to unity the element of largest magnitude in the currenteigenvector approximation and adjusting the other elements accordingly. Thisadjustment is made at each iteration. However, Geradin (1971) has tackled theproblem more directly, that is, by examining the Hessian itself and attempting toconstruct search directions which are mutually conjugate with respect to it. Thistreatment, though surely not the last word on the subject, is essentially repeatedhere.
The implementation details are my own.Firstly, consider the linear search subproblem, which in the current case canbe solved analytically. It is desired to minimiseR= (x+k t) TA( x+k t) / (x+k t) TB( x+k t)(19.33)with respect to k. For convenience this will be rewrittenR=N( k) /D( k )(19.34)with N and D used to denote the numerator and denominator, respectively.Differentiating with respect to k givesd R/ d k=0=(DdN/ d k- NdD/ d k) /D 2 .(19.35)Because of the positive definiteness of B, D can never be zero unlessx + kt = 0 .(19.36)Therefore, ignoring this last possibility, we set the numerator of expression(19.35) to zero to obtain the quadratic equationuk2 +v k + w =0(19.37)u = (t T At)(x T Bt)-(x T At)(t T Bt)(19.38)whereConjugate gradients method in linear algebrav= (t T At) (x T Bx) - (x T Ax) (t T Bt)TTTT245(19.39)w = (x At) (x Bx) - (x Ax) (x Bt).(19.40)x T At = t T A x(19.41)x T Bt=t T Bx.(19.42)Note that by symmetryandTherefore, only six inner products are needed in the evaluation of u, v and w.These areand(x T Ax)(x T At)and(tT At)(x T Bx)(x T Bt)and(tT Bt).The quadratic equation (19.37) has two roots, of which only one will correspond toa minimum.
Sincey (k )=0·5D 2 ( dR/ dk ) = u k 2 + v k + w(19.43)we get(19.44)At either of the roots of (19.37), this reduces to0 · 5D 2 ( d 2 R/ dk 2 )=2u k + v(19.45)so that a minimum has been found if d2R/dk 2 is positive, or2u k + v > 0 .(19.46)Substitution of both roots from the quadratic equation formula shows that thedesired root isk = [ -v + (v 2 - 4u w) ½ ] / ( 2u ) .(19.47)If v is negative, (19.47) can be used to evaluate the root. However, to avoid digitcancellation when v is positive,k = - 2w / [v + (v 2 - 4u w) ½ ](19.48)should be used. The linear search subproblem has therefore been resolved in astraightforward manner in this particular case.The second aspect of particularising the conjugate gradients algorithm to theminimisation of the Rayleigh quotient is the generation of the next conjugatedirection.
Note that the gradient of the Rayleigh quotient at x is given byg= 2 (Ax- RBx) / (x T Bx)(19.49)H= 2 (A- R B - Bxg T - gx T B) / (x T Bx) .(19.50)and the local Hessian by246Compact numerical methods for computersSubstitutingq= -gand the Hessian (19.50) into the expressionz=g T Ht/ t T Ht(16.5)(16.4)for the parameter in the two-term conjugate gradient recurrence, and noting thatgTt = 0(19.51)by virtue of the ‘exact’ linear searches, givesz=[g T (A-RB) t- (g T g) (x T Bt)]/[t T (A-RB) t].(19.52)The work of Geradin (1971) implies that z should be evaluated by computing theinner products within the square brackets and subtracting. I prefer to perform thesubtraction within each element of the inner product to reduce the effect of digitcancellation.Finally, condition (19.51) permits the calculation of the new value of theRayleigh quotient to be simplified.
Instead of the expression which results fromexpanding (19.33), we have from (19.49) evaluted at (x+kt), with (19.51), theexpressionR= (t T Ax+k t T At) / (t T Bx+k t T Bt).(19.53)This expression is not used in the algorithm. Fried (1972) has suggested severalother formulae for the recurrence parameter z of equation (19.52). At the time ofwriting, too little comparative testing has been carried out to suggest that onesuch formula is superior to any other.Algorithm 2.5. Rayleigh quotient minimisation by conjugate gradientsprocedure rqmcg( n : integer; {order of matrices}A, B : rmatrix; {matrices defining eigenproblem}var X : rvector; {eigenvector approximation, on bothinput and output to this procedure}var ipr : integer; {on input, a limit to the number ofmatrix products allowed, on output, the number ofmatrix products used}var rq : real); {Rayleigh quotient = eigenvalue approx.}{alg25.pas == Rayleigh quotient minimization by conjugate gradientsMinimize Rayleigh quotientX-transpose A X / X-transpose B Xthereby solving generalized symmetric matrix eigenproblemAX=rqBXfor minimal eigenvalue rq and its associated eigenvector.A and B are assumed symmetric, with B positive definite.While we supply explicit matrices here, only matrix productsare needed of the form v = A u, w = B u.Copyright 1988 J.C.Nash}varcount i, itn, itlimit : integer;Conjugate gradients method in linear algebraAlgorithm 25.
Rayleigh quotient minimisation by conjugate gradients (cont.)avec, bvec, yvec, zvec, g, t : rvector;beta, d, eps, g2, gg, oldg2, pa, pn, s2, step : real;t2, ta, tabt, tat, tb, tbt, tol, u, v, w, xat, xax, xbt, xbx : real;conv, fail : boolean;beginwriteln(‘alg25.pas -- Rayleigh quotient minimisation’);itlimit := ipr; {to save the iteration limit} {STEP 0}fail := false; {Algorithm has yet to fail.}conv := false; {Algorithm has yet to converge.}ipr := 0; {to initialize the iteration count}eps := calceps;tol := n*n*eps*eps; {a convergence tolerance}{The convergence tolerance, tol, should ideally be chosen relativeto the norm of the B matrix when used for deciding if matrix B issingular.
It needs a different scaling when deciding if the gradientis “small”. Here we have chosen a compromise value. Performance ofthis method could be improved in production codes by judicious choiceof the tolerances.}pa := big; {a very large initial value for the ‘minimum’ eigenvalue}while (ipr<=itlimit) and (not conv) dobegin {Main body of algorithm}matmul(n, A, X, avec); {STEP 1}matmul(n, B, X, bvec); {initial matrix multiplication}ipr := ipr+l ; {to count the number of products used}{STEP 2: Now form the starting Rayleigh quotient}xax := 0.0; xbx := 0.0; {accumulators for numerator and denominatorof the Rayleigh quotient}for i := l to n dobeginxax := xax+X[i]*avec[i]; xbx := xbx+X[i]*bvec[i];end; {loop on i for Rayleigh quotient}if xbx<=tol then halt; {STEP 3: safety check to avoid zero divide.This may imply a singular matrix B, or an inappropriatestarting vector X i.e.
one which is in the null space of B(implying B is singular!) or which is itself null.}rq := xax/xbx; {the Rayleigh quotient -- STEP 4}write(ipr,’ products -- ev approx. =’,rq:18);if rqcpa then {Principal convergence check, since if the Rayleighquotient has not been decreased in a major cycle, and we mustpresume that the minimum has been found. Note that this testrequires that we initialize pa to a large number.)begin {body of algorithm -- STEP 5}pa := rq; {to save the lowest value so far}gg := 0.0; {to initialize gradient norm.
Now calculate gradient.}for i := 1 to n do {STEP 6}beging[i] := 2.0*(avec[i]-rq*bvec[i])/xbx; gg := gg+g[i]*g[i];end; {gradient calculation}writeln(‘ squared gradient norm =’,gg:8);if gg>tol then {STEP 7}{Test to see if algorithm has converged. This test is unscaledand in some problems it is necessary to scale the tolerance tol247248Compact numerical methods for computersAlgorithm 25. Rayleigh quotient minimisation by conjugate gradients (cont.)to the size of the numbers in the problem.}begin {conjugate gradients search for improved eigenvector}{Now generate the first search direction.}for i := 1 to n do t[i] := -g[i]; {STEP 8}itn:= 0; {STEP 9}repeat {Major cg loop}itn := itn+l; {to count the conjugate gradient iterations}matmul(n.
A, t, yvec); {STEP 10}matmul(n, B, t, zvec); ipr := ipr+l;tat := 0.0; tbt := 0.0; xat := 0.0; xbt := 0.0; {STEP 11}for i := 1 to n dobeginxat := xat+X[i]*yvec[i]; tat := tat+t[i]*yvec[i];xbt := xbt+X[i]*zvec[i]; tbt := tbt+t[i]*zvec[i];end;{STEP 12 -- formation and solution of quadratic equation}u := tat*xbt-xat*tbt; v := tat*xbx-xax*tbt;w := xat*xbx-xax*xbt; d := v*v-4.0*U*W; {the discriminant}if d<0.0 then halt; {STEP 13 -- safety check}{Note: we may want a more imaginative response to this resultof the computations. Here we will assume imaginary roots ofthe quadradic cannot arise, but perform the check.}d := sqrt(d); {Now compute the roots in a stable manner -- STEP 14}if w0.0 then step := -2.0*w/(v+d) else step := 0.5*(d-v)/u;{STEP 15 -- update vectors}count := 0; {to count the number of unchanged vector components}xax := 0.0; xbx := 0.0;for i := l to n dobeginavec[i] := avec[i]+step*yvec[i];bvec[i] := bvec[i]+step*zvec[i];w := X[i]; X[i] := w+step*t[i];if (reltest+w)=(reltest+X[i]) then count := count+l;xax := xax+X[i]*avec[i]; xbx := xbx+X[i]*bvec[i];end; {loop on i}if xbx<=tol then halt {to avoid zero divide if B singular}else pn := xax/xbx; {STEP 16}if (count<n) and (pn<rq) then {STEPS 17 & 18}beginrq := pn gg := 0.0; {STEP 19}for i:= 1 to n dobeging[i] := 2.0*(avec[i]-pn*bvec[i])/xbx; gg := gg+g[i]*g[i];end; {loop on i}if gg>tol then {STEP 20}begin {STEP 21}xbt := 0.0; for i := 1 to n do xbt := xbt+X[i]*zvec[i];{STEP 22 -- compute formula (19.52)}tabt := 0.0; beta := 0.0;for i := l to n dobeginw := yvec[i]-pn*zvec[i]; tabt := tabt+t[i]*w;Conjugate gradients method in linear algebra249Algorithm 25.
Rayleigh quotient minimisation by conjugate gradients (cont.)beta := beta+g[i]*(w-g[i]*xbt); {to form the numerator}end; {loop on i}beta := beta/tabt; {STEP 23}{Now perform the recurrence for the next search vector.}for i := 1 to n do t[i] := beta*t[i]-g[i];end; {if gg>tol}end {if (count<n) and (pn<rq)}{Note: pn is computed from update information and may not beprecise.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
