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

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$?

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.

• it turns out that this is not the case. one possible answer may be that my original algorithm was no throwing out singular values due to the ratio of the smallest to the largest, but rather by the a given tolerance altogether. however, even after correcting that mistake i am still having some problems but i do not believe it is because of the linear dependence of a monomial basis set. – drjrm3 Dec 10 '11 at 16:54
• "i do not believe it is because of the linear dependence of a monomial basis set." - maybe not entirely, but I'd say it's a contributing factor. – J. M. Dec 12 '11 at 4:03

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.

• (I might elaborate further if there's interest...) – J. M. Dec 4 '11 at 3:11
• The Forsythe link does not resolve for me. – Matt Knepley Dec 4 '11 at 4:37
• @Matt: Strange; I already fired off an e-mail to SIAM. I've replaced the link in the meantime. – J. M. Dec 4 '11 at 4:45
• i have used legendre polynomials and the same effect occurs. i have a new question out about when to eliminate singular values. if i throw out too many, my answer is bad, if i don't through out enough, my answer is equally bad. there seems to be a magical number that works - how do i determine the best way to do this? – drjrm3 Dec 12 '11 at 3:59
• I wasn't even thinking about Legendre, @Laurbert. Forsythe's procedure generates orthogonal polynomials that are custom-made for your abscissas. In any event, the usual criterion for singular value removal is to dispose of singular values that are less than the largest singular value times machine epsilon. – J. M. Dec 12 '11 at 4:02