From: <i7...@us...> - 2010-05-20 06:00:20
|
Revision: 7323 http://octave.svn.sourceforge.net/octave/?rev=7323&view=rev Author: i7tiol Date: 2010-05-20 06:00:14 +0000 (Thu, 20 May 2010) Log Message: ----------- Deal with possible inaccuracy in bounds. - leasqr.m: Since projection algorithm might keep constraints with some inaccuracy, correct this for bounds, as long as bounds are the only constraints. - __dfdp__.m: Behave adequately if such an inaccuracy is present. Modified Paths: -------------- trunk/octave-forge/main/optim/inst/leasqr.m trunk/octave-forge/main/optim/inst/private/__dfdp__.m Modified: trunk/octave-forge/main/optim/inst/leasqr.m =================================================================== --- trunk/octave-forge/main/optim/inst/leasqr.m 2010-05-18 13:04:17 UTC (rev 7322) +++ trunk/octave-forge/main/optim/inst/leasqr.m 2010-05-20 06:00:14 UTC (rev 7323) @@ -85,16 +85,18 @@ %% since then some sanity tests are performed, and since the %% function 'dfdp.m' is guarantied not to violate constraints during %% determination of the numeric gradient only for those constraints - %% specified as 'bounds'. The first entry for general constraints - %% must be a differentiable vector valued function (say h), - %% specifying general inequality constraints of the form `h (p[, - %% idx]) >= 0'; p is the column vector of optimized paraters and the - %% optional argument idx is a logical index. h has to return the - %% values of all constraints if idx is not given, and has to return - %% only the indexed constraints if idx is given (so computation of - %% the other constraints can be spared). If a second entry for - %% general constraints is given, it must be a function (say dh) - %% which returnes a matrix whos rows contain the gradients of the + %% specified as 'bounds' (possibly with violations due to a certain + %% inaccuracy, however, except if no constraints except bounds are + %% specified). The first entry for general constraints must be a + %% differentiable vector valued function (say h), specifying general + %% inequality constraints of the form `h (p[, idx]) >= 0'; p is the + %% column vector of optimized paraters and the optional argument idx + %% is a logical index. h has to return the values of all constraints + %% if idx is not given, and has to return only the indexed + %% constraints if idx is given (so computation of the other + %% constraints can be spared). If a second entry for general + %% constraints is given, it must be a function (say dh) which + %% returnes a matrix whos rows contain the gradients of the %% constraint function h with respect to the optimized parameters. %% It has the form jac_h = dh (vh, p, dp, h, idx[, bounds]); p is %% the column vector of optimized parameters, and idx is a logical @@ -122,6 +124,8 @@ %% Field 'options.equc': equality constraints, specified the same %% way as inequality constraints (see field 'options.inequc'). %% Initial parameters must satisfy these constraints. + %% Note that there is possibly a certain inaccuracy in honoring + %% constraints, except if only bounds are specified. %% _Warning_: If constraints (or bounds) are set, returned guesses %% of corp, covp, and Z are generally invalid, even if no constraints %% are active for the final parameters. If equality constraints are @@ -425,6 +429,8 @@ error ('initial parameters violate equality constraints'); end end + have_constraints_except_bounds = ... + rv > 0 || erv > 0 || have_gencstr || have_genecstr; if (isfield (options, 'bounds')) bounds = options.bounds; if (rows (bounds) ~= n || columns (bounds) ~= 2) @@ -441,8 +447,7 @@ end lidx = pin < bounds(:, 1); uidx = pin > bounds(:, 2); - if (any (lidx | uidx) && (rv > 0 || have_gencstr || ... - erv > 0 || have_genecstr)) + if (any (lidx | uidx) && have_constraints_except_bounds) error ('initial parameters outside bounds, not corrected since other constraints are given'); end if (any (lidx)) @@ -808,6 +813,17 @@ % function if there is some non-miniscule % change p=chg+pprev; + %% since the projection method may have slightly violated + %% constraints due to inaccuracy, correct parameters to bounds + %% --- but only if no further constraints are given, otherwise + %% the inaccuracy in honoring them might increase by this + if (~have_constraints_except_bounds) + lidx = p < bounds(:, 1); + uidx = p > bounds(:, 2); + p(lidx) = bounds(lidx, 1); + p(uidx) = bounds(uidx, 2); + end + %% f = F (x, p); r = wt .* (y - f); r = r(:); Modified: trunk/octave-forge/main/optim/inst/private/__dfdp__.m =================================================================== --- trunk/octave-forge/main/optim/inst/private/__dfdp__.m 2010-05-18 13:04:17 UTC (rev 7322) +++ trunk/octave-forge/main/optim/inst/private/__dfdp__.m 2010-05-20 06:00:14 UTC (rev 7323) @@ -21,57 +21,73 @@ %% Meant to be called by interfaces 'dfdp.m' and 'dcdp.m', see there. m = length (f); - n=length(p); %dimensions - if (nargin < 6) + n = length (p); + if (nargin < 5) bounds = ones (n, 2); bounds(:, 1) = -Inf; bounds(:, 2) = Inf; end - prt=zeros(m,n); % initialise Jacobian to Zero - del = dp .* p; %cal delx=fract(dp)*param value(p) - idx = p == 0; - del(idx) = dp(idx); %if param=0 delx=fraction - idx = dp > 0; - del(idx) = abs (del(idx)); % not for one-sided intervals, changed - % direction of intervals could change - % behavior of optimization without bounds - min_del = min (abs (del), bounds(:, 2) - bounds(:, 1)); - for j=1:n - ps = p; - if (dp(j)~=0) - if (dp(j) < 0) - ps(j) = p(j) + del(j); - if (ps(j) < bounds(j, 1) || ps(j) > bounds(j, 2)) - t_del1 = max (bounds(j, 1) - p(j), - abs (del(j))); % - %non-positive - t_del2 = min (bounds(j, 2) - p(j), abs (del(j))); % - %non-negative - if (- t_del1 > t_del2) - del(j) = t_del1; - else - del(j) = t_del2; - end - ps(j) = p(j) + del(j); - end - tpf = func (ps); - prt(:, j) = (tpf(:) - f) / del(j); + prt = zeros (m, n); % initialise Jacobian to Zero + del = dp .* p; + absdel = abs (del); + idxa = p == 0; + del(idxa) = dp(idxa); + idxd = dp > 0; % double sided interval + idxz = dp == 0; + p1 = zeros (n, 1); + p2 = p1; + idxvs = false (n, 1); + idx1g2w = idxvs; + idx1le2w = idxvs; + + %% p may be slightly out of bounds due to inaccuracy, or exactly at + %% the bound -> single sided interval + idxvl = p <= bounds(:, 1); + idxvg = p >= bounds(:, 2); + p1(idxvl) = min (p(idxvl) + absdel(idxvl), bounds(idxvl, 2)); + idxd(idxvl) = false; + p1(idxvg) = max (p(idxvg) - absdel(idxvg), bounds(idxvg, 1)); + idxd(idxvg) = false; + idxs = ~(idxz | idxd); % single sided interval + + idxnv = ~(idxvl | idxvg); % current paramters within bounds + idxnvs = idxs & idxnv; % within bounds, single sided interval + idxnvd = idxd & idxnv; % within bounds, double sided interval + %% remaining single sided intervals + p1(idxnvs) = p(idxnvs) + del(idxnvs); % don't take absdel, this could + % change course of optimization without + % bounds with respect to previous + % versions + %% remaining single sided intervals, violating a bound -> take largest + %% possible direction of single sided interval + idxvs(idxnvs) = p1(idxnvs) < bounds(idxnvs, 1) | ... + p1(idxnvs) > bounds(idxnvs, 2); + del1 = p(idxvs) - bounds(idxvs, 1); + del2 = bounds(idxvs, 2) - p(idxvs); + idx1g2 = del1 > del2; + idx1g2w(idxvs) = idx1g2; + idx1le2w(idxvs) = ~idx1g2; + p1(idx1g2w) = max (p(idx1g2w) - absdel(idx1g2w), bounds(idx1g2w, 1)); + p1(idx1le2w) = min (p(idx1le2w) + absdel(idx1le2w), ... + bounds(idx1le2w, 2)); + %% double sided interval + p1(idxnvd) = min (p(idxnvd) + absdel(idxnvd), bounds(idxnvd, 2)); + p2(idxnvd) = max (p(idxnvd) - absdel(idxnvd), bounds(idxnvd, 1)); + + del(idxs) = p1(idxs) - p(idxs); + del(idxd) = p1(idxd) - p2(idxd); + + for j = 1:n + if (~idxz(j)) + ps = p; + ps(j) = p1(j); + tp1 = func (ps); + if (idxs(j)) + prt(:, j) = (tp1(:) - f) / del(j); else - if (p(j) - del(j) < bounds(j, 1)) - tp = bounds(j, 1); - ps(j) = tp + min_del(j); - elseif (p(j) + del(j) > bounds(j, 2)) - ps(j) = bounds(j, 2); - tp = ps(j) - min_del(j); - else - ps(j) = p(j) + del(j); - tp = p(j) - del(j); - min_del(j) = 2 * del(j); - end - f1 = func (ps); - f1 = f1(:); - ps(j) = tp; - tpf = func (ps); - prt(:, j) = (f1 - tpf(:)) / min_del(j); + ps(j) = p2(j); + tp2 = func (ps); + prt(:, j) = (tp1(:) - tp2(:)) / del(j); end end end This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |