Illustration of linear least-squares polynomial fitting
2017-09-27 Jeff Fessler, University of Michigan
2020-08-05 Julia 1.5.0
2021-08-23 Julia 1.6.2
2021-09-14 Pluto
xxxxxxxxxx10
1
md"""2
## Illustration of linear least-squares polynomial fitting 3
2017-09-27 Jeff Fessler, University of Michigan 4
​5
2020-08-05 Julia 1.5.0 6
​7
2021-08-23 Julia 1.6.28
​9
2021-09-14 Pluto10
"""xxxxxxxxxx5
1
begin2
using LinearAlgebra: Diagonal, svd3
using Random: seed!4
using Plots; default(markerstrokecolor=:auto)5
endxxxxxxxxxx2
1
sfun = (t) -> atan(4*(t-0.5)); # nonlinear function2
#s = (t) -> exp(-3t); # nonlinear function-1.04566
-1.09664
-0.903365
-0.867551
-0.932229
-0.934571
-0.78795
-0.798576
-0.936379
-0.725642
-0.58178
-0.586202
-0.474251
-0.254632
-0.311014
-0.111982
0.0391115
0.382789
0.75449
0.8368
0.997353
1.15451
1.19549
1.11689
1.00619
xxxxxxxxxx6
1
begin2
seed!(1) # seed rng3
M = 254
tm = sort(rand(M)) # M random sample locations5
y = sfun.(tm) + 0.1 * randn(M); # noisy samples6
endxxxxxxxxxx7
1
begin2
t0 = LinRange(0, 1, 101) # fine sampling for showing curve3
scatter(tm, y, color=:blue,4
xtick = 0:0.5:1, ytick = -1:1,5
label="y", xlabel="t", ylabel="y", ylim=(-1.3, 1.3))6
plot!(t0, sfun.(t0), color=:blue, label="s(t)", legend=:bottomright)7
end#11 (generic function with 1 method)xxxxxxxxxx1
1
Afun = (tt, deg) -> [t.^i for t in tt, i in 0:deg] # matrix of monomials25×4 Matrix{Float64}:
1.0 0.00790928 6.25568e-5 4.94779e-7
1.0 0.0203749 0.000415135 8.45833e-6
1.0 0.0769509 0.00592144 0.00045566
1.0 0.209472 0.0438787 0.00919137
1.0 0.210968 0.0445076 0.00938968
1.0 0.236033 0.0557117 0.0131498
1.0 0.251379 0.0631915 0.015885
â‹®
1.0 0.773223 0.597874 0.46229
1.0 0.859512 0.738761 0.634974
1.0 0.873544 0.763079 0.666583
1.0 0.951916 0.906145 0.862574
1.0 0.986666 0.973511 0.96053
1.0 0.999905 0.999809 0.999714xxxxxxxxxx7
1
begin2
# deg1 = 1 # polynomial degree3
A1 = Afun(tm, 1) # M × 2 matrix4
A2 = Afun(tm, 2) # M × 3 matrix5
# deg3 = 3 # polynomial degree6
A3 = Afun(tm, 3) # M × 3 matrix7
end-1.02238
-3.05827
14.7572
-9.58799
xxxxxxxxxx5
1
begin2
m4 = Int64.(round.(LinRange(1, M-1, 4))) # pick 4 points well separated3
A4 = A3[m4,:] # 4 × 4 matrix4
x4 = inv(A4) * y[m4] # inverse of 4×4 matrix to solve "y = A x"5
endxxxxxxxxxx1
1
plot!(t0, Afun(t0,3)*x4, color=:black, label="from 4 points")-1.01139
-0.996159
8.52431
-5.39276
xxxxxxxxxx1
1
xh = A3 \ y # initial mysterious (?) backslash solution using all M samplesxxxxxxxxxx1
1
plot!(t0, Afun(t0,3)*xh, color=:red, label="cubic from all M=$(M) points")xxxxxxxxxx1
1
#savefig("demo_ls_fit1.pdf")6.01282
2.31056
0.410402
0.0594208
xxxxxxxxxx1
1
U, s, V = svd(A3); s4×2 Matrix{Float64}:
-1.01139 -1.01139
-0.996159 -0.996159
8.52431 8.52431
-5.39276 -5.39276xxxxxxxxxx5
1
begin2
xh2 = V * Diagonal(1 ./ s) * (U' * y) # SVD-based solution3
xh3 = V * ( (1 ./ s) .* (U' * y) ) # mathematically equivalent alternate expression4
[xh2 xh3]5
end4×2 Matrix{Float64}:
-8.88178e-16 2.22045e-16
4.44089e-15 6.66134e-16
-1.42109e-14 -1.77636e-15
1.15463e-14 1.77636e-15xxxxxxxxxx1
1
[xh - xh2 xh2 - xh3]
