1

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)
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214
Alan M
  • 13
  • 3
  • What do mean by 'does not'? You should show the error! Show the prints as well. Not everyone can run your code to see these things for themselves. – hpaulj May 27 '20 at 06:50
  • Naively passing `sympy` expressions to `numpy/scipy` often gives problems. You have to examine the resulting array carefully before proceeding. `lambdify` is the safeist interface. – hpaulj May 27 '20 at 06:55

1 Answers1

2
>>> Hk
array([[4, -2],
       [-2, 2]], dtype=object)


>>> Hknp.dtype
dtype('int64')

The solution is Hk = np.array(Hks, dtype=np.int64).

orlp
  • 112,504
  • 36
  • 218
  • 315
  • Very useful to know now that dtype is important. In this case using float would be best. If the matrix calculation happens to result in a float but I declare the dtype as int64, Python (as one would expect) truncates to an integer. – Alan M May 28 '20 at 16:33