Algorithm for polynomial least squares fitting

Given the data points for ( is large), we want to find a degree polynomial that fits the data in the least squares sense. This polynomial is

To determine the coefficients we will use the method of normal equations. For this we build the matrix

and the vector of measurements

has a special form and is called a Vandermonde matrix.

The system of normal equations is

To determine the matrix we apply the matrix multiplication formula:

Clearly, so we need to determine only the elements of the upper triangle of .

For the right hand side vector, we have

! Compute P
DO i=1,n+1
   DO j=1,i
      P(i,j) = 0.0
      ipower=i+j-2
      DO k=1,m
         P(i,j) = P(i,j) + t(k)**ipower
      END DO
      P(j,i) = P(i,j)
   END DO
END DO

! Compute b
DO i=1,n+1
    b(i) = 0.0
    DO k=1,m
       b(i) = b(i) + t(k)**(i-1)*y(k)
    END DO
END DO

! Solve the system Pa=b etc.