Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[10622] trunk/octave-forge/main/optim
|
From: <i7...@us...> - 2012-06-12 14:09:04
|
Revision: 10622
http://octave.svn.sourceforge.net/octave/?rev=10622&view=rev
Author: i7tiol
Date: 2012-06-12 14:08:57 +0000 (Tue, 12 Jun 2012)
Log Message:
-----------
Completed optimization suite with new backend lm_feasible.
Modified Paths:
--------------
trunk/octave-forge/main/optim/DESCRIPTION
trunk/octave-forge/main/optim/NEWS
trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
trunk/octave-forge/main/optim/inst/nonlin_min.m
trunk/octave-forge/main/optim/inst/nonlin_residmin.m
trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
Added Paths:
-----------
trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m
Modified: trunk/octave-forge/main/optim/DESCRIPTION
===================================================================
--- trunk/octave-forge/main/optim/DESCRIPTION 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/DESCRIPTION 2012-06-12 14:08:57 UTC (rev 10622)
@@ -1,6 +1,6 @@
Name: Optim
-Version: 1.1.1
-Date: 2012-06-07
+Version: 1.2.0
+Date: 2012-06-12
Author: various authors
Maintainer: Octave-Forge community <oct...@li...>
Title: Optimization.
Modified: trunk/octave-forge/main/optim/NEWS
===================================================================
--- trunk/octave-forge/main/optim/NEWS 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/NEWS 2012-06-12 14:08:57 UTC (rev 10622)
@@ -1,3 +1,26 @@
+Summary of important user-visible changes for optim 1.2.0:
+-------------------------------------------------------------------
+
+ ** Together with the new backend "lm_feasible" there is now a
+ complete suite of backends for optimization with linear and
+ general equality and inequality constraints, for scalar valued
+ objective functions and for array valued model function, which
+ features, a.o., honouring of constraints throughout optimization
+ and handling of structure-based parameters. The respective user
+ functions (frontends) are
+
+ nonlin_min nonlin_residmin nonlin_curvefit
+
+ together with a user function for statistics
+
+ residmin_stat
+
+ ** The requirement of nonlin_min, nonlin_residmin, and
+ nonlin_curvefit for the general constraint functions being able to
+ honour an index of constraints has been removed, the respective
+ feature is still available by setting some options.
+
+
Summary of important user-visible changes for optim 1.1.0:
-------------------------------------------------------------------
Modified: trunk/octave-forge/main/optim/inst/nonlin_curvefit.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_curvefit.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -81,3 +81,21 @@
endif
endfunction
+
+%!demo
+%! ## Example for linear inequality constraints
+%! ## (see also the same example in 'demo nonlin_residmin')
+%!
+%! ## independents and observations
+%! indep = 1:5;
+%! obs = [1, 2, 4, 7, 14];
+%! ## model function:
+%! f = @ (p, x) p(1) * exp (p(2) * x);
+%! ## initial values:
+%! init = [.25; .25];
+%! ## linear constraints, A.' * parametervector + B >= 0
+%! A = [1; -1]; B = 0; # p(1) >= p(2);
+%! settings = optimset ("inequc", {A, B});
+%!
+%! ## start optimization
+%! [p, model_values, cvg, outp] = nonlin_curvefit (f, init, indep, obs, settings)
Modified: trunk/octave-forge/main/optim/inst/nonlin_min.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_min.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -22,10 +22,6 @@
## interface; any additionally needed constants can be supplied by
## wrapping the user functions into anonymous functions.
##
-## At the moment there is only a stochastic backend implemented, so that
-## some settings in the frontend have no use. A backend for
-## feasible-path sequential quadratic programming is being worked on.
-##
## The following description applies to usage with vector-based
## parameter handling. Differences in usage for structure-based
## parameter handling will be explained in a separate section below.
@@ -52,14 +48,14 @@
## currently can be @code{0} (maximum number of iterations exceeded),
## @code{1} (fixed number of iterations completed, e.g. in stochastic
## optimizers), @code{2} (parameter change less than specified precision
-## in two consecutive iterations), or @code{3} (improvement in objective
-## function less than specified).
+## in two consecutive iterations), @code{3} (improvement in objective
+## function less than specified), or @code{-4} (algorithm got stuck).
##
## @var{settings}:
##
-## @code{Algorithm}: String specifying the backend. Currently there is
-## no default, since only a stochastic backend is implemented
-## (@code{"siman"}). It is described in a separate section below.
+## @code{Algorithm}: String specifying the backend. Currently available
+## are @{"lm_feasible"} (default) and @code{"siman"}. They are described
+## in separate sections below.
##
## @code{objf_grad}: Function computing the gradient of the objective
## function with respect to the parameters, assuming residuals are
@@ -138,24 +134,29 @@
## @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, and has to return only the indexed constraints if
+## 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 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, except that
-## 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, and except that the
-## default 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.
+## 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}).
+## inequality constraints (see @code{inequc}). The respective additional
+## settings are named @code{equc_f_idx} and @code{equc_df_idx}.
##
## @code{cpiv}: Function for complementary pivoting, usable in
## algorithms for constraints. Default: @ cpiv_bard. Only the default
@@ -242,8 +243,30 @@
## parameters must contain arrays with dimensions reshapable to those of
## the respective parameters.
##
-## Description of backends (currently only one)
+## Description of backends
##
+## "lm_feasible"
+##
+## A Levenberg/Marquardt-like optimizer, attempting to honour
+## constraints throughout the course of optimization. This means that
+## the initial parameters must not violate constraints (to find an
+## initial feasible set of parameters, e.g. Octaves @code{sqp} can be
+## used, by specifying an objective function which is constant or which
+## returns the quadratic distance to the initial values). If the
+## constraints need only be honoured in the result of the optimization,
+## Octaves @code{sqp} may be preferable. The Hessian is either supplied
+## by the user or is approximated by the BFGS algorithm.
+##
+## Returned value @var{cvg} will be @code{2} or @code{3} for success and
+## @code{0} or @code{-4} for failure (see above for meaning).
+##
+## Backend-specific defaults are: @code{MaxIter}: 20, @code{fract_prec}:
+## @code{zeros (size (parameters))}, @code{max_fract_change}: @code{Inf}
+## for all parameters.
+##
+## Interpretation of @code{Display}: if set to @code{"iter"}, currently
+## only information on applying @code{max_fract_change} is printed.
+##
## "siman"
##
## A simulated annealing (stochastic) optimizer, changing all parameters
@@ -335,10 +358,14 @@
"fixed", [], \
"inequc", [], \
"equc", [], \
+ "inequc_f_idx", false, \
+ "inequc_df_idx", false, \
+ "equc_f_idx", false, \
+ "equc_df_idx", false, \
"TolFun", stol_default, \
"MaxIter", [], \
"Display", "off", \
- "Algorithm", [], \ # set reasonable default, once available
+ "Algorithm", "lm_feasible", \
"T_init", .01, \
"T_min", 1.0e-5, \
"mu_T", 1.005, \
@@ -351,10 +378,12 @@
return;
endif
- assign = @ assign; # Is this faster in repeated calls?
+ if (nargin < 2 || nargin > 3)
+ print_usage ();
+ endif
- if (nargin != 3)
- error ("incorrect number of arguments");
+ if (nargin == 2)
+ settings = struct ();
endif
if (ischar (f))
@@ -418,13 +447,36 @@
mc_struct = isstruct (mc);
emc_struct = isstruct (emc);
- ## correct "_pstruct" settings if functions are not supplied
+ ## correct "_pstruct" settings if functions are not supplied, handle
+ ## constraint functions not honoring indices
if (isempty (dfdp)) dfdp_pstruct = false; endif
if (isempty (hessian)) hessian_pstruct = false; endif
- if (isempty (f_genicstr)) f_inequc_pstruct = false; endif
- if (isempty (f_genecstr)) f_equc_pstruct = false; endif
- if (! user_df_gencstr) df_inequc_pstruct = false; endif
- if (! user_df_genecstr) df_equc_pstruct = false; endif
+ if (isempty (f_genicstr))
+ f_inequc_pstruct = false;
+ elseif (! optimget (settings, "inequc_f_idx", false))
+ f_genicstr = @ (p, varargin) apply_idx_if_given \
+ (f_genicstr (p, varargin{:}), varargin{:});
+ endif
+ if (isempty (f_genecstr))
+ f_equc_pstruct = false;
+ elseif (! optimget (settings, "equc_f_idx", false))
+ f_genecstr = @ (p, varargin) apply_idx_if_given \
+ (f_genecstr (p, varargin{:}), varargin{:});
+ endif
+ if (user_df_gencstr)
+ if (! optimget (settings, "inequc_df_idx", false))
+ df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_inequc_pstruct = false;
+ endif
+ if (user_df_genecstr)
+ if (! optimget (settings, "equc_df_idx", false))
+ df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_equc_pstruct = false;
+ endif
## some settings require a parameter order
if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
@@ -988,7 +1040,6 @@
endif
hook.Display = optimget (settings, "Display", "off");
hook.testing = optimget (settings, "debug", false);
- hook.new_s = optimget (settings, "lm_svd_feasible_alt_s", false);
hook.siman.T_init = optimget (settings, "T_init", .01);
hook.siman.T_min = optimget (settings, "T_min", 1.0e-5);
hook.siman.mu_T = optimget (settings, "mu_T", 1.005);
@@ -999,10 +1050,7 @@
optimget (settings, "trace_steps", false);
hook.siman_log = \
optimget (settings, "siman_log", false);
- backend = optimget (settings, "Algorithm");
- if (isempty (backend))
- error ("At the moment, there is no default algorithm, it must be specified.");
- endif
+ backend = optimget (settings, "Algorithm", "lm_feasible");
backend = map_matlab_algorithm_names (backend);
backend = map_backend (backend);
@@ -1246,7 +1294,7 @@
## passed function for complementary pivoting
## hook.cpiv = cpiv; # set before
- ## passed value of residual function for initial parameters
+ ## passed value of objective function for initial parameters
hook.f_pin = f_pin;
## passed options
@@ -1283,12 +1331,12 @@
function backend = map_backend (backend)
switch (backend)
- case "sqp_infeasible"
- backend = "__sqp__";
- case "sqp"
- backend = "__sqp__";
- case "sqp_feasible"
- backend = "__sqp_feasible__";
+ ## case "sqp_infeasible"
+ ## backend = "__sqp__";
+ ## case "sqp"
+ ## backend = "__sqp__";
+ case "lm_feasible"
+ backend = "__lm_feasible__";
case "siman"
backend = "__siman__";
otherwise
@@ -1336,7 +1384,24 @@
endfunction
+function ret = apply_idx_if_given (ret, varargin)
+
+ if (nargin > 1)
+ ret = ret(varargin{1});
+ endif
+
+endfunction
+
%!demo
+%! ## Example for default optimization (Levenberg/Marquardt with
+%! ## BFGS), one non-linear equality constraint. Constrained optimum is
+%! ## at p = [0; 1].
+%! objective_function = @ (p) p(1)^2 + p(2)^2;
+%! pin = [-2; 5];
+%! constraint_function = @ (p) p(1)^2 + 1 - p(2);
+%! [p, objf, cvg, outp] = nonlin_min (objective_function, pin, optimset ("equc", {constraint_function}))
+
+%!demo
%! ## Example for simulated annealing, two parameters, "trace_steps"
%! ## is true;
%! t_init = .2;
Modified: trunk/octave-forge/main/optim/inst/nonlin_residmin.m
===================================================================
--- trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/nonlin_residmin.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -129,21 +129,25 @@
## @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, and has to return only the indexed constraints if
+## 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 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, except that
-## 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, and except that the
-## default 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.
+## 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}).
@@ -281,3 +285,20 @@
[p, resid, cvg, outp] = __nonlin_residmin__ (varargin{:});
endfunction
+
+%!demo
+%! ## Example for linear inequality constraints
+%! ## (see also the same example in 'demo nonlin_curvefit')
+%!
+%! ## independents
+%! indep = 1:5;
+%! ## residual function:
+%! f = @ (p) p(1) * exp (p(2) * indep) - [1, 2, 4, 7, 14];
+%! ## initial values:
+%! init = [.25; .25];
+%! ## linear constraints, A.' * parametervector + B >= 0
+%! A = [1; -1]; B = 0; # p(1) >= p(2);
+%! settings = optimset ("inequc", {A, B});
+%!
+%! ## start optimization
+%! [p, residuals, cvg, outp] = nonlin_residmin (f, init, settings)
Added: trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m (rev 0)
+++ trunk/octave-forge/main/optim/inst/private/__lm_feasible__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -0,0 +1,505 @@
+## Copyright (C) 2012 Olaf Till <i7...@t-...>
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+function [p_res, objf, cvg, outp] = __lm_feasible__ (f, pin, hook)
+
+ ## some backend specific defaults
+ fract_prec_default = 0;
+ max_fract_step_default = Inf;
+
+ ## needed for some anonymous functions
+ if (exist ("ifelse") != 5)
+ ifelse = @ scalar_ifelse;
+ endif
+
+ n = length (pin);
+
+ ## passed constraints
+ mc = hook.mc; # matrix of linear constraints
+ vc = hook.vc; # vector of linear constraints
+ f_cstr = hook.f_cstr; # function of all constraints
+ df_cstr = hook.df_cstr; # function of derivatives of all constraints
+ n_gencstr = hook.n_gencstr; # number of non-linear constraints
+ eq_idx = hook.eq_idx; # logical index of equality constraints in all
+ # constraints
+ lbound = hook.lbound; # bounds, subset of linear inequality
+ ubound = hook.ubound; # constraints in mc and vc
+
+ ## passed values of constraints for initial parameters
+ pin_cstr = hook.pin_cstr;
+
+ ## passed return value of f for initial parameters
+ f_pin = hook.f_pin;
+
+ ## passed function for gradient of objective function
+ grad_f = hook.dfdp;
+
+ ## passed function for hessian of objective function
+ if (isempty (hessian = hook.hessian))
+ user_hessian = false;
+ A = eye (n);
+ else
+ user_hessian = true;
+ endif
+
+ ## passed function for complementary pivoting
+ cpiv = hook.cpiv;
+
+ ## passed options
+ ftol = hook.TolFun;
+ if (isempty (niter = hook.MaxIter)) niter = 20; endif
+ fixed = hook.fixed;
+ maxstep = hook.max_fract_change;
+ maxstep(isna (maxstep)) = max_fract_step_default;
+ pprec = hook.fract_prec;
+ pprec(isna (pprec)) = fract_prec_default;
+ verbose = strcmp (hook.Display, "iter");
+
+ ## some useful variables derived from passed variables
+ n_lcstr = size (vc, 1);
+ have_constraints_except_bounds = \
+ n_lcstr + n_gencstr > \
+ sum (lbound != -Inf) + sum (ubound != Inf);
+ ac_idx = true (n_lcstr + n_gencstr, 1); # index of all constraints
+ nc_idx = false (n_lcstr + n_gencstr, 1); # none of all constraints
+ gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); # gen. constr.
+
+ nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
+
+ ## backend-specific checking of options and constraints
+ ##
+ if (any (pin < lbound | pin > ubound) ||
+ any (pin_cstr.inequ.lin_except_bounds < 0) ||
+ any (pin_cstr.inequ.gen < 0) ||
+ any (abs (pin_cstr.equ.lin)) >= nz ||
+ any (abs (pin_cstr.equ.gen)) >= nz)
+ error ("Initial parameters violate constraints.");
+ endif
+ ##
+ idx = lbound == ubound;
+ if (any (idx))
+ warning ("lower and upper bounds identical for some parameters, fixing the respective parameters");
+ fixed(idx) = true;
+ endif
+ if (all (fixed))
+ error ("no free parameters");
+ endif
+ if (n_gencstr > 0 && any (! isinf (maxstep)))
+ warning ("setting both a maximum fractional step change of parameters and general constraints may result in inefficiency and failure");
+ endif
+
+ ## fill constant fields of hook for derivative-functions; some fields
+ ## may be backend-specific
+ dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
+ # the backend still may add to it
+
+ ## set up for iterations
+ p = pbest = pin;
+ vf = fbest = f_pin;
+ iter = 0;
+ done = false;
+ ll = 1;
+ ltab = [.1, 1, 1e2, 1e4, 1e6];
+ chgprev = Inf (n, 1);
+ df = [];
+ c_act = false (n, 1);
+ dca = zeros (n, 0);
+
+ while (! done)
+
+ iter++;
+
+ ## gradient of objective function
+ old_df = df;
+ df = grad_f (p, setfield (dfdp_hook, "f", f))(:);
+
+ ## constraints, preparation of some constants
+ v_cstr = f_cstr (p, ac_idx);
+ old_c_act = c_act;
+ old_dca = dca;
+ c_act = v_cstr < nz | eq_idx; # index of active constraints
+ if (any (c_act))
+
+ if (n_gencstr)
+ ## full gradient is needed later
+ dct = df_cstr (p, ac_idx, setfield (dfdp_hook, "f", v_cstr));
+ dct(:, fixed) = 0; # for user supplied dfdp; necessary?
+ dcat = dct(c_act, :);
+ else
+ dcat = df_cstr (p, c_act, setfield (dfdp_hook, "f", v_cstr));
+ dcat(:, fixed) = 0; # for user supplied dfdp; necessary?
+ endif
+
+ dca = dcat.';
+
+ a_eq_idx = eq_idx(c_act);
+
+ else
+
+ dca = zeros (n, 0);
+
+ endif
+
+ ## hessian of objectiv function
+ if (user_hessian)
+
+ A = hessian (p);
+ idx = isnan (A);
+ A(idx) = A.'(idx);
+ if (any (isnan (A(:))))
+ error ("some second derivatives undefined by user function");
+ endif
+ if (! isreal (A))
+ error ("second derivatives given by user function not real");
+ endif
+ if (! issymmetric (A))
+ error ("Hessian returned by user function not symmetric");
+ endif
+
+ elseif (iter > 1)
+
+ if (any (chg))
+
+ ## approximate Hessian of Lagrangian
+
+ ## I wonder if this hassle here and above with accounting for
+ ## changing active sets is indeed better than just approximating
+ ## the Hessian only of the objective function.
+ ##
+ ## index, over all constraints, of constraints active both
+ ## previously and currently
+ s_c_act = old_c_act & c_act;
+ ## index, over currently active constraints, of constraints
+ ## active both previously and currently
+ id_new = s_c_act(c_act);
+ ## index, over previously active constraints, of constraints
+ ## active both previously and currently
+ id_old = s_c_act(old_c_act);
+ ## gradients of currently active constraints which were also
+ ## active previously
+ dca_new_id = dca(:, id_new);
+ ## gradients of previously active constraints which are also
+ ## active currently
+ dca_old_id = old_dca(:, id_old);
+ ## index, over constraints active both previously and currently,
+ ## of (old) non-zero multipliers (bidx set below previously)
+ bidx_old_id = bidx(id_old);
+ ## index, over (old) non-zero multipliers, of constraints active
+ ## both previously and currently (bidx set below previously)
+ old_l_idx = id_old(bidx);
+
+ ## difference of derivatives of new and old active constraints,
+ ## multiplied by multipliers, as used for BFGS update (lb set
+ ## below previously)
+ dch = (dca_new_id(:, bidx_old_id) - \
+ dca_old_id(:, bidx_old_id)) * \
+ lb(old_l_idx);
+
+ y = df - old_df - dch;
+
+ ## Damped BFGS according to Nocedal & Wright, 2nd edition,
+ ## procedure 18.2.
+ chgt = chg.';
+ sAs = chgt * A * chg;
+ cy = chgt * y;
+ if (cy >= .2 * sAs)
+ th = 1;
+ else
+ if ((den1 = sAs - cy) == 0)
+ cvg = -4;
+ break;
+ endif
+ th = .8 * sAs / den1;
+ endif
+ Ac = A * chg;
+ r = th * y + (1 - th) * Ac;
+
+ if ((den2 = chgt * r) == 0 || sAs == 0)
+ cvg = -4;
+ break;
+ endif
+ A += r * r.' / den2 - Ac * Ac.' / sAs;
+
+ endif
+
+ endif
+
+ ## Inverse scaled decomposition A = G * (1 ./ L) * G.'
+ ##
+ ## make matrix Binv for scaling
+ Binv = diag (A);
+ nidx = ! (idx = Binv == 0);
+ Binv(nidx) = 1 ./ sqrt (abs (Binv(nidx)));
+ Binv(idx) = 1;
+ Binv = diag (Binv);
+ ## eigendecomposition of scaled A
+ [V, L] = eig (Binv * A * Binv);
+ L = diag (L);
+ ## A is symmetric, so V and L are real, delete any imaginary parts,
+ ## which might occur due to inaccuracy
+ V = real (V);
+ L = real (L);
+ ##
+ nminL = - min (L) * 1.1 / ltab(1);
+ G = Binv * V;
+
+ ## Levenberg/Marquardt
+ fgoal = (1 - ftol) * vf;
+ for l = ltab
+
+ ll = max (ll, nminL);
+ l = max (1e-7, ll * l);
+
+ R = G * diag (1 ./ (L + l)) * G.';
+
+ ## step computation
+ if (any (c_act))
+
+ ## some constraints are active, quadratic programming
+
+ tp = dcat * R;
+ [lb, bidx, ridx, tbl] = cpiv (- tp * df, tp * dca, a_eq_idx);
+ chg = R * (dca(:, bidx) * lb - df); # step direction
+
+ ## indices for different types of constraints
+ c_inact = ! c_act; # inactive constraints
+ c_binding = c_unbinding = nc_idx;
+ c_binding(c_act) = bidx; # constraints selected binding
+ c_unbinding(c_act) = ridx; # constraints unselected binding
+ c_nonbinding = c_act & ! (c_binding | c_unbinding); #
+ #constraints selected non-binding
+
+ else
+
+ ## no constraints are active, chg is the Levenberg/Marquardt step
+
+ chg = - R * df; # step direction
+
+ lb = zeros (0, 1);
+ bidx = false (0, 1);
+
+ ## indices for different types of constraints (meaning see above)
+ c_inact = ac_idx;
+ c_binding = nc_idx;
+ c_unbinding = nc_idx;
+ c_nonbinding = nc_idx;
+
+ endif
+
+ ## apply inactive and non-binding constraints to step width
+ ##
+ ## linear constraints
+ k = 1;
+ c_tp = c_inact(1:n_lcstr);
+ mcit = mc(:, c_tp).';
+ vci = vc(c_tp);
+ hstep = mcit * chg;
+ idx = hstep < 0;
+ if (any (idx))
+ k = min (1, min (- (vci(idx) + mcit(idx, :) * p) ./ \
+ hstep(idx)));
+ endif
+ ##
+ ## general constraints
+ if (n_gencstr)
+ c_tp = gc_idx & (c_nonbinding | c_inact);
+ if (any (c_tp) && any (f_cstr (p + k * chg, c_tp) < 0))
+ [k, fval, info] = \
+ fzero (@ (x) min (cat (1, \
+ f_cstr (p + x * chg, c_tp), \
+ k - x, \
+ ifelse (x < 0, -Inf, Inf))), \
+ 0);
+ if (info != 1 || abs (fval) >= nz)
+ error ("could not find stepwidth to satisfy inactive and non-binding general inequality constraints");
+ endif
+ endif
+ endif
+ ##
+ chg = k * chg;
+
+ ## if necessary, regain binding constraints and one of the
+ ## possibly active previously inactive or non-binding constraints
+ if (any (gc_idx & c_binding)) # none selected binding => none
+ # unselected binding
+ ptp1 = p + chg;
+
+ tp = true;
+ nt_nosuc = true;
+ lim = 20;
+ while (nt_nosuc && lim >= 0)
+ ## we keep d_p.' * inv (R) * d_p minimal in each step of the
+ ## inner loop
+ c_tp0 = c_inact | c_nonbinding;
+ c_tp1 = c_inact | (gc_idx & c_nonbinding);
+ btbl = tbl(bidx, bidx);
+ c_tp2 = c_binding;
+ ## once (any(tp)==false), it would not get true again even
+ ## with the following assignment
+ if (any (tp) && \
+ any (tp = f_cstr (ptp1, c_tp1) < nz))
+ ## keep only the first true entry in tp
+ tp(tp) = logical (cat (1, 1, zeros (sum (tp) - 1, 1)));
+ ## supplement binding index with one (the first) getting
+ ## binding in c_tp1
+ c_tp2(c_tp1) = tp;
+ ## gradient of this added constraint
+ caddt = dct(c_tp2 & ! c_binding, :);
+ cadd = caddt.';
+ C = dct(c_binding, :) * R * cadd;
+ Ct = C.';
+ T = [btbl, btbl * C; \
+ -Ct * btbl, caddt * R * cadd - Ct * btbl * C];
+ btbl = gjp (T, size (T, 1));
+ endif
+ dcbt = dct(c_tp2, :);
+ mfc = - R * dcbt.' * btbl;
+
+ ptp2 = ptp1;
+ nt_niter = nt_niter_start = 100;
+ while (nt_nosuc && nt_niter >= 0)
+ hv = f_cstr (ptp2, c_tp2);
+ if (all (abs (hv) < nz))
+ nt_nosuc = false;
+ chg = ptp2 - p;
+ else
+ ptp2 = ptp2 + mfc * hv; # step should be zero for each
+ # component for which the parameter is
+ # "fixed"
+ endif
+ nt_niter--;
+ endwhile
+
+ if (nt_nosuc || \
+ any (abs (chg) > abs (p .* maxstep)) || \
+ any (f_cstr (ptp2, c_tp0) < -nz))
+ ## if (nt_nosuc), regaining did not converge, else,
+ ## regaining violated type 3 and 4.
+ nt_nosuc = true;
+ ptp1 = (p + ptp1) / 2;
+ endif
+ if (! nt_nosuc && \
+ any ((tp = f_cstr (ptp2, c_unbinding)) < 0))
+ [discarded, id] = min(tp);
+ tid = find (ridx);
+ id = tid(id); # index within active constraints
+ unsuccessful_exchange = false;
+ if (abs (tbl(id, id)) < nz) # Bard: not absolute value
+ ## exchange this unselected binding constraint against a
+ ## binding constraint, but not against an equality
+ ## constraint
+ tbidx = bidx & ! a_eq_idx;
+ if (! any (tbidx))
+ unsuccessful_exchange = true;
+ else
+ [discarded, idm] = max (abs (tbl(tbidx, id)));
+ tid = find (tbidx);
+ idm = tid(idm); # -> index within active constraints
+ tbl = gjp (tbl, idm);
+ bidx(idm) = false;
+ ridx(idm) = true;
+ endif
+ endif
+ if (unsuccessful_exchange)
+ ## It probably doesn't look good now; this desperate last
+ ## attempt is not in the original algortithm, since that
+ ## didn't account for equality constraints.
+ ptp1 = (p + ptp1) / 2;
+ else
+ tbl = gjp (tbl, id);
+ bidx(id) = true;
+ ridx(id) = false;
+ c_binding = nc_idx;
+ c_binding(c_act) = bidx;
+ c_unbinding = nc_idx;
+ c_unbinding(c_act) = ridx;
+ endif
+ ## regaining violated type 2 constraints
+ nt_nosuc = true;
+ endif
+ lim--;
+ endwhile
+ if (nt_nosuc)
+ error ("could not regain binding constraints");
+ endif
+ else
+ ## check the maximal stepwidth and apply as necessary
+ ochg = chg;
+ idx = ! isinf (maxstep);
+ limit = abs (maxstep(idx) .* p(idx));
+ chg(idx) = min (max (chg(idx), - limit), limit);
+ if (verbose && any (ochg != chg))
+ printf ("Change in parameter(s): %s:maximal fractional stepwidth enforced", \
+ sprintf ("%d ", find (ochg != chg)));
+ endif
+ endif # regaining
+
+ aprec = abs (pprec .* pbest);
+ if (any (abs (chg) > 0.1 * aprec)) # only worth evaluating
+ # function if there is some
+ # non-miniscule change
+ p_chg = p + chg;
+ ## 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_chg < lbound;
+ uidx = p_chg > ubound;
+ p_chg(lidx, 1) = lbound(lidx, 1);
+ p_chg(uidx, 1) = ubound(uidx, 1);
+ chg(lidx, 1) = p_chg(lidx, 1) - p(lidx, 1);
+ chg(uidx, 1) = p_chg(uidx, 1) - p(uidx, 1);
+ endif
+ ##
+ if (! isreal (vf_chg = f (p_chg)))
+ error ("objective function not real");
+ endif
+ if (vf_chg < fbest)
+ pbest = p_chg;
+ fbest = vf_chg;
+ endif
+ if (vf_chg <= fgoal)
+ p = p_chg;
+ vf = vf_chg;
+ break;
+ endif
+ endif
+ endfor
+
+ ll = l;
+
+ aprec = abs (pprec .* pbest);
+ if (vf_chg < eps || vf_chg > fgoal)
+ cvg = 3;
+ done = true;
+ elseif (all (abs (chg) <= aprec) && all (abs (chgprev) <= aprec))
+ cvg = 2;
+ done = true;
+ elseif (iter == niter)
+ cvg = 0;
+ done = true;
+ else
+ chgprev = chg;
+ endif
+
+ endwhile
+
+ ## return result
+ p_res = pbest;
+ objf = fbest;
+ outp.niter = iter;
+
+endfunction
Modified: trunk/octave-forge/main/optim/inst/private/__lm_svd__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/private/__lm_svd__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -145,7 +145,7 @@
epsLlast=1;
epstab=[.1, 1, 1e2, 1e4, 1e6];
ac_idx = true (n_lcstr + n_gencstr, 1); % all constraints
- nc_idx = false (n_lcstr + n_gencstr, 1); % non of all constraints
+ nc_idx = false (n_lcstr + n_gencstr, 1); % none of all constraints
gc_idx = cat (1, false (n_lcstr, 1), true (n_gencstr, 1)); % gen. constr.
lc_idx = ~gc_idx;
@@ -159,10 +159,10 @@
# violated at start
if (any (c_act))
if (n_gencstr > 0)
+ %% full gradient is needed later
dct = df_cstr (p, ac_idx, ...
setfield (dfdp_hook, 'f', v_cstr));
dct(:, fixed) = 0; % for user supplied dfdp; necessary?
- dc = dct.';
dcat = dct(c_act, :);
else
dcat = df_cstr (p, c_act, ...
@@ -420,6 +420,12 @@
if (any(abs(chg) > 0.1*aprec))%--- % only worth evaluating
% function if there is some non-miniscule
% change
+ %% In the code of the outer loop before the inner loop pbest is
+ %% actually identical to p, since once they deviate, the outer
+ %% loop will not be repeated. Though the inner loop can still be
+ %% repeated in this case, pbest is not used in it. Since pprev
+ %% is set from pbest in the outer loop before the inner loop, it
+ %% is also identical to p up to here.
p=chg+pprev;
%% since the projection method may have slightly violated
%% constraints due to inaccuracy, correct parameters to bounds
Modified: trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m
===================================================================
--- trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-06-12 13:42:55 UTC (rev 10621)
+++ trunk/octave-forge/main/optim/inst/private/__nonlin_residmin__.m 2012-06-12 14:08:57 UTC (rev 10622)
@@ -69,6 +69,10 @@
"fixed", [], \
"inequc", [], \
"equc", [], \
+ "inequc_f_idx", false, \
+ "inequc_df_idx", false, \
+ "equc_f_idx", false, \
+ "equc_df_idx", false, \
"weights", [], \
"TolFun", stol_default, \
"MaxIter", [], \
@@ -142,12 +146,35 @@
mc_struct = isstruct (mc);
emc_struct = isstruct (emc);
- ## correct "_pstruct" settings if functions are not supplied
+ ## correct "_pstruct" settings if functions are not supplied, handle
+ ## constraint functions not honoring indices
if (isempty (dfdp)) dfdp_pstruct = false; endif
- if (isempty (f_genicstr)) f_inequc_pstruct = false; endif
- if (isempty (f_genecstr)) f_equc_pstruct = false; endif
- if (! user_df_gencstr) df_inequc_pstruct = false; endif
- if (! user_df_genecstr) df_equc_pstruct = false; endif
+ if (isempty (f_genicstr))
+ f_inequc_pstruct = false;
+ elseif (! optimget (settings, "inequc_f_idx", false))
+ f_genicstr = @ (p, varargin) apply_idx_if_given \
+ (f_genicstr (p, varargin{:}), varargin{:});
+ endif
+ if (isempty (f_genecstr))
+ f_equc_pstruct = false;
+ elseif (! optimget (settings, "equc_f_idx", false))
+ f_genecstr = @ (p, varargin) apply_idx_if_given \
+ (f_genecstr (p, varargin{:}), varargin{:});
+ endif
+ if (user_df_gencstr)
+ if (! optimget (settings, "inequc_df_idx", false))
+ df_gencstr = @ (varargin) df_gencstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_inequc_pstruct = false;
+ endif
+ if (user_df_genecstr)
+ if (! optimget (settings, "equc_df_idx", false))
+ df_genecstr = @ (varargin) df_genecstr (varargin{:})(varargin{2}, :);
+ endif
+ else
+ df_equc_pstruct = false;
+ endif
## some settings require a parameter order
if (pin_struct || ! isempty (pconf) || f_inequc_pstruct || \
@@ -995,3 +1022,11 @@
endif
endfunction
+
+function ret = apply_idx_if_given (ret, varargin)
+
+ if (nargin > 1)
+ ret = ret(varargin{1});
+ endif
+
+endfunction
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
