# Non-negative Least Squares to perform Inverse Laplace with weights

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]])


and

 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!

• 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. – smh Jun 17 '20 at 15:28