octave_packages/optim-1.2.0/jacobs.m

   1 ## Copyright (C) 2011 Fotios Kasolis <fotios.kasolis@gmail.com>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {Df =} jacobs (@var{x}, @var{f})
  18 ## @deftypefnx {Function File} {Df =} jacobs (@var{x}, @var{f}, @var{hook})
  19 ## Calculate the jacobian of a function using the complex step method.
  20 ##
  21 ## Let @var{f} be a user-supplied function. Given a point @var{x} at
  22 ## which we seek for the Jacobian, the function @command{jacobs} returns
  23 ## the Jacobian matrix @code{d(f(1), @dots{}, df(end))/d(x(1), @dots{},
  24 ## x(n))}. The function uses the complex step method and thus can be
  25 ## applied to real analytic functions.
  26 ##
  27 ## The optional argument @var{hook} is a structure with additional options. @var{hook}
  28 ## can have the following fields:
  29 ## @itemize @bullet
  30 ## @item
  31 ## @code{h} - can be used to define the magnitude of the complex step and defaults
  32 ## to 1e-20; steps larger than 1e-3 are not allowed.
  33 ## @item
  34 ## @code{fixed} - is a logical vector internally usable by some optimization
  35 ## functions; it indicates for which elements of @var{x} no gradient should be
  36 ## computed, but zero should be returned.
  37 ## @end itemize
  38 ##
  39 ## For example:
  40 ##
  41 ## @example
  42 ## @group
  43 ## f = @@(x) [x(1)^2 + x(2); x(2)*exp(x(1))];
  44 ## Df = jacobs ([1, 2], f)
  45 ## @end group
  46 ## @end example
  47 ## @end deftypefn
  48
  49 function Df = jacobs (x, f, hook)
  50
  51   if ( (nargin < 2) || (nargin > 3) )
  52     print_usage ();
  53   endif
  54
  55   if (ischar (f))
  56     f = str2func (f, "global");
  57   endif
  58
  59   n  = numel (x);
  60
  61   default_h = 1e-20;
  62   max_h = 1e-3; # must be positive
  63
  64   if (nargin > 2)
  65
  66     if (isfield (hook, "h"))
  67       if (! (isscalar (hook.h)))
  68         error ("complex step magnitude must be a scalar");
  69       endif
  70       if (abs (hook.h) > max_h)
  71         warning ("complex step magnitude larger than allowed, set to %e", ...
  72                  max_h)
  73         h = max_h;
  74       else
  75         h = hook.h;
  76       endif
  77     else
  78       h = default_h;
  79     endif
  80
  81     if (isfield (hook, "fixed"))
  82       if (numel (hook.fixed) != n)
  83         error ("index of fixed parameters has wrong dimensions");
  84       endif
  85       fixed = hook.fixed;
  86     else
  87       fixed = false (n, 1);
  88     endif
  89
  90   else
  91
  92     h = default_h;
  93     fixed = false (n, 1);
  94
  95   endif
  96
  97   if (all (fixed))
  98     error ("all elements of 'x' are fixed");
  99   endif
 100
 101   x = repmat (x(:), 1, n) + h * 1i * eye (n);
 102
 103   idx = find (! fixed);
 104
 105   ## after first evaluation, dimensionness of 'f' is known
 106   t_Df = imag (f (x(:, idx(1)))(:));
 107   dim = numel (t_Df);
 108
 109   Df = zeros (dim, n);
 110
 111   Df(:, idx(1)) = t_Df;
 112
 113   for count = idx(2:end)
 114     Df(:, count) = imag (f (x(:, count))(:));
 115   endfor
 116
 117   Df /=  h;
 118
 119 endfunction
 120
 121 %!assert (jacobs (1, @(x) x), 1)
 122 %!assert (jacobs (6, @(x) x^2), 12)
 123 %!assert (jacobs ([1; 1], @(x) [x(1)^2; x(1)*x(2)]), [2, 0; 1, 1])
 124 %!assert (jacobs ([1; 2], @(x) [x(1)^2 + x(2); x(2)*exp(x(1))]), [2, 1; 2*exp(1), exp(1)])
 125
 126 %% Test input validation
 127 %!error jacobs ()
 128 %!error jacobs (1)
 129 %!error jacobs (1, 2, 3, 4)
 130 %!error jacobs (@sin, 1, [1, 1])
 131 %!error jacobs (@sin, 1, ones(2, 2))
 132
 133 %!demo
 134 %! # Relative error against several h-values
 135 %! k = 3:20; h = 10 .^ (-k); x = 0.3*pi;
 136 %! err = zeros (1, numel (k));
 137 %! for count = 1 : numel (k)
 138 %!   err(count) = abs (jacobs (x, @sin, struct ("h", h(count))) - cos (x)) / abs (cos (x)) + eps;
 139 %! endfor
 140 %! loglog (h, err); grid minor;
 141 %! xlabel ("h"); ylabel ("|Df(x) - cos(x)| / |cos(x)|")
 142 %! title ("f(x)=sin(x), f'(x)=cos(x) at x = 0.3pi")