Approximating functions on ${\mathbb R}^+$ by exponential sums

Alexey Kuznetsov; Armin Mohammadioroojeh

arxiv: 2508.19095 · v2 · submitted 2025-08-26 · 🧮 math.NA · cs.NA

Approximating functions on {mathbb R}^+ by exponential sums

Alexey Kuznetsov , Armin Mohammadioroojeh This is my paper

Pith reviewed 2026-05-18 21:05 UTC · model grok-4.3

classification 🧮 math.NA cs.NA

keywords exponential sumsLaplace transformPadé approximationcontinued fractionsfunction approximationnumerical analysisGaussian densityprobability densities

0 comments

The pith

The multi-point Padé approximation of the Laplace transform using continued fractions produces accurate exponential sum approximations to functions on the positive reals.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper develops a technique to approximate functions on the positive real line as sums of exponentials with complex coefficients. It relies on building a rational approximation to the Laplace transform at multiple points using continued fractions, then inverting that rational function to get the exponential sum. The method is tested on the Gaussian function, lognormal and Gompertz-Makeham probability densities, the hockey-stick and unit step functions, and functions from gamma and Barnes G-function approximations. Readers in numerical analysis might value this because exponential sums often allow closed-form solutions for integrals or faster computations in applications involving positive quantities.

Core claim

We present a new method for approximating real-valued functions on R^+ by linear combinations of exponential functions with complex coefficients. The approach is based on a multi-point Padé approximation of the Laplace transform and employs a highly efficient continued fraction technique to construct the corresponding rational approximant. We demonstrate the accuracy of this method through a variety of examples, including the Gaussian function, probability density functions of the lognormal and Gompertz-Makeham distributions, the hockey stick and unit step functions, as well as a function arising in the approximation of the gamma and Barnes G-functions.

What carries the argument

Multi-point Padé approximation to the Laplace transform constructed via continued fractions, which generates a rational function whose inverse Laplace transform is the exponential sum.

Load-bearing premise

A rational approximant to the Laplace transform obtained via the continued-fraction multi-point Padé procedure will, upon inversion, yield an exponential sum whose error on R+ remains small for the listed target functions.

What would settle it

Direct comparison showing that the pointwise or integrated error of the resulting exponential sum exceeds the reported levels for the Gaussian function over a long interval on R+ would disprove the central claim.

Figures

Figures reproduced from arXiv: 2508.19095 by Alexey Kuznetsov, Armin Mohammadioroojeh.

**Figure 1.** Figure 1: Configuration of points {zj}1≤j≤p for even/odd values of p. The points are indexed so that Im(zj ) < Im(zj+1), and they satisfy the symmetry relation zp+1−j = ¯zj . (iv) Solve the multi-point Pad´e approximation problem: find a rational function R satisfying conditions (6) and (7); (v) If step (iv) is successful, compute the poles {−λj} of R. If all poles are simple, write R in partial fraction form (3) an… view at source ↗

**Figure 2.** Figure 2: The graphs of exp(−x 2 )−ϕ(x) for different values of M, A, and B. In each case we set n∞ = 2. We now turn to the results of various numerical experiments. For each example considered in the remainder of this section, the parameters cj and λj of all exponential sum approximations can be downloaded at https://kuznetsovmath.ca/. 3.1 Approximating the Gaussian function As our first example, we consider the G… view at source ↗

**Figure 3.** Figure 3: The graphs of η1(x) (black) and η2(x) (red) for different values of M, n∞, A, and B. functional equations and reflection formulas, these approximations can then be extended to compute logarithms of Γ(z) and G(z) throughout the cut plane C \ (−∞, 0]; see [16] for details. The paper [16] employed a two-point Pad´e approximation (interpolating at 0 and ∞) to produce two exponential sum approximations to f. Th… view at source ↗

**Figure 4.** Figure 4: Approximating the Gompertz–Makeham probability density function. Subplots (b) and (c) [PITH_FULL_IMAGE:figures/full_fig_p015_4.png] view at source ↗

**Figure 5.** Figure 5: The lognormal pdf gσ(x) with σ = 0.5 (black), σ = 1 (blue), and σ = 1.5 (red). where G(z) := Z ∞ 0 e −zxgσ(x) dx is the Laplace transform of gσ. Formula (45) follows directly by integration by parts. Now we can summarize the algorithm for approximating the lognormal pdf gσ(x) by exponential sums. 1. Compute the values of G(zj ) numerically using the double exponential quadrature (37). 2. Use (45) to obtain… view at source ↗

**Figure 6.** Figure 6: The graphs of gσ(x) − ϕe(x) with respect to ln(x). We set M = 30 in all cases, with n∞ = 6 in (a), (b), and n∞ = 4 in (c). 0 1 2 3 4 5 0 0.2 0.4 0.6 0.8 1 1.2 (a) h(x) (black) and ϕ(x) (red) 0 1 2 3 4 5 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 (b) h(x) − ϕ(x) [PITH_FULL_IMAGE:figures/full_fig_p017_6.png] view at source ↗

**Figure 7.** Figure 7: The graphs of h(x), ϕ(x), and h(x) − ϕ(x) for M = 5, n∞ = 2, A = 0.5, and B = 8.5. h(x) = 1 − x for 0 ≤ x ≤ 1, giving ξ0 = 1, ξ1 = −1, and ξj = 0 for all j ≥ 2. The Laplace transform of h is also available in closed form: F(z) = Z ∞ 0 e −zxh(x) dx = Z 1 0 e −zx(1 − x) dx = z −2 [PITH_FULL_IMAGE:figures/full_fig_p017_7.png] view at source ↗

**Figure 8.** Figure 8: The graphs of h(x) − ϕ(x) for different values of M, n∞, and B (A = 0 in all cases). 0 1 2 3 4 5 -6 -4 -2 0 2 4 6 10 -3 (a) M = 30, n∞ = 4, B = 83 0 1 2 3 4 5 -3 -2 -1 0 1 2 3 10 -3 (b) M = 60, n∞ = 6, B = 175 [PITH_FULL_IMAGE:figures/full_fig_p018_8.png] view at source ↗

**Figure 9.** Figure 9: The graphs of h(x) − ϕ(x) for different values of M, n∞, and B (A = 0 in all cases). 900, while the approximation shown on Figure 8c had the largest value of |cj | close to 1.5 × 106 . Such large values are undesirable since they may cause loss of precision in applications of these approximations. By experimenting with the parameters, we found that increasing B significantly reduces |cj | without substanti… view at source ↗

**Figure 10.** Figure 10: (a) The unit step function H(x) (black) and its exponential sum approximation ϕ(x) (red). (b)–(d) The errors H(x) − ϕ(x) for different values of M, B, and n∞. In all cases we set A = 0. Section 3.4). However, since H(0+) ̸= 0, our algorithm can also be applied directly to H, which we do in the following experiments. Figure 10a shows H(x) together with a 15-term exponential sum approximation ϕ(x), and Figu… view at source ↗

read the original abstract

We present a new method for approximating real-valued functions on ${\mathbb R}^+$ by linear combinations of exponential functions with complex coefficients. The approach is based on a multi-point Pad\'e approximation of the Laplace transform and employs a highly efficient continued fraction technique to construct the corresponding rational approximant. We demonstrate the accuracy of this method through a variety of examples, including the Gaussian function, probability density functions of the lognormal and Gompertz-Makeham distributions, the hockey stick and unit step functions, as well as a function arising in the approximation of the gamma and Barnes $G$-functions.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Desk Editor's Note private letter to a colleague

The continued-fraction multi-point Padé route to Laplace-transform approximants produces usable exponential-sum fits on the reported test cases, but the paper stays almost entirely empirical.

read the letter

The punchline is that this gives a practical tool for building exponential sum representations of functions on the positive line, and the examples suggest it works reliably for the cases they tested. The new element is the way they use a continued fraction to construct the multi-point Padé approximant to the Laplace transform before inverting it to get the sum of exponentials. This seems more efficient than solving large systems directly, and it is applied here to get good fits for the Gaussian, lognormal and Gompertz-Makeham densities, hockey-stick and unit step functions, plus some special function approximations. The paper does well in presenting these as concrete demonstrations rather than abstract claims. On the soft side, the work is almost entirely empirical. There are no theorems on convergence rates or error estimates that would tell us how the approximation behaves in general or for other functions. The accuracy is shown on specific instances, which is useful but leaves open questions about stability when the number of poles increases or for functions with different singularity structures. The weakest assumption is that the inverted rational approximant stays accurate on R+ without additional controls. This paper is for numerical analysts and practitioners who need exponential sum models for computation, such as in stochastic processes or integral equations. A reader interested in approximation methods for positive real functions would find the technique worth trying. I recommend sending it for peer review. The construction is internally consistent and the numerical evidence is there to evaluate, even if more analysis would help.

Referee Report

0 major / 4 minor

Summary. The manuscript presents a method for approximating real-valued functions on the positive reals by linear combinations of exponentials with complex coefficients. The construction proceeds by forming a multi-point Padé approximant to the Laplace transform of the target function via a continued-fraction algorithm, followed by partial-fraction decomposition to recover the exponential sum. Accuracy is demonstrated numerically on the Gaussian, lognormal and Gompertz-Makeham densities, the hockey-stick and unit-step functions, and functions appearing in gamma and Barnes G-function approximations.

Significance. If the reported numerical accuracies hold, the technique supplies a practical, transform-based route to exponential-sum representations that avoids direct nonlinear fitting. The continued-fraction construction of the rational approximant is efficient and leverages standard machinery, making the method reproducible for the listed classes of functions. This could be useful in probability, statistics, and special-function computations where exponential sums enable closed-form operations or fast evaluation.

minor comments (4)

The description of the continued-fraction algorithm for the multi-point Padé approximant (likely §2) would benefit from an explicit step-by-step outline or pseudocode to facilitate independent implementation.
In the numerical examples, the number of poles (or terms in the exponential sum) is chosen per function; a brief discussion of how this choice is made or a simple a-posteriori error indicator would improve clarity.
Figure captions for the error plots should state the precise norm or pointwise measure used and the range of the independent variable displayed.
A short comparison, even qualitative, with existing exponential-sum methods (e.g., Prony-type or nonlinear least-squares fits) on one or two examples would help situate the new approach.

Simulated Author's Rebuttal

0 responses · 0 unresolved

We thank the referee for the careful reading, positive assessment of the method, and recommendation for minor revision. The report correctly identifies the core construction via continued-fraction multi-point Padé approximation of the Laplace transform and the subsequent partial-fraction step. No major comments were raised, so we have no specific points to rebut or revise at this stage.

Circularity Check

0 steps flagged

No significant circularity in derivation chain

full rationale

The paper describes a standard construction: multi-point Padé approximation to the Laplace transform via continued fractions, followed by partial-fraction decomposition to obtain an exponential sum. This sequence follows classical approximation theory without any self-definitional steps, fitted parameters renamed as predictions, or load-bearing self-citations that reduce the central claim to a tautology. Concrete numerical examples for the listed target functions (Gaussian, lognormal density, etc.) serve as empirical validation rather than a forced equivalence. The derivation remains independent of its own outputs and rests on externally verifiable Laplace/Padé machinery.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The approach rests on the existence of the Laplace transform and on the convergence properties of continued-fraction Padé approximants; no free parameters, new entities, or ad-hoc axioms are introduced in the abstract.

axioms (1)

domain assumption The Laplace transform of each target function exists in a suitable half-plane and admits a rational approximation whose inverse Laplace transform is an exponential sum.
This background fact is required for the Padé step to produce the desired exponential representation.

pith-pipeline@v0.9.0 · 5626 in / 1263 out tokens · 38687 ms · 2026-05-18T21:05:32.556080+00:00 · methodology

discussion (0)

Lean theorems connected to this paper

Citations machine-checked in the Pith Canon. Every link opens the source theorem in the public Lean library.

IndisputableMonolith/Cost/FunctionalEquation.lean washburn_uniqueness_aczel unclear

?

unclear
Relation between the paper passage and the cited Recognition theorem.

The approach is based on a multi-point Padé approximation of the Laplace transform and employs a highly efficient continued fraction technique to construct the corresponding rational approximant.
IndisputableMonolith/Foundation/ArithmeticFromLogic.lean embed_eq_pow unclear

?

unclear
Relation between the paper passage and the cited Recognition theorem.

R(z) = ∑ c_j / (z + λ_j) … partial fraction decomposition

What do these tags mean?

matches: The paper's claim is directly supported by a theorem in the formal canon.
supports: The theorem supports part of the paper's argument, but the paper may add assumptions or extra steps.
extends: The paper goes beyond the formal theorem; the theorem is a base layer rather than the whole result.
uses: The paper appears to rely on the theorem as machinery.
contradicts: The paper's claim conflicts with a theorem or certificate in the canon.
unclear: Pith found a possible connection, but the passage is too broad, indirect, or ambiguous to say the theorem truly supports the claim.

Reference graph

Works this paper leans on

21 extracted references · 21 canonical work pages

[1]

O. Aberth. Iteration methods for finding all zeros of a polynomial simultaneously. Mathematics of Computation, 27(122):339–344, 1973. https://doi.org/10.2307/2005621

work page doi:10.2307/2005621 1973
[2]

D. H. Bailey. MPFUN2020: A thread-safe arbitrary precision package with special functions. 2020. https://www.davidhbailey.com/dhbsoftware/

work page 2020
[3]

G. A. Baker and P. Graves-Morris. Pad´ e Approximants, Part I: Basic Theory. Encyclopedia of Mathematics and Its Applications. Addison-Wesley Publishing, Reading, MA, 1981

work page 1981
[4]

Beylkin and L

G. Beylkin and L. Monz´ on. On approximation of functions by exponential sums. Applied and Computational Harmonic Analysis , 19(1):17–48, 2005. https://doi.org/10.1016/j.acha.2005. 01.003

work page doi:10.1016/j.acha.2005 2005
[5]

Applied and Computational Harmonic Analysis , series =

G. Beylkin and L. Monz´ on. Approximation by exponential sums revisited. Applied and Computa- tional Harmonic Analysis, 28(2):131–149, 2010. https://doi.org/10.1016/j.acha.2009.08.011

work page doi:10.1016/j.acha.2009.08.011 2010
[6]

Braess and W

D. Braess and W. Hackbusch. On the efficient computation of high-dimensional integrals and the approximation by exponential sums. In R. DeVore and A. Kunoth, editors, Multiscale, Nonlinear and Adaptive Approximation, pages 39–74, Berlin, Heidelberg, 2009. Springer Berlin Heidelberg

work page 2009
[7]

L. W. Ehrlich. A modified Newton method for polynomials. Communications of the ACM, 10(2):107– 108, 1967. https://doi.org/10.1145/363067.363115

work page doi:10.1145/363067.363115 1967
[8]

Feng and X

R. Feng and X. Jing. Analytical valuation and hedging of variable annuity guaranteed lifetime withdrawal benefits. Insurance: Mathematics and Economics , 72:36–48, 2017. https://doi.org/ 10.1016/j.insmatheco.2016.10.011

work page doi:10.1016/j.insmatheco.2016.10.011 2017
[9]

R. Feng, A. Kuznetsov, and F. Yang. Exponential functionals of L´ evy processes and variable annuity guaranteed benefits. Stochastic Processes and their Applications , 129(2):604–625, 2019. https://doi.org/10.1016/j.spa.2018.03.011

work page doi:10.1016/j.spa.2018.03.011 2019
[10]

Fernandez Rodriguez, L

A. Fernandez Rodriguez, L. de Santiago Rodrigo, E. Lopez Guillen, J. Rodriguez Ascariz, J. M. Miguel Jimenez, and L. Boquete. Coding Prony’s method in MATLAB and applying it to biomedical signal filtering. BMC Bioinformatics , 19:451, 2018. https://doi.org/10.1186/ s12859-018-2473-y. 20

work page 2018
[11]

Hackbusch

W. Hackbusch. Computation of best L∞ exponential sums for 1/x by Remez’ algorithm. Computing and Visualization in Science , 20:1–11, 2019. https://doi.org/10.1007/s00791-018-00308-4

work page doi:10.1007/s00791-018-00308-4 2019
[12]

Iscoe, K

I. Iscoe, K. Jackson, A. Kreinin, and X. Ma. An exponential approximation to the hockey stick function. preprint, 2010. http://defaultrisk.com/pp_quant_18.htm

work page 2010
[13]

D. W. Kammler. Characterization of best approximations by sums of exponentials. Journal of Approximation Theory, 9(2):173–191, 1973. https://doi.org/10.1016/0021-9045(73)90105-6

work page doi:10.1016/0021-9045(73)90105-6 1973
[14]

D. W. Kammler. Least squares approximation of completely monotonic functions by sums of ex- ponentials. SIAM Journal on Numerical Analysis , 16(5):801–818, 1979. http://www.jstor.org/ stable/2156635

work page arXiv 1979
[15]

Kunis, D

S. Kunis, D. Potts, and G. Steidl. Fast Gauss transforms with complex parameters using NFFTs. Journal of Numerical Mathematics , 14(4):295–303, 2006. https://doi.org/10.1515/ 156939506779874626

work page 2006
[16]

Kuznetsov

A. Kuznetsov. Computing the Barnes G-function and the gamma function in the entire complex plane. Journal of Computational and Applied Mathematics , 411:114270, 2022. https://doi.org/ 10.1016/j.cam.2022.114270

work page doi:10.1016/j.cam.2022.114270 2022
[17]

W. McLean. Exponential sum approximations for t−β. In J. Dick, F. Y. Kuo, and H. Wo´ zniakowski, editors, Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan, pages 911–930, Cham, 2018. Springer International Publishing

work page 2018
[18]

M. R. Osborne and G. K. Smyth. A modified Prony algorithm for exponential function fitting. SIAM Journal on Scientific Computing , 16(1):119–138, 1995. https://doi.org/10.1137/0916008

work page doi:10.1137/0916008 1995
[19]

J. R. Rice. Chebyshev approximation by exponentials. Journal of the Society for Industrial and Applied Mathematics, 10(1):149–161, 1962. https://doi.org/10.1137/0110013

work page doi:10.1137/0110013 1962
[20]

Takahasi and M

H. Takahasi and M. Mori. Double exponential formula for numerical integration. Publ. RIMS, Kyoto Univ., 9:721–741, 1974. https://doi.org/10.2977/prims/1195192451

work page doi:10.2977/prims/1195192451 1974
[21]

Wong and M

R. Wong and M. Wyman. Generalization of Watson’s Lemma. Canadian Journal of Mathematics , 24(2):185–208, 1972. 21

work page 1972

[1] [1]

O. Aberth. Iteration methods for finding all zeros of a polynomial simultaneously. Mathematics of Computation, 27(122):339–344, 1973. https://doi.org/10.2307/2005621

work page doi:10.2307/2005621 1973

[2] [2]

D. H. Bailey. MPFUN2020: A thread-safe arbitrary precision package with special functions. 2020. https://www.davidhbailey.com/dhbsoftware/

work page 2020

[3] [3]

G. A. Baker and P. Graves-Morris. Pad´ e Approximants, Part I: Basic Theory. Encyclopedia of Mathematics and Its Applications. Addison-Wesley Publishing, Reading, MA, 1981

work page 1981

[4] [4]

Beylkin and L

G. Beylkin and L. Monz´ on. On approximation of functions by exponential sums. Applied and Computational Harmonic Analysis , 19(1):17–48, 2005. https://doi.org/10.1016/j.acha.2005. 01.003

work page doi:10.1016/j.acha.2005 2005

[5] [5]

Applied and Computational Harmonic Analysis , series =

G. Beylkin and L. Monz´ on. Approximation by exponential sums revisited. Applied and Computa- tional Harmonic Analysis, 28(2):131–149, 2010. https://doi.org/10.1016/j.acha.2009.08.011

work page doi:10.1016/j.acha.2009.08.011 2010

[6] [6]

Braess and W

D. Braess and W. Hackbusch. On the efficient computation of high-dimensional integrals and the approximation by exponential sums. In R. DeVore and A. Kunoth, editors, Multiscale, Nonlinear and Adaptive Approximation, pages 39–74, Berlin, Heidelberg, 2009. Springer Berlin Heidelberg

work page 2009

[7] [7]

L. W. Ehrlich. A modified Newton method for polynomials. Communications of the ACM, 10(2):107– 108, 1967. https://doi.org/10.1145/363067.363115

work page doi:10.1145/363067.363115 1967

[8] [8]

Feng and X

R. Feng and X. Jing. Analytical valuation and hedging of variable annuity guaranteed lifetime withdrawal benefits. Insurance: Mathematics and Economics , 72:36–48, 2017. https://doi.org/ 10.1016/j.insmatheco.2016.10.011

work page doi:10.1016/j.insmatheco.2016.10.011 2017

[9] [9]

R. Feng, A. Kuznetsov, and F. Yang. Exponential functionals of L´ evy processes and variable annuity guaranteed benefits. Stochastic Processes and their Applications , 129(2):604–625, 2019. https://doi.org/10.1016/j.spa.2018.03.011

work page doi:10.1016/j.spa.2018.03.011 2019

[10] [10]

Fernandez Rodriguez, L

A. Fernandez Rodriguez, L. de Santiago Rodrigo, E. Lopez Guillen, J. Rodriguez Ascariz, J. M. Miguel Jimenez, and L. Boquete. Coding Prony’s method in MATLAB and applying it to biomedical signal filtering. BMC Bioinformatics , 19:451, 2018. https://doi.org/10.1186/ s12859-018-2473-y. 20

work page 2018

[11] [11]

Hackbusch

W. Hackbusch. Computation of best L∞ exponential sums for 1/x by Remez’ algorithm. Computing and Visualization in Science , 20:1–11, 2019. https://doi.org/10.1007/s00791-018-00308-4

work page doi:10.1007/s00791-018-00308-4 2019

[12] [12]

Iscoe, K

I. Iscoe, K. Jackson, A. Kreinin, and X. Ma. An exponential approximation to the hockey stick function. preprint, 2010. http://defaultrisk.com/pp_quant_18.htm

work page 2010

[13] [13]

D. W. Kammler. Characterization of best approximations by sums of exponentials. Journal of Approximation Theory, 9(2):173–191, 1973. https://doi.org/10.1016/0021-9045(73)90105-6

work page doi:10.1016/0021-9045(73)90105-6 1973

[14] [14]

D. W. Kammler. Least squares approximation of completely monotonic functions by sums of ex- ponentials. SIAM Journal on Numerical Analysis , 16(5):801–818, 1979. http://www.jstor.org/ stable/2156635

work page arXiv 1979

[15] [15]

Kunis, D

S. Kunis, D. Potts, and G. Steidl. Fast Gauss transforms with complex parameters using NFFTs. Journal of Numerical Mathematics , 14(4):295–303, 2006. https://doi.org/10.1515/ 156939506779874626

work page 2006

[16] [16]

Kuznetsov

A. Kuznetsov. Computing the Barnes G-function and the gamma function in the entire complex plane. Journal of Computational and Applied Mathematics , 411:114270, 2022. https://doi.org/ 10.1016/j.cam.2022.114270

work page doi:10.1016/j.cam.2022.114270 2022

[17] [17]

W. McLean. Exponential sum approximations for t−β. In J. Dick, F. Y. Kuo, and H. Wo´ zniakowski, editors, Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan, pages 911–930, Cham, 2018. Springer International Publishing

work page 2018

[18] [18]

M. R. Osborne and G. K. Smyth. A modified Prony algorithm for exponential function fitting. SIAM Journal on Scientific Computing , 16(1):119–138, 1995. https://doi.org/10.1137/0916008

work page doi:10.1137/0916008 1995

[19] [19]

J. R. Rice. Chebyshev approximation by exponentials. Journal of the Society for Industrial and Applied Mathematics, 10(1):149–161, 1962. https://doi.org/10.1137/0110013

work page doi:10.1137/0110013 1962

[20] [20]

Takahasi and M

H. Takahasi and M. Mori. Double exponential formula for numerical integration. Publ. RIMS, Kyoto Univ., 9:721–741, 1974. https://doi.org/10.2977/prims/1195192451

work page doi:10.2977/prims/1195192451 1974

[21] [21]

Wong and M

R. Wong and M. Wyman. Generalization of Watson’s Lemma. Canadian Journal of Mathematics , 24(2):185–208, 1972. 21

work page 1972