On Tue, Oct 04, 2011 at 03:01:08PM +0300, Fotios wrote:
> Στις 2011-10-04 14:31, ο/η Olaf Till έγραψε:
> >On Tue, Oct 04, 2011 at 01:18:11AM +0100, Carnë Draug wrote:
> >>On 4 October 2011 00:48, Fotios<fotios.kaso...@gmail.com> wrote:
> >>>Here is the deriv function with an added demo and minor correction.
> >>>
> >>>/Fotios
> >>This was supposed to be sent to the octave-forge. After talking to
> >>Fotios #octave, this function was renamed jacobs (there's already a
> >>deriv function) and added to the optim package.
> >>
> >>Carnë
> >I think such a function is really useful for the optim package, even
> >if it is intended to be used elsewhere. But in the current form, the
> >usefulness seems to be limited by the restriction that the passed
> >function 'f' is expected to compute exactly as many elements as there
> >are parameters. In most standard optimizations we either deal with
> >gradients of a scalar valued 'objective function' or with Jacobians of
> >a model function which usually returns more elements than there are
> >parameters. So I suggest that this restriction should be removed.
> >
> >If the above should be agreed, there are some minor issues, some of
> >them 'cosmetic':
> >
> >1. Since I currently try to standardize the interface to some of the
> >optimization function, could we change the order of argumets to
> >
> >'jacobs (x, f, ...)'
> >
> >?
> >
> >2. I think the "jacobs: ..." in the error messages is not necessary,
> >since the m-file in which the error originates is reported by Octave.
> >
> >3.
> >
> > if (abs (h)>= 1e-2)
> > warning ("jacobs: H is too big and the result should not be trusted");
> > endif
> >
> >How do you come to this value (1e-2) ? If there is a reason, there
> >should be probably a comment there explaining it.
> >
> >4. For the same reason as in 1., and for interfacing with some
> >optimization functions, can we
> >
> >- make the third argument a structure with a field 'h',
> >
> >- accept a further field 'fixed': a logical vector indicating for
> > which elements of 'x' no gradients are to be computed, but zero is
> > returned instead?
> >
> >
> >If there should be agreement about the main issue (dimension of 'f'),
> >I could suggest a patch for all the above, if you want me to.
> >
> >Olaf
>
> I absolutely agree with your general comment regarding optimization
> since the design variables are typically more than the objective
> functions - keep in mind though that i have not tried the function
> for mid-scale optimization problems and that i would argue about how
> useful it could be in that case.
Well, some of the optimization functions use some routine for finite
differencing, so using 'jacobs.m' instead has the potential to be
better. But I agree that the used objective or model functions would
have to be suitably programmed for this. May be it would be useful for
that if someone wrote a class for overloading "<", ">", "<=", ">=",
"max", "min" and maybe others ...
Anyway I have made a respective change, allowing arbitrary
dimensionness for 'f', as you can see in the patch.
> Regarding the interfacing comments
> you can apply all of your suggestions. There is no good reason for
> the 1e-2 warning since it always depends both on h and on the value
> of the 3rd derivative, BUT since log(|Df(h)-exact|) depends linearly
> on log(h) i though i had to prevent naive users from increasing h or
> at least warn them about the accuracy of the approximation (it is
> just a sloppy educational guess - after some testing - coded there
> with students in mind that can be removed as well). I could even
> suggest removing the user control over the step size h (demo needs
> to be removed) since this approach costs more (with respect to
> memory) compared to fd and the only gain from using it comes from
> the fact that you can get eps exact derivatives.
So I just changed the warning text to 'not allowed', sinces as you say
one cannot be sure that it is really 'too big'. Consequently, larger
values are reduced to the maximum value (1e-3 instead of 1e-2, to be
able to use ">" instead of ">=").
>
> /Fotis
I attach the patch, but as it is not very readable, also the new
version of the function.
All tests still pass and the demo works (after a slight change).
Do you think I can submit it so?
Olaf
Index: jacobs.m
===================================================================
--- jacobs.m (revision 8673)
+++ jacobs.m (working copy)
@@ -14,28 +14,33 @@
## along with this program; if not, see <http://www.gnu.org/licenses/>.
## -*- texinfo -*-
-## @deftypefn {Function File} {Df =} jacobs (@var{f}, @var{x})
-## @deftypefnx {Function File} {Df =} jacobs (@var{f}, @var{x}, @var{h})
+## @deftypefn {Function File} {Df =} jacobs (@var{x}, @var{f})
+## @deftypefnx {Function File} {Df =} jacobs (@var{x}, @var{f}, @var{hook})
## Calculate the jacobian of a function.
##
-## Let @var{f} be a user-supplied function. Given a point @var{x} at which
-## we seek for the Jacobian, the function @command{jacobs} returns the Jacobian
-## matrix @code{d(f(1), @dots{}, df(n))/d(x(1), @dots{}, x(n))}. The
-## function uses the complex step method and thus can be applied to real
-## analytic functions. The optional argument @var{h} can be used to define
-## the magnitude of the complex step and defaults to 1e-20.
+## Let @var{f} be a user-supplied function. Given a point @var{x} at
+## which we seek for the Jacobian, the function @command{jacobs} returns
+## the Jacobian matrix @code{d(f(1), @dots{}, df(end))/d(x(1), @dots{},
+## x(n))}. The function uses the complex step method and thus can be
+## applied to real analytic functions. The optional argument @var{hook}
+## is a structure with additional options. @var{hook} can have the
+## following fields: @code{h} can be used to define the magnitude of the
+## complex step and defaults to 1e-20; steps larger than 1e-3 are not
+## allowed. @code{fixed} is a logical vector internally usable by some
+## optimization functions; it indicates for which elements of @var{x} no
+## gradient should be computed, but zero should be returned.
##
## For example:
##
## @example
## @group
## f = @@(x) [x(1)^2 + x(2); x(2)*exp(x(1))];
-## Df = jacobs (f, [1,2])
+## Df = jacobs ([1, 2], f)
## @end group
## @end example
## @end deftypefn
-function Df = jacobs (f, x, h = 1e-20)
+function Df = jacobs (x, f, hook)
if ( (nargin < 2) || (nargin > 3) )
print_usage ();
@@ -45,20 +50,61 @@
f = str2func (f, "global");
endif
- if (! (isscalar (h)))
- error ("jacobs: H must be a scalar");
+ n = numel (x);
+
+ default_h = 1e-20;
+ max_h = 1e-3; # must be positive
+
+ if (nargin > 2)
+
+ if (isfield (hook, "h"))
+ if (! (isscalar (hook.h)))
+ error ("complex step magnitude must be a scalar");
+ endif
+ if (abs (hook.h) > max_h)
+ warning ("complex step magnitude larger than allowed, set to %e", \
+ max_h)
+ h = max_h;
+ else
+ h = hook.h;
+ endif
+ else
+ h = default_h;
+ endif
+
+ if (isfield (hook, "fixed"))
+ if (numel (hook.fixed) != n)
+ error ("index of fixed parameters has wrong dimensions");
+ endif
+ fixed = hook.fixed;
+ else
+ fixed = false (n, 1);
+ endif
+
+ else
+
+ h = default_h;
+ fixed = false (n, 1);
+
endif
- if (abs (h) >= 1e-2)
- warning ("jacobs: H is too big and the result should not be trusted");
+ if (all (fixed))
+ error ("all elements of 'x' are fixed");
endif
- n = numel (x);
- Df = zeros (n, n);
-
x = repmat (x(:), 1, n) + h * 1i * eye (n);
- for count = 1:n
+ idx = find (! fixed);
+
+ ## after first evaluation, dimensionness of 'f' is known
+ t_Df = imag (f (x(:, idx(1)))(:));
+ dim = numel (t_Df);
+
+ Df = zeros (dim, n);
+
+ Df(:, idx(1)) = t_Df;
+
+ for count = idx(2:end)
Df(:, count) = imag (f (x(:, count))(:));
endfor
@@ -66,10 +112,10 @@
endfunction
-%!assert (jacobs (@(x) x, 1), 1)
-%!assert (jacobs (@(x) x^2, 6), 12)
-%!assert (jacobs (@(x) [x(1)^2; x(1)*x(2)], [1; 1]), [2, 0; 1, 1])
-%!assert (jacobs (@(x) [x(1)^2 + x(2); x(2)*exp(x(1))], [1; 2]), [2, 1; 2*exp(1), exp(1)])
+%!assert (jacobs (1, @(x) x), 1)
+%!assert (jacobs (6, @(x) x^2), 12)
+%!assert (jacobs ([1; 1], @(x) [x(1)^2; x(1)*x(2)]), [2, 0; 1, 1])
+%!assert (jacobs ([1; 2], @(x) [x(1)^2 + x(2); x(2)*exp(x(1))]), [2, 1; 2*exp(1), exp(1)])
%% Test input validation
%!error jacobs ()
@@ -80,10 +126,10 @@
%!demo
%! # Relative error against several h-values
-%! k = 1:20; h = 10 .^ (-k); x = 0.3*pi;
+%! k = 3:20; h = 10 .^ (-k); x = 0.3*pi;
%! err = zeros (1, numel (k));
-%! for count = k
-%! err(count) = abs (jacobs (@sin, x, h(count)) - cos (x)) / abs (cos (x)) + eps;
+%! for count = 1 : numel (k)
+%! err(count) = abs (jacobs (x, @sin, struct ("h", h(count))) - cos (x)) / abs (cos (x)) + eps;
%! endfor
%! loglog (h, err); grid minor;
%! xlabel ("h"); ylabel ("|Df(x) - cos(x)| / |cos(x)|")
## Copyright (C) 2011 Fotios Kasolis <fotios.kaso...@gmail.com>
##
## 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} {Df =} jacobs (@var{x}, @var{f})
## @deftypefnx {Function File} {Df =} jacobs (@var{x}, @var{f}, @var{hook})
## Calculate the jacobian of a function.
##
## Let @var{f} be a user-supplied function. Given a point @var{x} at
## which we seek for the Jacobian, the function @command{jacobs} returns
## the Jacobian matrix @code{d(f(1), @dots{}, df(end))/d(x(1), @dots{},
## x(n))}. The function uses the complex step method and thus can be
## applied to real analytic functions. The optional argument @var{hook}
## is a structure with additional options. @var{hook} can have the
## following fields: @code{h} can be used to define the magnitude of the
## complex step and defaults to 1e-20; steps larger than 1e-3 are not
## allowed. @code{fixed} is a logical vector internally usable by some
## optimization functions; it indicates for which elements of @var{x} no
## gradient should be computed, but zero should be returned.
##
## For example:
##
## @example
## @group
## f = @@(x) [x(1)^2 + x(2); x(2)*exp(x(1))];
## Df = jacobs ([1, 2], f)
## @end group
## @end example
## @end deftypefn
function Df = jacobs (x, f, hook)
if ( (nargin < 2) || (nargin > 3) )
print_usage ();
endif
if (ischar (f))
f = str2func (f, "global");
endif
n = numel (x);
default_h = 1e-20;
max_h = 1e-3; # must be positive
if (nargin > 2)
if (isfield (hook, "h"))
if (! (isscalar (hook.h)))
error ("complex step magnitude must be a scalar");
endif
if (abs (hook.h) > max_h)
warning ("complex step magnitude larger than allowed, set to %e", \
max_h)
h = max_h;
else
h = hook.h;
endif
else
h = default_h;
endif
if (isfield (hook, "fixed"))
if (numel (hook.fixed) != n)
error ("index of fixed parameters has wrong dimensions");
endif
fixed = hook.fixed;
else
fixed = false (n, 1);
endif
else
h = default_h;
fixed = false (n, 1);
endif
if (all (fixed))
error ("all elements of 'x' are fixed");
endif
x = repmat (x(:), 1, n) + h * 1i * eye (n);
idx = find (! fixed);
## after first evaluation, dimensionness of 'f' is known
t_Df = imag (f (x(:, idx(1)))(:));
dim = numel (t_Df);
Df = zeros (dim, n);
Df(:, idx(1)) = t_Df;
for count = idx(2:end)
Df(:, count) = imag (f (x(:, count))(:));
endfor
Df /= h;
endfunction
%!assert (jacobs (1, @(x) x), 1)
%!assert (jacobs (6, @(x) x^2), 12)
%!assert (jacobs ([1; 1], @(x) [x(1)^2; x(1)*x(2)]), [2, 0; 1, 1])
%!assert (jacobs ([1; 2], @(x) [x(1)^2 + x(2); x(2)*exp(x(1))]), [2, 1;
2*exp(1), exp(1)])
%% Test input validation
%!error jacobs ()
%!error jacobs (1)
%!error jacobs (1, 2, 3, 4)
%!error jacobs (@sin, 1, [1, 1])
%!error jacobs (@sin, 1, ones(2, 2))
%!demo
%! # Relative error against several h-values
%! k = 3:20; h = 10 .^ (-k); x = 0.3*pi;
%! err = zeros (1, numel (k));
%! for count = 1 : numel (k)
%! err(count) = abs (jacobs (x, @sin, struct ("h", h(count))) - cos (x)) /
abs (cos (x)) + eps;
%! endfor
%! loglog (h, err); grid minor;
%! xlabel ("h"); ylabel ("|Df(x) - cos(x)| / |cos(x)|")
%! title ("f(x)=sin(x), f'(x)=cos(x) at x = 0.3pi")
------------------------------------------------------------------------------
All the data continuously generated in your IT infrastructure contains a
definitive record of customers, application performance, security
threats, fraudulent activity and more. Splunk takes this data and makes
sense of it. Business sense. IT sense. Common sense.
http://p.sf.net/sfu/splunk-d2dcopy1
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
