5
$\begingroup$

I am computing the Hessian of a scalar field, and tried using numdifftools. This seems to work, but was quite slow so I wrote my own approach using finite differences.

Here is my code for the Hessian:

def hessianComp ( func, x0, epsilon=1.e-5):
    f1 = scipy.optimize.approx_fprime( x0, func, epsilon=epsilon) 

    # Allocate space for the hessian
    n = x0.shape[0]
    hessian = np.zeros ( ( n, n ) )
    # The next loop fill in the matrix
    xx = x0
    for j in range( n ):
        xx0 = xx[j] # Store old value
        xx[j] = xx0 + epsilon # Perturb with finite difference
        # Recalculate the partial derivatives for this new point
        f2 = scipy.optimize.approx_fprime( xx, func, epsilon=epsilon) 
        hessian[:, j] = (f2 - f1)/epsilon # scale...
        xx[j] = xx0 # Restore initial value of x0        
    return hessian

I tried both on a test function using the following code:

def testfunction(x):
    return(x[0]**2 + x[1]**2)

out1 = hessianComp(testfunction, np.array([2.,2.]))
out2 = numdifftools.Hessian(testfunction)([2., 2.])

which returns 

out1 = array([[2.00000017, 0.        ],
             [0.        , 2.00000017]])

out2 = array([[2.00000000e+00, 1.04776726e-14],
             [1.04776726e-14, 2.00000000e+00]])

So for my test function, it seems to give the correct result. If, however, I try doing it on the actual function for which I was computing the Hessian I am not getting the same results. The function for which I am computing the Hessian is a scalar field:

def lambda2Field(x):
    out = solve_ivp(doubleGyreVar, t_span=(0, T), y0=[x[0], x[1], 1, 0, 0, 1],t_eval=[0, T], rtol = 1e-10, atol=1e-10)
    output = out.y
    J = output[2:,-1].reshape(2,2,order="F")
    CG = np.matmul(J.T , J)
    lambdas, xis = np.linalg.eig(CG)
    lambda_2 = np.max(np.abs(lambdas))
    return lambda_2

out1 = hessianComp(lambda2Field, np.array([2.,2.]))
out2 = numdifftools.Hessian(lambda2Field)([2., 2.])

and gives the following:

out1 = array([[-1.53769104e+18, -2.20719407e+21],
             [-2.20719407e+21,  1.39720111e+27]])

out2 = array([[-1.53767292e+18,  2.27457455e-07],
             [ 2.27457455e-07, -9.43198781e+16]])

Interestingly, only the first element (1,1) of my Hessian matrix is the same. The other elements are off by quite a bit. Could somebody help me understand where this could be coming from?

Thanks!

$\endgroup$
4
  • $\begingroup$ Remark probably unrelated to error: You could compute the singular values also directly from J with the svd method. // A little more relevant: How sure are you that the maximum singular value is twice differentiable at this point? $\endgroup$
    – Lutz Lehmann
    Jun 11, 2019 at 9:30
  • $\begingroup$ Why don't you compute the Hessian symbolically? $\endgroup$
    – Rodrigo de Azevedo
    Apr 6, 2020 at 7:24
  • $\begingroup$ @BarkingCat Since you are using finite differences to approximate the second derivatives, using the first order derivatives from python, why not approximate the second derivatives directly from function values? For instance, $$ f''_{xx}(x,y) \approx \dfrac{f(x+h,y)-2f(x,y)+f(x-h,y)}{h^2}$$ $$ f''_{yy}(x,y) \approx \dfrac{f(x,y+h)-2f(x,y)+f(x,y-h)}{h^2}$$ $\endgroup$
    – PierreCarre
    Sep 27, 2021 at 10:18
  • $\begingroup$ In my experience, these numerical methods tend to be horribly unstable. The automatic differentiation approach tends to be vastly superior. $\endgroup$
    – Igor Rivin
    Jun 11, 2022 at 1:30

2 Answers 2

Reset to default
3
$\begingroup$

May I recommend using a slightly different approach? It is still a numerical approach, although not based on finite differences.

You can use automatic differentiation to calculate Hessians. Check out the autograd package in Python. The Hessian can be computed as the Jacobian of the gradient using the following snippet:

from autograd import elementwise_grad as egrad
from autograd import jacobian
import autograd.numpy as np


def func(x):
    return np.sin(x[0]) * np.sin(x[1])


x_value = np.array([0.0, 0.0])  # note inputs have to be floats
H_f = jacobian(egrad(func))  # returns a function
print(H_f(x_value))

This will give the output

$$ \begin{bmatrix} 0 & 1 \\ 1 & 0 \\ \end{bmatrix} $$

$\endgroup$
0
$\begingroup$

Use a two-point formula instead to compute the derivative, it’s much more accurate. See numerical differentiation Wikipedia page.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .