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 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**2 + x**2) out1 = hessianComp(testfunction, np.array([2.,2.])) out2 = numdifftools.Hessian(testfunction)([2., 2.])
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, x, 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?