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}
   219   220   221   222   223   224   225   226   227   228   229