EXPLANATION FILE OF PROGRAM LSQPLY ================================== Least-Squares Fitting With Orthogonal Polynomials ------------------------------------------------- The major problem associated with applying the one-dimensional least-squares polynomial algorithm given before (see file Lstsqr.txt) is round-off error. This error is due to the large number of mathematical operations associated with implementing equation (1.4.4): t -1 t D = (X X) X Y By using Forsythe polynomials instead of simple powers of X, however, you can still apply the same general matrix equation while drastically reducing the number of calculations. (For a discussion of this technique, see Applied Numerical Methods, Ref. 10. It is interesting to note that the very signifi- cant round-off error implications of this alternative method are fairly obvious, but appear to have escaped forceful attention in the common litera- ture). We will denote the Forsythe polynomial of degree j as Pj(x). The associated approximation equation which is equivalent to equation (1.4.1) is y(X) = D + D P (x) + D P (X) + ... + D P (X) (1.6.1) 0 1 1 2 2 n n We will define the matrix P as | 1 P1(x1) P2(x1) ... Pn(x1) | | 1 | P = | 1 | | 1 | | 1 | | 1 P1(xm) P2(xm) ... Pn(xm) | This has the sa me appearance as the X matrix shown in section 1.4, and we can thereby write the approximation equation as Y = PD (1.6.2) Again, following the analysis of that section, we will define the error func- tion for the fit to be t Error = (Y - PD) (Y - PD) (1.6.3) The corresponding solution obtained for the coefficient vector, 0, is t -1 t D = (P P) P Y (1.6.4) The savings in computation (and thus in round-off error) is associated with the form of P^tP. If all the off-diagonal terms of that matrix could somehow be forced to be zero, the total number of calculations would be greatly re- duced. (P^tP)^-1 would then not require a full matrix inversion that is very prone to round-off error. ln fact, if that could be accomplished, the coef- ficients would simply be m Sum P (xi) yi i=1 j D = -------------- (1.6.5) j m 2 Sum P (xi) i=1 j The ability to force this condition rests in the choice of the Pj(x) poly- nomials. Forsythe found just a set. The recursion relation for calculating them is P (x) = (x - g ) P (x) - d P (x) (1.6.6) j j j-1 j j-2 with P (x) = 0 and P (x) = 1 -1 1 m 2 Sum xi |P (xi)| i=1 j-1 g = ----------------- (1.6.7) j m 2 sum |P (xi)| i=1 j-1 m Sum xi P (xi) P (xi) i=1 j-1 j-2 d = ------------------------ (1.6.8) j m 2 sum |P (xi)| i=1 j-2 Besides having the nice round-off error property, this formulation also has the advantage that the order of the fit can be increased by one degree without extensive recalculation. The calculation for the additional degree can simply be added on. Program LSQPLY performs these calculations. The user inputs the degree of fit (N), the number of data points (M), and the data pairs [X(I),Y(I)]. Also, an optional error reduction factor (E) can be included. If E = 0, the subrou- tine will provide a fit to the degree N. If E > 0, an increasing sequence of fits will be tried until the standard deviation either is reduced by less than that factor or increases, or until the Nth-degree fit is reached. ln the case of termination before the Nth degree, the actual degree of fit is retur- ned in L. It is apparent from the examples given that this algorithm is fairly resis- tant to round-off error. The precautions involved in using this subroutine are simple. The number of distinct data points must exceed the degree of fit (M >= N + 1). E must be nonnegative. This method is highly recommended for polynomial regression. From [BIBLI 01]. ----------------------------------------------------- End of file lsqply.txt