In this notebook, we'll do some simple examples of least-square fitting to polynomials, also called polynomial regression.
We'll start by creating some sample data: a hyperbola $y = \sqrt{x^2 + 1}$ plus some (gaussian) noise.
using PyPlot
x = linspace(0, 3, 200)
y = sqrt(x.^2 + 1) # hyperbola
yn = y + randn(length(y))*0.1 # + noise
plot(x, y, "k--")
plot(x, yn, "r.")
legend(["hyperbola", "hyperbola + noise"], loc="upper left")
PyObject <matplotlib.legend.Legend object at 0x31df9b5d0>
Our least-square solution minimizes the error norm $\Vert Ac - b\Vert_2$ where $c$ are the coefficients of a polynomial fit $c[1] + c[2] x + c[3] x^2 + \cdots + c[n+1] x^n$, $b$ is the vector of data points $yn[i]$ that we are fitting to, and $A$ is the matrix whose columns are $x^0, x^1, \ldots$. There are several ways to construct the matrix $A$ in Julia, for example for $n=2$:
A = [x.^0 x.^1 x.^2] # .^n is element-wise exponentiation and [a b c] concatenates columns into a matrix
200x3 Array{Float64,2}: 1.0 0.0 0.0 1.0 0.0150754 0.000227267 1.0 0.0301508 0.000909068 1.0 0.0452261 0.0020454 1.0 0.0603015 0.00363627 1.0 0.0753769 0.00568167 1.0 0.0904523 0.00818161 1.0 0.105528 0.0111361 1.0 0.120603 0.0145451 1.0 0.135678 0.0184086 1.0 0.150754 0.0227267 1.0 0.165829 0.0274993 1.0 0.180905 0.0327264 ⋮ 1.0 2.83417 8.03252 1.0 2.84925 8.1182 1.0 2.86432 8.20434 1.0 2.8794 8.29093 1.0 2.89447 8.37797 1.0 2.90955 8.46547 1.0 2.92462 8.55342 1.0 2.9397 8.64183 1.0 2.95477 8.73069 1.0 2.96985 8.82 1.0 2.98492 8.90978 1.0 3.0 9.0
# More generally, .^ is a "broadcasting" operation that can combine a column vector x and a row vector [0 1 2]
# into a matrix:
A = x .^ [0 1 2]
200x3 Array{Float64,2}: 1.0 0.0 0.0 1.0 0.0150754 0.000227267 1.0 0.0301508 0.000909068 1.0 0.0452261 0.0020454 1.0 0.0603015 0.00363627 1.0 0.0753769 0.00568167 1.0 0.0904523 0.00818161 1.0 0.105528 0.0111361 1.0 0.120603 0.0145451 1.0 0.135678 0.0184086 1.0 0.150754 0.0227267 1.0 0.165829 0.0274993 1.0 0.180905 0.0327264 ⋮ 1.0 2.83417 8.03252 1.0 2.84925 8.1182 1.0 2.86432 8.20434 1.0 2.8794 8.29093 1.0 2.89447 8.37797 1.0 2.90955 8.46547 1.0 2.92462 8.55342 1.0 2.9397 8.64183 1.0 2.95477 8.73069 1.0 2.96985 8.82 1.0 2.98492 8.90978 1.0 3.0 9.0
In Julia (similar to Matlab), if you do A \ b
with a square matrix $A$ it computes $A^{-1} b$ (albeit implicitly: it doesn't actually compute the inverse matrix, but instead performs Gaussian elimination to find the LU factorization of $A$). However, if you do A \ b
with a non-square matrix $A$ then it returns the least-square solution (via a QR factorization, but we won't worry about the algorithm for now).
Let's use the Interact
package to give us a slider control for the polynomial degree $n$, so that we can fit to polynomials of various degrees:
using Interact
f = figure()
@manipulate for n = 1:40
withfig(f) do
A = x .^ [0:n]'
c = A \ yn
plot(x, yn, "r.")
plot(x, A * c, "b-")
legend(["data", "fit to degree $n"], loc="upper left")
end
end
One of the things to notice is that, once the degree $n$ goes above 29 or so, the fit goes crazy. That is because the matrix $A$ is becoming so ill-conditioned that roundoff errors are killing us.
(In practice, one never fits to polynomials of such high degree; the same conditioning problems mean that the fit parameters become increasingly sensitive to errors in the data as you increase the polynomial degree.)
Let's use cond(A)
to compute and plot the condition number as a function of the degree $n$:
N = 1:20
semilogy(N, [cond(x .^ [0:n]') for n in N], "bo-")
xlabel(L"polynomial degree $n$")
ylabel(L"condition number of $A$")
title("condition number of least-square polynomial fitting")
PyObject <matplotlib.text.Text object at 0x31df252d0>
As you can see, the condition number increases exponentially with $n$. This is a general problem of with polynomial fits, especially with equispaced $x$ values.
(There are special polynomial-fitting techniques based on Chebyshev polynomials which allow high-degree polynomial interpolation without this kind of problem, but they don't use equispaced points.)