Inverse Laplace Transform with Concentrated Matrix-Exponential Functions

Gábor Horváth, Illés Horváth, Miklós Telek, Salah Al-Deen Almousa, Zsófia Talyigás

Features & Comparison

There are many inverse Laplace transform procedures published in the scientific literature. These procedures evaluate the known Laplace transform expression at certain points to approximate the inverse of the Laplace transform at a given point.

There is a class of methods, the Abate-Whitt framework, where the Laplace transform is evaluated always at the same points, independent of the form of the Laplace transform. The Abate-Whitt methods, in particular the Euler and the Gaver method, are very popular and widely used. Our new method, called CME, is far superior to them.

The CME method is the ideal choice when working with probability distributions. It always results valid probabilities due to the over- and under-shooting free behavior.

How does it work?

Assume we have a function \(h\). Its value at \(T\) can be expressed as $$ h(T) = \int_0^\infty \delta(T-t) h(t)\,dt,$$ where \(\delta(t)\) is the Dirac delta function. According to the Abate-Whitt framework the \(\delta(T-t)\) function is approximated by a weighted sum of exponentials consisting of \(n\) terms, that is $$\delta(T-t)\approx f_T^n(t)= \frac{1}{T}\sum_{k=1}^n \eta_k e^{-\frac{\beta_k}{T}t}.$$ With this approximation we have that \(h(T)\) can be expressed as the linear combination of its Laplace transform at certain points, thus $$ h(T) = \int_0^\infty \delta(T-t) h(t)\,dt \approx \int_0^\infty \frac{1}{T}\sum_{k=1}^n \eta_k e^{-\frac{\beta_k}{T}t} h(t)\,dt = \frac{1}{T}\sum_{k=1}^n \eta_k h^*\big(\frac{\beta_k}{T}\big),$$ where \(h^*(s)\) is the Laplace transform of \(h\).

We can observe that \(f_T^n(t)\) is just the re-scaled version of \(f_1^n(t)=\sum_{k=1}^n \eta_k e^{-\beta_k t}\). The exact form of this function plays a fundamental role in the Laplace transform inversion, and this is the point where the different ILT methods differ: the Euler, the Gaver and the CME methods provide different \(\eta_k\) coefficients and \(\beta_k\) exponents. In the Euler and the Gaver methods the algorithm to obtain these constants is relatively simple, but the resulting \(f_1^n(t)\) functions can go down to the negative range, too, which can be a source of over- or under-shooting when inverting the Laplace transform. In case of the CME approach determining the coefficients is a complex procedure, with a numerical optimization step also involved. However, the resulting \(f_1^n(t)\) (that is the density function of a concentrated matrix-exponential distribution) is always non-negative, and is very steep (a good approximation of the Dirac delta function). Furthermore, the coefficients can be pre-computed and stored in a file, the complexity of determining them does not make the ILT process slow neither does it make the implementation difficult.

Contact & Cite

The results leading to the CME-based ILT method have been achieved incrementally.

In the following paper we have published a conjecture that in the class of matrix-exponential distributions the cosine-squared form provides the minimal coefficient of variation. The parameters of these concentrated ME distributions where determined up to order 31:

In the next paper, 10 years later, a numerical optimizazion-based procedure is presented to obtain the parameters of CME distributions up to order 95. Based on that many results it was possible to examine the tendencies. The paper has shown (experimentally) that the squared coefficient of variation or CME distributions or order \(n\) is below \(2/n^2\), while for Erlang distributions it is only \(1/n\). It became clear that very concentrated distributions can be created as a mixture of exponential functions with a relatively low number of terms:

We recognized that these concentrated distributions can be used to develop an efficient inverse Laplace transform procedure. The corresponding publication is:

The following paper of the follow-up of this research, a procedure using only real eigenvalues:

In case of any questions about the method or the implementation do not hesitate to contact us:

• Gabor Horvath: ghorvath (at) hit (dot) bme (dot) hu
• Illes Horvath: horvath (dot) illes (dot) antal (at) gmail (dot) com
• Miklos Telek: telek (at) hit (dot) bme (dot) hu

Download

The implementation of the procedure in Mathematica, Matlab and IPython environments is available to download from this github repository. There is a separate file (iltcme.json) that contains the pre-computed parameters of the CME distributions. All three kinds of implementations load this file at the first time and use it for inverse Laplace transformation.

The Mathematica implementation also supports arbitrary precision computations. The demo included in the notebook demonstrates how much the Gaver method benefits from the increased precision, but the CME at machine precision is still the best in case of all functions we tested.

Demo

In the demo below you can play around with the CME method. You can inverse Laplace transform one of our test functions, or your own function, too (matrices and complex numbers are also allowed in the expression, for the syntax see mathjs.org).

The implementation is standard Javascript, the results are computed by your browser, on-the-fly, by using machine precision arithmetic. If you want to exclude a given method from the comparison, click the legend of the plot.