Here's the code. 2a works but 2b does not, even though both matrices have the same shape and the same type. Can anybody explain?
import sympy as sym
import numpy as np
x1, x2 = sym.symbols('x1 x2')
f = 2*x1**2 - 2*x1*x2 + x2**2 + 2*x1 - 2*x2
xk = np.array([[0 , 1]])
print("xk shape = ", xk.shape)
fmat = sym.Matrix([f])
H = sym.hessian(fmat, (x1, x2))
print("Symbolic H")
sym.pprint(H)
# 2a: Convert the SymPy Matrix to Numpy array using lambdify and then substitute values
Hnp = sym.lambdify((x1, x2), H, 'numpy')
Hknp = Hnp(xk[0][0], xk[0][1])
print("Hknp")
sym.pprint(Hknp)
print("Hknp type:", type(Hknp))
print("Hknp shape = ", Hknp.shape)
v, Q = np.linalg.eigh(Hknp)
print("v:", v)
print("Q:")
sym.pprint(Q)
# 2b: Substitute values into SymPy Matrix then convert to Numpy array
Hks = H.subs([(x1, xk[0][0]), (x2, xk[0][1])])
Hk = np.array(Hks)
print("Hk")
sym.pprint(Hk)
print("Hk type:", type(Hk))
print("Hk shape = ", Hk.shape)
v, Q = np.linalg.eigh(Hk)
print("v:", v)
print("Q:")
sym.pprint(Q)
>>> Hk
array([[4, -2],
[-2, 2]], dtype=object)
>>> Hknp.dtype
dtype('int64')
The solution is Hk = np.array(Hks, dtype=np.int64)
.