I'm trying to perform the inverse Laplace transform of a (noisy) dataset $y_i$ using Tikhonov regularization:

$$\min \sum_{i=1}^{N} \left(\int_0^\infty e^{-s_i t} f(t) \, dt - y_i \right)^2 - \lambda^2 \int_0^\infty f(t)^2 dt$$

where I know that $f(t) \geq 0$ for $t\geq 0$. I therefore convert that to a Non-negative least squares problem following the answer here with

$$\begin{aligned} A_{ij} &= e^{-s_i t_j} \sqrt{\Delta t} \\ b_i&=y_i \end{aligned} $$ where $\mathbf{A}$ is an $N \times L$ matrix.

Analytically we can calculate that $$ M(s_i, s_j) = \int_0^\infty e^{-s_i t}e^{-s_j t} dt = \frac{1}{s_i+s_j} \iff \mathbf{M} = \mathbf{A} \mathbf{A}^T$$ and numerically I get this up to a precision of $10^{-3}$ if I choose my $t$-axis with $L=1000$ points evenly spaced in log-space between $10^{-10}$ and $10^{10}$ (which is much better than with linear spacing).

My problem is this: I would like to include a diagonal weights matrix $W_{ii} = 1/\sigma_{y_i}^2$ as

$$\min \sum_{ij} \left(\int_0^\infty e^{-s_i t} f(t) \, dt - y_i \right) W_{ij} \left(\int_0^\infty e^{-s_j t} f(t) \, dt - y_j \right) - \lambda^2 \int_0^\infty f(t)^2 dt$$

and thought that I could just update $\mathbf{A}$ and $\mathbf{b}$ to $$\begin{aligned} \mathbf{A}' \to \mathbf{W}^{1/2} \mathbf{A} \\ \mathbf{b}' \to \mathbf{W}^{1/2} \mathbf{b} \end{aligned} $$ as described e.g. here.

However, if I now calculate $W_{ij} M(s_i, s_j)\longleftrightarrow \mathbf{A}' (\mathbf{A}')^T$ analytically and numerically I get different answers. This is the upper 4x4 part for analytical and numerical respectively:

array([[2.82777399e-06, 3.92720333e-06, 4.59359056e-06, 5.08745371e-06],
   [3.92720333e-06, 6.13584813e-06, 7.65547638e-06, 8.83179899e-06],
   [4.59359056e-06, 7.65547638e-06, 9.94943930e-06, 1.18061968e-05],
   [5.08745371e-06, 8.83179899e-06, 1.18061968e-05, 1.43013249e-05]])


 array([[2.82752363e-06, 1.88501575e-06, 1.41376182e-06, 1.13100945e-06],
   [8.18040652e-06, 6.13530489e-06, 4.90824391e-06, 4.09020326e-06],
   [1.49228376e-05, 1.19382701e-05, 9.94855841e-06, 8.52733578e-06],
   [2.28800939e-05, 1.90667449e-05, 1.63429242e-05, 1.43000587e-05]])

So only the diagonal is the same, but the off-diagonal terms are wrong and what's worse, the result is no longer symmetric. So what am I doing wrong here? Any help is much appreciated!

  • $\begingroup$ Just want to point out that there are efficient quadrature rules for accurately evaluating the Laplace integrals numerically. Look into Gauss-Laguerre, or alternatively, generalized Gaussian quadratures. There should be some implementations freely available on Netlib. $\endgroup$ – smh Jun 17 '20 at 15:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.