Elementary linear operations on two polynomials

Hypothetical scenario: An engineer (or anyone who uses this code for that matter!) fits a response surface to some data on day 1. They repeat this exercise on day 2, where the data is moderately different. Can we identify sensitivities, moments, and subspaces on the difference between the two responses.

Rationale: The engineer wants to understand what changed between day 1 and day 2. Now the simplest solution would be to create another response surface on the difference between the two datasets. But what if the input samples where the same, the basis was the same, but only the output samples were different?

The above, and a few other musings, motivate this post.

Consider two polynomials p_1 and p_2. Our goal is to be able to carry out operations of the form

p1 = Poly()
p2 = Poly()
p3 = p1 - p2

where p3 is an instant of the Poly class. There are two important ideas to take note of here. First, any linear operation — whether integration, differentiation, addition, subtraction – on a polynomial will yield another polynomial, which makes this easy. Second, if the support of the two polynomials is different then we most certainly may not have a polynomial. Naturally, one can develop such methods in PolyTrees() more aptly, but for starters let’s consider the simplest scenario.

We add the following code to Poly()


    def __add__(self, *args):
        """ Adds polynomials.

        Returns
        -------

        poly
            An instance of the Poly class.

        """
        polynew = deepcopy(self)
        for other in args:
            polynew.coefficients += other.coefficients
            if np.all( polynew._quadrature_points - other._quadrature_points ) == 0:
                polynew._model_evaluations += other._model_evaluations
        return polynew

    def __sub__(self, *args):
        """ Subtracts polynomials.

        Returns
        -------

        poly
            An instance of the Poly class.

        """
        polynew = deepcopy(self)
        for other in args:
            polynew.coefficients -= other.coefficients
            if np.all( polynew._quadrature_points - other._quadrature_points ) == 0:
                polynew._model_evaluations -= other._model_evaluations
        return polynew

which permits us to do the following:

def fun_1(x):
    return x**2 - 2*x + 3.2

def fun_2(x):
    return x**2

def fun_3(x):
    return -2*x + 10

p = Parameter(lower=-1., upper=1., order=2)
b = Basis('univariate')

# Poly 1
poly1 = Poly(p, b, method='numerical-integration')
poly1.set_model(fun_1)

# Poly 2
poly2 = Poly(p, b, method='numerical-integration')
poly2.set_model(fun_2)

# Poly 3
poly3 = Poly(p, b, method='numerical-integration')
poly3.set_model(fun_3)

# Test!
poly4 = poly1 + poly2
poly5 = poly1 - poly2 - poly3

# Plots!
poly4.plot_polyfit_1D(uncertainty=False)
poly5.plot_polyfit_1D(uncertainty=False)

Note that in practice, one can add or subtract as many polynomials as one desires!

index1

index2

Code commit can be found below. The next step would be to:

  • Have __add__, __sub__, and __mul__ operations on polynomials with different input sample locations.
  • Permit operations over different multi-index sets, provided the support for each parameter is the same.

Addition and subtraction for arbitrary input / output data pairs now complete. See: