octave_packages/optim-1.2.0/private/__lm_svd__.m

   1 %% Copyright (C) 1992-1994 Richard Shrager
   2 %% Copyright (C) 1992-1994 Arthur Jutan
   3 %% Copyright (C) 1992-1994 Ray Muzic
   4 %% Copyright (C) 2010, 2011 Olaf Till <olaf.till@uni-jena.de>
   5 %%
   6 %% This program is free software; you can redistribute it and/or modify it under
   7 %% the terms of the GNU General Public License as published by the Free Software
   8 %% Foundation; either version 3 of the License, or (at your option) any later
   9 %% version.
  10 %%
  11 %% This program is distributed in the hope that it will be useful, but WITHOUT
  12 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  13 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  14 %% details.
  15 %%
  16 %% You should have received a copy of the GNU General Public License along with
  17 %% this program; if not, see <http://www.gnu.org/licenses/>.
  18
  19 function [p, resid, cvg, outp] = __lm_svd__ (F, pin, hook)
  20
  21   %% This is a backend for optimization. This code was originally
  22   %% contained in leasqr.m, which is now a frontend.
  23
  24   %% some backend specific defaults
  25   fract_prec_default = 0;
  26   max_fract_step_default = Inf;
  27
  28   %% needed for some anonymous functions
  29   if (exist ('ifelse') ~= 5)
  30     ifelse = @ scalar_ifelse;
  31   end
  32
  33   %% passed constraints
  34   mc = hook.mc; % matrix of linear constraints
  35   vc = hook.vc; % vector of linear constraints
  36   f_cstr = hook.f_cstr; % function of all constraints
  37   df_cstr = hook.df_cstr; % function of derivatives of all constraints
  38   n_gencstr = hook.n_gencstr; % number of non-linear constraints
  39   eq_idx = hook.eq_idx; % logical index of equality constraints in all
  40                                 % constraints
  41   lbound = hook.lbound; % bounds, subset of linear inequality
  42   ubound = hook.ubound; % constraints in mc and vc
  43
  44   %% passed values of constraints for initial parameters
  45   pin_cstr = hook.pin_cstr;
  46
  47   %% passed return value of F for initial parameters
  48   f_pin = hook.f_pin;
  49
  50   %% passed derivative of residual function
  51   dfdp = hook.dfdp;
  52
  53   %% passed function for complementary pivoting
  54   cpiv = hook.cpiv;
  55
  56   %% passed options
  57   maxstep = hook.max_fract_change;
  58   maxstep(isna (maxstep)) = max_fract_step_default;
  59   pprec = hook.fract_prec;
  60   pprec(isna (pprec)) = fract_prec_default;
  61   stol = hook.TolFun;
  62   niter = hook.MaxIter;
  63   if (isempty (niter)) niter = 20; end
  64   wt = hook.weights;
  65   fixed = hook.fixed;
  66   verbose = strcmp (hook.Display, 'iter');
  67
  68   %% only preliminary, for testing
  69   if (isfield (hook, 'testing'))
  70     testing = hook.testing;
  71   else
  72     testing = false;
  73   end
  74   if (isfield (hook, 'new_s'))
  75     new_s = hook.new_s;
  76   else
  77     new_s = false;
  78   end
  79
  80   %% some useful variables derived from passed variables
  81   n_lcstr = size (vc, 1);
  82   have_constraints_except_bounds = ...
  83       n_lcstr + n_gencstr > ...
  84       sum (lbound ~= -Inf) + sum (ubound ~= Inf);
  85   n = length (pin);
  86   wtl = wt(:);
  87
  88   nz = 20 * eps; % This is arbitrary. Constraint function will be
  89                                 % regarded as <= zero if less than nz.
  90
  91   %% backend-specific checking of options and constraints
  92   if (have_constraints_except_bounds)
  93     if (any (pin_cstr.inequ.lin_except_bounds < 0) || ...
  94         (n_gencstr > 0 && any (pin_cstr.inequ.gen < 0)))
  95       warning ('initial parameters violate inequality constraints');
  96     end
  97     if (any (abs (pin_cstr.equ.lin) >= nz) || ...
  98         (n_gencstr > 0 && any (abs (pin_cstr.equ.gen) >= nz)))
  99       warning ('initial parameters violate equality constraints');
 100     end
 101   end
 102   idx = lbound == ubound;
 103   if (any (idx))
 104     warning ('lower and upper bounds identical for some parameters, fixing the respective parameters');
 105     fixed(idx) = true;
 106   end
 107   if (all (fixed))
 108     error ('no free parameters');
 109   end
 110   lidx = pin < lbound;
 111   uidx = pin > ubound;
 112   if (any (lidx | uidx) && have_constraints_except_bounds)
 113     warning ('initial parameters outside bounds, not corrected since other constraints are given');
 114   else
 115     if (any (lidx))
 116       warning ('some initial parameters set to lower bound');
 117       pin(lidx, 1) = lbound(lidx, 1);
 118     end
 119     if (any (uidx))
 120       warning ('some initial parameters set to upper bound');
 121       pin(uidx, 1) = ubound(uidx, 1);
 122     end
 123   end
 124   if (n_gencstr > 0 && any (~isinf (maxstep)))
 125     warning ('setting both a maximum fractional step change of parameters and general constraints may result in inefficiency and failure');
 126   end
 127
 128   %% fill constant fields of hook for derivative-functions; some fields
 129   %% may be backend-specific
 130   dfdp_hook.fixed = fixed; % this may be handled by the frontend, but
 131                                 % the backend still may add to it
 132
 133   %% set up for iterations
 134   %%
 135   p = pin;
 136   f = f_pin; fbest=f; pbest=p;
 137   m = prod (size (f));
 138   r = wt .* f;
 139   r = r(:);
 140   if (~isreal (r)) error ('weighted residuals are not real'); end
 141   ss = r.' * r;
 142   sbest=ss;
 143   chgprev=Inf*ones(n,1);
 144   cvg=0;
 145   epsLlast=1;
 146   epstab=[.1, 1, 1e2, 1e4, 1e6];
 147   ac_idx = true (n_lcstr + n_gencstr, 1); % all constraints
 148   nc_idx = false (n_lcstr + n_gencstr, 1); % none of all constraints
 149   gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); % gen. constr.
 150   lc_idx = ~gc_idx;
 151
 152   %% do iterations
 153   %%
 154   for iter = 1:niter
 155     deb_printf (testing, '\nstart outer iteration\n');
 156     v_cstr = f_cstr (p, ac_idx);
 157     %% index of active constraints
 158     c_act =  v_cstr < nz | eq_idx; # equality constraints might be
 159                                 # violated at start
 160     if (any (c_act))
 161       if (n_gencstr > 0)
 162         %% full gradient is needed later
 163         dct = df_cstr (p, ac_idx, ...
 164                        setfield (dfdp_hook, 'f', v_cstr));
 165         dct(:, fixed) = 0; % for user supplied dfdp; necessary?
 166         dcat = dct(c_act, :);
 167       else
 168         dcat = df_cstr (p, c_act, ...
 169                         setfield (dfdp_hook, 'f', v_cstr));
 170         dcat(:, fixed) = 0; % for user supplied dfdp; necessary?
 171       end
 172       dca = dcat.';
 173     end
 174     nrm = zeros (1, n);
 175     pprev=pbest;
 176     prt = dfdp (p, setfield (dfdp_hook, 'f', fbest(:)));
 177     prt(:, fixed) = 0; % for user supplied dfdp; necessary?
 178     r = wt .* -fbest;
 179     r = r(:);
 180     if (~isreal (r)) error ('weighted residuals are not real'); end
 181     sprev=sbest;
 182     sgoal=(1-stol)*sprev;
 183     msk = ~fixed;
 184     prt(:, msk) = prt(:, msk) .* wtl(:, ones (1, sum (msk)));
 185     nrm(msk) = sumsq (prt(:, msk), 1);
 186     msk = nrm > 0;
 187     nrm(msk) = 1 ./ sqrt (nrm(msk));
 188     prt = prt .* nrm(ones (1, m), :);
 189     nrm = nrm.';
 190     [prt,s,v]=svd(prt,0);
 191     s=diag(s);
 192     g = prt.' * r;
 193     for jjj=1:length(epstab)
 194       deb_printf (testing, '\nstart inner iteration\n');
 195       epsL = max(epsLlast*epstab(jjj),1e-7);
 196       %% printf ('epsL: %e\n', epsL); % for testing
 197
 198       %% Usage of this 'ser' later is equivalent to pre-multiplying the
 199       %% gradient with a positive-definit matrix, but not with a
 200       %% diagonal matrix, at epsL -> Inf; so there is a fallback to
 201       %% gradient descent, but not in general to descent for each
 202       %% gradient component. Using the commented-out 'ser' ((1 / (1 +
 203       %% epsL^2)) * (1 ./ se + epsL * s)) would be equivalent to using
 204       %% Marquardts diagonal of the Hessian-approximation for epsL ->
 205       %% Inf, but currently this gives no advantages in tests, even with
 206       %% constraints.
 207 %%% ser = 1 ./ sqrt((s.*s)+epsL);
 208       se = sqrt ((s.*s) + epsL);
 209       if (new_s)
 210         %% for testing
 211         ser = (1 / (1 + epsL^2)) * (1 ./ se + epsL * s);
 212       else
 213         ser = 1 ./ se;
 214       end
 215       tp1 = (v * (g .* ser)) .* nrm;
 216       if (any (c_act))
 217         deb_printf (testing, 'constraints are active:\n');
 218         deb_printf (testing, '%i\n', c_act);
 219         %% calculate chg by 'quadratic programming'
 220         nrme= diag (nrm);
 221         ser2 = diag (ser .* ser);
 222         mfc1 = nrme * v * ser2 * v.' * nrme;
 223         tp2 = mfc1 * dca;
 224         a_eq_idx = eq_idx(c_act);
 225         [lb, bidx, ridx, tbl] = cpiv (dcat * tp1, dcat * tp2, a_eq_idx);
 226         chg = tp1 + tp2(:, bidx) * lb; % if a parameter is 'fixed',
 227                                 % the respective component of chg should
 228                                 % be zero too, even here (with active
 229                                 % constraints)
 230         deb_printf (testing, 'change:\n');
 231         deb_printf (testing, '%e\n', chg);
 232         deb_printf (testing, '\n');
 233         %% indices for different types of constraints
 234         c_inact = ~c_act; % inactive constraints
 235         c_binding = nc_idx;
 236         c_binding(c_act) = bidx; % constraints selected binding
 237         c_unbinding = nc_idx;
 238         c_unbinding(c_act) = ridx; % constraints unselected binding
 239         c_nonbinding = c_act & ~(c_binding | c_unbinding); % constraints
 240                                 % selected non-binding
 241       else
 242         %% chg is the Levenberg/Marquardt step
 243         chg = tp1;
 244         %% indices for different types of constraints
 245         c_inact = ac_idx; % inactive constraints consist of all
 246                                 % constraints
 247         c_binding = nc_idx;
 248         c_unbinding = nc_idx;
 249         c_nonbinding = nc_idx;
 250       end
 251       %% apply constraints to step width (since this is a
 252       %% Levenberg/Marquardt algorithm, no line-search is performed
 253       %% here)
 254       k = 1;
 255       c_tp = c_inact(1:n_lcstr);
 256       mcit = mc(:, c_tp).';
 257       vci = vc(c_tp);
 258       hstep = mcit * chg;
 259       idx = hstep < 0;
 260       if (any (idx))
 261         k = min (1, min (- (vci(idx) + mcit(idx, :) * pprev) ./ ...
 262                          hstep(idx)));
 263       end
 264       if (k < 1)
 265         deb_printf (testing, 'stepwidth: linear constraints\n');
 266       end
 267       if (n_gencstr > 0)
 268         c_tp = gc_idx & (c_nonbinding | c_inact);
 269         if (any (c_tp) && any (f_cstr (pprev + k * chg, c_tp) < 0))
 270           [k, fval, info] = ...
 271               fzero (@ (x) min (cat (1, ...
 272                                      f_cstr (pprev + x * chg, c_tp), ...
 273                                      k - x, ...
 274                                      ifelse (x < 0, -Inf, Inf))), ...
 275                      0);
 276           if (info ~= 1 || abs (fval) >= nz)
 277             error ('could not find stepwidth to satisfy inactive and non-binding general inequality constraints');
 278           end
 279           deb_printf (testing, 'general constraints limit stepwidth\n');
 280         end
 281       end
 282       chg = k * chg;
 283
 284       if (any (gc_idx & c_binding)) % none selected binding =>
 285                                 % none unselected binding
 286         deb_printf (testing, 'general binding constraints must be regained:\n');
 287         %% regain binding constraints and one of the possibly active
 288         %% previously inactive or non-binding constraints
 289         ptp1 = pprev + chg;
 290
 291         tp = true;
 292         nt_nosuc = true;
 293         lim = 20;
 294         while (nt_nosuc && lim >= 0)
 295           deb_printf (testing, 'starting from new value of p in regaining:\n');
 296           deb_printf (testing, '%e\n', ptp1);
 297           %% we keep d_p.' * inv (mfc1) * d_p minimal in each step of
 298           %% the inner loop; this is both sensible (this metric
 299           %% considers a guess of curvature of sum of squared residuals)
 300           %% and convenient (we have useful matrices available for it)
 301           c_tp0 = c_inact | c_nonbinding;
 302           c_tp1 = c_inact | (gc_idx & c_nonbinding);
 303           btbl = tbl(bidx, bidx);
 304           c_tp2 = c_binding;
 305           if (any (tp)) % if none before, does not get true again
 306             tp = f_cstr (ptp1, c_tp1) < nz;
 307             if (any (tp)) % could be less clumsy, but ml-compatibility..
 308               %% keep only the first true entry in tp
 309               tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
 310               %% supplement binding index with one (the first) getting
 311               %% binding in c_tp1
 312               c_tp2(c_tp1) = tp;
 313               %% gradient of this added constraint
 314               caddt = dct(c_tp2 & ~c_binding, :);
 315               cadd = caddt.';
 316               C = dct(c_binding, :) * mfc1 * cadd;
 317               Ct = C.';
 318               G = [btbl, btbl * C; ...
 319                    -Ct * btbl, caddt * mfc1 * cadd - Ct * btbl * C];
 320               btbl = gjp (G, size (G, 1));
 321             end
 322           end
 323           dcbt = dct(c_tp2, :);
 324           mfc = - mfc1 * dcbt.' * btbl;
 325           deb_printf (testing, 'constraints to regain:\n');
 326           deb_printf (testing, '%i\n', c_tp2);
 327
 328           ptp2 = ptp1;
 329           nt_niter_start = 100;
 330           nt_niter = nt_niter_start;
 331           while (nt_nosuc && nt_niter >= 0)
 332             hv = f_cstr (ptp2, c_tp2);
 333             if (all (abs (hv) < nz))
 334               nt_nosuc = false;
 335               chg = ptp2 - pprev;
 336             else
 337               ptp2 = ptp2 + mfc * hv; % step should be zero for each
 338                                 % component for which the parameter is
 339                                 % 'fixed'
 340             end
 341             nt_niter = nt_niter - 1;
 342           end
 343           deb_printf (testing, 'constraints after regaining:\n');
 344           deb_printf (testing, '%e\n', hv);
 345           if (nt_nosuc || ...
 346               any (abs (chg) > abs (pprev .* maxstep)) || ...
 347               any (f_cstr (ptp2, c_tp0) < -nz))
 348             if (nt_nosuc)
 349               deb_printf (testing, 'regaining did not converge\n');
 350             else
 351               deb_printf (testing, 'regaining violated type 3 and 4\n');
 352             end
 353             nt_nosuc = true;
 354             ptp1 = (pprev + ptp1) / 2;
 355           end
 356           if (~nt_nosuc)
 357             tp = f_cstr (ptp2, c_unbinding);
 358             if (any (tp) < 0) % again ml-compatibility clumsyness..
 359               [discarded, id] = min(tp);
 360               tid = find (ridx);
 361               id = tid(id); % index within active constraints
 362               unsuccessful_exchange = false;
 363               if (abs (tbl(id, id)) < nz) % Bard: not absolute value
 364                 %% exchange this unselected binding constraint against a
 365                 %% binding constraint, but not against an equality
 366                 %% constraint
 367                 tbidx = bidx & ~a_eq_idx;
 368                 if (~any (tbidx))
 369                   unsuccessful_exchange = true;
 370                 else
 371                   [discarded, idm] = max (abs (tbl(tbidx, id)));
 372                   tid = find (tbidx);
 373                   idm = tid(idm); % -> index within active constraints
 374                   tbl = gjp (tbl, idm);
 375                   bidx(idm) = false;
 376                   ridx(idm) = true;
 377                 end
 378               end
 379               if (unsuccessful_exchange)
 380                 %% It probably doesn't look good now; this desperate
 381                 %% last attempt is not in the original algortithm, since
 382                 %% that didn't account for equality constraints.
 383                 ptp1 = (pprev + ptp1) / 2;
 384               else
 385                 tbl = gjp (tbl, id);
 386                 bidx(id) = true;
 387                 ridx(id) = false;
 388                 c_binding = nc_idx;
 389                 c_binding(c_act) = bidx;
 390                 c_unbinding = nc_idx;
 391                 c_unbinding(c_act) = ridx;
 392               end
 393               nt_nosuc = true;
 394               deb_printf (testing, 'regaining violated type 2\n');
 395             end
 396           end
 397           if (~nt_nosuc)
 398             deb_printf (testing, 'regaining successful, converged with %i iterations:\n', ...
 399             nt_niter_start - nt_niter);
 400             deb_printf (testing, '%e\n', ptp2);
 401           end
 402           lim = lim - 1;
 403         end
 404         if (nt_nosuc)
 405           error ('could not regain binding constraints');
 406         end
 407       else
 408         %% check the maximal stepwidth and apply as necessary
 409         ochg=chg;
 410         idx = ~isinf(maxstep);
 411         limit = abs(maxstep(idx).*pprev(idx));
 412         chg(idx) = min(max(chg(idx),-limit),limit);
 413         if (verbose && any(ochg ~= chg))
 414           disp(['Change in parameter(s): ', ...
 415                 sprintf('%d ',find(ochg ~= chg)), 'maximal fractional stepwidth enforced']);
 416         end
 417       end
 418       aprec=abs(pprec.*pbest);       %---
 419       %% ss=scalar sum of squares=sum((wt.*f)^2).
 420       if (any(abs(chg) > 0.1*aprec))%---  % only worth evaluating
 421                                 % function if there is some non-miniscule
 422                                 % change
 423         %% In the code of the outer loop before the inner loop pbest is
 424         %% actually identical to p, since once they deviate, the outer
 425         %% loop will not be repeated. Though the inner loop can still be
 426         %% repeated in this case, pbest is not used in it. Since pprev
 427         %% is set from pbest in the outer loop before the inner loop, it
 428         %% is also identical to p up to here.
 429         p=chg+pprev;
 430         %% since the projection method may have slightly violated
 431         %% constraints due to inaccuracy, correct parameters to bounds
 432         %% --- but only if no further constraints are given, otherwise
 433         %% the inaccuracy in honoring them might increase by this
 434         if (~have_constraints_except_bounds)
 435           lidx = p < lbound;
 436           uidx = p > ubound;
 437           p(lidx, 1) = lbound(lidx, 1);
 438           p(uidx, 1) = ubound(uidx, 1);
 439           chg(lidx, 1) = p(lidx, 1) - pprev(lidx, 1);
 440           chg(uidx, 1) = p(uidx, 1) - pprev(uidx, 1);
 441         end
 442         %%
 443         f = F (p);
 444         r = wt .* f;
 445         r = r(:);
 446         if (~isreal (r))
 447           error ('weighted residuals are not real');
 448         end
 449         ss = r.' * r;
 450         deb_printf (testing, 'sbest: %.16e\n', sbest);
 451         deb_printf (testing, 'sgoal: %.16e\n', sgoal);
 452         deb_printf (testing, '   ss: %.16e\n', ss);
 453         if (ss<sbest)
 454           pbest=p;
 455           fbest=f;
 456           sbest=ss;
 457         end
 458         if (ss<=sgoal)
 459           break;
 460         end
 461       end                          %---
 462     end
 463     %% printf ('epsL no.: %i\n', jjj); % for testing
 464     epsLlast = epsL;
 465     if (verbose)
 466       hook.plot_cmd (f);
 467     end
 468     if (ss < eps) % in this case ss == sbest
 469       cvg = 3; % there is no more suitable flag for this
 470       break;
 471     end
 472     if (ss>sgoal)
 473       cvg = 3;
 474       break;
 475     end
 476     aprec=abs(pprec.*pbest);
 477     %% [aprec, chg, chgprev]
 478     if (all(abs(chg) <= aprec) && all(abs(chgprev) <= aprec))
 479       cvg = 2;
 480       if (verbose)
 481         fprintf('Parameter changes converged to specified precision\n');
 482       end
 483       break;
 484     else
 485       chgprev=chg;
 486     end
 487   end
 488
 489   %% set further return values
 490   %%
 491   p = pbest;
 492   resid = fbest;
 493   outp.niter = iter;
 494
 495 function deb_printf (do_printf, varargin)
 496
 497   %% for testing
 498
 499   if (do_printf)
 500     printf (varargin{:})
 501   end
 502
 503 function fval = scalar_ifelse (cond, tval, fval)
 504
 505   %% needed for some anonymous functions, builtin ifelse only available
 506   %% in Octave > 3.2; we need only the scalar case here
 507
 508   if (cond)
 509     fval = tval;
 510   end