Computing the unitary best approximant to the exponential function

TL;DR

This work advances computational methods for the unitary best approximation of the exponential on the imaginary axis by $(n,n)$-rational functions. It presents two robust algorithms—the interpolation-based Maehly/BRASIL approach and the AAA-Lawson method—along with a posteriori bounds on the uniform error and two a priori estimates for selecting the frequency $\omega$ to achieve a target accuracy. The interpolation-based method, especially with the combined node-correction strategy, demonstrates strong convergence and scalability to high degrees, while the AAA-Lawson approach provides an alternative that achieves unitary approximants in many settings but with less consistent uniformity control. These results have practical impact for stable time integration and matrix exponential computations where unitarity is advantageous.

Abstract

Unitary best approximation to the exponential function on an interval on the imaginary axis has been introduced recently. In the present work two algorithms are considered to compute this best approximant: an algorithm based on rational interpolation in successively corrected interpolation nodes and the AAA-Lawson method. Moreover, a posteriori bounds are introduced to evaluate the quality of a computed approximant and to show convergence to the unitary best approximant in practice. Two a priori estimates -- one based on experimental data, and one based on an asymptotic error estimate -- are introduced to determine the underlying frequency for which the unitary best approximant achieves a given accuracy. Performance of algorithms and estimates is verified by numerical experiments. In particular, the interpolation-based algorithm converges to the unitary best approximant within a small number of iterations in practice.

Figures (7)

• Figure 1: These plots show computed data for the scaling factor $\xi=\xi(n,\varepsilon)$ with $\omega = (n+1)\pi\xi$ s.t. $\omega$ satisfies the problem \ref{['eq:problemfindw']}. The lines in these plots illustrate $\xi(n,\varepsilon)$ (left) and $-\log (\xi(n,\varepsilon))$ (right) over $n$ for different values of $\varepsilon$. In particular, the lines marked by symbols '$+$, $\times$, $\Delta$, $\Diamond$, $\circ$, $*$' and '$\square$' correspond to $\varepsilon = 10^{0}, 10^{-2}, \ldots, 10^{-12}$, respectively.
• Figure 2: The solid lines in these plots show the polynomials $p_a$ (left) and $p_b$ (right) corresponding to \ref{['eq:pabformulae']}. The experimental values of $\widetilde{a}_\varepsilon$ (left) and $\widetilde{b}_\varepsilon$ (right) for $\varepsilon = 10^{-14}, \sqrt{10}\cdot 10^{-14}, 10^{-13}, \ldots, 1$ are marked by '$\circ$' symbols. The values of $p_a(\log \varepsilon)$ and $p_b(\log \varepsilon)$ for $\varepsilon=2$ are marked by '$\times$' symbols. The dashed line shows extrapolated data for $\varepsilon<10^{-14}$.
• Figure 3: The symbols '$\circ$' mark a set of error objectives $\varepsilon$ over degrees $n$ which are used to construct estimates for $\omega$ with the aim that the respective unitary best approximant attains an error of $\varepsilon$. The errors of the computed unitary best approximants using the estimates $\omega_e(n,\varepsilon)$\ref{['eq:est2w']} and $\omega_a(n,\varepsilon)$\ref{['eq:asymestw']} for $\omega$ are marked by '$+$' and '$\times$', respectively. Marks close to each other correspond to the same value of $\varepsilon$, and marks are neglected in case the error computed for $\omega_e$ or $\omega_a$ is too far from the error objective. The dashed line shows $10^{-2(n-4)/3}$ over $n$.
• Figure 4: These plots illustrate the performance of the interpolation-based algorithm for $n=32$ and different frequencies $\omega$ over the number of iterations. The top and bottom rows show the error $\|r-\exp(\omega\cdot)\|$ and the error in uniformity $\delta$, respectively, for the interpolation-based algorithm using the combined strategy (left column) and the strategy of the BRASIL algorithm (right column) for interpolation nodes correction. The combined strategy is introduced in Subsection \ref{['subsec:bestalgorithm']}, and for the present example, only consists of the strategy of Maehly's second method. In each plot different lines show results for the different frequencies $\omega$ provided in Table \ref{['table:wanderror']}, i.e., $\omega = 95.48$ '$\times$', $91.35$ '$+$', $84.16$ '$\Delta$', $77.86$ '$\Diamond$', $72.19$ '$\circ$', $67.03$ '$*$', and $62.29$ '$\square$'. The reference errors from this table are also illustrated in the plots in the top row as dotted horizontal lines.
• Figure 5: These plots show the error $\|r-\exp(\omega\cdot)\|$ and the error in uniformity $\delta$ of the interpolation-based algorithm over the number of iterations using different strategies for interpolation nodes correction similar to \ref{['fig:n32brib']}. The present figure shows results for $n=256$ and $\omega = 797.18$ '$\times$', $791.45$ '$+$', $780.93$ '$\Delta$', $771.16$ '$\Diamond$', $761.89$ '$\circ$', $753.01$ '$*$', and $744.44$ '$\square$' as provided in Table \ref{['table:wanderror']}. For $\omega=761.89$, $753.01$ and $744.44$ the combined strategy applies BRASIL's strategy in an initial phase, and for these iterations the errors are illustrated by dashed lines without marks in the plots. Otherwise, the combined strategy applies the strategy of Maehly's second method for the present examples.

Theorems & Definitions (13)

• Remark 1.1
• Remark 1.2
• Remark 2.1
• Proposition 2.2
• proof
• Proposition 2.3
• proof
• Proposition 2.4
• proof
• Proposition 2.5