Exponential divided differences via Chebyshev polynomials

TL;DR

This work tackles the numerical instability and high cost of evaluating high-order exponential divided differences for dynamically evolving node sets. It introduces a Chebyshev-polynomial–based algorithm that combines a Chebyshev--Bessel expansion with a direct recurrence for Chebyshev divided differences, achieving a per-evaluation cost of $O(qN)$ and a truncation length $N=\Theta(c)$ for fixed interval half-width $c$. An incremental update scheme enables ${O}(N)$ updates when inserting or removing a single node, independent of the order $q$, which is crucial for Monte Carlo trajectories. The approach also yields a streamlined representation for scaled exponentials $e^{-eta[x_0,...,x_q]}$, with $eta$-dependence entering only through modified Bessel functions, allowing efficient reuse across multiple $\beta$ values. A public C++ reference implementation accompanies the paper, and numerical experiments validate stability, accuracy, and significant practical speedups over traditional methods in large-scale, dynamic-node applications such as PMR-based quantum Monte Carlo.

Abstract

Exponential divided differences arise in numerical linear algebra, matrix-function evaluation, and quantum Monte Carlo simulations, where they serve as kernel weights for time evolution and observable estimation. Efficient and numerically stable evaluation of high-order exponential divided differences for dynamically evolving node sets remains a significant computational challenge. We present a Chebyshev-polynomial-based algorithm that addresses this problem by combining the Chebyshev-Bessel expansion of the exponential function with a direct recurrence for Chebyshev divided differences. The method achieves a computational cost of ${\cal O}(qN)$, where $q$ is the divided-difference order and $N$ is the Chebyshev truncation length. We show that $N$ scales linearly with the spectral width through the decay of modified Bessel coefficients, while the dependence on $q$ enters only through structural polynomial constraints. We further develop an incremental update scheme for dynamic node sets that enables the insertion or removal of a single node in ${\cal O}(N)$ time when the affine mapping interval is held fixed. A full \texttt{C++} reference implementation of the algorithms described in this work is publicly available.

Figures (5)

• Figure 1: Modified-Bessel ratio coefficients $R_n(c)=I_n(c)/I_0(c)$ governing the Chebyshev--Bessel expansion of $e^{cy}$ on $y\in[-1,1]$. Left: the dependence on $c$ for representative orders $n$. Right: the decay in $n$ for representative $c$ values, illustrating the concentration of weight around $n=\Theta(c)$ that drives the truncation length $N=\Theta(c)$ at fixed tolerance.
• Figure 2: Median number of Chebyshev terms required for convergence as a function of interval half-width $c$ for various fixed divided difference orders $q$. Lines show the median (Q2) over 100 random trials per $(c,q)$ combination, with shaded regions indicating the interquartile range (Q1-Q3). The linear scale presentation reveals that convergence requirements grow sub-linearly with interval size, with high-order curves ($q \geq 50$) exhibiting particularly weak $c$-dependence. For small intervals ($c \lesssim 1$), all curves converge rapidly with $N_{\mathrm{terms}} \approx q + 10-15$ nearly independent of $c$. At large intervals ($c \gtrsim 50$), the curves fan out but maintain modest growth, with the highest tested order ($q=100$) requiring only 170 terms even at $c=200$. The narrow quartile ranges demonstrate excellent consistency across random node configurations.
• Figure 3: Convergence regime analysis demonstrating weak parameter dependencies in extreme limits. (a) Large-interval regime ($c \gg q$): For fixed large intervals ($c=50, 100, 200$), convergence exhibits mild $q$-dependence. A 30-fold increase in divided difference order ($q: 1 \to 30$) produces only 9--16% increases in required terms. The flat curves reflect a saturation effect where slow Bessel ratio decay at large $c$ dominates the convergence requirement, rendering the method nearly insensitive to derivative order in this regime. (b) High-order regime ($q \gg c$): For fixed high orders ($q=50, 70, 100$), convergence exhibits mild $c$-dependence. A 100-fold increase in interval size ($c: 0.1 \to 10$) produces only 10--22% increases in required terms, with weaker sensitivity at higher $q$. The near-horizontal curves confirm that when $q$ is large, the triangular structure of the Chebyshev recurrence (which ensures $D_n[q]=0$ for $n<q$) dominates convergence behavior independent of interval size.
• Figure 4: Median number of Chebyshev terms required for convergence as a function of divided-difference order $q$ for various fixed interval half-widths $c$. Lines show the median (Q2) over 100 random trials per $(c,q)$ combination, with shaded regions indicating the interquartile range (Q1--Q3). The linear-scale presentation reveals an approximately linear relationship. For larger intervals ($c \ge 50$), the curves exhibit approximately linear growth with a subleading square-root correction, reflecting weaker $q$-dependence in the regime where interval size dominates convergence behavior. The small quartile spread across all parameter combinations confirms that convergence is determined primarily by $(c,q)$ rather than specific node positions.
• Figure 5: Median relative difference between MNP and Chebyshev normalized methods for ratios of exponential divided differences. The heatmap shows $|\text{ratio}_{\text{Cheby}}-\text{ratio}_{\text{MNP}}| /|\text{ratio}_{\text{Cheby}}|$ across $(c,q)$ parameter space, where $c$ is the interval half-width and $q$ is the divided-difference order. Colors indicate agreement on a logarithmic scale, from green ($<10^{-12}$) to red ($>10^{-1}$). Red crosses denote MNP failures, defined as cases where the algorithm produces a negative or non-finite result. Each cell represents the median over 100 random trials.

Theorems & Definitions (6)

• Proposition 2.1: Basic properties of divided differences
• proof
• Theorem 2.2: Chebyshev divided-difference recurrence
• proof
• Theorem 6.1: Chebyshev--Bessel representation of scaled exponential divided differences
• proof