1

I faced an issue with Numba's njit tool in Python. I noticed that the function gives different results when run with @numba.njit and when run as plain Python code. Particularly, after debugging, I noticed that the discrepancy in calculation occurs when performing matrix inversion using numpy. Please see my test code below. The value for matrix A and vector b are in the following csv files that can be accessed via the following links: A.csv and b.csv

The results from the plain Python function is the correct one. Please help me solve this issue! Do I need use a Numba wrapper function around the numpy matrix inversion function to resolve what seems to be a numerical issue?

Kind regrads, and I look forward to hearing from you guys soon :)

Ahmad

@numba.njit
def cal_Test_jit(A,b):
    c = np.linalg.inv(A)@b
    return c, np.linalg.inv(A)

def cal_Test(A,b):
    c = np.linalg.inv(A)@b
    return c, np.linalg.inv(A)

A = np.loadtxt(open("A.csv", "rb"), delimiter=",")
b = np.loadtxt(open("b.csv", "rb"), delimiter=",")

c_jit, Ai_jit = cal_Test_jit(A,b)
c, Ai = cal_Test(A,b)
err_c = abs(c-c_jit)
err_A = abs(Ai_jit-Ai)

# ploting the error in the parameters
plt.figure()
plt.plot(err_c)

# only ploting the error in first three columns of A
fig, ax = plt.subplots(1,3)
ax[0].plot(err_A[:,0])
ax[1].plot(err_A[:,1])
ax[2].plot(err_A[:,2])

1 Answers1

2

One way to use numba in your problem is to add:

@numba.jit(forceobj=True)

which will get the true results but execution time comparison is as using njit (different results) > this method (exact results) == plain Python e.g. using colab TPU:

1000 loops, best of 5: 545 µs per loop     # using njit
1000 loops, best of 5: 505 µs per loop     # this method
1000 loops, best of 5: 500 µs per loop     # plain Python

but as it is recommended before on another SO question, Numba is very useful for optimizing a subset of pure Python, especially loops, close to the performance of optimized C code @BatWannaBe and there are almost always better approaches than inverting a matrix, e.g. np.linalg.solve @Humer512; which is well explained on another SO question @ali_m.

Ali_Sh
  • 2,667
  • 3
  • 43
  • 66
  • 1
    Lovely. Thank you so much, @Ali_Sh, for the prompt response. I used np.linalg.solve, and it works. This will get me going. Later, I will take a look at the application the flag forceobj=True and the additional post you referenced. It seems that np.linalg.solve is similar to MATLAB's A\b. Many thanks again, Ali_Sh :) – Ahmad Abuaish Nov 19 '21 at 19:15