Take the 2-minute tour ×

Is there a way, given a set of values (x,f(x)), to find the polynomial of a given degree that best fits the data?

I know polynomial interpolation, which is for finding a polynomial of degree n given n+1 data points, but here there are a large number of values and we want to find a low-degree polynomial (find best linear fit, best quadratic, best cubic, etc.). It might be related to least squares...

More generally, I would like to know the answer when we have a multivariate function -- points like (x,y,f(x,y)), say -- and want to find the best polynomial (p(x,y)) of a given degree in the variables. (Specifically a polynomial, not splines or Fourier series.)

Both theory and code/libraries (preferably in Python, but any language is okay) would be useful.

share|improve this question

14 Answers 14

up vote 37 down vote accepted

Thanks for everyone's replies. Here is another attempt at summarizing them. Pardon if I say too many "obvious" things: I knew nothing about least squares before, so everything was new to me.

NOT polynomial interpolation

Polynomial interpolation is fitting a polynomial of degree n given n+1 data points, e.g. finding a cubic that passes exactly through four given points. As said in the question, this was not want I wanted—I had a lot of points and wanted a small-degree polynomial (which will only approximately fit, unless we've been lucky)—but since some of the answers insisted on talking about it, I should mention them :) Lagrange polynomial, Vandermonde matrix, etc.

What is least-squares?

"Least squares" is a particular definition/criterion/"metric" of "how well" a polynomial fits. (There are others, but this is simplest.) Say you are trying to fit a polynomial p(x,y) = a + bx + cy + dx2 + ey2 + fxy to some given data points (xi,yi,Zi) (where "Zi" was "f(xi,yi)" in the question). With least-squares the problem is to find the "best" coefficients (a,b,c,d,e,f), such that what is minimized (kept "least") is the "sum of squared residuals", namely

S = ∑i (a + bxi + cyi + dxi2 + eyi2 + fxiyi - Zi)2


The important idea is that if you look at S as a function of (a,b,c,d,e,f), then S is minimized at a point at which its gradient is 0. This means that for example ∂S/∂f=0, i.e. that

i2(a + … + fxiyi - Zi)xiyi = 0

and similar equations for a, b, c, d, e. Note that these are just linear equations in a…f. So we can solve them with Gaussian elimination or any of the usual methods.

This is still called "linear least squares", because although the function we wanted was a quadratic polynomial, it is still linear in the parameters (a,b,c,d,e,f). Note that the same thing works when we want p(x,y) to be any "linear combination" of arbitrary functions fj, instead of just a polynomial (= "linear combination of monomials").


For the univariate case (when there is only variable x — the fj are monomials xj), there is Numpy's polyfit:

>>> import numpy
>>> xs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
>>> ys = [1.1, 3.9, 11.2, 21.5, 34.8, 51, 70.2, 92.3, 117.4, 145.5]
>>> p = numpy.poly1d(numpy.polyfit(xs, ys, deg=2))
>>> print p
1.517 x + 2.483 x + 0.4927

For the multivariate case, or linear least squares in general, there is SciPy. As explained in its documentation, it takes a matrix A of the values fj(xi). (The theory is that it finds the Moore-Penrose pseudoinverse of A.) With our above example involving (xi,yi,Zi), fitting a polynomial means the fj are the monomials x()y(). The following finds the best quadratic (or best polynomial of any other degree, if you change the "degree = 2" line):

from scipy import linalg
import random

n = 20
x = [100*random.random() for i in range(n)]
y = [100*random.random() for i in range(n)]
Z = [(x[i]+y[i])**2 + 0.01*random.random() for i in range(n)]

degree = 2
A = []
for i in range(n):
    for xd in range(degree+1):
        for yd in range(degree+1-xd):
            A[i].append((x[i]**xd)*(y[i]**yd)) #f_j(x_i)

c,_,_,_ = linalg.lstsq(A,Z)
j = 0
for xd in range(0,degree+1):
    for yd in range(0,degree+1-xd):
        print " + (%.2f)x^%dy^%d" % (c[j], xd, yd),
        j += 1


 + (0.01)x^0y^0  + (-0.00)x^0y^1  + (1.00)x^0y^2  + (-0.00)x^1y^0  + (2.00)x^1y^1  + (1.00)x^2y^0

so it has discovered that the polynomial is x2+2xy+y2+0.01. [The last term is sometimes -0.01 and sometimes 0, which is to be expected because of the random noise we added.]

Alternatives to Python+Numpy/Scipy are R and Computer Algebra Systems: Sage, Mathematica, Matlab, Maple. Even Excel might be able to do it. Numerical Recipes discusses methods to implement it ourselves (in C, Fortran).


  • It is strongly influenced by how the points are chosen. When I had x=y=range(20) instead of the random points, it always produced 1.33x2+1.33xy+1.33y2, which was puzzling... until I realised that because I always had x[i]=y[i], the polynomials were the same: x2+2xy+y2 = 4x2 = (4/3)(x2+xy+y2). So the moral is that it is important to choose the points carefully to get the "right" polynomial. (If you can chose, you should choose Chebyshev nodes for polynomial interpolation; not sure if the same is true for least squares as well.)
  • Overfitting: higher-degree polynomials can always fit the data better. If you change the degree to 3 or 4 or 5, it still mostly recognizes the same quadratic polynomial (coefficients are 0 for higher-degree terms) but for larger degrees, it starts fitting higher-degree polynomials. But even with degree 6, taking larger n (more data points instead of 20, say 200) still fits the quadratic polynomial. So the moral is to avoid overfitting, for which it might help to take as many data points as possible.
  • There might be issues of numerical stability I don't fully understand.
  • If you don't need a polynomial, you can obtain better fits with other kinds of functions, e.g. splines (piecewise polynomials).
share|improve this answer
Chebyshev nodes work for least squares. –  Jason S Dec 31 '08 at 23:57
@Jason:Are you sure that Chebyshev nodes are known to be the best points to choose even for least squares? [It seems there is the different question of picking Chebyshev polynomials themselves as the approximating polynomials, for a different kind of fit from least squares -- "minimax" polynomial.] –  ShreevatsaR Jan 1 '09 at 0:44
df/dx = 0 does not necessarily mean that f is minimized, it might also be maximized. –  quant_dev Feb 20 '09 at 21:11
True, but that's not what we're saying here. We're saying that (assuming the function has partial derivatives etc.), any minimum occurs (either on the boundary or) at a point where the gradient is 0. –  ShreevatsaR Feb 20 '09 at 21:34
about your concerns on numerical stability: defining a polynomial (= "linear combination of monomials") is a dangerous thing to do, because (saying it in non mathematical words) monomials (say) above degree 4 are terribly similar to each other in the area around 0 and then they just "grow crazy". a much better thing is to decide in which interval you are attempting to fit the data, redefine your independent variable so that you are actually fitting in (-1, 1), and look for a linear combination of good polynomials instead of monomials. I would use the Chebyshev set. –  mariotomo Jul 27 '09 at 17:44

Yes, the way this is typically done is by using least squares. There are other ways of specifying how well a polynomial fits, but the theory is simplest for least squares. The general theory is called linear regression.

Your best bet is probably to start with Numerical Recipes.

R is free and will do everything you want and more, but it has a big learning curve.

If you have access to Mathematica, you can use the Fit function to do a least squares fit. I imagine Matlab and its open source counterpart Octave have a similar function.

share|improve this answer
This is helpful, but do you know which of these do multivariate fitting? –  ShreevatsaR Dec 19 '08 at 22:11
least squares can be multivariate. gil strang's "intro to applied mathematics" has a very nice, readable discussion. –  duffymo Dec 28 '08 at 15:55
Yes, thanks... when I asked that comment about multivariate fitting I didn't know enough about least squares :-) –  ShreevatsaR Dec 31 '08 at 23:14

For (x, f(x)) case:

import numpy

x = numpy.arange(10)
y = x**2

coeffs = numpy.polyfit(x, y, deg=2)
poly = numpy.poly1d(coeffs)
print poly
yp = numpy.polyval(poly, x)
print (yp-y)
share|improve this answer

Lagrange polynomials (as @j w posted) give you an exact fit at the points you specify, but with polynomials of degree more than say 5 or 6 you can run into numerical instability.

Least squares gives you the "best fit" polynomial with error defined as the sum of squares of the individual errors. (take the distance along the y-axis between the points you have and the function that results, square them, and sum them up) The MATLAB polyfit function does this, and with multiple return arguments, you can have it automatically take care of scaling/offset issues (e.g. if you have 100 points all between x=312.1 and 312.3, and you want a 6th degree polynomial, you're going to want to calculate u = (x-312.2)/0.1 so the u-values are distributed between -1 and +=).

NOTE that the results of least-squares fits are strongly influenced by the distribution of x-axis values. If the x-values are equally spaced, then you'll get larger errors at the ends. If you have a case where you can choose the x values and you care about the maximum deviation from your known function and an interpolating polynomial, then the use of Chebyshev polynomials will give you something that is close to the perfect minimax polynomial (which is very hard to calculate). This is discussed at some length in Numerical Recipes.

Edit: From what I gather, this all works well for functions of one variable. For multivariate functions it is likely to be much more difficult if the degree is more than, say, 2. I did find a reference on Google Books.

share|improve this answer
Thanks for the reference, BTW. The relevant stuff is a few pages later; 4.10.4 on page 231. The same thing works for multivariate polynomials of higher degree too without great difficulty, although there is the concern of "overfitting". –  ShreevatsaR Dec 21 '08 at 16:17

Bare in mind that a polynomial of higher degree ALWAYS fits the data better. Polynomials of higher degree typically leads to highly improbable functions (see Occam's Razor), though (overfitting). You want to find a balance between simplicity (degree of polynomial) and fit (e.g. least square error). Quantitatively, there are tests for this, the Akaike Information Criterion or the Bayesian Information Criterion. These tests give a score which model is to be prefered.

share|improve this answer
Right, I realised after a while that there was some trade-off between fit and simplicity, as you said. Thanks for the information about the criteria. –  ShreevatsaR Dec 20 '08 at 16:39

at college we had this book which I still find extremely useful: Conte, de Boor; elementary numerical analysis; Mc Grow Hill. The relevant paragraph is 6.2: Data Fitting.
example code comes in FORTRAN, and the listings are not very readable either, but the explanations are deep and clear at the same time. you end up understanding what you are doing, not just doing it (as is my experience of Numerical Recipes).
I usually start with Numerical Recipes but for things like this I quickly have to grab Conte-de Boor.

maybe better posting some code... it's a bit stripped down, but the most relevant parts are there. it relies on numpy, obviously!

def Tn(n, x):
  if n==0:
    return 1.0
  elif n==1:
    return float(x)
    return (2.0 * x * Tn(n - 1, x)) - Tn(n - 2, x)

class ChebyshevFit:

  def __init__(self):
    self.Tn = Memoize(Tn)

  def fit(self, data, degree=None):
    """fit the data by a 'minimal squares' linear combination of chebyshev polinomials.

    cfr: Conte, de Boor; elementary numerical analysis; Mc Grow Hill (6.2: Data Fitting)

    if degree is None:
      degree = 5

    data = sorted(data)
    self.range = start, end = (min(data)[0], max(data)[0])
    self.halfwidth = (end - start) / 2.0
    vec_x = [(x - start - self.halfwidth)/self.halfwidth for (x, y) in data]
    vec_f = [y for (x, y) in data]

    mat_phi = [numpy.array([self.Tn(i, x) for x in vec_x]) for i in range(degree+1)]
    mat_A = numpy.inner(mat_phi, mat_phi)
    vec_b = numpy.inner(vec_f, mat_phi)

    self.coefficients = numpy.linalg.solve(mat_A, vec_b)
    self.degree = degree

  def evaluate(self, x):
    """use Clenshaw algorithm


    x = (x-self.range[0]-self.halfwidth) / self.halfwidth

    b_2 = float(self.coefficients[self.degree])
    b_1 = 2 * x * b_2 + float(self.coefficients[self.degree - 1])

    for i in range(2, self.degree):
      b_1, b_2 = 2.0 * x * b_1 + self.coefficients[self.degree - i] - b_2, b_1
      b_0 = x*b_1 + self.coefficients[0] - b_2

    return b_0
share|improve this answer
Thanks again; that's quite clear. Why is it good to normalise the range to (-1,1), BTW? –  ShreevatsaR Jul 28 '09 at 12:52
because that's the range where the chebyshev polynomials are good behaving. in fact inside that range you can characterize them this way: T_n(x) = cos(n*acos(x)). this formula has no meaning for x not in (-1, 1). –  mariotomo Jul 28 '09 at 14:39
I have tested my module against numpy.polyfit (just that single example from the page you point to) and I was a bit surprised seeing that my fit matches the numpy.polyfit (even for extrapolation) up to the 15th digit. I should try more bad conditioned cases... if they still match then maybe numpy uses chebyshev polynomials behind the scenes and returns the corresponding monomial coefficients... –  mariotomo Jul 28 '09 at 14:45

If you want to fit the (xi, f(xi)) to an polynomial of degree n then you would set up a linear least squares problem with the data (1, xi, xi, xi^2, ..., xi^n, f(xi) ). This will return a set of coefficients (c0, c1, ..., cn) so that the best fitting polynomial is *y = c0 + c1 * x + c2 * x^2 + ... + cn * x^n.*

You can generalize this two more than one dependent variable by including powers of y and combinations of x and y in the problem.

share|improve this answer

The lagrange polynomial is in some sense the "simplest" interpolating polynomial that fits a given set of data points.

It is sometimes problematic because it can vary wildly between data points.

share|improve this answer
Also, the Lagrange polynomial has degree n-1 with n points -- that's what I wrote in the question about polynomial interpolation -- it doesn't give the best-fit polynomial for a given degree. –  ShreevatsaR Dec 19 '08 at 21:43

Not an answer--

I've always thought that for game development it would be really nice to have a program that let you draw some arbitrary (x, f(x)) style curve with the mouse, add numbers to the scales, and have it print out a formula to use for f.

Using a pure polynomial solution seems to me to have the problem that--well like I don't think it could come up with sin(x) if you input a sine wave. I could be wrong, I haven't done serious maths in decades.

share|improve this answer
For visual applications, this is usually done with splines -- piecewise polynomials that fit together smoothly -- rather than with polynomials. Polynomials wiggle too much when they fit data. –  John D. Cook Dec 19 '08 at 21:54
WEll, technically you could. Transcendental functions can be represented to any degree of accuracy with Taylor Series. However, it takes a polynomial of infinite degree to exactly match the function, so I'd probably look at sinusoidal methods for this type of problem. –  stalepretzel Dec 21 '08 at 16:36

It's called approximation theory, a vast field in numerical analysis, well researched.

share|improve this answer
Sure it's probably a vast ocean, but it would be nice to know something from it. :-) –  ShreevatsaR Dec 20 '08 at 16:34

Remember, there's a big difference between approximating the polynomial and finding an exact one.

For example, if I give you 4 points, you could

  1. Approximate a line with a method like least squares
  2. Approximate a parabola with a method like least squares
  3. Find an exact cubic function through these four points.

Be sure to select the method that's right for you!

share|improve this answer
Yes, I know :-) That is why I mentioned "polynomial interpolation" in the question, which is for finding an exact cubic through four points, or an exact n-th degree polynomial through n+1 points. –  ShreevatsaR Dec 21 '08 at 17:17

Some theory:

Vandermonde matrix: http://en.wikipedia.org/wiki/Vandermonde_matrix

Moore-Penrose pseudoinverse: http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse

share|improve this answer

It's rather easy to scare up a quick fit using Excel's matrix functions if you know how to represent the least squares problem as a linear algebra problem. (That depends on how reliable you think Excel is as a linear algebra solver.)

share|improve this answer

Here's a VB6 code library that might be of use.

VB6 Polynomial Regression

share|improve this answer

Your Answer


