I am trying to write software that does some engineering calculations.
 I currently use Micrsoft Excel but now have to write the new program
independent of Excel. In Excel, I used the function called "LINEST"
 
I don't know the algorithm that LINEST uses and need to come up with
this somehow. I am trying to "derive the polynomial co-efficients from
a set of data points. Depending on the data spread the polynomial may
need to be anything from 3rd to 6th order. Is there a "standard" way
of doing this? As long as it is reasonably accurate we'll be happy.
i.e. given a set of typically 10 points x1, y1   x2,y2  
x3,y3...............x10,y10 we need to calculate the coefficients
a,b,c,d,e,f for each term in the corresponding polynomial ax5 + bx4  +
cx3  + dx2  + ex + f = y?

Ideally the solution should be given in VB.net. I don't want to have
to buy add-in software so need the source code.

---

Request for Question Clarification by mathtalk-ga on 18 Feb 2005 05:43 PST

Hi, rwallacej-ga:

The Excel LINEST function produces a "least squares" fit, i.e. the
polynomial f(x) of specified degree N such that:

  SUM (y_i - f(x_i))^2

is minimized.

Using basic calculus the minimum is determined (if all the x_i are
distinct) by a set of N+1 equations in the N+1 unknown coefficients of
f(x).

These are called the normal equations, and this linear system can be
formulated as a matrix equation with a positive definite symmetric
matrix of coefficients. The simplest case is a constant polynomial,
which gives one equation in one unknown (equivalent to averaging all
the y-values).

If you are limiting the scope of the calculation to polynomials of
degree 6 or less, then there will be a need to solve linear
systems/matrix equations of size 7x7 or smaller.

Since you are building a standalone software package that does
engineering calculations, I wonder if you already have some functions
for solving small linear systems.  As I've implied above, the
coefficient matrix is nonsingular in theory if the x-values are all
distinct.  However least-squares matrices are notorious for having
high "condition numbers", which means that the accuracy of the
solution is sensitive to rounding errors.

So it would be helpful to know a bit more about the libraries and
calculation environment you would already have at hand, in addition to
needing to solve the normal equations for a least squares
approximation.

regards, mathtalk-ga

If the number of data points (observations) is greater than the order
of the polynomial you wish to fit, then you presumably want to use a
least squares method to determine the best estimate of the parameters
in the polynomial equation.  There are zillions of websites out there
that explain the basic idea of least squares.  For example, see:
http://www.efunda.com/math/leastsquares/leastsquares.cfm
http://mathworld.wolfram.com/LeastSquaresFitting.html
http://hilltop.bradley.edu/~jhahn/note7.pdf
and many more

Of course, there are even more of those old standbys, "books" that
cover this topic in infinite detail.  A good, simple introduction is
"Data Reduction and Error ANalysis for the Physical Sciences", 1969 by
P.R. Bevington.

The specific application of least squares to fitting a polynomial is
also discussed in many web sites.  For example, see:
http://www.efunda.com/math/leastsquares/lstsqrmdcurve.cfm
http://mathworld.wolfram.com/LeastSquaresFittingPolynomial.html

Again, there are bazillions of books out there that cover this, many
of which include code for polynomial fitting routines (e.g., the
FORTRAN routine POLFIT in the book by Bevington mentioned above).

Finally, there are lots of commercial outfits out there who will
happily sell you code libraries to do least squares fitting (among
other things).  For instance, http://www.eobj.com/eolsfit-dll.html.

Hello, I solved and programmed this (and similar problems about ten
years ago) in QBasic. I am posting that LeastSquare solution below. I
would have produced a VB.NET solution if I were Google Reasearcher or
had the time. You should be able to do that if familiar with the BASIC
programming language.
--------------------------------------------------------------------------

REM TO RUN THIS PROGRAM, PRESS F5
15 CLS 0: PRINT CHR$(7), "VARIABLE PAIR LEAST SQUARE POLYNOMIAL
REPRESENTATION (MAY '95)"
20 PRINT "DESIGNED BY 'TIMOGOSE'"
22 INPUT "ARE YOU USING A PRINTER (Y/N) ", PTER$
24 IF PTER$ = "Y" THEN PE = 1: LPRINT "PROGRAM 'LEASTSQ' ; DESIGNER :
'TIMOGOSE' (M.Sc)"
30 PRINT "GIVEN P POINTS OF AN ORDERED PAIR (X,Y),THE PROGRAMM GENERATES ANY"
40 PRINT "DESIRED MTH DEGREE POLYNOMIAL Y=M(X)  [M<=P-1] THAT GIVES BEST FIT"
50 INPUT "NUMBER OF POINTS ", P%: P = P%: G = 180: IF P > 10 THEN G = 250
60 IF P < 2 OR P > 254 THEN PRINT "BE REALISTIC": GOTO 50
70 REDIM H(P), F(255, P), E(P, P), C(P), G$(P), XX$(P), Y$(P), A(P,
P), B(P, 1), X(G, P), V(P, P), Y(P), Z(P)
75 GOTO 100
80 INPUT "DESIRED DEGREE OF POLYNOMIAL (NUMERIC) ", M%: M = M%: N = M + 1
90 IF M < 1 OR M > P - 1 THEN PRINT "RIDICULOUS VALUE [ 1 < = M < =";
P - 1; "]": GOTO 80
95 GOTO 200
100 PRINT "THE PAIRS ARE (X,Y) "
110 FOR I = 1 TO P: PRINT "COORDINATE "; I: INPUT H(I): INPUT C(I): NEXT I
115 CLS 0: PRINT "X"
120 FOR I = 1 TO P: PRINT H(I); : NEXT I: PRINT
125 PRINT "Y"
130 FOR I = 1 TO P: PRINT C(I); : NEXT I: PRINT
140 INPUT "ARE THERE ERRORS IN THE ENTRIES (Y/N) ", ER$
150 IF ER$ = "N" THEN GOSUB 620: GOTO 80
160 INPUT "NUMBER OF POINT LOCATIONS TO CORRECT ", PP
170 FOR I = 1 TO PP: PRINT "LOCATION "; I;
172 INPUT "LOCATION COORDINATE ", PCO
173 IF PCO < 1 OR PCO > P THEN PRINT "OUTSIDE RANGE ": GOTO 172
174 PRINT "X,Y OF LOCATION COORDINATE "; PCO
175 INPUT H(PCO): INPUT C(PCO): NEXT I
176 GOTO 115
200 REM SHALL SOLVE LINEAR SYSTEM BY GAUSS-SEIDEL
250 IF M > 127 OR M = P - 1 THEN 350
260 F(0, P) = P
270 FOR Q = 1 TO 2 * M: FOR K = 1 TO P
280 F(Q, K) = H(K) ^ Q + F(Q, K - 1): NEXT K, Q
290 FOR I = 1 TO N: FOR J = 1 TO N
300 Q = I + J - 2: A(I, J) = F(Q, P)
310 NEXT J, I
320 FOR I = 1 TO N: FOR K = 1 TO P
330 E(I, K) = E(I, K - 1) + C(K) * H(K) ^ (I - 1)
340 NEXT K: B(I, 1) = E(I, P): NEXT I: GOTO 400
350 I = 2
360 IF H(I) = 0 THEN T = H(I): TT = C(I): H(I) = H(1): C(I) = C(1):
H(1) = T: C(1) = TT: GOTO 380
370 I = I + 1: IF I <= P THEN 360
380 FOR I = 1 TO N: FOR J = 1 TO N
390 A(I, J) = H(I) ^ (J - 1): NEXT J: B(I, 1) = C(I): NEXT I
400 REM PREVIOUS LINES BIULT UP LINEAR SYSTEM
410 FOR I = 1 TO N: FOR J = 1 TO N: A(I, J) = A(I, J) / 100: NEXT J:
B(I, 1) = B(I, 1) / 100: NEXT I: REM TO AVOID OVERFLOW
420 FOR J = 1 TO N
430 X(0, J) = 0: B(J, 1) = B(J, 1) / A(J, J)
440 FOR K = 1 TO N: V(J, K) = -1 * A(J, K) / A(J, J)
450 NEXT K, J
460 FOR MM = 0 TO G - 1: FOR J = 1 TO N
470 GOSUB 500
480 X(MM + 1, J) = Y(J - 1) + Z(N) + B(J, 1): NEXT J, MM
490 GOTO 590
500 Y(0) = 0: Z(J) = 0
510 FOR K = 1 TO J - 1
520 Y(K) = Y(K - 1) + V(J, K) * X(MM + 1, K)
530 NEXT K
540 FOR K = J + 1 TO N
550 Z(K) = Z(K - 1) + V(J, K) * X(MM, K)
560 NEXT K
570 REM EXIT SUB
580 RETURN
590 PRINT "A RELATION IS Y =": FOR J = 1 TO N: X(G, J) = INT(1000 *
X(G, J) + .5) / 1000
592 XX$(J) = STR$(X(G, J)): Y$(J) = STR$(J - 1): G$(J) = XX$(J) + "X^"
593 G$(J) = G$(J) + Y$(J): PRINT G$(J) + " + ";
594 NEXT J: PRINT
595 IF PTER$ = "N" THEN 598
596 LPRINT "A RELATION IS Y = ": FOR J = 1 TO N: LPRINT G$(J) + " + ";
597 NEXT J: LPRINT
598 INPUT "DO YOU WISH FOR ANOTHER APPROXIMATING POLYNOMIAL (Y/N) ", AP$
599 IF AP$ = "Y" THEN 80
600 INPUT "MORE (NEW) LEASTSQUARE PAIRS "; MO$
605 IF MO$ = "Y" THEN CLS 0: GOTO 50
610 END
620 REM SUB GIVES HARDCOPY OF POINT PAIRS
630 IF PE = 0 THEN 670
640 LPRINT "      X    ,     Y    "
650 FOR I = 1 TO P: LPRINT USING "######.####"; H(I); C(I)
660 NEXT I: LPRINT
670 RETURN

I've not had a reply from rwallacej-ga, so let me at least sketch how
I intended to approach the problem.

About twenty years ago I wrote a routine (in Apple II Basis!) which
used the conjugate gradient method to solve the normal equations (for
a polynomial least squares fit; the program used other methods to find
best fit polynomials by total absolute error and min-max criteria).

Unknown to me at the time, C. C. Paige and M. A. Saunders had written
an algorithm called LSQR which shares the low storage characteristics
of conjugate gradient methods and empirically more robust for
ill-conditioned matrices.  If there is a "standard" practice now for
these problems, this is it.

So I've been working on porting this LSQR algorithm to VB.Net, with
the intention of placing it (with Google Answer's approval) under a
Common Public License.  Stanford University hosts a site where CPL'd
implementations of the algorithm are available in Fortran, C, and
Matlab code:

[LSQR: Sparse Equations and Least Squares]
http://www.stanford.edu/group/SOL/software/lsqr.html


regards, mathtalk-ga