octave_packages/optim-1.2.0/private/__sqp__.m

   1 ## Copyright (C) 2012 Olaf Till <i7tiol@t-online.de>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 3 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## The simulated annealing code is translated and adapted from siman.c,
  17 ## written by Mark Galassi, of the GNU Scientific Library.
  18
  19 function [p_res, objf, cvg, outp] = __sqp__ (f, pin, hook)
  20
  21   ## needed for some anonymous functions
  22   if (exist ("ifelse") != 5)
  23     ifelse = @ scalar_ifelse;
  24   endif
  25
  26   n = length (pin);
  27
  28   ## passed constraints
  29   mc = hook.mc; # matrix of linear constraints
  30   vc = hook.vc; # vector of linear constraints
  31   f_cstr = hook.f_cstr; # function of all constraints
  32   df_cstr = hook.df_cstr; # function of derivatives of all constraints
  33   n_gencstr = hook.n_gencstr; # number of non-linear constraints
  34   eq_idx = hook.eq_idx; # logical index of equality constraints in all
  35                                 # constraints
  36   lbound = hook.lbound; # bounds, subset of linear inequality
  37   ubound = hook.ubound; # constraints in mc and vc
  38
  39   ## passed values of constraints for initial parameters
  40   pin_cstr = hook.pin_cstr;
  41
  42   ## passed return value of f for initial parameters
  43   f_pin = hook.f_pin;
  44
  45   ## passed function for gradient of objective function
  46   grad_f = hook.dfdp;
  47
  48   ## passed function for hessian of objective function
  49   if (isempty (hessian = hook.hessian))
  50     user_hessian = false;
  51     R = eye (n);
  52   else
  53     user_hessian = true;
  54   endif
  55
  56   ## passed function for complementary pivoting
  57   cpiv = hook.cpiv;
  58
  59   ## passed options
  60   ftol = hook.TolFun;
  61   niter = hook.MaxIter;
  62   if (isempty (niter)) niter = 20; endif
  63   fixed = hook.fixed;
  64
  65   ## some useful variables derived from passed variables
  66   ##
  67   ## ...
  68
  69   nz = 20 * eps; # This is arbitrary. Accuracy of equality constraints.
  70
  71   ## backend-specific checking of options and constraints
  72   ##
  73   if (any (pin < lbound | pin > ubound) ||
  74       any (pin_cstr.inequ.lin_except_bounds < 0) ||
  75       any (pin_cstr.inequ.gen < 0) ||
  76       any (abs (pin_cstr.equ.lin)) >= nz ||
  77       any (abs (pin_cstr.equ.gen)) >= nz)
  78     error ("Initial parameters violate constraints.");
  79   endif
  80   ##
  81   if (all (fixed))
  82     error ("no free parameters");
  83   endif
  84
  85   ## fill constant fields of hook for derivative-functions; some fields
  86   ## may be backend-specific
  87   dfdp_hook.fixed = fixed; # this may be handled by the frontend, but
  88                                 # the backend still may add to it
  89
  90   ## set up for iterations
  91   p = pin;
  92   f = f_pin;
  93   n_iter = 0;
  94   done = false;
  95
  96   while (! done)
  97
  98     niter++;
  99
 100     if (user_hessian)
 101
 102       H = hessian (p);
 103       idx = isnan (H);
 104       H(idx) = H.'(idx);
 105       if (any (isnan (H(:))))
 106         error ("some second derivatives undefined by user function");
 107       endif
 108       if (! isreal (H))
 109         error ("second derivatives given by user function not real");
 110       endif
 111       if (! issymmetric (H))
 112         error ("Hessian returned by user function not symmetric");
 113       endif
 114
 115       R = directional_discrimination (H);
 116
 117     endif
 118
 119
 120   endwhile
 121
 122   ## return result
 123
 124 endfunction
 125
 126 function R = directional_discrimination (A)
 127
 128   ## A is expected to be real and symmetric without checking
 129   ##
 130   ## "Directional discrimination" (Bard, Nonlinear Parameter Estimation,
 131   ## Academic Press, 1974). Compute R, which is "similar" to computing
 132   ## inv(H), but succeeds even if H is singular; R is positive definite
 133   ## and is tuned to avoid large steps of one parameter with respect to
 134   ## the steps of the others.
 135
 136   ## make matrix Binv for scaling
 137   Binv = diag (A);
 138   nidx = ! (idx = Binv == 0);
 139   Binv(nidx) = 1 ./ sqrt (abs (Binv(nidx)));
 140   Binv(idx) = 1;
 141   Binv = diag (Binv);
 142
 143   ## eigendecomposition of scaled hessian
 144   [V, L] = eig (Binv * A * Binv);
 145
 146   ## A is symmetric, so V and L are real, delete any imaginary parts,
 147   ## which might occur due to inaccuracy
 148   V = real (V);
 149   L = real (L);
 150
 151   ## actual directional discrimination, does not exactly follow Bard
 152   L = abs (diag (L)); # R should get positive definite
 153   L = max (L, .001 * max (L)); # avoids relatively large steps of
 154                                # parameters
 155
 156   G = Binv * V;
 157
 158   R = G * diag (1 ./ L) * G.';
 159
 160 endfunction