Page 225 - Compact Numerical Methods For Computers

214 Compact numerical methods for computers

Algorithm 23. Modified Marquardt method for minimising a nonlinear sum-of-squares
function (cont.)

begin
nljac(i, n, Bvec, X); {puts i’th row of Jacobian in X}
res:=nlres(i, n, Bvec, notcomp); {This calculation is not really
necessary. The results of the sum of squares calculation
can be saved in a residual vector to avoid the
recalculation. However, this way saves storing a possibly
large vector. NOTE: we ignore notcomp here, since the
computation of the residuals has already been proven
possible at STEP 1.}
for j:=1 to n do
begin
v[j]:=vu[j]+X[j]*res; {to accumulate the gradient}
q:=(j*(j-l)) div 2; {to store the correct position in the
row-order vector form of the lower triangle of a symmetric
matrix}
for k:=1 to j do a[q+k]:=a[q+k]+X[j]*X[k];
end; {loop on j}
end; {loop on i}
for j:=1 to nn2 do c[j]:=a[j]; {STEP 5 -- copy a and b}
for j:=1 to n do X[j]:=Bvec[j]; {to save the best parameters}
end; {if calcmat}
writeln(‘LAMDA=’,lambda:8); {STEP 6}
for j:=1 to n do
begin
q:=(i*(j+l)) div 2;
a[q]:=c[q]*(1.0+lambda)+phi*lambda; {damping added}
delta[j]:=-v[j]; {to set RHS in equations 17.24}
if j>1 then
for i:=1 to (j-1) do a[q-i]:=c[q-i];
end; {loop on j}
notcomp:=false; {to avoid undefined value}
Choldcmp(n, a, singmat); {STEP 7 -- Choleski factorization}
if (not singmat) then {matrix successfully decomposed}
begin {STEP 8 -- Choleski back-substitution}
Cholback(n, a, delta); {delta is the change in the parameters}
count:=0; {to count parameters unchanged in update}
for i:=1 to n do {STEP 9}
begin
Bvec[i]:=X[i]+delta[i]; {new = old + update}
if (reltest+Bvec[i])=(reltest+X[i]) then count:=count+1;
{Here the global constant reltest is used to test for the
equality of the new and old parameters.}
end; {loop on i over parameters)
if count<n then {STEP 10: parameters have been changed}
begin {compute the sum of squares, checking computability}
p:=0.0; i:=0; {initialization}
repeat
i:=i+1; res:=nlres(i,n,Bvec,notcomp);
if (not notcomp) then p:=p+res*res;
until notcomp or (i>=n);
ifn:=ifn+1; {count the function evaluation}
end; {if count<n}