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);