Discrete inverse Laplace transform (constrains and regularisation)

Specific Domains Numerics
, , ,
8

Turns out that Ipopt.jl does the trick. This way I could apply the constrains via the JuMP interface, whenn additionally using the second derivative operator for regularization it actually works pretty well on discrete data even with lots of noise:reg_constrained :

using LinearAlgebra
using JuMP, Ipopt 
using QuadGK
using Plots

function rilt(t, y, smin, smax, N;α=1)

    # Log spaced sampling
    s = smin * (smax / smin).^range(0, 1, length=N + 1)
    ds = diff(s)
    # central differences (input noise will make this actually unimportant)
    sc = (s[1:end - 1] + s[2:end]) / 2

    # zero rate constant for y-offset in input data
    push!(ds, 1)
    push!(sc, 0)

    # kernel
    a(i, j) = exp(-t[i] * sc[j]) * ds[j]

    # convert a(i,j) to  matrix
    A = [a(i, j) for i = eachindex(t), j = eachindex(sc)]

    # add rows for regularization 
        # Identity
        # L= Matrix{Float32}(I(length(sc)))   
    # Second derivative regularization with smooth transitions to 0 at the edges
    L= zeros(N+2,N+1)   
    δ(i,j)= i==j ? 1 : 0    
    L= [δ(i,j)/ds[j]^2+δ(i,j+2)/ds[j]^2-2*δ(i+1,j)/ds[j]^2 for i=1:N+2, j=1:N+1]
    L[:,end] .=0  #Do not regularize on the y offset   
    
    AR = vcat(A, α * L)
    # add  entry to y to store the regularization 
    yr = vcat(y, zeros(size(L,1 )))

    model = Model(Ipopt.Optimizer)
    set_silent(model)
    @variable(model,x[1:N+1])
    @constraint(model, [i=1:N], x[i] >=0);
    @objective(model, Min, sum((AR*x-yr).^2))
    optimize!(model)
    return (sc[1:N],value.(x)[1:N])
end

dl(s;d=1,σ=1)=  1/sqrt(2*π*σ)*exp(-(s-d)^2/(2*σ^2))
l(t;d=1,σ=1) = quadgk(s-> (dl(s;d=4*d,σ=σ)+dl(s;d=d,σ=σ))*exp(-t*s),0,Inf,rtol=1e-8  )[1] 

#generate test data:
t=0:0.01:10
t=10.0 .^(-4:0.01:2)
y=l.(t;d=1,σ=0.1)

noise = (rand(length(y)) .-0.5) * 0.01 
s=10.0 .^(-1:0.001:1)
plot(s,dl.(s;d=1,σ=0.1)+dl.(s;d=4,σ=0.1),label="truth",xaxis=:log,ylims=(-0.2,2))
plot!(rilt(t, y, 0.1,   10, 200, α=2e-7)...,label="ILT w/o noise",xlabel="s")
plot!(rilt(t, y .+ noise.+1.1, 0.1, 10, 200 , α=2e-7)...,xaxis=:log,label="ILT w noise ")

The regularization parameter still needs to be set manually, but it shouldn’t be too hard to automate this
@rafael.guerra: If you don’t mind pushing all these extra dependencies into your package, I would try prepare a pull request with all this and more neatly implemented.

1 Like