Page 249 - Compact Numerical Methods For Computers

236 Compact numerical methods for computers
Substitution of
b =b +k t (19.7)
j + 1 j j j
into (15.11) gives an equation for k , permitting the optimal step length to be
j
computed. For convenience in explaining this, the subscript j will be omitted.
Thus from substituting, we obtain
T
f(k)=½(b+k ) H(b+kt) -c (b+k t)+(anyscalar). (19.8)
T
t
The derivative of this with respect to k is
T T
f' (k) =t H(b+k t) -c t
T T
=t (Hb-c) +k t Ht
T T
=t g+kt Ht (19.9)
so that f'(k)=0 implies
T
T
k=-t g/t Ht. (19.10)
But from the same line of reasoning as that which follows equation (16.6), the
accurate line searches imply that the function has been minimised on a hyperplane
spanned by the gradient directions which go to make up t except for g itself: thus
T T
k=+g g/t Ht. (19.11)
Note the sign has changed since t will have component in direction g with a
coefficient of -1.
The recurrence formula for this problem can be generated by one of the
formulae (16.8), (16.10) or (16.11). The last of these is the simplest and is the one
that has been chosen. Because the problem has a genuine quadratic form it may
be taken to have converged after n steps if the gradient does not become small
earlier. However, the algorithm which follows restarts the iteration in case
rounding errors during the conjugate gradients steps have caused the solution to
be missed.
Algorithm 24. Solution of a consistent set of linear equations by conjugate gradients

procedure lecg( n : integer; {order of the problem}
H : rmatrix; {coefficient matrix in linear equations}
C : rvector; (right hand side vector}
var Bvec : rvector; {the solution vector -- on input, we
must supply some guess to this vector}
var itcount : integer; {on input, a limit to the number of
conjugate gradients iterations allowed; on output,
the number of iterations used (negative if the
limit is exceeded)}
var ssmin : real); {the approximate minimum sum of squared
residuals for the solution}
{alg24.pas == linear equations by conjugate gradients
This implementation uses an explicit matrix H. However,
we only need the result of multiplying H times a vector,
that is, where the procedure call
matmul( n, H, t, v);