# [SOLVED] Python code to calculate base Lagrange Polynomial is taking too long for large numbers of data

## Issue

So, I need help minimizing the time it takes to run the code with large numbers of data only by using NumPy. I think the for loops made my code inefficient.. But I do not know how to make the for loop into a list comprehension, which might help it run faster..

``````def lagrange(p,node,n,x):
m=[]

#base lagrange polynomial
for i in range(n):
for j in range(p+1):
L=1
for k in range(p+1):
if k!=j:
L= L*(x[i] - node[k])/(node[j] - node[k])

m.append(L)
lagrange= np.array(m).reshape(n,p+1)
return lagrange

def interpolant(a,b,p,n,x,f):
m=[]
node=np.linspace(a,b,p+1)
for j in range(n):
polynomial=0
for i in range(p+1):
polynomial += f(node[i]) * lagrange(p,node,n,x)
m.append(polynomial)

interpolant = np.array(inter)

return interpolant

``````

## Solution

It appears the value of `lagrange_poly(...)` is recomputed `n*(p+1)` times for no reason which is very very expensive! You can compute it once before the loop, store it in a variable and reuse the variable later.

Here is the fixed code:

``````def uniform_poly_interpolation(a,b,p,n,x,f,produce_fig):
inter=[]
xhat=np.linspace(a,b,p+1)
#use for loop to iterate interpolant.
mat = lagrange_poly(p,xhat,n,x,1e-10)
for j in range(n):
po=0
for i in range(p+1):
po += f(xhat[i]) * mat[i,j]
inter.append(po)

interpolant = np.array(inter)
return interpolant
``````

This should be much much faster.

Moreover, the execution is slow because accessing scalar values of Numpy arrays from CPython is very slow. Numpy is designed to work with array and not to extract scalar values in loops. Additionally, the loop CPython interpreter are relatively slow. You can solve this problem efficiently with Numba that compile your code to a very fast native code using a JIT-compiler.

Here is the Numba code:

``````import numba as nb

@nb.njit
def lagrange_poly(p, xhat, n, x, tol):
error_flag = 0
er = 1
lagrange_matrix = np.empty((n, p+1), dtype=np.float64)

for l in range(p):
if abs(xhat[l] - xhat[l+1]) < tol:
error_flag = er

# Base lagrange polynomial
for i in range(n):
for j in range(p+1):
L = 1.0
for k in range(p+1):
if k!=j:
L = L * (x[i] - xhat[k]) / (xhat[j] - xhat[k])
lagrange_matrix[i, j] = L

return lagrange_matrix, error_flag
``````

Overall, this should be several order of magnitude faster.