Cálculo del hessiano en Python mediante diferencias finitas

Question

Estoy calculando el hessiano de un campo escalar y lo he intentado con numdifftools. Esto parece funcionar, pero era bastante lento, así que escribí mi propio enfoque utilizando diferencias finitas.

Aquí está mi código para el 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

He probado ambos en una función de prueba utilizando el siguiente código:

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

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

que devuelve

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

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

Así que para mi función de prueba, parece dar el resultado correcto. Sin embargo, si intento hacerlo con la función real para la que estaba calculando el hessiano, no obtengo los mismos resultados. La función para la que estoy calculando el hessiano es un campo escalar:

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.])

y da lo siguiente:

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]])

Curiosamente, sólo el primer elemento (1,1) de mi matriz hessiana es el mismo. Los demás elementos difieren bastante. ¿Podría alguien ayudarme a entender de dónde puede venir esto?

Gracias.

Answer 1

¿Puedo recomendar un enfoque ligeramente diferente? Sigue siendo un enfoque numérico, aunque no basado en diferencias finitas.

Puede utilizar la diferenciación automática para calcular los hessianes. Compruebe el autograd en Python. El hessiano puede calcularse como el jacobiano del gradiente utilizando el siguiente fragmento:

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))

El resultado será

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