/usr/share/octave/packages/optim-1.3.0/nonlin_residmin.m is in octave-optim 1.3.0-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 | ## Copyright (C) 2010-2013 Olaf Till <i7tiol@t-online.de>
##
## 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/>.
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{p}, @var{resid}, @var{cvg}, @var{outp}] =} nonlin_residmin (@var{f}, @var{pin})
## @deftypefnx {Function File} {[@var{p}, @var{resid}, @var{cvg}, @var{outp}] =} nonlin_residmin (@var{f}, @var{pin}, @var{settings})
## Frontend for nonlinear minimization of residuals returned by a model
## function.
##
## The functions supplied by the user have a minimal
## interface; any additionally needed constants (e.g. observed values)
## can be supplied by wrapping the user functions into anonymous
## functions.
##
## 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.
##
## @var{f}: function returning the array of residuals. It gets a column
## vector of real parameters as argument. In gradient determination,
## this function may be called with an informational second argument,
## whose content depends on the function for gradient determination.
##
## @var{pin}: real column vector of initial parameters.
##
## @var{settings}: structure whose fields stand for optional settings
## referred to below. The fields can be set by @code{optimset()} with
## Octave versions 3.3.55 or greater; with older Octave versions, the
## fields must be set directly as structure-fields in the correct case.
##
## The returned values are the column vector of final parameters
## @var{p}, the final array of residuals @var{resid}, an integer
## @var{cvg} indicating if and how optimization succeeded or failed, and
## a structure @var{outp} with additional information, curently with
## only one field: @var{niter}, the number of iterations. @var{cvg} is
## greater than zero for success and less than or equal to zero for
## 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).
##
## @var{settings}:
##
## @code{Algorithm}: String specifying the backend. Default:
## @code{"lm_svd_feasible"}. The latter is currently the only backend
## distributed with this package. It is described in a separate section
## below.
##
## @code{dfdp}: Function computing the jacobian of the residuals with
## respect to the parameters, assuming residuals are reshaped to a
## vector. Default: finite differences. Will be called with the column
## vector of parameters and an informational structure as arguments. The
## structure has the fields @code{f}: value of residuals for current
## parameters, reshaped to a column vector, @code{fixed}: logical vector
## indicating which parameters are not optimized, so these partial
## derivatives need not be computed and can be set to zero,
## @code{diffp}, @code{diff_onesided}, @code{lbound}, @code{ubound}:
## identical to the user settings of this name, @code{plabels}:
## 1-dimensional cell-array of column-cell-arrays, each column with
## labels for all parameters, the first column contains the numerical
## indices of the parameters. The default jacobian function will call
## the model function with the second argument set with fields @code{f}:
## as the @code{f} passed to the jacobian function, @code{plabels}:
## cell-array of 1x1 cell-arrays with the entries of the
## column-cell-arrays of @code{plabels} as passed to the jacobian
## function corresponding to current parameter, @code{side}: @code{0}
## for one-sided interval, @code{1} or @code{2}, respectively, for the
## sides of a two-sided interval, and @code{parallel}: logical scalar
## indicating parallel computation of partial derivatives.
##
## @code{diffp}: column vector of fractional intervals (doubled for
## central intervals) supposed to be used by jacobian functions
## performing finite differencing. Default: @code{.001 * ones (size (parameters))}. The default jacobian function will use these as
## absolute intervals for parameters with value zero.
##
## @code{diff_onesided}: logical column vector indicating that one-sided
## intervals should be used by jacobian functions performing finite
## differencing. Default: @code{false (size (parameters))}.
##
## @code{complex_step_derivative_f},
## @code{complex_step_derivative_inequc},
## @code{complex_step_derivative_equc}: logical scalars, default: false.
## Estimate Jacobian of model function, general inequality constraints,
## and general equality constraints, respectively, with complex step
## derivative approximation. Use only if you know that your model
## function, function of general inequality constraints, or function of
## general equality constraints, respectively, is suitable for this. No
## user function for the respective Jacobian must be specified.
##
## @code{cstep}: scalar step size for complex step derivative
## approximation. Default: 1e-20.
##
## @code{parallel_local}: logical scalar, default: false. Estimate
## Jacobians of model function and of constraints in parallel processes.
## Works for default finite difference Jacobian function and for complex
## step derivatives. Due to overhead, a speed advantage can only be
## expected if model function or constraint functions are time consuming
## enough.
##
## @code{fixed}: logical column vector indicating which parameters
## should not be optimized, but kept to their inital value. Fixing is
## done independently of the backend, but the backend may choose to fix
## additional parameters under certain conditions.
##
## @code{lbound}, @code{ubound}: column vectors of lower and upper
## bounds for parameters. Default: @code{-Inf} and @code{+Inf},
## respectively. The bounds are non-strict, i.e. parameters are allowed
## to be exactly equal to a bound. The default jacobian function will
## respect bounds (but no further inequality constraints) in finite
## differencing.
##
## @code{inequc}: Further inequality constraints. Cell-array containing
## up to four entries, two entries for linear inequality constraints
## and/or one or two entries for general inequality constraints. Either
## linear or general constraints may be the first entries, but the two
## entries for linear constraints must be adjacent and, if two entries
## are given for general constraints, they also must be adjacent. The
## 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 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{f_inequc_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{df_inequc_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}).
##
## @code{cpiv}: Function for complementary pivoting, usable in
## algorithms for constraints. Default: @ cpiv_bard. Only the default
## function is supplied with the package.
##
## @code{weights}: Array of weights for the residuals. Dimensions must
## match.
##
## @code{TolFun}: Minimum fractional improvement in objective function
## (e.g. sum of squares) in an iteration (termination criterium). Default:
## .0001.
##
## @code{MaxIter}: Maximum number of iterations (termination criterium).
## Default: backend-specific.
##
## @code{fract_prec}: Column Vector, minimum fractional change of
## parameters in an iteration (termination criterium if violated in two
## consecutive iterations). Default: backend-specific.
##
## @code{max_fract_change}: Column Vector, enforced maximum fractional
## change in parameters in an iteration. Default: backend-specific.
##
## @code{Display}: String indicating the degree of verbosity. Default:
## @code{"off"}. Possible values are currently @code{"off"} (no
## messages) and @code{"iter"} (some messages after each iteration).
## Support of this setting and its exact interpretation are
## backend-specific.
##
## @code{plot_cmd}: Function enabling backend to plot results or
## intermediate results. Will be called with current computed
## residuals. Default: plot nothing.
##
## @code{debug}: Logical scalar, default: @code{false}. Will be passed
## to the backend, which might print debugging information if true.
##
## Structure-based parameter handling
##
## The setting @code{param_order} is a cell-array with names of the
## optimized parameters. If not given, and initial parameters are a
## structure, all parameters in the structure are optimized. If initial
## parameters are a structure, it is an error if @code{param_order} is
## not given and there are any non-structure-based configuration items
## or functions.
##
## The initial parameters @var{pin} can be given as a structure
## containing at least all fields named in @code{param_order}. In this
## case the returned parameters @var{p} will also be a structure.
##
## Each user-supplied function can be called with the argument
## containing the current parameters being a structure instead of a
## column vector. For this, a corresponding setting must be set to
## @code{true}: @code{f_pstruct} (model function), @code{df_pstruct}
## (jacobian of model function), @code{f_inequc_pstruct} (general
## inequality constraints), @code{df_inequc_pstruct} (jacobian of
## general inequality constraints), @code{f_equc_pstruct} (general
## equality constraints), and @code{df_equc_pstruct} (jacobian of
## general equality constraints). If a jacobian-function is configured
## in such a way, it must return the columns of the jacobian as fields
## of a structure under the respective parameter names.
##
## Similarly, for specifying linear constraints, instead of the matrix
## (called @code{m} above), a structure containing the rows of the
## matrix in fields under the respective parameter names can be given.
## In this case, rows containing only zeros need not be given.
##
## The vector-based settings @code{lbound}, @code{ubound},
## @code{fixed}, @code{diffp}, @code{diff_onesided}, @code{fract_prec},
## and @code{max_fract_change} can be replaced by the setting
## @code{param_config}. It is a structure that can contain fields named
## in @code{param_order}. For each such field, there may be subfields
## with the same names as the above vector-based settings, but
## containing a scalar value for the respective parameter. If
## @code{param_config} is specified, none of the above
## vector/matrix-based settings may be used.
##
## Additionally, named parameters are allowed to be non-scalar real
## arrays. In this case, their dimensions are given by the setting
## @code{param_dims}, a cell-array of dimension vectors, each containing
## at least two dimensions; if not given, dimensions are taken from the
## initial parameters, if these are given in a structure. Any
## vector-based settings or not structure-based linear constraints then
## must correspond to an order of parameters with all parameters
## reshaped to vectors and concatenated in the user-given order of
## parameter names. Structure-based settings or structure-based initial
## parameters must contain arrays with dimensions reshapable to those of
## the respective parameters.
##
## Description of backends (currently only one)
##
## "lm_svd_feasible"
##
## A Levenberg/Marquardt algorithm using singular value decomposition
## and featuring constraints which must be met by the initial parameters
## and are attempted to be kept met throughout the optimization.
##
## Parameters with identical lower and upper bounds will be fixed.
##
## Returned value @var{cvg} will be @code{0}, @code{2}, or @code{3}.
##
## 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
## @code{plot_cmd} is evaluated for each iteration, and some further
## diagnostics may be printed.
##
## Specific option: @code{lm_svd_feasible_alt_s}: if falling back to
## nearly gradient descent, do it more like original Levenberg/Marquardt
## method, with descent in each gradient component; for testing only.
##
## @seealso {nonlin_curvefit}
## @end deftypefn
function [p, resid, cvg, outp] = nonlin_residmin (varargin)
if (nargin == 1)
p = __nonlin_residmin__ (varargin{1});
return;
endif
if (nargin < 2 || nargin > 3)
print_usage ();
endif
if (nargin == 2)
varargin{3} = struct ();
endif
varargin{4} = struct ();
[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)
|