octave_packages/optim-1.2.0/private/__lm_feasible__.m

   1 ## Copyright (C) 2012 Olaf Till <i7tiol@t-online.de>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 function [p_res, objf, cvg, outp] = __lm_feasible__ (f, pin, hook)
  17
  18   ## some backend specific defaults
  19   fract_prec_default = 0;
  20   max_fract_step_default = Inf;
  21
  22   ## needed for some anonymous functions
  23   if (exist ("ifelse") != 5)
  24     ifelse = @ scalar_ifelse;
  25   endif
  26
  27   n = length (pin);
  28
  29   ## passed constraints
  30   mc = hook.mc; # matrix of linear constraints
  31   vc = hook.vc; # vector of linear constraints
  32   f_cstr = hook.f_cstr; # function of all constraints
  33   df_cstr = hook.df_cstr; # function of derivatives of all constraints
  34   n_gencstr = hook.n_gencstr; # number of non-linear constraints
  35   eq_idx = hook.eq_idx; # logical index of equality constraints in all
  36                                 # constraints
  37   lbound = hook.lbound; # bounds, subset of linear inequality
  38   ubound = hook.ubound; # constraints in mc and vc
  39
  40   ## passed values of constraints for initial parameters
  41   pin_cstr = hook.pin_cstr;
  42
  43   ## passed return value of f for initial parameters
  44   f_pin = hook.f_pin;
  45
  46   ## passed function for gradient of objective function
  47   grad_f = hook.dfdp;
  48
  49   ## passed function for hessian of objective function
  50   if (isempty (hessian = hook.hessian))
  51     user_hessian = false;
  52     A = eye (n);
  53   else
  54     user_hessian = true;
  55   endif
  56
  57   ## passed function for complementary pivoting
  58   cpiv = hook.cpiv;
  59
  60   ## passed options
  61   ftol = hook.TolFun;
  62   if (isempty (niter = hook.MaxIter)) niter = 20; endif
  63   fixed = hook.fixed;
  64   maxstep = hook.max_fract_change;
  65   maxstep(isna (maxstep)) = max_fract_step_default;
  66   pprec = hook.fract_prec;
  67   pprec(isna (pprec)) = fract_prec_default;
  68   verbose = strcmp (hook.Display, "iter");
  69
  70   ## some useful variables derived from passed variables
  71   n_lcstr = size (vc, 1);
  72   have_constraints_except_bounds = \
  73       n_lcstr + n_gencstr > \
  74       sum (lbound != -Inf) + sum (ubound != Inf);
  75   ac_idx = true (n_lcstr + n_gencstr, 1); # index of all constraints
  76   nc_idx = false (n_lcstr + n_gencstr, 1); # none of all constraints
  77   gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); # gen. constr.
  78
  79   nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
  80
  81   ## backend-specific checking of options and constraints
  82   ##
  83   if (any (pin < lbound | pin > ubound) ||
  84       any (pin_cstr.inequ.lin_except_bounds < 0) ||
  85       any (pin_cstr.inequ.gen < 0) ||
  86       any (abs (pin_cstr.equ.lin)) >= nz ||
  87       any (abs (pin_cstr.equ.gen)) >= nz)
  88     error ("Initial parameters violate constraints.");
  89   endif
  90   ##
  91   idx = lbound == ubound;
  92   if (any (idx))
  93     warning ("lower and upper bounds identical for some parameters, fixing the respective parameters");
  94     fixed(idx) = true;
  95   endif
  96   if (all (fixed))
  97     error ("no free parameters");
  98   endif
  99   if (n_gencstr > 0 && any (! isinf (maxstep)))
 100     warning ("setting both a maximum fractional step change of parameters and general constraints may result in inefficiency and failure");
 101   endif
 102
 103   ## fill constant fields of hook for derivative-functions; some fields
 104   ## may be backend-specific
 105   dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
 106                                 # the backend still may add to it
 107
 108   ## set up for iterations
 109   p = pbest = pin;
 110   vf = fbest = f_pin;
 111   iter = 0;
 112   done = false;
 113   ll = 1;
 114   ltab = [.1, 1, 1e2, 1e4, 1e6];
 115   chgprev = Inf (n, 1);
 116   df = [];
 117   c_act = false (n, 1);
 118   dca = zeros (n, 0);
 119
 120   while (! done)
 121
 122     iter++;
 123
 124     ## gradient of objective function
 125     old_df = df;
 126     df = grad_f (p, setfield (dfdp_hook, "f", f))(:);
 127
 128     ## constraints, preparation of some constants
 129     v_cstr = f_cstr (p, ac_idx);
 130     old_c_act = c_act;
 131     old_dca = dca;
 132     c_act =  v_cstr < nz | eq_idx; # index of active constraints
 133     if (any (c_act))
 134
 135       if (n_gencstr)
 136         ## full gradient is needed later
 137         dct = df_cstr (p, ac_idx, setfield (dfdp_hook, "f", v_cstr));
 138         dct(:, fixed) = 0; # for user supplied dfdp; necessary?
 139         dcat = dct(c_act, :);
 140       else
 141         dcat = df_cstr (p, c_act, setfield (dfdp_hook, "f", v_cstr));
 142         dcat(:, fixed) = 0; # for user supplied dfdp; necessary?
 143       endif
 144
 145       dca = dcat.';
 146
 147       a_eq_idx = eq_idx(c_act);
 148
 149     else
 150
 151       dca = zeros (n, 0);
 152
 153     endif
 154
 155     ## hessian of objectiv function
 156     if (user_hessian)
 157
 158       A = hessian (p);
 159       idx = isnan (A);
 160       A(idx) = A.'(idx);
 161       if (any (isnan (A(:))))
 162         error ("some second derivatives undefined by user function");
 163       endif
 164       if (! isreal (A))
 165         error ("second derivatives given by user function not real");
 166       endif
 167       if (! issymmetric (A))
 168         error ("Hessian returned by user function not symmetric");
 169       endif
 170
 171     elseif (iter > 1)
 172
 173       if (any (chg))
 174
 175         ## approximate Hessian of Lagrangian
 176
 177         ## I wonder if this hassle here and above with accounting for
 178         ## changing active sets is indeed better than just approximating
 179         ## the Hessian only of the objective function.
 180         ##
 181         ## index, over all constraints, of constraints active both
 182         ## previously and currently
 183         s_c_act = old_c_act & c_act;
 184         ## index, over currently active constraints, of constraints
 185         ## active both previously and currently
 186         id_new = s_c_act(c_act);
 187         ## index, over previously active constraints, of constraints
 188         ## active both previously and currently
 189         id_old = s_c_act(old_c_act);
 190         ## gradients of currently active constraints which were also
 191         ## active previously
 192         dca_new_id = dca(:, id_new);
 193         ## gradients of previously active constraints which are also
 194         ## active currently
 195         dca_old_id = old_dca(:, id_old);
 196         ## index, over constraints active both previously and currently,
 197         ## of (old) non-zero multipliers (bidx set below previously)
 198         bidx_old_id = bidx(id_old);
 199         ## index, over (old) non-zero multipliers, of constraints active
 200         ## both previously and currently (bidx set below previously)
 201         old_l_idx = id_old(bidx);
 202
 203         ## difference of derivatives of new and old active constraints,
 204         ## multiplied by multipliers, as used for BFGS update (lb set
 205         ## below previously)
 206         dch = (dca_new_id(:, bidx_old_id) - \
 207                dca_old_id(:, bidx_old_id)) * \
 208             lb(old_l_idx);
 209
 210         y = df - old_df - dch;
 211
 212         ## Damped BFGS according to Nocedal & Wright, 2nd edition,
 213         ## procedure 18.2.
 214         chgt = chg.';
 215         sAs = chgt * A * chg;
 216         cy = chgt * y;
 217         if (cy >= .2 * sAs)
 218           th = 1;
 219         else
 220           if ((den1 = sAs - cy) == 0)
 221             cvg = -4;
 222             break;
 223           endif
 224           th = .8 * sAs / den1;
 225         endif
 226         Ac = A * chg;
 227         r = th * y + (1 - th) * Ac;
 228
 229         if ((den2 = chgt * r) == 0 || sAs == 0)
 230           cvg = -4;
 231           break;
 232         endif
 233         A += r * r.' / den2 - Ac * Ac.' / sAs;
 234
 235       endif
 236
 237     endif
 238
 239     ## Inverse scaled decomposition A = G * (1 ./ L) * G.'
 240     ##
 241     ## make matrix Binv for scaling
 242     Binv = diag (A);
 243     nidx = ! (idx = Binv == 0);
 244     Binv(nidx) = 1 ./ sqrt (abs (Binv(nidx)));
 245     Binv(idx) = 1;
 246     Binv = diag (Binv);
 247     ## eigendecomposition of scaled A
 248     [V, L] = eig (Binv * A * Binv);
 249     L = diag (L);
 250     ## A is symmetric, so V and L are real, delete any imaginary parts,
 251     ## which might occur due to inaccuracy
 252     V = real (V);
 253     L = real (L);
 254     ##
 255     nminL = - min (L) * 1.1 / ltab(1);
 256     G = Binv * V;
 257
 258     ## Levenberg/Marquardt
 259     fgoal = (1 - ftol) * vf;
 260     for l = ltab
 261
 262       ll = max (ll, nminL);
 263       l = max (1e-7, ll * l);
 264
 265       R = G * diag (1 ./ (L + l)) * G.';
 266
 267       ## step computation
 268       if (any (c_act))
 269
 270         ## some constraints are active, quadratic programming
 271
 272         tp = dcat * R;
 273         [lb, bidx, ridx, tbl] = cpiv (- tp * df, tp * dca, a_eq_idx);
 274         chg = R * (dca(:, bidx) * lb - df); # step direction
 275
 276         ## indices for different types of constraints
 277         c_inact = ! c_act; # inactive constraints
 278         c_binding = c_unbinding = nc_idx;
 279         c_binding(c_act) = bidx; # constraints selected binding
 280         c_unbinding(c_act) = ridx; # constraints unselected binding
 281         c_nonbinding = c_act & ! (c_binding | c_unbinding); #
 282                                 #constraints selected non-binding
 283
 284       else
 285
 286         ## no constraints are active, chg is the Levenberg/Marquardt step
 287
 288         chg = - R * df; # step direction
 289
 290         lb = zeros (0, 1);
 291         bidx = false (0, 1);
 292
 293         ## indices for different types of constraints (meaning see above)
 294         c_inact = ac_idx;
 295         c_binding = nc_idx;
 296         c_unbinding = nc_idx;
 297         c_nonbinding = nc_idx;
 298
 299       endif
 300
 301       ## apply inactive and non-binding constraints to step width
 302       ##
 303       ## linear constraints
 304       k = 1;
 305       c_tp = c_inact(1:n_lcstr);
 306       mcit = mc(:, c_tp).';
 307       vci = vc(c_tp);
 308       hstep = mcit * chg;
 309       idx = hstep < 0;
 310       if (any (idx))
 311         k = min (1, min (- (vci(idx) + mcit(idx, :) * p) ./ \
 312                          hstep(idx)));
 313       endif
 314       ##
 315       ## general constraints
 316       if (n_gencstr)
 317         c_tp = gc_idx & (c_nonbinding | c_inact);
 318         if (any (c_tp) && any (f_cstr (p + k * chg, c_tp) < 0))
 319           [k, fval, info] = \
 320               fzero (@ (x) min (cat (1, \
 321                                      f_cstr (p + x * chg, c_tp), \
 322                                      k - x, \
 323                                      ifelse (x < 0, -Inf, Inf))), \
 324                      0);
 325           if (info != 1 || abs (fval) >= nz)
 326             error ("could not find stepwidth to satisfy inactive and non-binding general inequality constraints");
 327           endif
 328         endif
 329       endif
 330       ##
 331       chg = k * chg;
 332
 333       ## if necessary, regain binding constraints and one of the
 334       ## possibly active previously inactive or non-binding constraints
 335       if (any (gc_idx & c_binding)) # none selected binding => none
 336                                 # unselected binding
 337         ptp1 = p + chg;
 338
 339         tp = true;
 340         nt_nosuc = true;
 341         lim = 20;
 342         while (nt_nosuc && lim >= 0)
 343           ## we keep d_p.' * inv (R) * d_p minimal in each step of the
 344           ## inner loop
 345           c_tp0 = c_inact | c_nonbinding;
 346           c_tp1 = c_inact | (gc_idx & c_nonbinding);
 347           btbl = tbl(bidx, bidx);
 348           c_tp2 = c_binding;
 349           ## once (any(tp)==false), it would not get true again even
 350           ## with the following assignment
 351           if (any (tp) && \
 352               any (tp = f_cstr (ptp1, c_tp1) < nz))
 353             ## keep only the first true entry in tp
 354             tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
 355             ## supplement binding index with one (the first) getting
 356             ## binding in c_tp1
 357             c_tp2(c_tp1) = tp;
 358             ## gradient of this added constraint
 359             caddt = dct(c_tp2 & ! c_binding, :);
 360             cadd = caddt.';
 361             C = dct(c_binding, :) * R * cadd;
 362             Ct = C.';
 363             T = [btbl, btbl * C; \
 364                  -Ct * btbl, caddt * R * cadd - Ct * btbl * C];
 365             btbl = gjp (T, size (T, 1));
 366           endif
 367           dcbt = dct(c_tp2, :);
 368           mfc = - R * dcbt.' * btbl;
 369
 370           ptp2 = ptp1;
 371           nt_niter = nt_niter_start = 100;
 372           while (nt_nosuc && nt_niter >= 0)
 373             hv = f_cstr (ptp2, c_tp2);
 374             if (all (abs (hv) < nz))
 375               nt_nosuc = false;
 376               chg = ptp2 - p;
 377             else
 378               ptp2 = ptp2 + mfc * hv; # step should be zero for each
 379                                 # component for which the parameter is
 380                                 # "fixed"
 381             endif
 382             nt_niter--;
 383           endwhile
 384
 385           if (nt_nosuc || \
 386               any (abs (chg) > abs (p .* maxstep)) || \
 387               any (f_cstr (ptp2, c_tp0) < -nz))
 388             ## if (nt_nosuc), regaining did not converge, else,
 389             ## regaining violated type 3 and 4.
 390             nt_nosuc = true;
 391             ptp1 = (p + ptp1) / 2;
 392           endif
 393           if (! nt_nosuc && \
 394               any ((tp = f_cstr (ptp2, c_unbinding)) < 0))
 395             [discarded, id] = min(tp);
 396             tid = find (ridx);
 397             id = tid(id); # index within active constraints
 398             unsuccessful_exchange = false;
 399             if (abs (tbl(id, id)) < nz) # Bard: not absolute value
 400               ## exchange this unselected binding constraint against a
 401               ## binding constraint, but not against an equality
 402               ## constraint
 403               tbidx = bidx & ! a_eq_idx;
 404               if (! any (tbidx))
 405                 unsuccessful_exchange = true;
 406               else
 407                 [discarded, idm] = max (abs (tbl(tbidx, id)));
 408                 tid = find (tbidx);
 409                 idm = tid(idm); # -> index within active constraints
 410                 tbl = gjp (tbl, idm);
 411                 bidx(idm) = false;
 412                 ridx(idm) = true;
 413               endif
 414             endif
 415             if (unsuccessful_exchange)
 416               ## It probably doesn't look good now; this desperate last
 417               ## attempt is not in the original algortithm, since that
 418               ## didn't account for equality constraints.
 419               ptp1 = (p + ptp1) / 2;
 420             else
 421               tbl = gjp (tbl, id);
 422               bidx(id) = true;
 423               ridx(id) = false;
 424               c_binding = nc_idx;
 425               c_binding(c_act) = bidx;
 426               c_unbinding = nc_idx;
 427               c_unbinding(c_act) = ridx;
 428             endif
 429             ## regaining violated type 2 constraints
 430             nt_nosuc = true;
 431           endif
 432           lim--;
 433         endwhile
 434         if (nt_nosuc)
 435           error ("could not regain binding constraints");
 436         endif
 437       else
 438         ## check the maximal stepwidth and apply as necessary
 439         ochg = chg;
 440         idx = ! isinf (maxstep);
 441         limit = abs (maxstep(idx) .* p(idx));
 442         chg(idx) = min (max (chg(idx), - limit), limit);
 443         if (verbose && any (ochg != chg))
 444           printf ("Change in parameter(s): %s:maximal fractional stepwidth enforced", \
 445                   sprintf ("%d ", find (ochg != chg)));
 446         endif
 447       endif # regaining
 448
 449       aprec = abs (pprec .* pbest);
 450       if (any (abs (chg) > 0.1 * aprec)) # only worth evaluating
 451                                 # function if there is some
 452                                 # non-miniscule change
 453         p_chg = p + chg;
 454         ## since the projection method may have slightly violated
 455         ## constraints due to inaccuracy, correct parameters to bounds
 456         ## --- but only if no further constraints are given, otherwise
 457         ## the inaccuracy in honoring them might increase by this
 458         if (! have_constraints_except_bounds)
 459           lidx = p_chg < lbound;
 460           uidx = p_chg > ubound;
 461           p_chg(lidx, 1) = lbound(lidx, 1);
 462           p_chg(uidx, 1) = ubound(uidx, 1);
 463           chg(lidx, 1) = p_chg(lidx, 1) - p(lidx, 1);
 464           chg(uidx, 1) = p_chg(uidx, 1) - p(uidx, 1);
 465         endif
 466         ##
 467         if (! isreal (vf_chg = f (p_chg)))
 468           error ("objective function not real");
 469         endif
 470         if (vf_chg < fbest)
 471           pbest = p_chg;
 472           fbest = vf_chg;
 473         endif
 474         if (vf_chg <= fgoal)
 475           p = p_chg;
 476           vf = vf_chg;
 477           break;
 478         endif
 479       endif
 480     endfor
 481
 482     ll = l;
 483
 484     aprec = abs (pprec .* pbest);
 485     if (vf_chg < eps || vf_chg > fgoal)
 486       cvg = 3;
 487       done = true;
 488     elseif (all (abs (chg) <= aprec) && all (abs (chgprev) <= aprec))
 489       cvg = 2;
 490       done = true;
 491     elseif (iter == niter)
 492       cvg = 0;
 493       done = true;
 494     else
 495       chgprev = chg;
 496     endif
 497
 498   endwhile
 499
 500   ## return result
 501   p_res = pbest;
 502   objf = fbest;
 503   outp.niter = iter;
 504
 505 endfunction