Example 5: Hessian Approximation
This example uses the same data as in the One-Sided Differences example. In this example numerical differentiation is used to approximate the Hessian matrix of 
. This symmetric Hessian matrix is ![Hf = \left[ {\begin{array}{*{20}c} {\frac{{\partial ^2 f}} {{\partial y_1^2 }}} & {\frac{{\partial ^2 f}} {{\partial y_1 \partial y_2 }}} \\ {\frac{{\partial ^2 f}} {{\partial y_1 \partial y_2 }}} & {\frac{{\partial ^2 f}} {{\partial y_2^2 }}} \\ \end{array} } \right]](eqn_0207.png)
Our method is based on casting the matrix 
as the 2 by 2 Jacobian matrix of the gradient function. Each inner evaluation of the gradient function is itself computed using NumericalDerivatives
. Central differences are used for both the inner and outer numerical differentiation. Because of the inherent error in both processes, the expected accuracy is about the 
power of machine precision. Note that the approximation obtained is not symmetric. However, the difference between the off-diagonal elements provides an error estimate of that term.
import com.imsl.math.*;
import java.text.*;
public class NumericalDerivativesEx5 extends NumericalDerivatives {
static int m = 1, n = 2;
static double a, b, c, v = 0.0;
class InnerNumericalDerivatives extends NumericalDerivatives {
public InnerNumericalDerivatives(NumericalDerivatives.Function fcn) {
super(fcn);
}
// Override evaluateF.
public double[] evaluateF(int varIndex, double[] y) {
double[] valueF = new double[m];
valueF[0] = a * Math.exp(b * y[0]) + c * y[0] * y[1] * y[1];
return valueF;
}
}
public NumericalDerivativesEx5(NumericalDerivatives.Function fcn) {
super(fcn);
}
// Override evaluateF.
public double[] evaluateF(int varIndex, double[] y) {
int[] iopt = new int[n];
double[] valueF = new double[m];
double[] scale = new double[n];
double[] fac = new double[n];
// This is the analytic gradient. for comparison only.
// stateh(1)=a*b*exp(b*y(1))+c*y(2)**2
// stateh(2)=2*c*y(1)*y(2)
//
// Each request for a gradient evaluation uses the
// functionality of numerical evaluation, but with
// the same numerical code.
iopt[0] = NumericalDerivatives.CENTRAL;
iopt[1] = NumericalDerivatives.CENTRAL;
// Set the increment used at the default value.
// Set defaults for increments and scaling:
fac[0] = 1.4901161193847656E-8;
fac[1] = 1.4901161193847656E-8;
// Change scale to account for large value of a.
switch (varIndex) {
case 1:
scale[0] = 1.e0;
scale[1] = 8.e8;
break;
case 2:
scale[0] = 1.e4;
scale[1] = 8.e8;
break;
}
NumericalDerivatives.Function fcn = new NumericalDerivatives.Function() {
public double[] f(int varIndex, double[] y) {
return new double[m];
}
};
InnerNumericalDerivatives derv = new InnerNumericalDerivatives(fcn);
derv.setPercentageFactor(fac);
derv.setDifferencingMethods(iopt);
derv.setScalingFactors(scale);
derv.setInitialF(valueF);
double[][] fjac = derv.evaluateJ(y);
// Since the function is never evaluated at the
// initial point, hold back until the request is made.
// Copy gradient value into array expected by
// outer loop computing the Hessian matrix.
double[] tmp = new double[n];
for (int i = 0; i < n; i++) {
tmp[i] = fjac[0][i];
}
return tmp;
}
public static void main(String args[]) {
int[] iopth = new int[n];
double u;
double[] fach = new double[n];
double[] stateh = new double[n];
double[] scaleh = new double[n];
double[][] actual = new double[n][n];
double[] y = new double[n];
// Define data and point of evaluation:
a = 2.5e6;
b = 3.4e0;
c = 4.5e0;
y[0] = 2.1e0;
y[1] = 3.2e0;
// Machine precision, for measuring errors
u = 2.220446049250313e-016;
// Compute expected relative error using two applications
// of central differences.
v = Math.pow(3.e0 * u, 2.e0 / 3.e0);
v = Math.pow(3 * v, 2. / 3.);
// Set increments and scaling:
fach[0] = 1.4901161193847656E-8;
fach[1] = 1.4901161193847656E-8;
iopth[0] = NumericalDerivatives.CENTRAL;
iopth[1] = NumericalDerivatives.CENTRAL;
// Compute true values of partials.
actual[0][0] = a * b * b * Math.exp(b * y[0]);
actual[1][0] = 2 * c * y[1];
actual[0][1] = 2 * c * y[1];
actual[1][1] = 2 * c * y[0];
// Set the increment used at the default value.
scaleh[0] = 1;
scaleh[1] = 8.e5;
NumericalDerivatives.Function fcn = new NumericalDerivatives.Function() {
public double[] f(int varIndex, double[] y) {
return new double[n];
}
};
NumericalDerivativesEx5 derv2 = new NumericalDerivativesEx5(fcn);
derv2.setPercentageFactor(fach);
derv2.setDifferencingMethods(iopth);
derv2.setScalingFactors(scaleh);
derv2.setInitialF(stateh);
double[][] h = derv2.evaluateJ(y);
NumberFormat nf = NumberFormat.getInstance();
nf.setMaximumFractionDigits(2);
nf.setMinimumFractionDigits(2);
PrintMatrixFormat pmf = new PrintMatrixFormat();
pmf.setNumberFormat(nf);
new PrintMatrix("Numerical Hessian:").print(pmf, h);
new PrintMatrix("Analytic Hessian:").print(pmf, actual);
// Since the function is never evaluated at the
// initial point, hold back until the request is made.
// Subtract the actual hessian matrix values and check.
h[0][0] = (h[0][0] - actual[0][0]) / h[0][0] / v;
h[1][0] = (h[1][0] - actual[1][0]) / h[1][0] / v;
h[0][1] = (h[0][1] - actual[0][1]) / h[0][1] / v;
h[1][1] = (h[1][1] - actual[1][1]) / h[1][1] / v;
new PrintMatrix("Hessian Matrix, Expected Normalized " +
"Relative Error, |all entries|").print(h);
}
}
Output
Numerical Hessian:
0 1
0 36,455,292,905.82 28.80
1 28.80 18.90
Analytic Hessian:
0 1
0 36,455,280,444.94 28.80
1 28.80 18.90
Hessian Matrix, Expected Normalized Relative Error, |all entries|
0 1
0 0.914 0.037
1 0.036 0
Link to Java source.