In [1]:
using Random, LinearAlgebra, Plots
# --- Generate data ---
nx = 80
f(x) = x * exp(-x)
xs_raw = range(0, 6, length=nx)
xs = @. xs_raw / 6
Random.seed!(3)
b_true = f.(xs_raw)
noise = 0.08 * randn(nx)
b = b_true + noise .* b_true.^0.5
# --- Biharmonic regularization matrix ---
function biharmonic_matrix(n)
B = zeros(n, n)
for i in 1:n, j in 1:n
if i > 1 && j > 1
B[i, j] = i * (i - 1) * j * (j - 1) / (i + j - 3)
end
end
println("--- Biharmonic regularization matrix (B matrix) ---")
show(IOContext(stdout, :limit => true), "text/plain", B)
println()
println()
return Symmetric(B)
end
# --- Polynomial fit with curvature regularization ---
function poly_fit_biharmonic(n, λ)
A = [x^j for x in xs, j in 1:n]
println("--- Polynomial fit with curvature regularization (A matrix) ---")
show(IOContext(stdout, :limit => true), "text/plain", A)
println()
println()
if λ == 0
# --- pure least-squares solution ---
c = A \ b
else
# --- regularized least squares ---
B = biharmonic_matrix(n)
c = (A' * A + λ * B) \ (A' * b)
end
# --- Evaluate fitted polynomial ---
p_scaled = [sum(c[j] * x^j for j in 1:n) for x in xs]
return p_scaled
end
# --- Example: compare regularized vs unregularized fits ---
n = 40
λ = 3e-4
p_reg = poly_fit_biharmonic(n, λ)
p_unreg = poly_fit_biharmonic(n, 0)
plot(xs_raw, b, label="Noisy data", lw=1)
plot!(xs_raw, b_true, label="True f(x)", lw=2)
plot!(xs_raw, p_reg, label="Curvature reg. n=$n, λ=$(λ)", lw=2)
plot!(xs_raw, p_unreg, label="No reg. (λ=0)", lw=2)
xlabel!("x")
ylabel!("f(x)")
title!("Curvature (biharmonic) regularization on polynomial fit")
--- Polynomial fit with curvature regularization (A matrix) ---
80×40 Matrix{Float64}:
0.0 0.0 0.0 … 0.0 0.0
0.0126582 0.000160231 2.02824e-6 9.82977e-75 1.24427e-76
0.0253165 0.000640923 1.62259e-5 5.40397e-63 1.36809e-64
0.0379747 0.00144208 5.47624e-5 3.98357e-56 1.51275e-57
0.0506329 0.00256369 0.000129807 2.97087e-51 1.50424e-52
0.0632911 0.00400577 0.00025353 … 1.78803e-47 1.13166e-48
0.0759494 0.00576831 0.000438099 2.18999e-44 1.66328e-45
0.0886076 0.00785131 0.000695685 8.94061e-42 7.92206e-43
0.101266 0.0102548 0.00103846 1.63325e-39 1.65393e-40
0.113924 0.0129787 0.00147858 1.61436e-37 1.83915e-38
0.126582 0.0160231 0.00202824 … 9.82977e-36 1.24427e-36
0.139241 0.0193879 0.00269958 4.04444e-34 5.6315e-35
0.151899 0.0230732 0.00350479 1.20396e-32 1.8288e-33
⋮ ⋱
0.860759 0.740907 0.637743 0.00288663 0.00248469
0.873418 0.762859 0.666294 0.00510103 0.00445533
0.886076 0.785131 0.695685 … 0.00894061 0.00792206
0.898734 0.807723 0.725928 0.015546 0.0139717
0.911392 0.830636 0.757035 0.0268232 0.0244465
0.924051 0.85387 0.789019 0.0459341 0.0424454
0.936709 0.877423 0.82189 0.0780873 0.0731451
0.949367 0.901298 0.855663 … 0.131805 0.125132
0.962025 0.925493 0.890347 0.22094 0.212549
0.974684 0.950008 0.925957 0.367859 0.358546
0.987342 0.974844 0.962504 0.608461 0.600759
1.0 1.0 1.0 1.0 1.0
--- Biharmonic regularization matrix (B matrix) ---
40×40 Matrix{Float64}:
0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0
0.0 4.0 6.0 8.0 10.0 76.0 78.0 80.0
0.0 6.0 12.0 18.0 24.0 222.0 228.0 234.0
0.0 8.0 18.0 28.8 40.0 432.615 444.6 456.585
0.0 10.0 24.0 40.0 57.1429 703.0 722.927 742.857
0.0 12.0 30.0 51.4286 75.0 … 1028.78 1058.57 1088.37
0.0 14.0 36.0 63.0 93.3333 1406.0 1447.53 1489.09
0.0 16.0 42.0 74.6667 112.0 1831.07 1886.18 1941.33
0.0 18.0 48.0 86.4 130.909 2300.73 2371.2 2441.74
0.0 20.0 54.0 98.1818 150.0 2812.0 2899.57 2987.23
0.0 22.0 60.0 110.0 169.231 … 3362.17 3468.51 3575.0
0.0 24.0 66.0 121.846 188.571 3948.77 4075.5 4202.45
0.0 26.0 72.0 133.714 208.0 4569.5 4718.2 4867.2
⋮ ⋱
0.0 58.0 168.0 324.8 523.871 17838.6 18513.6 19192.7
0.0 60.0 174.0 336.774 543.75 18818.8 19535.5 20256.7
0.0 62.0 180.0 348.75 563.636 … 19811.8 20571.0 21335.3
0.0 64.0 186.0 360.727 583.529 20817.2 21619.8 22427.8
0.0 66.0 192.0 372.706 603.429 21834.4 22681.0 23533.7
0.0 68.0 198.0 384.686 623.333 22862.8 23754.3 24652.4
0.0 70.0 204.0 396.667 643.243 23902.0 24839.2 25783.3
0.0 72.0 210.0 408.649 663.158 … 24951.5 25935.0 26926.0
0.0 74.0 216.0 420.632 683.077 26011.0 27041.4 28080.0
0.0 76.0 222.0 432.615 703.0 27079.9 28158.0 29244.8
0.0 78.0 228.0 444.6 722.927 28158.0 29284.3 30420.0
0.0 80.0 234.0 456.585 742.857 29244.8 30420.0 31605.2
--- Polynomial fit with curvature regularization (A matrix) ---
80×40 Matrix{Float64}:
0.0 0.0 0.0 … 0.0 0.0
0.0126582 0.000160231 2.02824e-6 9.82977e-75 1.24427e-76
0.0253165 0.000640923 1.62259e-5 5.40397e-63 1.36809e-64
0.0379747 0.00144208 5.47624e-5 3.98357e-56 1.51275e-57
0.0506329 0.00256369 0.000129807 2.97087e-51 1.50424e-52
0.0632911 0.00400577 0.00025353 … 1.78803e-47 1.13166e-48
0.0759494 0.00576831 0.000438099 2.18999e-44 1.66328e-45
0.0886076 0.00785131 0.000695685 8.94061e-42 7.92206e-43
0.101266 0.0102548 0.00103846 1.63325e-39 1.65393e-40
0.113924 0.0129787 0.00147858 1.61436e-37 1.83915e-38
0.126582 0.0160231 0.00202824 … 9.82977e-36 1.24427e-36
0.139241 0.0193879 0.00269958 4.04444e-34 5.6315e-35
0.151899 0.0230732 0.00350479 1.20396e-32 1.8288e-33
⋮ ⋱
0.860759 0.740907 0.637743 0.00288663 0.00248469
0.873418 0.762859 0.666294 0.00510103 0.00445533
0.886076 0.785131 0.695685 … 0.00894061 0.00792206
0.898734 0.807723 0.725928 0.015546 0.0139717
0.911392 0.830636 0.757035 0.0268232 0.0244465
0.924051 0.85387 0.789019 0.0459341 0.0424454
0.936709 0.877423 0.82189 0.0780873 0.0731451
0.949367 0.901298 0.855663 … 0.131805 0.125132
0.962025 0.925493 0.890347 0.22094 0.212549
0.974684 0.950008 0.925957 0.367859 0.358546
0.987342 0.974844 0.962504 0.608461 0.600759
1.0 1.0 1.0 1.0 1.0
Out[1]:
