[Feature request] Smoothing factor

[getzze]

This is mentioned in #5 (comment).
The feature is provided by the Dierckx.jl package and in R as mentioned in the comment, but also by scipy https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html
It is a very useful parameter to avoid overfitting.

[tomasaschan]

@getzze Do you have any resources on the mathematics behind such a smoothing factor? Something one could look at to implement it?

[getzze]

From wikipedia I found this open-access paper: https://projecteuclid.org/euclid.ss/1038425655
The minimizer is given in equation (13).
Or more generalist, Chapter 5.4 Smoothing splines of The elements of statistical learning by Hastie, Tibshirani and Friedman.

Looking at the implementation in fortran in DIERCKX (http://www.netlib.org/dierckx/ the curfit.f file) , I wonder if both smoothing techniques are equivalent. Quoting from the scipy doc for the wrapper for the fortran lib:

s : Positive smoothing factor used to choose the number of knots. Number of knots will be increased until the smoothing condition is satisfied:

sum((w[i] * (y[i]-spl(x[i])))**2, axis=0) <= s

[tomasaschan]

Thanks - the open-access paper can hopefully be really useful if someone wants to tackle implementing this!

For the DIERCKX implementation I can't find a mention of a license anywhere, which means the implementation is copyrighted by the author. Because of that, unless you can somewhere find an open-source (MIT-compatible) license for the code we can't accept a contribution based on that implementation without written consent from the author (who does, conveniently, list some contact info in the source files).

[getzze]

Coming back to this issue.
I tried implementing the smoothing factor as in the open-access paper I found before. Because the structure of the package is really complex and is really optimized, I had to directly modify the prefiltering.jl file, just for a proof of concept.

using BandedMatrices
using SparseArrays
using LinearAlgebra: lu

function prefilter!(
    ::Type{TWeights}, ret::TCoefs, it::BSpline
    ) where {TWeights,TCoefs<:AbstractArray}
    local buf, shape, retrs
    sz = size(ret)
    first = true
    for dim in 1:ndims(ret)
        M, b = prefiltering_system(TWeights, eltype(TCoefs), sz[dim], degree(it))
        ### TEST REGULARIZATION
        if true
            λ = 100000
            n = sz[dim]
            k = 3
#            D1 = BandedMatrix( (0=>-ones(TWeights, n-1), 1=>ones(TWeights, n-1)), (n-1,n), (0,1))
            D1 = spdiagm(0 => -ones(TWeights,n-1), 1 => ones(TWeights,n-1))[1:end-1, :]
            D = D1
            for i in 1:k-1
                D = D1[1+i:end, 1+i:end] * D
#                lmul!(D1[1+i:end, 1+i:end], D)
            end
            Q = D' * D
            K = cholesky(Matrix(Matrix(M)' * Matrix(M) + λ/n * Q))
            B = Matrix(M)' * popwrapper(ret)
            ldiv!(popwrapper(ret), K, B)
        else
            A_ldiv_B_md!(popwrapper(ret), M, popwrapper(ret), dim, b)
        end
    end
    ret
end

Then some code for testing:

using Revise
using Interpolations
using Random

using Plots

Random.seed!(11)

xplot = range(0, 1, length=101)


f(x) = x*sin(x/0.3)

xi = range(0, 1, length=11)
yi = f.(xi) + 0.1*randn(size(xi))
xii = xi[2:end-1]

it = interpolate(yi, BSpline(Cubic(Line(OnGrid()))))

plot()
plot!(xplot, f.(xplot), label="f")
scatter!(xi, yi, label="xi")
plot!(xi[1:length(xi)], [it[i] for i in 1:length(xi)], markershape=:diamond, label="interp")

I changed the value of λ and k directly in prefiltering.jl (thanks Revise.jl !) and I get smoothed interpolation that looks as expected on this example:

• k=2, λ=0.01
  

• k=2, λ=1
  

• k=2, λ=100
  

• k=2, λ=10000
  

• k=3, λ=0.01
  

• k=3, λ=1
  

• k=3, λ=100
  

• k=3, λ=10000
  

[nluetts]

It is great that this problem was fixed with #348, but I think the docstring of Interpolations.interpolate should say something about "smoothing" when addressing lambda. I was searching the docs for "smoothing" and could not find the relevant info. Only digging into the issues revealed the solution. If you like I could add one or two sentences, but it is perhaps a too small change for its own PR.

[mkitti]

It is great that this problem was fixed with #348, but I think the docstring of Interpolations.interpolate should say something about "smoothing" when addressing lambda. I was searching the docs for "smoothing" and could not find the relevant info. Only digging into the issues revealed the solution. If you like I could add one or two sentences, but it is perhaps a too small change for its own PR.

A separate issue or PR would be appreciated.