Mercurial > forge
changeset 11208:0c64c9856133 octave-forge
Fix some bugs in termination criteria, argument testing, and documentation.
author | i7tiol |
---|---|
date | Sun, 04 Nov 2012 18:42:25 +0000 |
parents | 2bcda689a109 |
children | 47f4637b3ee1 |
files | main/optim/DESCRIPTION main/optim/inst/nonlin_min.m main/optim/inst/nonlin_residmin.m main/optim/inst/private/__collect_constraints__.m main/optim/inst/private/__lm_feasible__.m main/optim/inst/private/__lm_svd__.m main/optim/inst/private/__nonlin_residmin__.m |
diffstat | 7 files changed, 149 insertions(+), 70 deletions(-) [+] |
line wrap: on
line diff
--- a/main/optim/DESCRIPTION Sun Nov 04 17:27:46 2012 +0000 +++ b/main/optim/DESCRIPTION Sun Nov 04 18:42:25 2012 +0000 @@ -1,6 +1,6 @@ Name: Optim Version: 1.2.1 -Date: 2012-10-01 +Date: 2012-11-04 Author: various authors Maintainer: Octave-Forge community <octave-dev@lists.sourceforge.net> Title: Optimization.
--- a/main/optim/inst/nonlin_min.m Sun Nov 04 17:27:46 2012 +0000 +++ b/main/optim/inst/nonlin_min.m Sun Nov 04 18:42:25 2012 +0000 @@ -129,30 +129,30 @@ ## two entries for linear constraints are a matrix (say @code{m}) and a ## vector (say @code{v}), specifying linear inequality constraints of ## the form @code{m.' * parameters + v >= 0}. The first entry for -## general constraints must be a differentiable vector valued function -## (say @code{h}), specifying general inequality constraints of the form -## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized -## paraters and the optional argument @code{idx} is a logical index. -## @code{h} has to return the values of all constraints if @code{idx} is -## not given. It may choose to return only the indexed constraints if -## @code{idx} is given (so computation of the other constraints can be -## spared); in this case, the additional setting @code{inequc_f_idx} has -## to be set to @code{true}. In gradient determination, this function -## may be called with an informational third argument, whose content -## depends on the function for gradient determination. If a second entry -## for general inequality constraints is given, it must be a function -## computing the jacobian of the constraints with respect to the -## parameters. For this function, the description of @code{dfdp} above -## applies, with 2 exceptions: 1) it is called with 3 arguments since it -## has an additional argument @code{idx} --- a logical index --- at -## second position, indicating which rows of the jacobian must be -## returned (if the function chooses to return only indexed rows, the -## additional setting @code{inequc_df_idx} has to be set to -## @code{true}). 2) the default jacobian function calls @code{h} with 3 -## arguments, since the argument @code{idx} is also supplied. Note that -## specifying linear constraints as general constraints will generally -## waste performance, even if further, non-linear, general constraints -## are also specified. +## general constraints must be a differentiable column-vector valued +## function (say @code{h}), specifying general inequality constraints of +## the form @code{h (p[, idx]) >= 0}; @code{p} is the column vector of +## optimized parameters and the optional argument @code{idx} is a +## logical index. @code{h} has to return the values of all constraints +## if @code{idx} is not given. It may choose to return only the indexed +## constraints if @code{idx} is given (so computation of the other +## constraints can be spared); in this case, the additional setting +## @code{inequc_f_idx} has to be set to @code{true}. In gradient +## determination, this function may be called with an informational +## third argument, whose content depends on the function for gradient +## determination. If a second entry for general inequality constraints +## is given, it must be a function computing the jacobian of the +## constraints with respect to the parameters. For this function, the +## description of @code{dfdp} above applies, with 2 exceptions: 1) it is +## called with 3 arguments since it has an additional argument +## @code{idx}, a logical index, at second position, indicating which +## rows of the jacobian must be returned (if the function chooses to +## return only indexed rows, the additional setting @code{inequc_df_idx} +## has to be set to @code{true}). 2) the default jacobian function calls +## @code{h} with 3 arguments, since the argument @code{idx} is also +## supplied. Note that specifying linear constraints as general +## constraints will generally waste performance, even if further, +## non-linear, general constraints are also specified. ## ## @code{equc}: Equality constraints. Specified the same way as ## inequality constraints (see @code{inequc}). The respective additional @@ -163,13 +163,13 @@ ## function is supplied with the package. ## ## @code{TolFun}: Minimum fractional improvement in objective function -## in an iteration (abortion criterium). Default: .0001. +## in an iteration (termination criterium). Default: .0001. ## -## @code{MaxIter}: Maximum number of iterations (abortion criterium). +## @code{MaxIter}: Maximum number of iterations (termination criterium). ## Default: backend-specific. ## ## @code{fract_prec}: Column Vector, minimum fractional change of -## parameters in an iteration (abortion criterium if violated in two +## parameters in an iteration (termination criterium if violated in two ## consecutive iterations). Default: backend-specific. ## ## @code{max_fract_change}: Column Vector, enforced maximum fractional @@ -994,6 +994,26 @@ warning ("some fixed parameters outside bounds"); endif + if (any (diffp <= 0)) + error ("some elements of 'diffp' non-positive"); + endif + + if (cstep <= 0) + error ("'cstep' non-positive"); + endif + + if ((hook.TolFun = optimget (settings, "TolFun", stol_default)) < 0) + error ("'TolFun' negative"); + endif + + if (any (fract_prec < 0)) + error ("some elements of 'fract_prec' negative"); + endif + + if (any (max_fract_change < 0)) + error ("some elements of 'max_fract_change' negative"); + endif + ## dimensions of linear constraints if (isempty (mc)) mc = zeros (np, 0); @@ -1033,7 +1053,6 @@ endif #### collect remaining settings - hook.TolFun = optimget (settings, "TolFun", stol_default); hook.MaxIter = optimget (settings, "MaxIter"); if (ischar (hook.cpiv = optimget (settings, "cpiv", @ cpiv_bard))) hook.cpiv = str2func (hook.cpiv);
--- a/main/optim/inst/nonlin_residmin.m Sun Nov 04 17:27:46 2012 +0000 +++ b/main/optim/inst/nonlin_residmin.m Sun Nov 04 18:42:25 2012 +0000 @@ -49,8 +49,8 @@ ## failure; its possible values depend on the used backend and currently ## can be @code{0} (maximum number of iterations exceeded), @code{2} ## (parameter change less than specified precision in two consecutive -## iterations), or @code{3} (improvement in objective function --- e.g. -## sum of squares --- less than specified). +## iterations), or @code{3} (improvement in objective function -- e.g. +## sum of squares -- less than specified). ## ## @var{settings}: ## @@ -124,30 +124,30 @@ ## two entries for linear constraints are a matrix (say @code{m}) and a ## vector (say @code{v}), specifying linear inequality constraints of ## the form @code{m.' * parameters + v >= 0}. The first entry for -## general constraints must be a differentiable vector valued function -## (say @code{h}), specifying general inequality constraints of the form -## @code{h (p[, idx]) >= 0}; @code{p} is the column vector of optimized -## paraters and the optional argument @code{idx} is a logical index. -## @code{h} has to return the values of all constraints if @code{idx} is -## not given. It may choose to return only the indexed constraints if -## @code{idx} is given (so computation of the other constraints can be -## spared); in this case, the additional setting @code{inequc_f_idx} has -## to be set to @code{true}. In gradient determination, this function -## may be called with an informational third argument, whose content -## depends on the function for gradient determination. If a second entry -## for general inequality constraints is given, it must be a function -## computing the jacobian of the constraints with respect to the -## parameters. For this function, the description of @code{dfdp} above -## applies, with 2 exceptions: 1) it is called with 3 arguments since it -## has an additional argument @code{idx} --- a logical index --- at -## second position, indicating which rows of the jacobian must be -## returned (if the function chooses to return only indexed rows, the -## additional setting @code{inequc_df_idx} has to be set to -## @code{true}). 2) the default jacobian function calls @code{h} with 3 -## arguments, since the argument @code{idx} is also supplied. Note that -## specifying linear constraints as general constraints will generally -## waste performance, even if further, non-linear, general constraints -## are also specified. +## general constraints must be a differentiable column-vector valued +## function (say @code{h}), specifying general inequality constraints of +## the form @code{h (p[, idx]) >= 0}; @code{p} is the column vector of +## optimized paraters and the optional argument @code{idx} is a logical +## index. @code{h} has to return the values of all constraints if +## @code{idx} is not given. It may choose to return only the indexed +## constraints if @code{idx} is given (so computation of the other +## constraints can be spared); in this case, the additional setting +## @code{inequc_f_idx} has to be set to @code{true}. In gradient +## determination, this function may be called with an informational +## third argument, whose content depends on the function for gradient +## determination. If a second entry for general inequality constraints +## is given, it must be a function computing the jacobian of the +## constraints with respect to the parameters. For this function, the +## description of @code{dfdp} above applies, with 2 exceptions: 1) it is +## called with 3 arguments since it has an additional argument +## @code{idx}, a logical index, at second position, indicating which +## rows of the jacobian must be returned (if the function chooses to +## return only indexed rows, the additional setting @code{inequc_df_idx} +## has to be set to @code{true}). 2) the default jacobian function calls +## @code{h} with 3 arguments, since the argument @code{idx} is also +## supplied. Note that specifying linear constraints as general +## constraints will generally waste performance, even if further, +## non-linear, general constraints are also specified. ## ## @code{equc}: Equality constraints. Specified the same way as ## inequality constraints (see @code{inequc}). @@ -160,14 +160,14 @@ ## match. ## ## @code{TolFun}: Minimum fractional improvement in objective function -## (e.g. sum of squares) in an iteration (abortion criterium). Default: +## (e.g. sum of squares) in an iteration (termination criterium). Default: ## .0001. ## -## @code{MaxIter}: Maximum number of iterations (abortion criterium). +## @code{MaxIter}: Maximum number of iterations (termination criterium). ## Default: backend-specific. ## ## @code{fract_prec}: Column Vector, minimum fractional change of -## parameters in an iteration (abortion criterium if violated in two +## parameters in an iteration (termination criterium if violated in two ## consecutive iterations). Default: backend-specific. ## ## @code{max_fract_change}: Column Vector, enforced maximum fractional
--- a/main/optim/inst/private/__collect_constraints__.m Sun Nov 04 17:27:46 2012 +0000 +++ b/main/optim/inst/private/__collect_constraints__.m Sun Nov 04 18:42:25 2012 +0000 @@ -21,6 +21,12 @@ if (isempty (cstr)) return; endif + for id = 1 : length (cstr) + if (ischar (cstr{id})) + cstr{id} = str2func (cstr{id}); + endif + endfor + if (ismatrix (tp = cstr{1}) || isstruct (tp)) mc = tp; vc = cstr{2};
--- a/main/optim/inst/private/__lm_feasible__.m Sun Nov 04 17:27:46 2012 +0000 +++ b/main/optim/inst/private/__lm_feasible__.m Sun Nov 04 18:42:25 2012 +0000 @@ -65,6 +65,12 @@ maxstep(isna (maxstep)) = max_fract_step_default; pprec = hook.fract_prec; pprec(isna (pprec)) = fract_prec_default; + ## keep absolute precision positive for non-null relative precision; + ## arbitrary value, added to parameters before multiplying with + ## relative precision + add_pprec = zeros (n, 1); + add_pprec(pprec > 0) = sqrt (eps); + ## verbose = strcmp (hook.Display, "iter"); ## some useful variables derived from passed variables @@ -123,7 +129,7 @@ ## gradient of objective function old_df = df; - df = grad_f (p, setfield (dfdp_hook, "f", f))(:); + df = grad_f (p, setfield (dfdp_hook, "f", vf))(:); ## constraints, preparation of some constants v_cstr = f_cstr (p, ac_idx); @@ -256,7 +262,7 @@ G = Binv * V; ## Levenberg/Marquardt - fgoal = (1 - ftol) * vf; + fgoal = vf - (abs (vf) + sqrt (eps)) * ftol; for l = ltab ll = max (ll, nminL); @@ -446,10 +452,11 @@ endif endif # regaining - aprec = abs (pprec .* pbest); + aprec = pprec .* (abs (pbest) + add_pprec); if (any (abs (chg) > 0.1 * aprec)) # only worth evaluating # function if there is some # non-miniscule change + skipped = false; p_chg = p + chg; ## since the projection method may have slightly violated ## constraints due to inaccuracy, correct parameters to bounds @@ -471,18 +478,26 @@ pbest = p_chg; fbest = vf_chg; endif - if (vf_chg <= fgoal) + if (vf_chg < fgoal) # <, not <=, since fgoal can be equal to vf + # if TolFun <= eps p = p_chg; vf = vf_chg; break; endif + else + skipped = true; + break; endif endfor ll = l; - aprec = abs (pprec .* pbest); - if (vf_chg < eps || vf_chg > fgoal) + aprec = pprec .* (abs (pbest) + add_pprec); + if (skipped) + cvg = 2; + done = true; + elseif (vf_chg >= fgoal) # >=, not >, since fgoal can be equal to vf + # if TolFun <= eps cvg = 3; done = true; elseif (all (abs (chg) <= aprec) && all (abs (chgprev) <= aprec))
--- a/main/optim/inst/private/__lm_svd__.m Sun Nov 04 17:27:46 2012 +0000 +++ b/main/optim/inst/private/__lm_svd__.m Sun Nov 04 18:42:25 2012 +0000 @@ -30,6 +30,8 @@ ifelse = @ scalar_ifelse; end + n = length (pin); + %% passed constraints mc = hook.mc; % matrix of linear constraints vc = hook.vc; % vector of linear constraints @@ -58,6 +60,12 @@ maxstep(isna (maxstep)) = max_fract_step_default; pprec = hook.fract_prec; pprec(isna (pprec)) = fract_prec_default; + ## keep absolute precision positive for non-null relative precision; + ## arbitrary value, added to parameters before multiplying with + ## relative precision + add_pprec = zeros (n, 1); + add_pprec(pprec > 0) = sqrt (eps); + ## stol = hook.TolFun; niter = hook.MaxIter; if (isempty (niter)) niter = 20; end @@ -82,7 +90,6 @@ have_constraints_except_bounds = ... n_lcstr + n_gencstr > ... sum (lbound ~= -Inf) + sum (ubound ~= Inf); - n = length (pin); wtl = wt(:); nz = 20 * eps; % This is arbitrary. Constraint function will be @@ -415,7 +422,7 @@ sprintf('%d ',find(ochg ~= chg)), 'maximal fractional stepwidth enforced']); end end - aprec=abs(pprec.*pbest); %--- + aprec = pprec .* (abs (pbest) + add_pprec); %% ss=scalar sum of squares=sum((wt.*f)^2). if (any(abs(chg) > 0.1*aprec))%--- % only worth evaluating % function if there is some non-miniscule @@ -431,6 +438,7 @@ %% 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 + skipped = false; if (~have_constraints_except_bounds) lidx = p < lbound; uidx = p > ubound; @@ -455,9 +463,13 @@ fbest=f; sbest=ss; end - if (ss<=sgoal) + if (ss < sgoal) # <, not <=, since sgoal can be equal to sprev + # if TolFun <= eps break; end + else + skipped = true; + break; end %--- end %% printf ('epsL no.: %i\n', jjj); % for testing @@ -465,15 +477,20 @@ if (verbose) hook.plot_cmd (f); end + if (skipped) + cvg = 2; + break; + end if (ss < eps) % in this case ss == sbest cvg = 3; % there is no more suitable flag for this break; end - if (ss>sgoal) + if (ss >= sgoal) # >=, not >, since sgoal can be equal to sprev if + # TolFun <= eps cvg = 3; break; end - aprec=abs(pprec.*pbest); + aprec = pprec .* (abs (pbest) + add_pprec); %% [aprec, chg, chgprev] if (all(abs(chg) <= aprec) && all(abs(chgprev) <= aprec)) cvg = 2;
--- a/main/optim/inst/private/__nonlin_residmin__.m Sun Nov 04 17:27:46 2012 +0000 +++ b/main/optim/inst/private/__nonlin_residmin__.m Sun Nov 04 18:42:25 2012 +0000 @@ -658,6 +658,26 @@ warning ("some fixed parameters outside bounds"); endif + if (any (diffp <= 0)) + error ("some elements of 'diffp' non-positive"); + endif + + if (cstep <= 0) + error ("'cstep' non-positive"); + endif + + if ((hook.TolFun = optimget (settings, "TolFun", stol_default)) < 0) + error ("'TolFun' negative"); + endif + + if (any (fract_prec < 0)) + error ("some elements of 'fract_prec' negative"); + endif + + if (any (max_fract_change < 0)) + error ("some elements of 'max_fract_change' negative"); + endif + ## dimensions of linear constraints if (isempty (mc)) mc = zeros (np, 0); @@ -683,6 +703,9 @@ if (any (size (weights) != size (f_pin))) error ("dimension of weights and residuals must match"); endif + if (any (weights(:) < 0)) + error ("some weights negative") + endif ## note initial values of linear constraits pin_cstr.inequ.lin_except_bounds = mc.' * pin + vc; @@ -703,7 +726,6 @@ endif #### collect remaining settings - hook.TolFun = optimget (settings, "TolFun", stol_default); hook.MaxIter = optimget (settings, "MaxIter"); if (ischar (hook.cpiv = optimget (settings, "cpiv", @ cpiv_bard))) hook.cpiv = str2func (hook.cpiv);