I'm using NumPy to find out langrage polynomial interpolation. I'm using 2 for loop to find out langrage polynomial, but I want to reduce 2nd for loop so that my code time complexity can be less. Can you guys please suggest me how I can cut a for loop to increase my time complexity.

import numpy as np
import matplotlib.pyplot as plt
class LagrangePolynomial:
    def __init__(self, data_x, data_y):

        assert type(data_x) == np.ndarray, "data_x is not numpy array"
        assert type(data_y) == np.ndarray, "data_y is not numpy array"
        assert len(data_x) == len(data_y), "length of data_x and data_y must be equal"

        self.x = data_x
        self.y = data_y

        self.degree = len(data_x) - 1

    def __str__(self):
        strL = f"LagrangePolynomial of order {self.degree}\n"
        strL += "p(x) = "
        for i in self.y:
            if i == 0:
            elif i >= 0:
                strL += f"+ {i}*l_{count + 1}(x) "
                strL += f"- {-i}*l_{count + 1}(x) "
        return strL
    def __call__(self,x):
        y_interp = np.zeros(len(x))
        for i in range(len(x)):
            #remove this for loop.
            for j in range(self.degree + 1):
               upper=np.prod(np.delete(self.x - x[i],j))
               lower=np.prod(np.delete(self.x[j] - self.x, j))
               y_interp[i] +=(upper/lower)*self.y[j]
        return y_interp
data_x = np.array([-3.,-2.,-1.,0.,1.,3.,4.])
data_y = np.array([-60.,-80.,6.,1.,45.,30.,16.])

p = LagrangePolynomial(data_x, data_y)

#generating 100 points from -3 to 4 in order to create a smooth line
x = np.linspace(-3, 4, 100)
y_interp = p(x)
plt.plot(x, y_interp, 'g--')
plt.plot(data_x, data_y, 'ro')
plt.legend(['interpolated', 'node points'], loc = 'lower right')
plt.title('Lagrange Polynomial')


The second loop can be eliminated by creating a 2D array by tiling the self.x array with the np.tile method. The removal of elements in de 1D array with np.delete can than be replaced with a boolean mask (unfortunately this results in a 1D array which needs to be reshaped).

The calculation of the 2D array can be done outside the first loop, as can be the calculation of lower as they do not depend on x.

See my rewrite of the __call__ method:

def __call__(self,x):

    y_interp = np.zeros(len(x))
    n = self.degree + 1
    mask =  ~np.eye(n, dtype=np.bool)

    tiled = np.tile(self.x, (n,1))
    masked = np.reshape(tiled[mask], (n, n-1))

    lower = np.prod(self.x[:,np.newaxis] - masked, axis=1)
    for i in range(len(x)):

        upper = np.prod(masked -x[i], axis=1)
        y_interp[i] = np.sum((upper/lower) * self.y)
    return y_interp
  • \$\begingroup\$ thanks a lot for your code. \$\endgroup\$ – Oronno Akash Aug 18 '20 at 7:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.