I want to conduct an Inverse Laplace Transform (ILT) on sampled data.
I am aware of InverseLaplace.jl, but sadly it doesn’t work well on discrete (sampled) data due to a missing interface and regularization. So I decided to give my terribly limited skills in numerics a try and do this on my own (I appreciate any advice).
The problem is that I want to solve
for F(s). A discetized version wold be:
with \Delta s_j = s_{j+1}s_j. So that we can write this as matrix equation
And I can append a regularization to the columns of K an y:
and solve this by a simple s=K\y
This is what I came up with so far:
#Regularizet inverse laplace transform
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 yoffset 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
L= Matrix{Float64}(I(length(sc)))
AR = vcat(A, α * L)
# add entry to y to store the regularization
yr = vcat(y, zeros(size(L,1 )))
# solve
sy = (AR \ yr)
(s[1:N], sy[1:N])
end
A short test shows that it actually works somehow:
using QuadGK
using Plots
dl(s;d=1,σ=1)= 1/sqrt(2*π*σ)*exp((sd)^2/(2*σ^2))
l(t;d=1,σ=1) = quadgk(s> (dl(s;d=2*d,σ=σ)+dl(s;d=d,σ=σ))*exp(t*s),0,Inf,rtol=1e6 )[1]
#generate test data:
t=0:0.1:10
y=l.(t;d=1,σ=0.1)
noise = rand(length(y)) * 0.0001
plot(s,dl.(s;d=1,σ=0.1)+dl.(s;d=2,σ=0.1),label="truth",xaxis=:log)
plot!(rilt(t, y, 0.01, 100, 250, α=0)...,label="ILT w/o noise",xlabel="s")
plot!(rilt(t, y .+ noise, 0.01, 100, 250 , α=2e5)...,xaxis=:log,label="ILT w noise + reg")
tldr;
Now the actual questions (probably not the last ones):

How can I constrain the solution to positive values only (this would likely make the ILT much better for my use case).

I found a paper with the title Stabilization of the inverse Laplace transform of multiexponential decay through introduction of a second dimension . However, it is formulated too vaguely for my knowledge. Can anyone point me to a correct way of implementing it?

Finding a proper regularization parameter can be determined from the Lcurve by the point of maximum curvature. How to find it efficiently?
I’m pretty sure there are a lot of publications out there to my questions, but I have no overview in this area and am often completely lost in the face of all the topic related vocabulary .
Surely I would also be happy to be pointed to a package that does all this for me, and what I have overlooked so far… thanks so far…