🔢 mathematics

Approximation by short exponential sums with geometric error decay based on Gauss quadrature

This paper introduces a numerically stable method based on Gauss-Laguerre and Gauss-Hermite quadratures to construct short exponential sum approximations with geometric error decay for the functions $1/(a+x)$ and $e^{-x^2/2\sigma}$ , which are subsequently applied to achieve highly accurate approximations of $\log(x)$ and the error function.

Original authors: Gerlind Plonka, Yannick Riebe, Annie Cuyt

Published 2026-06-03

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Gerlind Plonka, Yannick Riebe, Annie Cuyt

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 are trying to describe a very smooth, curvy mountain range using only a few straight lines. It's impossible to get it perfect with just one line, but if you break the mountain into smaller and smaller sections, you can use a short, simple line for each section to get a picture that looks almost exactly like the real thing.

This paper is about a clever mathematical trick to do exactly that, but with specific types of "mountains" (mathematical functions) that are crucial for science and engineering. The authors, Gerlind Plonka, Yannick Riebe, and Annie Cuyt, have developed a method to approximate two very important shapes:

The Reciprocal Curve ( $1/(a+x)$ ): Think of this as a slide that starts high and slowly flattens out as it goes to the right.
The Bell Curve ( $e^{-x^2/2\sigma}$ ): This is the famous "Gaussian" or "Normal Distribution" shape, like a smooth hill that rises and falls symmetrically.

The Problem: One Size Doesn't Fit All

Usually, when mathematicians try to approximate these curves with "exponential sums" (which are just adding up a few simple exponential waves), they run into a problem. If you try to make the approximation perfect for the whole curve at once, you need a huge number of waves, which is slow and messy to compute. If you use too few, the error (the difference between the real curve and your approximation) gets huge in some places.

The Solution: The "Zoom-In" Strategy

The authors propose a "divide and conquer" strategy. Instead of trying to cover the entire infinite curve with one set of rules, they chop the curve into consecutive intervals (segments).

For the Reciprocal Curve: The segments get longer and longer as you move out, like zooming out on a camera.
For the Bell Curve: The segments are all the same length, like tiles on a floor.

On each specific segment, they use a very short list of exponential waves (just 8 or 10 of them) to create a near-perfect copy of the curve.

The Secret Sauce: Gauss Quadrature

How do they find the perfect waves for each segment? They use a tool called Gauss Quadrature.

Imagine you are trying to measure the area under a curve. You could take a million tiny slices and add them up, but that takes forever. Gauss Quadrature is like a magic ruler that tells you exactly where to place just a few "measurement points" (nodes) and how much weight to give each one so that the sum is incredibly accurate.

The authors realized that these specific curves ( $1/(a+x)$ and the Bell Curve) can be rewritten as integrals (areas under curves) that are perfectly suited for these magic rulers:

For the Reciprocal Curve, they use Gauss-Laguerre quadrature.
For the Bell Curve, they use Gauss-Hermite quadrature.

The "Geometric Decay" Promise

The most exciting part of their discovery is the error guarantee.

Usually, when you approximate something, the error might drop slowly. The authors prove that their method makes the error drop geometrically. This means if you double your effort (or just tweak your parameters), the error doesn't just get a little smaller; it gets exponentially smaller.

They can promise you that if you choose a specific setting (called $\rho$ ), the error will be smaller than $\rho^{-2N}$ . In plain English: "If you use a short list of 10 waves, the mistake will be so tiny it's practically invisible, even on a computer."

Why This Matters (According to the Paper)

The authors show that this method is not just theoretically sound but also computationally cheap.

Stability: Once you calculate the "magic numbers" (nodes and weights) for the first segment, you can easily transform them to work for all the other segments. You don't have to start from scratch every time.
Precision: They tested this with double-precision arithmetic (the standard for high-end computing) and achieved errors as small as $10^{-15}$ . That is roughly the difference between the width of a human hair and the width of the Earth.

Real-World Applications Mentioned

The paper explicitly demonstrates two major applications of this technique:

Approximating the Logarithm ( $\log(x)$ ): By integrating their approximation of the reciprocal curve, they create a highly accurate, short exponential sum for the logarithm function.
Approximating the Error Function ( $\text{erf}(x)$ ): By integrating their approximation of the Bell Curve, they create a highly accurate way to calculate the error function, which is vital in statistics and signal processing.

Summary

Think of this paper as inventing a new type of Lego set. Instead of trying to build a whole castle with one giant, complex piece, they give you a set of small, simple, standard bricks. They show you exactly how to snap these bricks together on different sections of the castle so that the final result is indistinguishable from the real thing, and they prove mathematically that the gaps between the bricks are so small they don't even matter.

This allows scientists to use these complex functions in calculations (like simulating quantum chemistry or analyzing signals) much faster and with incredible precision, without needing supercomputers to handle the heavy lifting.

Technical Summary: Approximation by Short Exponential Sums with Geometric Error Decay Based on Gauss Quadrature

Problem Statement
The paper addresses the challenge of approximating two specific analytic functions using short exponential sums of a fixed, user-defined length $N$ :

The reciprocal function $f_1(x) = \frac{1}{a+x}$ for $a > 0$ on the domain $[0, \infty)$ .
The Gaussian function $f_2(x) = e^{-x^2/2\sigma}$ for $\sigma > 0$ on the domain $\mathbb{R}$ .

The primary objective is to achieve a geometric error decay of the form $\rho^{-2N}$ (where $\rho > 1$ is a user-defined parameter) while maintaining numerical stability. The authors note that achieving this decay rate simultaneously for all $x \in [0, \infty)$ is impossible with a single fixed exponential sum. Consequently, the problem is reformulated to construct approximations over consecutive intervals $[b_j, b_{j+1})$ , where the interval lengths depend on the desired decay parameter $\rho$ .

Methodology
The proposed method relies on representing the target functions via parametric integral representations that allow for the direct application of Gaussian quadrature formulas.

For $f_1(x)$ : The authors utilize the integral representation $\frac{1}{x+a} = \frac{1}{ab} \int_0^\infty e^{-s(\frac{x+a}{ab}-1)} e^{-s} ds$ . They apply Gauss-Laguerre quadrature to this integral. A key theoretical contribution is the derivation of a new, sharp error estimate (Theorem 2.1) tailored to the specific structure of the integrand, which improves upon standard Peano kernel estimates.
For $f_2(x)$ : The Gaussian is represented using the Fourier transform and the identity $e^{-x^2/2\sigma} = \frac{1}{\sqrt{\pi}} \int_{-\infty}^\infty e^{-u^2} \cos(\frac{x}{\sqrt{\sigma}}\sqrt{2}u) du$ . The authors apply Gauss-Hermite quadrature to this integral. Similar to the reciprocal case, they derive specific error bounds based on the expansion of the integrand into Hermite polynomials.

Key Contributions and Construction

Interval Strategy: The method partitions the domain into consecutive intervals $[b_j, b_{j+1})$ $[b_{j}, b_{j + 1})$ .
- For $f_1$ , the intervals grow exponentially: $b_{j+1} - b_j = \frac{2a}{\rho-1} (\frac{\rho+1}{\rho-1})^j$ .
- For $f_2$ , the intervals are equidistant: $b_{j+1} - b_j = \frac{2}{\rho} \sqrt{\frac{N\sigma}{e}}$ .
Parameter Recurrence: A significant feature of the construction is that the parameters (weights and exponents) for the exponential sum on any interval $j > 0$ are easily computed from the parameters on the initial interval $[b_0, b_1)$ via simple coordinate transformations. This ensures the method is computationally efficient and numerically stable without requiring the re-computation of quadrature nodes and weights for every interval.
Error Control: The technique guarantees a consistent relative error structure across all intervals. The absolute error is bounded by $c^{(j)} \rho^{-2N}$ , where the sequence $c^{(j)}$ decays exponentially, ensuring the approximation remains highly accurate even as $x$ increases.
Applications:
- Logarithm: By integrating the approximation of $\frac{1}{t+1}$ , the authors derive a short exponential sum approximation for $\log(x)$ on $[1, \infty)$ with predictable geometric error decay.
- Error Function: By integrating the approximation of the Gaussian, they construct a highly accurate approximation for the error function $\text{erf}(x)$ .

Results and Numerical Experiments
The paper provides numerical examples for sum lengths $N=8$ and $N=10$ using double-precision arithmetic.

Reciprocal Function: For $N=8$ and $\rho=4$ , the method achieves absolute errors below $10^{-13}$ on the first interval, with improved bounds reaching $10^{-15}$ when using precomputed correction terms (Remark 2.2).
Gaussian Function: For $N=8$ and $\rho=4$ , the approximation of $e^{-x^2/2}$ yields errors in the range of $10^{-15}$ across multiple intervals.
Logarithm and Error Function: The derived approximations for $\log(x)$ and $\text{erf}(x)$ demonstrate errors reaching the range of $10^{-15}$ , effectively utilizing the full precision of double-precision arithmetic.
The authors emphasize that for $N \le 10$ , the method is sufficient to achieve machine-precision accuracy, and all required parameters can be computed with high accuracy using standard Gauss quadrature routines.

Significance and Claims
The authors claim their approach is distinct from existing methods (such as Prony's method, Remez algorithms, or differential equation approaches) because it allows the user to fix the desired error bound and the sum length $N$ in advance. Unlike methods that optimize for specific norms on compact intervals, this technique provides a systematic way to extend approximations to semi-infinite or infinite domains with guaranteed geometric error decay.

The paper asserts that the method is:

Computationally Efficient: Parameters are derived via simple transformations rather than solving nonlinear optimization problems.
Numerically Stable: It avoids the instability often associated with Prony-type methods.
Versatile: It successfully handles both the reciprocal and Gaussian functions, and by extension, their integrals (logarithm and error function).
Theoretically Grounded: The error estimates are derived analytically and are shown to be tight in numerical experiments, particularly when refined using the precomputed sums mentioned in Remark 2.2.

The work concludes that this quadrature-based strategy offers a robust alternative for representing kernels and Green's functions in applications requiring high-precision, short exponential sums.