How to fit a polynomial with some of the coefficients constrained?

Question

Using NumPy's polyfit (or something similar) is there an easy way to get a solution where one or more of the coefficients are constrained to a specific value?

For example, we could find the ordinary polynomial fitting using:

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])
z = np.polyfit(x, y, 3)

yielding

array([ 0.08703704, -0.81349206,  1.69312169, -0.03968254])

But what if I wanted the best fit polynomial where the third coefficient (in the above case z[2]) was required to be 1? Or will I need to write the fitting from scratch?

Answer 1

In this case, I would use curve_fit or lmfit; I quickly show it for the first one.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def func(x, a, b, c, d):
  return a + b * x + c * x ** 2 + d * x ** 3

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

print(np.polyfit(x, y, 3))

popt, _ = curve_fit(func, x, y)
print(popt)

popt_cons, _ = curve_fit(func, x, y, bounds=([-np.inf, 2, -np.inf, -np.inf], [np.inf, 2.001, np.inf, np.inf]))
print(popt_cons)

xnew = np.linspace(x[0], x[-1], 1000)

plt.plot(x, y, 'bo')
plt.plot(xnew, func(xnew, *popt), 'k-')
plt.plot(xnew, func(xnew, *popt_cons), 'r-')
plt.show()

This will print:

[ 0.08703704 -0.81349206  1.69312169 -0.03968254]
[-0.03968254  1.69312169 -0.81349206  0.08703704]
[-0.14331349  2.         -0.95913556  0.10494372]

So in the unconstrained case, polyfit and curve_fit give identical results (just the order is different), in the constrained case, the fixed parameter is 2, as desired.

The plot looks then as follows:

In lmfit you can also choose whether a parameter should be fitted or not, so you can then also just set it to a desired value (check this answer).

Answer 2

For completeness, with lmfit the solution would look like this:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

def func(x, a, b, c, d):
    return a + b * x + c * x ** 2 + d * x ** 3

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

pmodel = Model(func)
params = pmodel.make_params(a=1, b=2, c=1, d=1)

params['b'].vary = False 

result = pmodel.fit(y, params, x=x)

print(result.fit_report())

xnew = np.linspace(x[0], x[-1], 1000)
ynew = result.eval(x=xnew)

plt.plot(x, y, 'bo')
plt.plot(x, result.best_fit, 'k-')
plt.plot(xnew, ynew, 'r-')
plt.show()

which would print a comprehensive report, including uncertainties, correlations and fit statistics as:

[[Model]]
    Model(func)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 10
    # data points      = 6
    # variables        = 3
    chi-square         = 0.066
    reduced chi-square = 0.022
    Akaike info crit   = -21.089
    Bayesian info crit = -21.714
[[Variables]]
    a:  -0.14331348 +/- 0.109441 (76.37%) (init= 1)
    b:   2 (fixed)
    c:  -0.95913555 +/- 0.041516 (4.33%) (init= 1)
    d:   0.10494371 +/- 0.008231 (7.84%) (init= 1)
[[Correlations]] (unreported correlations are <  0.100)
    C(c, d)                      = -0.987
    C(a, c)                      = -0.695
    C(a, d)                      =  0.610

and produce a plot of

Note that lmfit.Model has many improvements over curve_fit, including automatically naming parameters based on function arguments, allowing any parameter to have bounds or simply be fixed without requiring nonsense like having upper and lower bounds that are almost equal. The key is that lmfit uses Parameter objects that have attributes instead of plain arrays of fitting variables. lmfit also supports mathematical constraints, composite models (eg, adding or multiplying models), and has superior reports.

Answer 3

Here is a way to do this using scipy.optimize.curve_fit:

First, let's recreate your example (as a sanity check):

import numpy as np
from scipy.optimize import curve_fit
​
def f(x, x3, x2, x1, x0):
    """this is the polynomial function"""
    return x0 + x1*x + x2*(x*x) + x3*(x*x*x)
​
popt, pcov = curve_fit(f, x, y)

print(popt)
#array([ 0.08703704, -0.81349206,  1.69312169, -0.03968254])

Which matches the values you get from np.polyfit().

Now adding the constraints for x1:

popt, pcov = curve_fit(
    f, 
    x,
    y,
    bounds = ([-np.inf, -np.inf, .999999999, -np.inf], [np.inf, np.inf, 1.0, np.inf])
)
print(popt)
#array([ 0.04659264, -0.48453866,  1.        ,  0.19438046])

I had to use .999999999 because the lower bound must be strictly less than the upper bound.

Alternatively, you could define your function with the constrained coefficient as a constant, and get the values for the other 3:

def f_new(x, x3, x2, x0):
    x1 = 1
    return x0 + x1*x + x2*(x*x) + x3*(x*x*x)
popt, pcov = curve_fit(f_new, x, y)
print(popt)
#array([ 0.04659264, -0.48453866,  0.19438046])

Answer 4

Sorry for the resurrection

..but I felt that this answer was missing.

To fit a polynomial we solve the following system of equations:

a0*x0^n + a1*x0^(n-1) .. + an*x0^0 = y0
a0*x1^n + a1*x1^(n-1) .. + an*x1^0 = y1
                 ...
a0*xm^n + a1*xm^(n-1) .. + an*xm^0 = ym

Which is a problem of the form V @ a = y

where "V" is a Vandermonde matrix:

[[x0^n  x0^(n-1)  1],
 [x1^n  x1^(n-1)  1],
        ...
 [xm^n  xm^(n-1)  1]]

"y" is a column vector holding the y-values:

[[y0],
 [y1],
  ...
 [ym]]

..and "a" is the column vector of coefficients that we are solving for:

[[a0],
 [a1],
  ...
 [an]]

This problem can be solved using linear least squares as follows:

import numpy as np

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

deg = 3
V = np.vander(x, deg + 1)
z, *_ = np.linalg.lstsq(V, y, rcond=None)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

..which produces the same solution as the polyfit method:

z = np.polyfit(x, y, deg)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

Instead we want a solution where a2 = 1

substituting a2 = 1 into the system of equations from the beginning of the answer, and then moving the corresponding term from the lhs to the rhs we get:

a0*x0^n + a1*x0^(n-1) + 1*x0^(n-2) .. + an*x0^0 = y0
a0*x1^n + a1*x1^(n-1) + 1*x0^(n-2) .. + an*x1^0 = y1
                 ...
a0*xm^n + a1*xm^(n-1) + 1*x0^(n-2) .. + an*xm^0 = ym

=>

a0*x0^n + a1*x0^(n-1) .. + an*x0^0 = y0 - 1*x0^(n-2)
a0*x1^n + a1*x1^(n-1) .. + an*x1^0 = y1 - 1*x0^(n-2)
                 ...
a0*xm^n + a1*xm^(n-1) .. + an*xm^0 = ym - 1*x0^(n-2)

This corresponds to removing column 2 from the Vandermonde matrix and subtracting it from the y-vector as follows:

y_ = y - V[:, 2]
V_ = np.delete(V, 2, axis=1)
z_, *_ = np.linalg.lstsq(V_, y_, rcond=None)
z_ = np.insert(z_, 2, 1)

print(z_)
# [ 0.04659264 -0.48453866  1.          0.19438046]

Notice that I inserted the 1 in the coefficient vector after solving the linear least-squares problem, we are no longer solving for a2 since we set it to 1 and removed it from the problem.

For completeness this is what the solution looks like when plotted:

and the complete code that I used:

import numpy as np

x = np.array([0.0, 1.0, 2.0, 3.0,  4.0,  5.0])
y = np.array([0.0, 0.8, 0.9, 0.1, -0.8, -1.0])

deg = 3
V = np.vander(x, deg + 1)
z, *_ = np.linalg.lstsq(V, y, rcond=None)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

z = np.polyfit(x, y, deg)

print(z)
# [ 0.08703704 -0.81349206  1.69312169 -0.03968254]

y_ = y - V[:, 2]
V_ = np.delete(V, 2, axis=1)
z_, *_ = np.linalg.lstsq(V_, y_, rcond=None)
z_ = np.insert(z_, 2, 1)

print(z_)
# [ 0.04659264 -0.48453866  1.          0.19438046]

from matplotlib import pyplot as plt

plt.plot(x, y, 'o', label='data')
plt.plot(x, V @ z, label='polyfit')
plt.plot(x, V @ z_, label='constrained (a2 = 0)')

plt.legend()

plt.show()