spCompHessian

Parameters

g

a Scilab function computing the gradient

sp

a sparse matrix

options

a sequence of optional named arguments allowing to customize the engine (see available options below)

hess

the returned engine, an MList of _spCompHessian type (fields are described below)

Description

First and foremost, you need a Scilab function g which computes the gradient (preferably exact) of the function for wich the Hessian is to be approximated.

The spCompHessian function consructs a compressed Hessian engine from the given sparsity pattern using column intersection graph coloring. Internally, after using a heuristic coloring algorithm, a "seed matrix" is generated, available as hess.seed, and later used by the engine to compute (approximate) directional derivatives of the gradient function g in each seed column direction. Then the uncompressed Hessian is recovered. The typical call sequence starts by calling spCompHessian then continues by invoking the returned MList with the vectors at which the Hessian is to be computed:

hess = spCompHessian(g,sp);
H0 = hess(x0);
H1 = hess(x1);

The options allows to change some parameters of the coloring algorithm and to choose the derivative approximation scheme as well as the step vector. The available options are the following:

• Coloring : the allowed values are "ACYCLIC_FOR_INDIRECT_RECOVERY" (the default), "STAR"","RESTRICTED_STAR" and "DISTANCE_TWO"

• Ordering : this field allows to change the way columns of the sparsity pattern are ordered before applying the coloring algorithm. The possible values are "NATURAL" (the default), "SMALLEST_LAST", "LARGEST_FIRST", "INCIDENCE_DEGREE", "RANDOM".

• FiniteDifferenceType : this field gives the chosen approximation scheme. The possible values are "CENTERED" (the default order 2 centered scheme), "FORWARD" (the usual forward order 1 scheme) and "COMPLEXSTEP" (the order 2 complex step scheme). If you choose the complex step scheme check that g accepts complex input and handles it correctly (for example transposition operator has to be the dot-prefixed non-conjugate transposition .').

• FiniteDifferenceStepSize : this is the relative stepsize (see also "TypicalX" below). For a given x0 the actual step is

  sign(x0).*FiniteDifferenceStepSize.*max(TypicalX,x0)

  It can be a scalar, or a vector with the same size as x0. Its default value depends on the chosen approximation scheme: sqrt(%eps) for the "FORWARD" scheme, %eps^(1/3) for the "CENTERED" scheme and 1e-100 for the "COMPLEXSTEP" scheme.

• TypicalX : a vector of the same size as x0, with typical magnitude of components. The default value is ones(x0).

• Vectorized : "on" or "off" : allows to inform the engine if g can be called with a matrix X as argument, with the convention that g(X)(:,i) == g(X(:,i)).

  for each column of X.

The hess MList can be invoked with more than one arguments. In this case the remaining argument after the first are used by the engine as complimentary arguments when calling the gradient g.

Examples

function f=fun(x)
    n = length(x);
    f = norm(x(1:n-2).*x(2:n-1).*x(3:n))^2;
endfunction

function g=grad(x)
    n = length(x);
    g = zeros(n,1);
    k = 3:n-2;
    g(1)   = 2*x(1)*x(2)^2*x(3)^2;
    g(2)   = 2*(x(2)*x(3)^2*x(4)^2+x(1)^2*x(2)*x(3)^2);
    g(k)   = 2*(x(k).*x(k+1).^2.*x(k+2).^2+x(k-1).^2.*x(k).*x(k+1).^2+x(k-2).^2.*x(k-1).^2.*x(k))
    g(n-1) = 2*(x(n-3)^2*x(n-2)^2*x(n-1)+x(n-2)^2*x(n-1)*x(n)^2); 
    g(n)   = 2*x(n-2)^2*x(n-1)^2*x(n);
endfunction

n = 1000;
// Hessian is pentadiagonal
i = (1:n).*.ones(n,1);
sp = sparse(abs(i-i') <= 2);

hess = spCompHessian(grad, sp, FiniteDifferenceType = "COMPLEXSTEP")
x0 = rand(n,1);
tic
H = hess(x0);
toc
tic
Hd = numderivative(grad,x0);
toc
disp(max(abs(H-Hd)))

See also

• spCompJacobian — Consruct a compressed Jacobian engine from its sparsity pattern using column intersection graph coloring
• numderivative — approximate derivatives of a function (Jacobian or Hessian)

Bibliography

A. H. Gebremedhin, D. C. Nguyen, Md. M. A. Patwary, A. Pothen}, ColPack: Software for graph coloring and related problems in scientific computing, ACM Trans. Math. Softw., N. 40(1), pp. 1-31, 2013, https://doi.org/10.1145/2513109.2513110.