Least-squares curve fitting

This notebook is demonstrates how to find a least-squares degree n-1 polynomial to m>n noisy (x,y) datapoints, using a Vandermonde matrix, QR decomposition, and backsubstitution.

Get some noisy datapoints

In [17]:
# define an underlying functional relation y=f(x)
f(x) = 1 - 3x + 2x.^2

# collect m noisy datapoints, uniformly spaced
m = 31
xdata = linspace(0,2,m)
ydata = f(xdata) + (rand(m)-0.5)

# plot noisy data and "true" curve
using PyPlot
x = linspace(-0.5, 2.25, 100)
plot(x, f(x), "b--")
plot(xdata, ydata, "bo")
legend(("f(x)", "noisy data"),loc="upper left")
WARNING: Method definition f(Any) in module Main at In[8]:2 overwritten at In[17]:2.
Out[17]:
PyObject <matplotlib.legend.Legend object at 0x7f6971259c10>

Fit the noisy datapoints with a polynomial interpolant

In [18]:
# fit m noisy datapoints with an m-1 degree polynomial 
# using Newton Divided Differences

include("polyinterp.jl")
(c,b) = polyinterp(xdata, ydata)

plot(xdata, ydata, "bo")
plot(x, f(x), "b--")
plot(x, polyeval(c,b,x), "r-")
legend(("noisy data", "f(x)", "polynomial interpolant"),loc="upper left")
ylim(-1,5)
WARNING: Method definition polyeval(Any, Any) in module Main at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:6 overwritten at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:6.
WARNING: replacing docs for 'polyeval :: Tuple{Any,Any}' in module 'Main'.
WARNING: Method definition polyeval(Any, Any, Any) in module Main at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:19 overwritten at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:19.
WARNING: replacing docs for 'polyeval :: Tuple{Any,Any,Any}' in module 'Main'.
WARNING: Method definition polyinterp(Any, Any) in module Main at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:41 overwritten at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:41.
WARNING: replacing docs for 'polyinterp :: Tuple{Any,Any}' in module 'Main'.
WARNING: Method definition backsolve(Any, Any) in module Main at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:59 overwritten at /home/gibson/professional/teaching/math753/lecture-10-26-leastsquares/polyinterp.jl:59.
Out[18]:
(-1,5)

Least-squares curve fit

In [19]:
# Construct m x n Vandermonde matrix for least-squares polynomial 
n = 4
X = zeros(m,n)
X[:,1] = ones(m)
for j = 2:n
    X[:,j] = X[:,j-1] .* xdata
end
In [20]:
X
Out[20]:
31×4 Array{Float64,2}:
 1.0  0.0        0.0         0.0        
 1.0  0.0666667  0.00444444  0.000296296
 1.0  0.133333   0.0177778   0.00237037 
 1.0  0.2        0.04        0.008      
 1.0  0.266667   0.0711111   0.018963   
 1.0  0.333333   0.111111    0.037037   
 1.0  0.4        0.16        0.064      
 1.0  0.466667   0.217778    0.10163    
 1.0  0.533333   0.284444    0.151704   
 1.0  0.6        0.36        0.216      
 1.0  0.666667   0.444444    0.296296   
 1.0  0.733333   0.537778    0.39437    
 1.0  0.8        0.64        0.512      
 ⋮                                      
 1.0  1.26667    1.60444     2.0323     
 1.0  1.33333    1.77778     2.37037    
 1.0  1.4        1.96        2.744      
 1.0  1.46667    2.15111     3.15496    
 1.0  1.53333    2.35111     3.60504    
 1.0  1.6        2.56        4.096      
 1.0  1.66667    2.77778     4.62963    
 1.0  1.73333    3.00444     5.2077     
 1.0  1.8        3.24        5.832      
 1.0  1.86667    3.48444     6.5043     
 1.0  1.93333    3.73778     7.22637    
 1.0  2.0        4.0         8.0        
In [21]:
# Solve least-squares X c = y problem with QR decomp and backsubstitution
(Q,R) = qr(X);
c = backsolve(R, Q'*ydata)

# plot data, f(x), and least-squares polynomial fit
plot(xdata, ydata, "bo")
plot(x, f(x), "b--")
plot(x, polyeval(c,x), "r-")
legend(("noisy data", "f(x)", "least-squares polynomial"),loc="upper left")
ylim(-1,5)
Out[21]:
(-1,5)
In [16]:
c
Out[16]:
4-element Array{Float64,1}:
  0.663366 
 -1.97982  
  1.3518   
  0.0628596
In [ ]:
norm(X*c - ydata)
In [ ]:
norm(X*c - ydata)/sqrt(m)
In [ ]: