Linear regression via SVD not producing best fit with escalating polynomial degree

Question

I am using a basic singular value decomposition (via LAPACK) routine in FORTRAN to solve an overdetermined system in the form of $A\cdot X = B$ where $\mathrm{size}(A) = [m,n]$ with $m > n$.

My sample data points come from a noisy sine function and I am trying to use linear regression with $x^i$ as my basis functions. I find that (with the noisy sine function) I get very good approximations when I keep my polynomial low. That is, when I fit a function of the form $$ a_0 + a_1 x + a_2 x^2 + \ldots + a_Nx^N $$ with $N \lesssim 10$, I get great results. When I allow $N$ to get higher (I have 1000 data points), say, up to 250, my fit goes to hell and my sum of squares is $\approx 5.6\times 10^{41}$ (as opposed to about 200 with $N \lesssim 10$).

Why does this happen? If a degree 6 polynomial, for example, provides the best fit, then shouldn't a degree $N$ polynomial produce that same fit with $0$ for $a_7 \ldots a_N$?

Answer 1

The problem here is that the condition number for your monomial basis becomes very large, meaning that the value of the monomial is very sensitive to the value of the coefficient. Thus, when you try to compute in this basis for high degree polynomials, your answer is very inaccurate. There are polynomial bases that are better behaved, like any of the families of orthogonal polynomials.

Answer 2

Sections of Vandermonde matrices do tend to be ill conditioned. If indeed you were using the singular value decomposition for your fitting, you should have already seen a rash of tiny singular values. Did you remember to zero out those tiny singular values before computing the least squares solution?

The ill-condition is easily visualized: if you plot successive members of the family $x^k$ for increasing $k$, you'll find that the members look almost indistinguishable for high enough $k$. This manifests itself in the columns becoming nearly linearly dependent.

Alternatively, like Matt says, you can use orthogonal polynomials. (I talked about them here, but I suppose I should say something here for convenience.) If your data points are equispaced in the abscissas, you can use the Gram polynomials as a basis. If they are irregularly spaced, a fair bit of work is needed; see Forsythe's paper for details.