by Tobias Jawecki

This package is still under development. Future versions might have different API. The current version is not fully documented.

Besides new code, this package also contains some routines of the baryrat package¹, and variants of the AAA-Lawson²³ method which is part of the Chebfun package⁴.

Install the package python -m pip install rexpi and run numerical examples from the /examples folder.

We consider rational approximations $$r(\mathrm{i} x) \approx \mathrm{e}^{\mathrm{i}\omega x},\quad\quad\textrm{for}\quad x\in[-1,1]\quad\textrm{and a frequency}\quad\omega>0,$$ where $r$ is a unitary $(n,n)$-rational function, i.e., $r=p/q$ where $p$ and $q$ are polynomials of degree $\leq n$ and $|r(\mathrm{i} x)|=1$ for $x\in\mathbb{R}$. Our focus is best approximations in a Chebyshev sense which minimize the uniform error $$\max_{x\in[-1,1]}| r(\mathrm{i} x) - \mathrm{e}^{\mathrm{i}\omega x} |.$$

Following a recent work⁵, the unitary best approximation to $\mathrm{e}^{\mathrm{i}\omega x}$ is uniquely characterized by an equioscillating phase error. In particular, unitary rational functions are of the form $r(\mathrm{i} x) = \mathrm{e}^{\mathrm{i}g(x)}$ for $x\in\mathbb{R}$ where $g$ is a phase function, and we introduce the phase error $g(x) - \omega x$ in addition to the approximation error $r(\mathrm{i} x) - \mathrm{e}^{\mathrm{i}\omega x}$. The phase and approximation error satisfy the relation $|r(\mathrm{i} x) - \mathrm{e}^{\mathrm{i}\omega x}| = 2|\sin((g(x) - \omega x)/2)|$, and for the case that the approximation error is non-maximal, i.e., $\max_{x\in[-1,1]}| r(\mathrm{i} x) - \mathrm{e}^{\mathrm{i}\omega x} |<2$, the best unitary rational approximation to $\mathrm{e}^{\mathrm{i}\omega x}$ corresponds to the phase function $g$ which minimizes the phase error $\max_{x\in[-1,1]} |g(x) - \omega x|$. Following⁵, the corresponding phase error $g(x)-\omega x$ equioscillates at $2n+2$ nodes $\eta_1< \eta_2< \ldots <\eta_{2n+2}$, namely, $$g(\eta_j) - \omega \eta_j = (-1)^{j+1} \max_{x\in[-1,1]}| g(x) - \omega x |,\quad\quad j=1,\ldots,2n+2.$$ In the present package we introduce numerical algorithms to compute unitary best approximations mostly based on consequences of this equioscillation property, and certainly relying on the result that the unitary best approximation is unique and non-degenerate.

As a consequence of the equioscillation property, the unitary best approximation attains $2n+1$ points of interpolation $x_1 <\ldots < x_{2n+1}$ in the sense of $r(\mathrm{i} x_j) = \mathrm{e}^{\mathrm{i}\omega x_j}$. The existence of such interpolation nodes motivates the brib algorithm introduced in the present package - a modified BRASIL algorithm which aims to compute a uniform best approximation by rational interpolation and by iteratively adapting the underlying interpolation nodes. In addition, the present package also provides a variant of the AAA-Lawson method to compute approximations to $\mathrm{e}^{\mathrm{i}\omega x}$. It has been shown recently in⁶ that rational interpolants and rational approximations constructed by the AAA-Lawson algorithm are unitary, and thus, these approaches fit well to compute unitary rational best approximations.

• The best rational interpolation based brib algorithm. This is a modified BRASIL algorithm to compute unitary best approximations to $\mathrm{e}^{\mathrm{i}\omega x}$. r,_,_ = brib(w, n) computes an $(n,n)$-rational approximation to $\mathrm{e}^{\mathrm{i}\omega x}$ for $x\in[-1,1]$ in barycentric rational form. For a start, we recommend using n=10 and w = buerrest_getw(n, tol) for the frequency $\omega$ with tol = 1e-6.

import numpy as np
import matplotlib.pyplot as plt
import rexpi
n=10
tol=1e-6
w = rexpi.buerrest_getw(n, tol)
r, _, _ = rexpi.brib(w,n)
xs = np.linspace(-1,1,5000)
err = r(1j*xs)-np.exp(1j*w*xs)
errmax = np.max(np.abs(err))
plt.plot(xs,np.abs(err),[-1,1],[errmax,errmax],':')

• linearizedLawson: A version of the AAA-Lawson method to approximate $\mathrm{e}^{\mathrm{i}\omega x}$ using Chebyshev nodes as pre-assigned support nodes. The approximation computed by the linearizedLawson algorithm is comparable with the approximation computed by the brib algorithm, however, the latter seems to perform better in various settings. A basic example for the linearizedLawson algorithm:

import numpy as np
import matplotlib.pyplot as plt
import rexpi
n=10
tol=1e-6
w = rexpi.buerrest_getw(n, tol)
r, _ = rexpi.linearizedLawson(w=w, n=n, nlawson=100,nx=1000)
xs = np.linspace(-1,1,5000)
err = r(1j*xs)-np.exp(1j*w*xs)
errmax = np.max(np.abs(err))
plt.semilogy(xs,np.abs(err),[-1,1],[errmax,errmax],':')

• Error estimation. This package contains implementations of a priori estimates buerrest, buerrest_getw, and buerrest_getn for the degree $n$ and the frequency $\omega$ s.t. the approximation error is bounded by a given tolerance.
• Stopping criteria for the brib algorithm. In case the maximal number of iterations in the birb algorithm is not reached, then this algorithm returns a rational approximation with relative deviation $\delta$ smaller than tolequi=1e-3. We use a slightly different definition of the relative deviation compared to the BRASIL algorithm. The relative deviation is defined as $$\delta = 1-\frac{\min_{j=1,\ldots,2n+2} \alpha_j}{\max_{j=1,\ldots,2n+2} \alpha_j},$$ where $\alpha_j>0$ refer to the magnitude of the phase error at non-uniform alternating points $\tau_j$ which are local maxima or minima, including boundary points $x=-1$ and $x=1$, with phase error $$g(\tau_j) - \omega \tau_j = (-1)^{j+\iota} \alpha_j,\quad\quad j=1,\ldots,2n+2,\quad \iota\in{0,1}.$$ Using the relative deviation as a stopping criteria shows to be mathematical justified. [A variant of lower bounds on the minimax error sometimes referred to de la Vallée-Poussin] Assuming we have $2n+2$ non-uniform alternating points, then the phase error of the unitary best approximation satisfies $$\varepsilon_\star = \min_{g~:~ r(\mathrm{i} x) = \mathrm{e}^{\mathrm{i}g(x)},~ r\in U_n}~ \max_{x\in[-1,1]} |g(x)-\omega x| \geq \min_{j=1,\ldots,2n+2} \alpha_j.$$ Let $\varepsilon>0$ denote the phase error of the computed approximation. Since the phase error $\varepsilon$ is either attained at a local maxima, minima or at the boundary, we have $\varepsilon = \max_{j} \alpha_j$. Thus, $$\varepsilon \geq \varepsilon_\star \geq \varepsilon (1-\delta).$$ Thus, a computed approximation with a small relative deviation $\delta$ has a phase error (and approximation error) close to the minimal achievable phase error (approximation error).
• The brib algorithm is based on rational interpolation and the Lawson method is based on minimizing a linearized error. In addition to ideas in⁶ to preserve unitarity of such approximations in computer arithmetic, we introduce algorithms to also preserve symmetry in computer arithmetic, i.e., $r(-\mathrm{i} x) = r(\mathrm{i} x)^{-1}$ for $x\in\mathbb{R}$, a property which holds true for the unitary best approximation. Symmetry and unitarity are preserved when computing weights in barycentric rational form and when computing the poles of the constructed barycentric rational function, i.e., poles are computed as real poles or complex conjugate pairs of poles by design.
• Applications of these algorithms are approximations to the scalar exponential function and approximations to exponentials of skew-Hermitian matrices or the action of such matrix exponentials.