Python code to calculate base Lagrange Polynomial is taking too long for large numbers of data-CodePudding

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_poly(p,xhat,n,x,tol):
    matrix=[]
    error_flag=0
    er=1
    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
            for k in range(p 1):
                if k!=j:
                    L= L*(x[i] - xhat[k])/(xhat[j] - xhat[k])
                
            matrix.append(L)
    lagrange_matrix= np.transpose(np.array(matrix).reshape(n,p 1))
    return lagrange_matrix, error_flag


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.
        for j in range(n):
                po=0 
                for i in range(p 1):
                    po  = f(xhat[i]) * lagrange_poly(p,xhat,n,x,1e-10)[0][i,j] 
        inter.append(po) 

        interpolant = np.array(inter) 

                return interpolant

CodePudding user response：

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)[0]
        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.