🔢 mathematics

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

This contribution presents a new, highly efficient method for approximating real-valued functions on the positive real axis using complex exponential sums, which is derived from a multi-point Padé approximation of the Laplace transform via continued fractions and validated by various numerical examples.

Original authors: Alexey Kuznetsov, Armin Mohammadioroojeh

Published 2026-05-05

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Alexey Kuznetsov, Armin Mohammadioroojeh

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you have drawn a complex, wavy curve on a sheet of paper. This curve represents a real function, such as the movement of a stock price, the growth of a population, or the shape of a bell curve. Now imagine you must describe this curve but are allowed to use only a very specific set of building blocks: exponential decay curves. These are the smooth, falling shapes you see when a cup of coffee cools down or when a radioactive substance decays.

The challenge is: How do you stack a few of these simple "cooling coffee" curves together to perfectly mimic your complex, wavy target curve?

This is exactly what the work of Alexey Kuznetsov and Armin Mohammadioroojeh addresses. They have developed a new, highly efficient "recipe" for creating these approximations. Here is how their method works, broken down into simple concepts.

1. The Problem: The "Translation" Difficulty

The authors begin with a difficult problem. Trying to fit a complex curve directly with simple exponential pieces is like trying to solve a 3D puzzle blindfolded. It is chaotic and computationally intensive.

However, there is a magical tool in mathematics: the Laplace transform. You can think of this as a translator or a magic lens.

If you view your complex curve through this lens, it transforms into a different form called a "rational function" (essentially one polynomial divided by another).
The authors recognized an important secret: If you can approximate the "translated" version well, the original version will also be approximated well.

Instead of trying to fit the complex curve directly, they decided to fit the "translated" version first.

2. The Solution: The "Multi-Point" Strategy

To fit the translated version, they use a technique called Padé approximation. Imagine trying to guess the shape of a hidden object by touching it at specific points.

The old way: You might touch it at only one or two places (such as at the origin and at infinity).
The new way: The authors say, "Let us touch it at many specific places simultaneously."

They select a collection of points in a specific pattern (a V-shape or a fan shape in the complex plane) and force their mathematical model to match the target exactly at all these points simultaneously. This is called Multi-Point Padé approximation.

3. The Engine: The "Continued Fraction" Machine

Once they have decided where to touch the curve, they must actually build the model. They use a mathematical tool called a continued fraction.

Think of a continued fraction like a Russian Matryoshka doll or a stack of pancakes. You start with a base layer, then add a layer on top, then another, and so on.
The authors use a specific, fast algorithm (based on this nesting structure) to calculate the required exact layers.
Why is this cool? Most other methods require solving a huge, chaotic system of equations (like untangling a large ball of cables). The continued fraction method is like a conveyor belt; it builds the solution step by step, very quickly and without getting stuck.

4. The Result: Back to Reality

Once the "Russian Matryoshka" (the rational function) is built, the authors reverse the process:

They decompose the rational function (like taking apart the Matryoshka dolls).
This provides them with the exact list of exponential curves (the "cooling coffee" shapes) and their specific weights.
When you add these weighted curves together, you get a near-perfect copy of the original complex function.

5. What They Tested It On

The authors did not just talk about theory; they tested their recipe on a variety of "stubborn cases":

The Gaussian (Bell Curve): A fundamental shape in statistics.
Financial "Hockey Stick" functions: Shapes used to model financial payouts that run flat and then rise steeply.
Probability distributions: Complex shapes used to model human mortality (Gompertz–Makeham) and income distributions (Lognormal).
Step functions: Sudden jumps from zero to one, like flipping a light switch.

In almost every case, their method delivered a more accurate result with fewer building blocks (fewer exponential terms) than previous famous methods (such as the Beylkin-Monzon method).

6. The "Secret Sauce" Parameters

The method has a few "knobs" that the user can adjust (parameters $A$ , $B$ , $M$ , and $n_\infty$ ).

$M$ (Number of terms): If you increase this, you add more building blocks, making the fit tighter.
$n_\infty$ (Behavior near zero): This knob controls how well the approximation fits the curve right at the beginning (near zero).
$A$ and $B$ (The shape of the touch points): These determine the V-shaped pattern where the authors "touch" the translated curve. The authors found that adjusting these knobs allowed them to fix a common problem: sometimes the mathematics produces coefficients (weights) that are astronomically large, leading to computer errors. By adjusting $B$ , they could keep the weights small and stable without losing accuracy.

Summary

In short, the work presents a new, fast, and clever method for breaking down complicated mathematical shapes into simple, manageable parts. By translating the problem into a different mathematical "language" (the Laplace domain), using a multi-point "touching" strategy, and building the solution with a fast "Matryoshka" algorithm, they can reconstruct complex functions with high precision and with fewer parts than before.

Note: The work focuses exclusively on the mathematical construction and numerical testing of these approximations. It does not claim to solve specific medical problems or predict future market crashes, but rather offers a powerful new tool for anyone who needs to efficiently compute or simulate this type of function.

Technical Summary: Approximation of Functions on $\mathbb{R}_+$ by Exponential Sums

Problem Statement
The paper addresses the problem of approximating real-valued functions $f$ , defined on the interval $[0, \infty)$ , by linear combinations of exponential functions with complex coefficients:
$\phi(x) = \sum_{j=1}^{M} c_j e^{-\lambda_j x}, \quad c_j, \lambda_j \in \mathbb{C}.$
Such approximations are crucial for fast algorithms, separated representations of kernels and Green's functions, as individual exponential terms are analytically tractable (easy to integrate, convolve, and separate). While historical methods like the Prony method and modern approaches such as the Beylkin-Monzón algorithm (based on Hankel matrices and SVD) exist, the authors propose a new algorithm to improve efficiency and accuracy, particularly for functions whose Laplace transforms are not known in closed form.

Methodology
The core of the proposed method is a multi-point Padé approximation of the Laplace transform of the target function $f$ . The approach relies on the theoretical relationship between the $L_p$ approximation of a function $f$ and the $L_\infty$ approximation of its Laplace transform $F(z)$ .

The algorithm proceeds as follows:

Formulation: The goal is the construction of a rational function $R(z)$ that approximates the Laplace transform $F(z)$ . If $R(z)$ is a rational function with simple poles at $-\lambda_j$ , its partial fraction decomposition yields the coefficients $c_j$ and exponents $\lambda_j$ for the exponential sum $\phi(x)$ .
Constraints: The rational approximant $R(z)$ $R (z)$ is constructed to satisfy two types of conditions:
- Interpolation at finite points: $R(z_j) = F(z_j)$ for a set of $p$ points $\{z_j\}$ in the right half-plane ( $\text{Re}(z) \ge 0$ ).
- Asymptotic expansion at infinity: $R(z)$ agrees with the first $n_\infty$ terms of the asymptotic expansion of $F(z)$ as $z \to \infty$ . This expansion is derived from the Taylor series of $f(x)$ at $x=0$ via Watson's lemma.
Parameter Selection: The interpolation points $\{z_j\}$ are chosen on a piecewise linear curve $\gamma_{A,B}$ in the complex plane, defined by two parameters $A \ge 0$ and $B > 0$ . This curve consists of two segments connecting $A-Bi$ to $0$ and $0$ to $A+Bi$. This specific geometry is chosen to control the error in a sector $S_\theta$ of the complex plane, exploiting the maximum modulus principle.
Solution via Continued Fractions: To solve the resulting multi-point Padé approximation problem, the authors apply a continued fraction algorithm. This recursive method constructs the rational function $R(z)$ by iteratively reducing the number of interpolation points and updating the coefficients of the asymptotic expansion. This approach avoids solving large systems of linear equations and offers speed and simplicity, though it requires the assumption that the leading coefficient of the expansion at infinity ( $\xi_0$ ) is non-zero (a workaround is provided for the case $\xi_0 = 0$ ).
Recovery: Once $R(z)$ is found, its poles are computed and a partial fraction decomposition is performed to obtain the parameters $c_j$ and $\lambda_j$ for the final exponential sum $\phi(x)$ .

Main Contributions

New Algorithm: A novel algorithm for constructing exponential sum approximations based on the multi-point Padé approximation of the Laplace transform, solved using a continued fraction technique.
Handling Non-Closed Forms: The method effectively handles functions whose Laplace transforms are not available in closed form by utilizing high-precision numerical integration (double-exponential quadrature) to evaluate $F(z_j)$ .
Parameter Control: The introduction of parameters $A$ , $B$ , $M$ (number of terms), and $n_\infty$ (order of expansion at zero) allows for fine-tuning the approximation to minimize errors in specific norms ( $L_1$ or $L_\infty$ ) and for specific ranges of $x$ .
High-Precision Implementation: The authors provide a Fortran implementation using arbitrary-precision arithmetic (100 decimal places) to manage the numerical instability inherent in root-finding for high-degree polynomials and in iterative continued fractions.

Numerical Results
The paper demonstrates the effectiveness of the method through several numerical experiments, comparing the results with existing literature (particularly the Beylkin-Monzón method) and theoretical bounds:

Gaussian Function: Achieved errors on the order of $10^{-23}$ with $M=24$ terms.
Gamma and Barnes-G Functions: The method provided approximations for the logarithms of these functions with accuracy in double and quadruple precision, using 20% and 33% fewer terms, respectively, compared to previous two-point Padé approaches.
Gompertz-Makeham Distribution: Achieved an $L_\infty$ error of $\approx 2.2 \times 10^{-12}$ with 28 terms, outperforming a previous approximation with 15 terms that yielded an error of $\approx 10^{-6}$ .
Lognormal Distribution: Successful approximation of the density function by first approximating its cumulative distribution function (to handle the vanishing derivative at zero) and subsequently differentiating the result.
Hockey-Stick and Unit Step Functions: The method achieved localized errors and smaller $L_\infty$ errors than the Beylkin-Monzón method. The authors noted that increasing the parameter $B$ effectively reduced the magnitude of the coefficients $c_j$ , thereby mitigating potential accuracy losses in applications.

Significance and Claims
The authors claim that their method represents a competitive and often superior alternative to existing techniques, particularly regarding the number of terms required to achieve a certain accuracy. The continued fraction approach is highlighted for its computational efficiency and recursive nature.

The paper maintains a modest tone regarding theoretical guarantees. The authors explicitly state that they provide no rigorous proofs concerning the convergence or accuracy of the method. Instead, the significance is supported by "numerous numerical experiments" indicating that the algorithm works well and justifies further study. The work is presented as a practical, efficient tool for constructing exponential sum approximations, with code and results made publicly available for verification and reuse. Future work is suggested to investigate various geometries for interpolation points and to extend the method to functions whose Laplace transforms exist only in a shifted half-plane.