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>
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
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
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/>.
19 function [p, resid, cvg, outp] = __lm_svd__ (F, pin, hook)
21 %% This is a backend for optimization. This code was originally
22 %% contained in leasqr.m, which is now a frontend.
24 %% some backend specific defaults
25 fract_prec_default = 0;
26 max_fract_step_default = Inf;
28 %% needed for some anonymous functions
29 if (exist ('ifelse') ~= 5)
30 ifelse = @ scalar_ifelse;
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
41 lbound = hook.lbound; % bounds, subset of linear inequality
42 ubound = hook.ubound; % constraints in mc and vc
44 %% passed values of constraints for initial parameters
45 pin_cstr = hook.pin_cstr;
47 %% passed return value of F for initial parameters
50 %% passed derivative of residual function
53 %% passed function for complementary pivoting
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;
63 if (isempty (niter)) niter = 20; end
66 verbose = strcmp (hook.Display, 'iter');
68 %% only preliminary, for testing
69 if (isfield (hook, 'testing'))
70 testing = hook.testing;
74 if (isfield (hook, 'new_s'))
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);
88 nz = 20 * eps; % This is arbitrary. Constraint function will be
89 % regarded as <= zero if less than nz.
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');
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');
102 idx = lbound == ubound;
104 warning ('lower and upper bounds identical for some parameters, fixing the respective parameters');
108 error ('no free parameters');
112 if (any (lidx | uidx) && have_constraints_except_bounds)
113 warning ('initial parameters outside bounds, not corrected since other constraints are given');
116 warning ('some initial parameters set to lower bound');
117 pin(lidx, 1) = lbound(lidx, 1);
120 warning ('some initial parameters set to upper bound');
121 pin(uidx, 1) = ubound(uidx, 1);
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');
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
133 %% set up for iterations
136 f = f_pin; fbest=f; pbest=p;
140 if (~isreal (r)) error ('weighted residuals are not real'); end
143 chgprev=Inf*ones(n,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.
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
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, :);
168 dcat = df_cstr (p, c_act, ...
169 setfield (dfdp_hook, 'f', v_cstr));
170 dcat(:, fixed) = 0; % for user supplied dfdp; necessary?
176 prt = dfdp (p, setfield (dfdp_hook, 'f', fbest(:)));
177 prt(:, fixed) = 0; % for user supplied dfdp; necessary?
180 if (~isreal (r)) error ('weighted residuals are not real'); end
182 sgoal=(1-stol)*sprev;
184 prt(:, msk) = prt(:, msk) .* wtl(:, ones (1, sum (msk)));
185 nrm(msk) = sumsq (prt(:, msk), 1);
187 nrm(msk) = 1 ./ sqrt (nrm(msk));
188 prt = prt .* nrm(ones (1, m), :);
190 [prt,s,v]=svd(prt,0);
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
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
207 %%% ser = 1 ./ sqrt((s.*s)+epsL);
208 se = sqrt ((s.*s) + epsL);
211 ser = (1 / (1 + epsL^2)) * (1 ./ se + epsL * s);
215 tp1 = (v * (g .* ser)) .* nrm;
217 deb_printf (testing, 'constraints are active:\n');
218 deb_printf (testing, '%i\n', c_act);
219 %% calculate chg by 'quadratic programming'
221 ser2 = diag (ser .* ser);
222 mfc1 = nrme * v * ser2 * v.' * nrme;
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
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
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
242 %% chg is the Levenberg/Marquardt step
244 %% indices for different types of constraints
245 c_inact = ac_idx; % inactive constraints consist of all
248 c_unbinding = nc_idx;
249 c_nonbinding = nc_idx;
251 %% apply constraints to step width (since this is a
252 %% Levenberg/Marquardt algorithm, no line-search is performed
255 c_tp = c_inact(1:n_lcstr);
256 mcit = mc(:, c_tp).';
261 k = min (1, min (- (vci(idx) + mcit(idx, :) * pprev) ./ ...
265 deb_printf (testing, 'stepwidth: linear constraints\n');
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), ...
274 ifelse (x < 0, -Inf, Inf))), ...
276 if (info ~= 1 || abs (fval) >= nz)
277 error ('could not find stepwidth to satisfy inactive and non-binding general inequality constraints');
279 deb_printf (testing, 'general constraints limit stepwidth\n');
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
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);
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
313 %% gradient of this added constraint
314 caddt = dct(c_tp2 & ~c_binding, :);
316 C = dct(c_binding, :) * mfc1 * cadd;
318 G = [btbl, btbl * C; ...
319 -Ct * btbl, caddt * mfc1 * cadd - Ct * btbl * C];
320 btbl = gjp (G, size (G, 1));
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);
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))
337 ptp2 = ptp2 + mfc * hv; % step should be zero for each
338 % component for which the parameter is
341 nt_niter = nt_niter - 1;
343 deb_printf (testing, 'constraints after regaining:\n');
344 deb_printf (testing, '%e\n', hv);
346 any (abs (chg) > abs (pprev .* maxstep)) || ...
347 any (f_cstr (ptp2, c_tp0) < -nz))
349 deb_printf (testing, 'regaining did not converge\n');
351 deb_printf (testing, 'regaining violated type 3 and 4\n');
354 ptp1 = (pprev + ptp1) / 2;
357 tp = f_cstr (ptp2, c_unbinding);
358 if (any (tp) < 0) % again ml-compatibility clumsyness..
359 [discarded, id] = min(tp);
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
367 tbidx = bidx & ~a_eq_idx;
369 unsuccessful_exchange = true;
371 [discarded, idm] = max (abs (tbl(tbidx, id)));
373 idm = tid(idm); % -> index within active constraints
374 tbl = gjp (tbl, idm);
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;
389 c_binding(c_act) = bidx;
390 c_unbinding = nc_idx;
391 c_unbinding(c_act) = ridx;
394 deb_printf (testing, 'regaining violated type 2\n');
398 deb_printf (testing, 'regaining successful, converged with %i iterations:\n', ...
399 nt_niter_start - nt_niter);
400 deb_printf (testing, '%e\n', ptp2);
405 error ('could not regain binding constraints');
408 %% check the maximal stepwidth and apply as necessary
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']);
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
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.
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)
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);
447 error ('weighted residuals are not real');
450 deb_printf (testing, 'sbest: %.16e\n', sbest);
451 deb_printf (testing, 'sgoal: %.16e\n', sgoal);
452 deb_printf (testing, ' ss: %.16e\n', ss);
463 %% printf ('epsL no.: %i\n', jjj); % for testing
468 if (ss < eps) % in this case ss == sbest
469 cvg = 3; % there is no more suitable flag for this
476 aprec=abs(pprec.*pbest);
477 %% [aprec, chg, chgprev]
478 if (all(abs(chg) <= aprec) && all(abs(chgprev) <= aprec))
481 fprintf('Parameter changes converged to specified precision\n');
489 %% set further return values
495 function deb_printf (do_printf, varargin)
503 function fval = scalar_ifelse (cond, tval, fval)
505 %% needed for some anonymous functions, builtin ifelse only available
506 %% in Octave > 3.2; we need only the scalar case here