Page 224 - Compact Numerical Methods For Computers
P. 224
Minimising a nonlinear sum of squares 213
Algorithm 23. Modified Marquardt method for minimising a nonlinear sum-of-squares
function (cont.)
a, c: smatvec;
delta, v : rvector;
dec, eps, inc, lambda, p, phi, res : real;
count, i, ifn, igrad, j, k, nn2, q : integer;
notcomp, singmat, calcmat: boolean;
begin
writeln(‘alg23.pas -- Nash Marquardt nonlinear least squares’);
with Workdata do
begin {STEP 0 partly in procedure call}
if nlls = false then halt; {cannot proceed if we do not have a nonlinear
least squares problem available}
Fmin:=big; {safety setting of the minimal function value}
inc:=10.0; {increase factor for damping coefficient lambda}
dec:=0.4; {decrease factor for damping coefficient lambda}
eps:=calceps; {machine precision}
lambda:=0.0001; {initialize damping factor}
phi:=1.0; {set the Nash damping factor}
ifn:=0; igrad:=0; {set the function and gradient evaluation counts}
calcmat:=true; {to force calculation of the J-transpose * J matrix and
J-transpose * residual on the first iteration}
nn2:=(n*(n+l)) div 2; {elements in the triangular form of the inner
product matrix -- ensure this is an integer}
p:=0.0; {STEP 1}
for i:=1 to m do
begin
res:-lllres(i, n, Bvec, notcomp); {the residual}
{writeln(‘res[‘,i,’]=’,res);
if notcomp then halt; {safety check on initial evaluation}
p:=p+res*res; {sum of squares accumulation}
end;
ifn:=ifn+1; {count the function evaluation}
Fmin:=p; {to save best sum of squares so far}
count:=0; {to avoid convergence immediately}
{STEP 2}
while count<n do
begin {main Marquardt iteration}
{NOTE: in this version we do not reduce the damping parameter here,
The structure of Pascal lends itself better to adjusting the damping
parameter below.}
if calcmat then
begin {Compute sum of squares and cross-products matrix}
writeln(igrad,’ ’,ifn,’ sum of squares=’,Fmin);
for i:=1 to n do
begin
write(Bvec[i]:10:fi,’ ’);
if (7 * (i div 7) = i) and (i<r) then writeln;
end; {loop on i}
writeln;
igrad:=igrad+1; {STEP 3}
for j:=1 to nn2 do a[j]:=0.0;
for j:=1 to n do v[j]:=0.0;
for i:=1 to m do {STEP 4}