Fastest way to fit a parabola to set of points?

Given a set of points, what's the fastest way to fit a parabola to them? Is it doing the 开发者_JS百科least squares calculation or is there an iterative way?

Thanks

Edit: I think gradient descent is the way to go. The least squares calculation would have been a little bit more taxing (having to do qr decomposition or something to keep things stable).

If the points have no error associated, you may interpolate by three points. Otherwise least squares or any equivalent formulation is the way to go.

I recently needed to find a parabola that passes through 3 points.

suppose you have (x1,y1), (x2,y2) and (x3,y3) and you want the parabola

y-y0 = a*(x-x0)^2

to pass through them: find y0, x0, and a.

You can do some algebra and get this solution (providing the points aren't all on a line) :

let c = (y1-y2) / (y2-y3)
x0    = ( -x1^2 + x2^2 + c*( x2^2 - x3^2 ) )  /  (2.0*( -x1+x2 + c*x2 - c*x3 ))
a     = (y1-y2)  /  ( (x1-x0)^2 - (x2-x0)^2 )
y0    = y1 - a*(x1-x0)^2

Note in the equation for c if y2==y3 then you've got a problem. So in my algorithm I check for this and swap say x1, y1 with x2, y2 and then proceed.

hope that helps!

Paul Probert

A calculated solution is almost always faster than an iterative solution. The "exception" would be for low iteration counts and complex calculations.

I would use the least squares method. I've only every coded it for linear regression fits but it can be used for parabolas (I had reason to look it up recently - sources included an old edition of "Numerical Recipes" Press et al; and "Engineering Mathematics" Kreyzig).

ALGORITHM FOR PARABOLA

1. Read no. of data points n and order of polynomial Mp .
2. Read data values .
3. If n< Mp [ Regression is not possible ] stop else continue ;
4. Set M=Mp + 1 ;
5. Compute co-efficient of C-matrix .
6. Compute co-efficient of B-matrix .
7. Solve for the co-efficients a1,a2,. . . . . . . an .
8. Write the co-efficient .
9. Estimate the function value at the glren of independents variables .

Using the free arbitrary accuracy math program "PARI" (for Mac or PC): Here is how I would fit a parabola to a set of 641 points, and I also show how to find the minimum of that parabola:

0. Set a high number of digits of precision:

   \p 300
   

1. Write the data points to a text file separated by one space for each data point
   (use ASCII characters in base ten, no space at file start or file end, and no returns, write extremely large or small floating points as for example "9.0E-23" but not "9.0D-23" ).

2. make a string to point to that file:

   fileone="./desktop/data.txt"
   

3. read that file into PARI using the following instructions:

   fileopen(fileone,r)
   readsplit(file) = my(cmd);cmd="perl -ne \"chomp; print '[' . join(',', split(/ +/)) . ']\n';\"";eval(externstr(Str(cmd," ",file)))
   readsplit(fileone) 
   

4. Label that data with a name:

   in = %
   V = in[1]
   

5. Define a least squares fit function:

   lsf(X,Y,n) = my(M=matrix(#X,n+1,i,j,X[i]^(j-1)));fit=Polrev(matsolve(M~*M,M~*Y~))
   

6. Apply that lsf function to your 641 data points:

   lsf([-320..320],V, 2)
   

7. Then if you want to show the minimum of that parabolic fit, enter:

   xextreme = solve (x=-1000,1000,eval(deriv(fit)));print (xextreme*(124.5678-123.5678)/640+(124.5678+123.5678)/2);x=xextreme;print(eval(fit))
   

(I had to adjust for my particular x-axis scaling before the "print" statement in that command line above).

(Note: A sacrifice made to simplify this algorithm causes it to work only when the data set has equally spaced x-axis coordinates.)

I was worried that my last post was too compact to follow and too hard to convert to other environments. I would like to show here how to solve the generalized problem of parabolic data fitting explicitly without specialized matrix math terminology; and so that each multiplication, division, subtraction and addition can be seen at once. To save ink this fit reparameterizes the x-axis as evenly spaced points centered on zero so that odd powered sums all get eliminated (saving a lot of space and time), so the x-coordinates of the N data points are effectively labeled by points of this vector: X=[-(N-1)/2..(N-1)/2]. For example "xextreme" will be returned versus those integer indices and so (if desired) a simple (consumes very little CPU time) linear transformation must be applied after the algorithm below to get it versus your problem's particular x-axis labels.

This is written in the language of the free program "PARI" but all the commands are simple to translate to any language.

Step 1: assign a label to the y-axis data:

? V=[5,2,1,2,5]

"PARI" confirms that entry:

%280 = [5, 2, 1, 2, 5]

Then type in the following processing algorithm which calculates a best fit parabola through any y-axis data set with constant x-axis separation:

? g=#V;h=(g-1)*g*(g+1)/3;i=h*(3*g*g-7)/5;\
a=sum(i=1,g,V[i]);b=sum(i=1,g,(2*i-1-g)*V[i]);c=sum(i=1,g,(2*i-1-g)*(2*i-1-g)*V[i]);\
A=matdet([a,c;h,i])/matdet([g,h;h,i]);B=b/h*2;C=matdet([g,h;a,c])/matdet([g,h;h,i])*4;\
xextreme=-B/(2*C);yextreme=-B*B/(4*C)+A;fit=Polrev([A,B,C]);\
print("\n","y of extreme is ",yextreme,"\n","which occurs this many data points from center of data: ",xextreme)                                                                                       
                                                                          
(Note for non-PARI users:
the command "matdet([a,c;h,i])" 
is just another way of entering "a*i-c*h") 

Those commands then produce the following screen output:

y of extreme is 1 which occurs this many data points from center of data: 0

The algorithm stores the polynomial of the fit in the variable "fit":

? fit
%282 = x^2 + 1
?

(Note that to make that algorithm short the x-axis labels are assigned as X=[-(N-1)/2..(N-1)/2], thus they are X=[-2,-1,0,1,2] To correct that for the same polynomial as parameterized by an x-axis coordinate data set of say X=[−1,0,1,2,3]: just apply a simple linear transform, in this case: "x^2 + 1" --> "(t - 1)^2 + 1".)