Code covered by the BSD License

Complex step Hessian

by

Yi Cao (view profile)

02 Jan 2008 (Updated )

Calculate Hessian using complex step differentiation

[A,g,z]=hessiancsd(fun,x)
```function [A,g,z]=hessiancsd(fun,x)
% HESSIANCSD    Complex Step Hessian
% H = hessiancsd(f,x) returns the numberical (m x n) Hessian matrix of a
% scalar function, f(x) at the refdrence point, x (n-vector).
% [H,g] = jacobiancsd(f,x) also returns the gradient function, g=f'(x).
% With all outputs
% [H,g,z] = jacobiancsd(f,x) returns the function value, z=f(x) as the third output.
%
% Example
% f=@(x)x(1)*x(2)^2+x(1)*x(3)^2+x(1)^2;
% x=1:3;
% [H,g,z]=hessiancsd(f,x)
% H =
%     2     4     6
%     4     2     0
%     6     0     2
% g =
%    15
%     4
%     6
% z =
%    14
%
% By Yi Cao at Cranfield University, 02/01/2008
%
if nargout>2
z=fun(x);                       % function value
end
n=numel(x);                         % size of independent
A=zeros(n);                         % allocate memory for the Hessian matrix
g=zeros(n,1);                       % allocate memory for the gradient matrix
h=sqrt(eps);                        % differentiation step size
h2=h*h;                             % square of step size
for k=1:n                           % loop for each independent variable
x1=x;                           % reference point
x1(k)=x1(k)+h*i;                % increment in kth independent variable
if nargout>1
v=fun(x1);                  % function call with a comlex step increment
end
for l=k:n                       % loop for off diagonal Hessian
x2=x1;                      % reference with kth increment
x2(l)=x2(l)+h;              % kth + lth (positive real) increment
u1=fun(x2);                 % function call with a double increment
x2(l)=x1(l)-h;              % kth + lth (negative real) increment
u2=fun(x2);                 % function call with a double increment
A(k,l)=imag(u1-u2)/h2/2;    % Hessian (central + complex step)
A(l,k)=A(k,l);              % symmetric
end
end
```