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 f(y_1 ,y_2 ). 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]
Our method is based on casting the matrix Hf 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 4/9 = \left( {2/3} \right)^2 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.