Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[Feature request] Smoothing factor #254

Closed
getzze opened this issue Oct 11, 2018 · 6 comments · Fixed by #348
Closed

[Feature request] Smoothing factor #254

getzze opened this issue Oct 11, 2018 · 6 comments · Fixed by #348
Labels

Comments

@getzze
Copy link
Contributor

getzze commented Oct 11, 2018

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
Copy link
Contributor

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

@getzze
Copy link
Contributor Author

getzze commented Oct 11, 2018

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
Copy link
Contributor

tomasaschan commented Oct 11, 2018

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
Copy link
Contributor Author

getzze commented Nov 22, 2019

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
    regularization_k=2_λ=0 01

  • k=2, λ=1
    regularization_k=2_λ=1

  • k=2, λ=100
    regularization_k=2_λ=100

  • k=2, λ=10000
    regularization_k=2_λ=10000

  • k=3, λ=0.01
    regularization_k=3_λ=0 01

  • k=3, λ=1
    regularization_k=3_λ=1

  • k=3, λ=100
    regularization_k=3_λ=100

  • k=3, λ=10000
    regularization_k=3_λ=10000

@nluetts
Copy link

nluetts commented Jul 5, 2021

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
Copy link
Collaborator

mkitti commented Jul 5, 2021

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

5 participants