Add smoothing tests using AD

Author: jipolanco

test/Project.toml

@@ -4,6 +4,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

test/smoothing.jl

@@ -2,8 +2,67 @@
 
 using BSplineKit
 using QuadGK: quadgk
+using ReverseDiff
 using Test
 
+# This is the objective function that `fit` is supposed to minimise.
+# We can verify this using automatic differentiation: the gradient wrt the spline
+# coefficients should be zero.
+function smoothing_objective(cs, R::AbstractBSplineBasis, xs, ys; weights = nothing, λ)
+    # Construct spline from coefficients and knots
+    S = Spline(R, cs)
+    S″ = Derivative(2) * S
+
+    # Compute first term of objective (loss) function
+    T = eltype(cs)
+    loss = zero(T)
+    for i in eachindex(xs, ys)
+        w = weights === nothing ? 1 : weights[i]
+        loss += w * abs2(ys[i] - S(xs[i]))
+    end
+
+    # Integrate roughness term interval by interval
+    for i in eachindex(xs)[1:end-1]
+        # Note: S″(x) is linear within each interval, and thus the integrand is quadratic.
+        # Therefore, a two-point GL quadrature is exact (weights = 1 and locations = ±1/√3).
+        a = xs[i]
+        b = xs[i + 1]
+        Δ = (b - a) / 2
+        xc = (a + b) / 2
+        gl_weight = 1
+        gl_ξ = 1 / sqrt(T(3))
+        for ξ in (-gl_ξ, +gl_ξ)
+            x = Δ * ξ + xc
+            loss += λ * Δ * gl_weight * abs2(S″(x))
+        end
+    end
+
+    loss
+end
+
+function check_zero_gradient(S::Spline, xs, ys; weights = nothing, λ)
+    cs = coefficients(S)
+    R = basis(S)  # usually a RecombinedBSplineBasis
+
+    # Not sure how useful this is...
+    ∇f = similar(cs)  # gradient wrt coefficients
+    inputs = (cs,)
+    results = (∇f,)
+    # all_results = map(DiffResults.GradientResult, results)
+    cfg = ReverseDiff.GradientConfig(inputs)
+
+    # Compute gradient
+    ReverseDiff.gradient!(results, cs -> smoothing_objective(cs, R, xs, ys; weights, λ), inputs, cfg)
+
+    # Verify that |∇f|² is negligible. Note that is has the same units as |y_i|² ≡ Y², since
+    # f ~ Y² and therefore ∂f/∂cⱼ ~ Y. So we compare it with the sum of |y_i|².
+    reference = sum(abs2, ys)
+    err = sum(abs2, ∇f)
+    @test err / reference < 1e-12
+
+    nothing
+end
+
 # Returns the integral of |S''(x)| (the "curvature") over the whole spline.
 function total_curvature(S::Spline)
     ts = knots(S)
@@ -41,10 +100,11 @@ end
         λs = [0.0, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2]
         curvatures = similar(λs)
         distances = similar(λs)
-        for i in eachindex(λs)
-            S = @inferred fit(xs, ys, λs[i])
+        for (i, λ) in pairs(λs)
+            S = @inferred fit(xs, ys, λ)
             curvatures[i] = total_curvature(S)
             distances[i] = distance_from_data(S, xs, ys)
+            check_zero_gradient(S, xs, ys; λ)
         end
         @test issorted(curvatures; rev = true)  # in decreasing order (small λ => large curvature)
         @test issorted(distances)  # in increasing order (large λ => large distance from data)
@@ -55,10 +115,12 @@ end
         weights = fill!(similar(xs), 1)
         S = fit(xs, ys, λ)
         Sw = @inferred fit(xs, ys, λ; weights)  # equivalent to the default (all weights are 1)
+        check_zero_gradient(Sw, xs, ys; λ, weights)
         @test coefficients(S) == coefficients(Sw)
         # Now give more weight to point i = 3
         weights[3] = 1000
         Sw = fit(xs, ys, λ; weights)
+        check_zero_gradient(Sw, xs, ys; λ, weights)
         @test abs(Sw(xs[3]) - ys[3]) < abs(S(xs[3]) - ys[3])  # the new curve is closer to the data point i = 3
         @test total_curvature(Sw) > total_curvature(S)  # since we give more importance to data fitting (basically, the sum of weights is larger)
         @test distance_from_data(Sw, xs, ys) < distance_from_data(S, xs, ys)