1 ## Copyright (C) 2012 Olaf Till <i7tiol@t-online.de>
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.
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.
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/>.
16 function [p_res, objf, cvg, outp] = __lm_feasible__ (f, pin, hook)
18 ## some backend specific defaults
19 fract_prec_default = 0;
20 max_fract_step_default = Inf;
22 ## needed for some anonymous functions
23 if (exist ("ifelse") != 5)
24 ifelse = @ scalar_ifelse;
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
37 lbound = hook.lbound; # bounds, subset of linear inequality
38 ubound = hook.ubound; # constraints in mc and vc
40 ## passed values of constraints for initial parameters
41 pin_cstr = hook.pin_cstr;
43 ## passed return value of f for initial parameters
46 ## passed function for gradient of objective function
49 ## passed function for hessian of objective function
50 if (isempty (hessian = hook.hessian))
57 ## passed function for complementary pivoting
62 if (isempty (niter = hook.MaxIter)) niter = 20; endif
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");
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.
79 nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
81 ## backend-specific checking of options and constraints
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.");
91 idx = lbound == ubound;
93 warning ("lower and upper bounds identical for some parameters, fixing the respective parameters");
97 error ("no free parameters");
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");
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
108 ## set up for iterations
114 ltab = [.1, 1, 1e2, 1e4, 1e6];
115 chgprev = Inf (n, 1);
117 c_act = false (n, 1);
124 ## gradient of objective function
126 df = grad_f (p, setfield (dfdp_hook, "f", f))(:);
128 ## constraints, preparation of some constants
129 v_cstr = f_cstr (p, ac_idx);
132 c_act = v_cstr < nz | eq_idx; # index of active constraints
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, :);
141 dcat = df_cstr (p, c_act, setfield (dfdp_hook, "f", v_cstr));
142 dcat(:, fixed) = 0; # for user supplied dfdp; necessary?
147 a_eq_idx = eq_idx(c_act);
155 ## hessian of objectiv function
161 if (any (isnan (A(:))))
162 error ("some second derivatives undefined by user function");
165 error ("second derivatives given by user function not real");
167 if (! issymmetric (A))
168 error ("Hessian returned by user function not symmetric");
175 ## approximate Hessian of Lagrangian
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.
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
192 dca_new_id = dca(:, id_new);
193 ## gradients of previously active constraints which are also
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);
203 ## difference of derivatives of new and old active constraints,
204 ## multiplied by multipliers, as used for BFGS update (lb set
206 dch = (dca_new_id(:, bidx_old_id) - \
207 dca_old_id(:, bidx_old_id)) * \
210 y = df - old_df - dch;
212 ## Damped BFGS according to Nocedal & Wright, 2nd edition,
215 sAs = chgt * A * chg;
220 if ((den1 = sAs - cy) == 0)
224 th = .8 * sAs / den1;
227 r = th * y + (1 - th) * Ac;
229 if ((den2 = chgt * r) == 0 || sAs == 0)
233 A += r * r.' / den2 - Ac * Ac.' / sAs;
239 ## Inverse scaled decomposition A = G * (1 ./ L) * G.'
241 ## make matrix Binv for scaling
243 nidx = ! (idx = Binv == 0);
244 Binv(nidx) = 1 ./ sqrt (abs (Binv(nidx)));
247 ## eigendecomposition of scaled A
248 [V, L] = eig (Binv * A * Binv);
250 ## A is symmetric, so V and L are real, delete any imaginary parts,
251 ## which might occur due to inaccuracy
255 nminL = - min (L) * 1.1 / ltab(1);
258 ## Levenberg/Marquardt
259 fgoal = (1 - ftol) * vf;
262 ll = max (ll, nminL);
263 l = max (1e-7, ll * l);
265 R = G * diag (1 ./ (L + l)) * G.';
270 ## some constraints are active, quadratic programming
273 [lb, bidx, ridx, tbl] = cpiv (- tp * df, tp * dca, a_eq_idx);
274 chg = R * (dca(:, bidx) * lb - df); # step direction
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
286 ## no constraints are active, chg is the Levenberg/Marquardt step
288 chg = - R * df; # step direction
293 ## indices for different types of constraints (meaning see above)
296 c_unbinding = nc_idx;
297 c_nonbinding = nc_idx;
301 ## apply inactive and non-binding constraints to step width
303 ## linear constraints
305 c_tp = c_inact(1:n_lcstr);
306 mcit = mc(:, c_tp).';
311 k = min (1, min (- (vci(idx) + mcit(idx, :) * p) ./ \
315 ## general constraints
317 c_tp = gc_idx & (c_nonbinding | c_inact);
318 if (any (c_tp) && any (f_cstr (p + k * chg, c_tp) < 0))
320 fzero (@ (x) min (cat (1, \
321 f_cstr (p + x * chg, c_tp), \
323 ifelse (x < 0, -Inf, Inf))), \
325 if (info != 1 || abs (fval) >= nz)
326 error ("could not find stepwidth to satisfy inactive and non-binding general inequality constraints");
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
342 while (nt_nosuc && lim >= 0)
343 ## we keep d_p.' * inv (R) * d_p minimal in each step of the
345 c_tp0 = c_inact | c_nonbinding;
346 c_tp1 = c_inact | (gc_idx & c_nonbinding);
347 btbl = tbl(bidx, bidx);
349 ## once (any(tp)==false), it would not get true again even
350 ## with the following assignment
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
358 ## gradient of this added constraint
359 caddt = dct(c_tp2 & ! c_binding, :);
361 C = dct(c_binding, :) * R * cadd;
363 T = [btbl, btbl * C; \
364 -Ct * btbl, caddt * R * cadd - Ct * btbl * C];
365 btbl = gjp (T, size (T, 1));
367 dcbt = dct(c_tp2, :);
368 mfc = - R * dcbt.' * btbl;
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))
378 ptp2 = ptp2 + mfc * hv; # step should be zero for each
379 # component for which the parameter is
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.
391 ptp1 = (p + ptp1) / 2;
394 any ((tp = f_cstr (ptp2, c_unbinding)) < 0))
395 [discarded, id] = min(tp);
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
403 tbidx = bidx & ! a_eq_idx;
405 unsuccessful_exchange = true;
407 [discarded, idm] = max (abs (tbl(tbidx, id)));
409 idm = tid(idm); # -> index within active constraints
410 tbl = gjp (tbl, idm);
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;
425 c_binding(c_act) = bidx;
426 c_unbinding = nc_idx;
427 c_unbinding(c_act) = ridx;
429 ## regaining violated type 2 constraints
435 error ("could not regain binding constraints");
438 ## check the maximal stepwidth and apply as necessary
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)));
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
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);
467 if (! isreal (vf_chg = f (p_chg)))
468 error ("objective function not real");
484 aprec = abs (pprec .* pbest);
485 if (vf_chg < eps || vf_chg > fgoal)
488 elseif (all (abs (chg) <= aprec) && all (abs (chgprev) <= aprec))
491 elseif (iter == niter)