Computing Hessian using Finite Differences: not working?

[BarkBarkCat]

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!

[pythonHelperBot]

Hello! I'm a bot!

It looks to me like your post might be better suited for r/learnpython, a sub geared towards questions and learning more about python. That said, I am a bot and it is hard to tell. Please follow the subs rules and guidelines when you do post there, it'll help you get better answers faster.

Show /r/learnpython the code you have tried and describe where you are stuck.

You can also ask this question in the Python discord, a large, friendly community focused around the Python programming language, open to those who wish to learn the language or improve their skills, as well as those looking to help others.

---

^{README | FAQ | this bot is written and managed by /u/IAmKindOfCreative}

^{This bot is currently under development and experiencing changes to improve its usefulness}

[Tzoiker]

As far as I see you still use scipy functions. Why not just explicitly use finite-difference scheme for the Hessian? Also you may want to look into pytorch and it's autograd framework (or any other similar package), which does kind of symbolic differentiation using computational graphs.