GitHub

Spline grid fitting

This page illustrates finite-series modeling with splines with the intention of comparing to INR.

URL.

This page was generated from a single Julia file: 01-spline.jl.

In any such Julia documentation, you can access the source code using the "Edit on GitHub" link in the top right.

The corresponding notebook can be viewed in nbviewer here: 01-spline.ipynb, and opened in binder here: 01-spline.ipynb.

Setup

Packages needed here.

using BSplineKit
using LinearAlgebra: I, norm
using MIRTjim: jim, prompt
using InteractiveUtils: versioninfo
using Plots
using Random: seed!
default(markerstrokecolor=:auto, label="")
seed!(0)
Random.TaskLocalRNG()

The following line is helpful when running this file as a script; this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

Data

N = 100
tmax = 3
t = tmax * sort(rand(N)) # random samples on [0,3]
f(t) = 1 < t < 2 # rect function
y = f.(t); # noiseless data for now

plot data and true signal (finely sampled)

yaxis = ((-0.2, 1.4), -0.2:0.2:1.4)
tf = range(-0.3, tmax+0.3, 1001) # fine grid
ytf = f.(tf) # true function on a fine grid
p0 = scatter(t, y; yaxis, color=:black, label="data")
plot!(p0, tf, ytf, label="true")
Example block output

B-spline settings

degree = 3
M = 40
knots = range(0, tmax, M) |> collect;

Make B-spline basis objects

basis = BSplineBasis(BSplineOrder(degree + 1), knots) # ; periodic=true
42-element BSplineBasis of order 4, domain [0.0, 3.0]
 knots: [0.0, 0.0, 0.0, 0.0, 0.0769231, 0.153846, 0.230769, 0.307692, 0.384615, 0.461538  …  2.53846, 2.61538, 2.69231, 2.76923, 2.84615, 2.92308, 3.0, 3.0, 3.0, 3.0]

Plot finely sampled B-spline basis functions

Bf = hcat([b.(tf) for b in basis]...)
pb1 = plot(tf, Bf; title="Basis functions, degree=$degree")
Example block output

Sample B-spline basis objects at measurement points

B = hcat([b.(t) for b in basis]...)
pb2 = deepcopy(pb1)
scatter!(pb2, t, B; title="Basis function samples")
Example block output

LS fit of B-splines to data

xh = B \ y
yhf = +([b.(tf) * xh[i] for (i,b) in enumerate(basis)]...);

Utility functions for NRMSE

do_nrmse(yr::AbstractVector) = norm(yr - ytf) / norm(ytf)
round3(x) = round(x, digits=3);

p1 = deepcopy(p0)
plot!(p1, tf, yhf; label="bspline, NRMSE=$(round3(do_nrmse(yhf)))")
Example block output

Tikhonov regularization

Evidently the unregularized B-spline fit here is over-fitting, even though $N < M$, so next we try simple Tikhonov regularization.

function do_fit(β::Real)
    A = [B; sqrt(β)*I]
    yz = [y; zeros(length(basis))]
    xr = A \ yz
    yr = +([b.(tf) * xr[i] for (i,b) in enumerate(basis)]...)
end
nrmse_fit(β::Real) = do_nrmse(do_fit(β));

β = 2e-1
yfr = do_fit(β);

p2 = deepcopy(p0)
plot!(p2, tf, yfr,
 label = "bspline β=$(round3(β)), NRMSE=$(round3(do_nrmse(yfr)))")
Example block output

Parameter tuning

What regularization parameter β to choose?

Cross-validation is one approach to choose.

Here we use the simpler "oracle" approach of finding the β value that leads to the best fit to the true function (which would be unknown in practice).

It turns out here that Tikhonov regularization only slightly reduces NRMSE, even in this very favorable case where we use the ground-truth "oracle" to select β.

Apparently Tikhonov regularization is ineffective for (uniformly spaced) B-spline fitting to non-uniformly data.

βlist = 10 .^ (-5:0.2:1)
nrmse = nrmse_fit.(βlist)
βbest = argmin(nrmse_fit, βlist)
pn = plot(log10.(βlist), 100*nrmse; title="NRMSE (%)", xlabel="log10(β)",
 yaxis = ("NRMSE (%)", (10, 30)),
)
plot!(log10(βbest) * [1, 1], [10, 40], color=:red)
Example block output
yfb = do_fit(βbest)
p3 = deepcopy(p0)
nrmse_best = nrmse_fit(βbest)
plot!(p3, tf, yfb,
 label="bspline β=$(round3(βbest)), NRMSE=$(round3(nrmse_best))")
Example block output

Extensions

  • 2D image instead of 1D signal
  • comparison with implicit neural representation (INR)

Reproducibility

This page was generated with the following version of Julia:

using InteractiveUtils: versioninfo
io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')
12-element Vector{SubString{String}}:
 "Julia Version 1.12.2"
 "Commit ca9b6662be4 (2025-11-20 16:25 UTC)"
 "Build Info:"
 "  Official https://julialang.org release"
 "Platform Info:"
 "  OS: Linux (x86_64-linux-gnu)"
 "  CPU: 4 × AMD EPYC 7763 64-Core Processor"
 "  WORD_SIZE: 64"
 "  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)"
 "  GC: Built with stock GC"
 "Threads: 1 default, 1 interactive, 1 GC (on 4 virtual cores)"
 ""

And with the following package versions

import Pkg; Pkg.status()
Status `~/work/Work/Work/docs/Project.toml`
  [093aae92] BSplineKit v0.19.1
  [e30172f5] Documenter v1.16.1
  [b964fa9f] LaTeXStrings v1.4.0
  [7031d0ef] LazyGrids v1.1.0
  [98b081ad] Literate v2.21.0
  [170b2178] MIRTjim v0.26.0
  [91a5bcdd] Plots v1.41.2
  [bf9f9b4b] Work v0.0.1 `~/work/Work/Work`
  [b77e0a4c] InteractiveUtils v1.11.0
  [37e2e46d] LinearAlgebra v1.12.0
  [44cfe95a] Pkg v1.12.0
  [9a3f8284] Random v1.11.0

This page was generated using Literate.jl.