Rational minimax approximation via adaptive barycentric representations

TL;DR

This paper introduces adaptive barycentric representations for rational minimax approximation, enabling robust algorithms that outperform traditional methods, especially near singularities, demonstrated by high-degree approximations in standard floating point.

Contribution

It develops a new adaptive barycentric approach for rational minimax approximation, combining classical and iterative methods for improved robustness and efficiency.

Findings

Able to compute high-degree rational approximations in standard precision

Outperforms previous methods requiring extended precision

Achieves quadratic convergence with combined algorithms

Abstract

Computing rational minimax approximations can be very challenging when there are singularities on or near the interval of approximation - precisely the case where rational functions outperform polynomials by a landslide. We show that far more robust algorithms than previously available can be developed by making use of rational barycentric representations whose support points are chosen in an adaptive fashion as the approximant is computed. Three variants of this barycentric strategy are all shown to be powerful: (1) a classical Remez algorithm, (2) a "AAA-Lawson" method of iteratively reweighted least-squares, and (3) a differential correction algorithm. Our preferred combination, implemented in the Chebfun MINIMAX code, is to use (2) in an initial phase and then switch to (1) for generically quadratic convergence. By such methods we can calculate approximations up to type (80, 80) of $|x|$ on $[-1, 1]$ in standard 16-digit floating point arithmetic, a problem for which Varga, Ruttan, and Carpenter required 200-digit extended precision.

Tables

Table 1. Table 1: Best approximation to five difficult functions by the barycentric rational Remez algorithm. f 1 ′′ superscript subscript 𝑓 1 ′′ f_{1}^{\prime\prime} is discontinuous at x = 1 / 2 𝑥 1 2 x=1/\sqrt{2} , f 2 ′ superscript subscript 𝑓 2 ′ f_{2}^{\prime} is discontinuous at x = 0 𝑥 0 x=0 , f 3 ′ superscript subscript 𝑓 3 ′ f_{3}^{\prime} is unbounded as x → 0 → 𝑥 0 x\to 0 , f 4 subscript 𝑓 4 f_{4} has two sharp peaks at x = ± 0.6 𝑥 plus-or-minus 0.6 x=\pm 0.6 , and f 5 subscript 𝑓 5 f_{5} has a logarithmic singularity at x = 0 𝑥 0 x=0 .

$i$	$f_i$	$[a,b]$	$(m,n)$	${\Vert f-r^{\ast }\Vert }_{\infty }$
1		$[0,1]$	$(22,22)$	$2.439\times {10}^{-9}$
2	$\left|x\right|\sqrt{\left|x\right|}$	$[-0.7,2]$	$(17,71)$	$4.371\times {10}^{-8}$
3	$x^3+\frac{\sqrt[3]{x}{e}^{-x^2}}{8}$	$[-0.2,0.5]$	$(45,23)$	$2.505\times {10}^{-5}$
4	$\frac{100\pi \left(x^2-0.36\right)}{\sinh \left(100\pi \left(x^2-0.36\right)\right)}$	$[-1,1]$	$(38,38)$	$1.780\times {10}^{-12}$
5	$-\frac{1}{\log \left|x\right|}$	$[-0.1,0.1]$	$(8,8)$	$1.52\times {10}^{-2}$

Table 2. Table 2: Best type ( 16 , 16 ) 16 16 (16,16) approximations to four functions using the barycentric DC algorithm. X 𝑋 X consists of 20000 20000 20000 equispaced points inside [ − 1 , 1 ] 1 1 [-1,1] .

$i$	$f_i$	${\Vert f_i-r^{\ast }\Vert }_{X,\infty }$
1	${\sum }_{k=0}^{\infty }{2}^{-k}\cos \left(3^{k}x\right)$	$0.1377$
2	$\min \{\text{sech}\left(3\sin \left(10x\right)\right),\sin \left(9x\right)\}$	$0.0610$
3	$\sqrt{\left|x^3\right|}+\left|x+0.5\right|$	$1.2057\cdot {10}^{-4}$
4	$\left(\frac{1}{2}\text{erf}\frac{x}{\sqrt{0.0002}}+\frac{3}{2}\right)e^{-x}$	$6.2045\cdot {10}^{-6}$

Equations

$${\mathcal{R}}_{m,n}=\left\{\frac{p}{q}:p\in {\mathbb{R}}_m[x],q\in {\mathbb{R}}_n[x]\right\}.$$

$$\underset{r\in {\mathcal{R}}_{m,n}}{\min }\Vert f-r{\Vert }_{\infty },$$

$$f(x_{\ell })-r^{\ast }(x_{\ell })=(-1{)}^{\ell +1}\lambda ,\ell =0,\dots ,m+n+1-d,$$

$$r(z)=\frac{N(z)}{D(z)}=\sum_{k=0}^{n}\dfrac{\alpha_{k}}{z-t_{k}}\bigg{/}\sum_{k=0}^{n}\dfrac{\beta_{k}}{z-t_{k}},$$

$${\omega }_t(z)=\prod _{k=0}^n(z-t_k),$$

$$r(z)=\sum_{k=0}^{n}\dfrac{f(t_{k})\beta_{k}}{z-t_{k}}\bigg{/}\sum_{k=0}^{n}\dfrac{\beta_{k}}{z-t_{k}},$$

$$r(z)=\frac{{\omega }_t(z)N(z)}{{\omega }_t(z)D(z)}=\frac{{\prod }_{k=0}^n(z-t_k){\sum }_{k=0}^n\frac{{\alpha }_k}{z-t_k}}{{\prod }_{k=0}^n(z-t_k){\sum }_{k=0}^n\frac{{\beta }_k}{z-t_k}}=:\frac{p(z)}{q(z)}.$$

$$V_m=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ t_0& t_1& \cdots & t_n\\ ⋮& ⋮& & ⋮\\ t_0^{n-1-m}& t_1^{n-1-m}& \cdots & t_n^{n-1-m}\end{array}\right].$$

$$V_n=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ t_0& t_1& \cdots & t_m\\ ⋮& ⋮& & ⋮\\ t_0^{m-1-n}& t_1^{m-1-n}& \cdots & t_m^{m-1-n}\end{array}\right],$$

$$\widehat{r}(x)=\sum_{k=0}^{n}\dfrac{\alpha_{k}(1+\epsilon_{\alpha_{k}})}{x-t_{k}}\bigg{/}\sum_{k=0}^{n}\dfrac{\beta_{k}(1+\epsilon_{\beta_{k}})}{x-t_{k}},$$

$${\hat{\alpha }}_k={\alpha }_k(1+{\delta }_{{\alpha }_k}),{\delta }_{{\alpha }_k}=O({\kappa }_{\alpha }u),{\hat{\beta }}_k={\beta }_k(1+{\delta }_{{\beta }_k}),{\delta }_{{\beta }_k}=O({\kappa }_{\beta }u),k=0,\dots ,n,$$

$$\left|\frac{r(x)-\hat{r}(x)}{r(x)}\right|\le u(n+3+O({\kappa }_{\alpha }))\frac{{\sum }_{k=0}^n\left|\frac{{\alpha }_k}{x-t_k}\right|}{\left|{\sum }_{k=0}^n\frac{{\alpha }_k}{x-t_k}\right|}+u(n+2+O({\kappa }_{\beta }))\frac{{\sum }_{k=0}^n\left|\frac{{\beta }_k}{x-t_k}\right|}{\left|{\sum }_{k=0}^n\frac{{\beta }_k}{x-t_k}\right|}+O(u^2).$$

$$a\le x_0^{(k)}<\cdots <x_{m+n+1}^{(k)}\le b.$$

$$f(x_{\ell }^{(k)})-r_k(x_{\ell }^{(k)})=(-1{)}^{\ell +1}{\lambda }_k,\ell =0,\dots ,m+n+1.$$

$$s(-1{)}^{\ell }\left(f(x_{\ell }^{(k+1)})-r_k(x_{\ell }^{(k+1)})\right)\ge \left|{\lambda }_k\right|,\ell =0,\dots ,m+n+1,$$

$$w(x_{\ell }^{(k)})\left(f(x_{\ell }^{(k)})-r_k(x_{\ell }^{(k)})\right)=(-1{)}^{\ell +1}{\lambda }_k,\ell =0,\dots ,m+n+1$$

$$s(-1{)}^{\ell }w(x_{\ell }^{(k+1)})\left(f(x_{\ell }^{(k+1)})-r_k(x_{\ell }^{(k+1)})\right)\ge \left|{\lambda }_k\right|,\ell =0,\dots ,m+n+1,$$

$$f(x_{\ell })-r(x_{\ell })=(-1{)}^{\ell +1}\lambda ,\ell =0,\dots ,2n+1$$

$$p(x)=\sum _{k=0}^n{c}_{p,k}{\phi }_k(x),q(x)=\sum _{k=0}^n{c}_{q,k}{\phi }_k(x).$$

$$p(x_{\ell })=q(x_{\ell })\left(f(x_{\ell })-(-1{)}^{\ell +1}\lambda \right),$$

$${\Phi }_x{c}_p=\left(\left[\begin{array}{cccc}f(x_0)& & & \\ & f(x_1)& & \\ & & \ddots & \\ & & & f(x_{2n+1})\end{array}\right]-\lambda \left[\begin{array}{cccc}-1& & & \\ & 1& & \\ & & -1& \\ & & & \ddots \end{array}\right]\right){\Phi }_x{c}_q,$$

$$\left[\begin{array}{cc}D{\Phi }_x& -FD{\Phi }_x\end{array}\right]\left[\begin{array}{c}{c}_p\\ c_q\end{array}\right]=\lambda \left[\begin{array}{cc}0& -SD{\Phi }_x\end{array}\right]\left[\begin{array}{c}{c}_p\\ c_q\end{array}\right],$$

$$D{\Phi }_x=\left[\begin{array}{cc}{Q}_1& Q_2\end{array}\right]\left[\begin{array}{c}R\\ 0\end{array}\right]=Q_{1}R.$$

$$Q_2^{T}F{Q}_{1}R{c}_q=\lambda Q_2^{T}S{Q}_{1}R{c}_q.$$

$$V_x^T{\Omega }_x{V}_x=0,$$

$${\left(V_x^T{\Omega }_x{V}_x\right)}_{i,j}=\sum _{\ell =0}^{2n+1}{x}_{\ell }^{i+j}\frac{1}{{\omega }_x^{\prime }(x_{\ell })}=(x^{i+j})[x_0,\dots ,x_{2n+1}]=0,i,j\in \left\{0,\dots ,n\right\},$$

$${\Phi }_x^T{\Omega }_x{\Phi }_x=0.$$

$${\Phi }_x^T{\Omega }_{x}F{\Phi }_x{c}_q=\lambda {\Phi }_x^T{\Omega }_{x}S{\Phi }_x{c}_q.$$

$$Q_1^{T}SF{Q}_{1}y=\lambda y.$$

$$R{c}_p=Q_1^T(F-\lambda S){\left|{\Omega }_x\right|}^{1/2}{\Phi }_x{c}_q=Q_1^{T}F{Q}_{1}y.$$

Full text

1 Introduction

The problem we are interested in is that of approximating functions $f\in\mathcal{C}([a,b])$ using type $(m,n)$ rational approximations with real coefficients, in the $L^{\infty}$ setting. The set of feasible approximations is

[TABLE]

Given $f$ and prescribed nonnegative integers $m,n$, the goal is to compute

[TABLE]

where $\|\cdot\|_{\infty}$ denotes the infinity norm over $[a,b]$, i.e., $\|f-r\|_{\infty}=\max_{x\in[a,b]}|f(x)-r(x)|$. The minimizer of (2) is known to exist and to be unique [58, Ch. 24].

Let the minimax (or best) approximation be written $r^{*}=p/q\in\mathcal{R}_{m,n}$, where $p$ and $q$ have no common factors. The number $d=\min\left\{m-\deg p,m-\deg q\right\}$ is called the defect of $r^{*}$. It is known that there exists a so-called alternant (or reference) set consisting of ordered nodes $a\leqslant x_{0}<x_{1}<\cdots<x_{m+n+1-d}\leqslant b$, where $f-r^{*}$ takes its global extremum over $[a,b]$ with alternating signs. In other words, we have the beautiful equioscillation property [58, Theorem 24.1]

[TABLE]

where $|\lambda|=\|f-r^{*}\|_{\infty}$. Minimax approximations with $d>0$ are called degenerate, and they can cause problems for computation. Accordingly, unless otherwise stated, we make the assumption that $d=0$ for (2). In practice, degeneracy most often arises due to symmetries in approximating even or odd functions, and we check for these cases explicitly to make sure they are treated properly. Other degeneracies can usually be detected by examining in succession the set of best approximations of types $(m-k,n-k),(m-k+1,n-k+1),\ldots,(m,n)$ with $k=\min\left\{m,n\right\}$ [11, p. 161].

In the approximation theory literature [15, 50, 63, 11, 40], two algorithms are usually considered for the numerical solution of (2), the rational Remez and differential correction (DC) algorithms. The various challenges that are inherent in rational approximations can, more often than not, make the use of such methods difficult. Finding the best polynomial approximation, by contrast, can usually be done robustly by a standard implementation of the linear version of the Remez algorithm [47]. This might explain why the current software landscape for minimax rational approximations is rather barren. Nevertheless, implementations of the rational Remez algorithm are available in some mathematical software packages: the Mathematica MiniMaxApproximation function, the Maple numapprox[minimax] routine and the MATLAB Chebfun [24] remez code. The Boost C++ libraries [1] also contain an implementation.

Over the years, the applications that have benefited most from minimax rational approximations come from recursive filter design in signal processing [23, 13] and the representation of special functions [18, 19]. Apart from such practical motivations, we believe it worthwhile to pursue robust numerical methods for computing these approximations because of their fundamental importance to approximation theory. A new development of this kind has already resulted from the algorithms described here: the discovery that type $(k,k)$ rational approximations to $x^{n}$, for $n\gg k$, converge geometrically at the rate $O(9.28903\cdots^{-k})$ [44].

In this paper we present elements that greatly improve the numerical robustness of algorithms for computing best rational approximations. The key idea is the use of barycentric representations with adaptively chosen basis functions, which can overcome the numerical difficulties frequently encountered when $f$ has nonsmooth points. For instance, when trying to approximate $f(x)=|x|$ on $[-1,1]$ using standard IEEE double precision arithmetic in MATLAB, our barycentric Remez algorithm can compute rational approximants of type up to $(82,82)$—higher than that obtained by Varga, Ruttan and Carpenter in [62] using $200$-digit arithmetic111Chebfun’s previous remez command (until version 5.6.0 in December 2016) could only go up to type $(8,8)$..

A similar Remez iteration using the barycentric representation was described by Ioni\textcommabelowtă [35, Sec. 3.2.3] in his PhD thesis. We adopt the same set of support points (see Section 4.3), and our analysis justifies its choice: we prove its optimality in a certain sense. A difference from Ioni\textcommabelowtă’s treatment is that we reduce the core computational task to a symmetric eigenvalue problem, rather than a generalized eigenproblem as in [35]. The bigger difference is that Ioni\textcommabelowtă treated just the core iteration for approximations of type $(n,n)$, whereas we generalize the approach to type $(m,n)$ and include the initialization strategies that are crucial for making the entire procedure into a fully practical algorithm.

This work is motivated by the recent AAA algorithm [43] for rational approximation, which uses adaptive barycentric representations with great success. A large part of the text is focused on introducing a robust version of the rational Remez algorithm, followed by a discussion of two other methods for discrete $\ell_{\infty}$ rational approximation: the AAA-Lawson algorithm (efficient at least in the early stages, but non-robust) and the DC algorithm (robust, but not very efficient). We shall see how all three algorithms benefit from an adaptive barycentric basis. In practice, we advocate using the Remez algorithm, mainly for its convergence properties (usually quadratic [21], unlike AAA-Lawson, which converges linearly at best), practical speed (an eigenvalue-based Remez implementation is usually much faster than a linear programming-based DC method), and its ability to work with the interval $[a,b]$ directly rather than requiring a discretization (unlike both AAA-Lawson and DC). AAA-Lawson is used mainly as an efficient approach to initialize the Remez algorithm.

The paper is organized as follows. In Section 2 we review the barycentric representation for rational functions. Sections 3 to 6 are the core of the paper; here we develop the barycentric rational Remez algorithm with adaptive basis functions. Numerical experiments are presented in Section 7. We describe the AAA-Lawson algorithm in Section 8, and in Section 9 we briefly present the barycentric version of the differential correction algorithm. Section 10 presents a flow chart of minimax and an example of how to compute a best approximation in Chebfun.

2 Barycentric rational functions

All of our methods are made possible by a barycentric representation of $r$, in which both the numerator and denominator are given as partial fraction expansions. Specifically, we consider

[TABLE]

where $n\in\mathbb{N}$, $\alpha_{0},\ldots,\alpha_{n}$ and $\beta_{0},\ldots,\beta_{n}$ are sets of real coefficients and $t_{0},\ldots,t_{n}$ is a set of distinct real support points. The names $N$ and $D$ stand for “numerator” and “denominator”.

If we denote by $\omega_{t}$ the node polynomial associated with $t_{0},\ldots,t_{n}$,

[TABLE]

then $p(z)=\omega_{t}(z)N(z)$ and $q(z)=\omega_{t}(z)D(z)$ are both polynomials in $\mathbb{R}_{n}[x]$. We thus get $r(z)=p(z)/q(z)$, meaning that $r$ is a type $(n,n)$ rational function. (This is not necessarily sharp; $r$ may also be of type $(\mu,\nu)$ with $\mu<n$ and/or $\nu<n$.) At each point $t_{k}$ with nonzero $\alpha_{k}$ or $\beta_{k}$, formula (4) is undefined, but this is a removable singularity with $\lim_{z\rightarrow t_{k}}r(z)=\alpha_{k}/\beta_{k}$ (or a simple pole in the case $\alpha_{k}\neq 0,\beta_{k}=0$), meaning $r$ is a rational interpolant to the values $\left\{\alpha_{k}/\beta_{k}\right\}$ at the support points $\left\{t_{k}\right\}$.

Much of the literature on barycentric representations exploits this interpolatory property [55, 7, 10, 8, 27, 12] by taking $\alpha_{k}=f(t_{k})\beta_{k}$, so that $r$ is an interpolant to some given function values $f(t_{0}),\ldots,f(t_{n})$ at the support points. In this case

[TABLE]

with the coefficients $\left\{\beta_{k}\right\}$ commonly known as barycentric weights; we have $r(t_{k})=f(t_{k})$ as long as $\beta_{k}\neq 0$. While such a property is useful and convenient when we want to compute good approximations to $f$ (see in particular the AAA algorithm), for a best rational approximation $r^{*}$ we do not know a priori where $r^{*}$ will intersect $f$, so enforcing interpolation is not always an option. (We use interpolation for Remez but not for AAA-Lawson or DC.) Formula (4), on the other hand, has $2n+1$ degrees of freedom and can be used to represent any rational function of type $(n,n)$ by appropriately choosing $\left\{\alpha_{k}\right\}$ and $\left\{\beta_{k}\right\}$ [43, Theorem 2.1]. We remark that variants of (4) also form the basis for the popular vector fitting [31, 30] method used to match frequency response measurements of dynamical systems. A crucial difference is that the support points $\{t_{k}\}$ in vector fitting are selected to approximate poles of $f$, whereas, as we shall describe in detail, we choose them so that our representation uses a numerically stable basis.

2.1 Representing rational functions of nondiagonal type

Functions $r$ expressed in the barycentric form (4) range precisely over the set of all rational functions of (not necessarily exact) type $(n,n)$. When one requires rational functions of type $(m,n)$ with $m\neq n$, additional steps are needed to enforce the type.

The approach we have followed, which we shall now describe, is a linear algebraic one based on previous work by Berrut and Mittelmann [9], where we make use of Vandermonde matrices to impose certain conditions that limit the numerator or denominator degree. An alternative might be to avoid such matrices and constrain the barycentric representation more directly to have a certain number of poles or zeros at $z=\infty$. This is a matter for future research.

To examine the situation, we first suppose $m<n$ and convert $r$ into the conventional polynomial quotient representation

[TABLE]

The numerator $p$ is a polynomial of degree at most $n$. Further, it can be seen (either via direct computation or from [9, eq. (1)]) that $p$ is of degree $m\ (<n)$ if and only if the vector $\alpha=[\alpha_{0},\ldots,\alpha_{n}]^{T}$ lies in a subspace spanned by the null space of the (transposed) Vandermonde matrix

[TABLE]

That is, to enforce $r\in\mathcal{R}_{m,n}$ with $m<n$, we require $\alpha\in\mbox{span}(P_{m})$, where $P_{m}\in\mathbb{R}^{(n+1)\times(m+1)}$ has orthonormal columns, obtained by taking the full QR factorization $V_{m}^{T}=\begin{bmatrix}P_{m}^{\perp}&P_{m}\end{bmatrix}\begin{bmatrix}R_{m}\\ 0\end{bmatrix}$, where $P_{m}^{\perp}\in\mathbb{R}^{(n+1)\times(n-m)}$, $R_{m}\in\mathbb{R}^{(n-m)\times(n-m)}$. Note that $R_{m}$ is nonsingular if the support points $\{t_{k}\}$ are distinct.

Similarly, for $m>n$, we need to take $m+1$ terms in (4), that is, r(z)=\sum_{k=0}^{m}\alpha_{k}(z-t_{k})^{-1}\big{/}\sum_{k=0}^{m}\beta_{k}(z-t_{k})^{-1}, and force $\beta\in\mbox{span}(P_{n})$, where $\mbox{span}(P_{n})$ is the null space of the matrix

[TABLE]

obtained by the QR factorization $V_{n}^{T}=\begin{bmatrix}P_{n}^{\perp}&P_{n}\end{bmatrix}\begin{bmatrix}R_{n}\\ 0\end{bmatrix}$, where $P_{n}^{\perp}\in\mathbb{R}^{(m+1)\times(m-n)}$, $R_{n}\in\mathbb{R}^{(m-n)\times(m-n)}$.

In Section 4.4 we describe how to use the matrices $P_{m},P_{n}$ in specific situations. Since these matrices are obtained via $V_{m},V_{n}$ in (7)–(8) and real-valued Vandermonde matrices are usually highly ill-conditioned [4, 5, 48], care is needed when computing their null spaces, as extracting the orthogonal factors in QR (or SVD) is susceptible to numerical errors. Berrut and Mittelmann [9] suggest a careful elimination process to remedy this (for a slightly different problem). Here, in view of the Krylov-type structure of the matrices $V_{m}^{T}$ and $V_{n}^{T}$, we propose the following simpler approach, based on an Arnoldi-style orthogonalization:

Let $Q=[1,\ldots,1]^{T}$ when $m>n$, and $Q=[f(t_{0}),\ldots,f(t_{n})]^{T}$ when $m<n$, and normalize to have Euclidean norm 1. 2. 2.

Let $q$ be the last column of $Q$. Take the projection of $\mbox{diag}(t_{0},\ldots,t_{\max(m,n)})q$ onto the orthogonal complement of $Q$, normalize, and append it to the right of $Q$. Repeat this $|m-n|$ times to obtain $Q\in\mathbb{C}^{(\max(m,n)+1)\times(|m-n|)}$. In MATLAB, this is q = Q(:,end); q = diag(t)q; for i = 1:size(Q,2), q = q-Q(:,i)(Q(:,i)’*q); end, q = q/norm(q); Q = [Q,q];. 3. 3.

Take the orthogonal complement $Q^{\perp}$ of $Q$ via computing the QR factorization of $Q$. $Q^{\perp}$ is the desired matrix, $P_{m}$ or $P_{n}$.

Note that the matrix $Q$ in the final step is well conditioned ($\kappa_{2}(Q)=1$ in exact arithmetic), so the final QR factorization is a stable computation.

2.2 Why does the barycentric representation help?

The choice of the support points $\left\{t_{k}\right\}$ is very important numerically, and indeed it is the flexibility of where to place these points that is the source of the power of barycentric representations. If the points are well chosen, the basis functions $1/(x-t_{k})$ lead to a representation of $r$ that is much better conditioned (often exponentially better) than the conventional representation as a ratio of polynomials. We motivate and explain our adaptive choice of $\left\{t_{k}\right\}$ for the Remez algorithm in Sections 4.3 and 4.5. The analogous choices for AAA-Lawson and DC are discussed in Sections 8.5 and 9.2.

To understand why a barycentric representation is preferable for rational approximation, we first consider the standard quotient representation $p/q$. It is well known that a polynomial will vary in size by exponentially large factors over an interval unless its roots are suitably distributed (approximating a minimal-energy configuration). If $p/q$ is a rational approximation, however, the zeros of $p$ and $q$ will be positioned by approximation considerations, and if $f$ has singularities or near-singularities they will be clustered near those points. In the clustering region, $p$ and $q$ will be exponentially smaller than in other parts of the interval and will lose much or all of their relative accuracy. Since the quotient $p/q$ depends on that relative accuracy, its accuracy too will be lost.

A barycentric quotient $N/D$, by contrast, is composed of terms that vary in size just algebraically across the interval, not exponentially, so this effect does not arise. If the support points are suitably clustered, $N$ and $D$ may have approximately uniform size across the interval (away from their poles, which cancel in the quotient), as illustrated in Figure 1.

2.3 Numerical stability of evaluation

Regarding the evaluation of $r$ in the barycentric representation, Higham’s analysis in [34, p. 551] (presented for barycentric polynomial interpolation, but equally valid for (4)) shows that evaluating $r(x)$ is backward stable in the sense that the computed value $\widehat{r}(x)$ satisfies

[TABLE]

where $\epsilon_{\alpha_{k}},\epsilon_{\beta_{k}}$ denote quantities of size $O(u)$, or more precisely, bounded by $(1+u)^{3n+4}$. In other words, $\widehat{r}(x)$ is an exact evaluation of (4) for slightly perturbed $\{\alpha_{k}\},\{\beta_{k}\}$. Note that when $r$ represents a polynomial (as assumed in [34]), (9) does not imply backward stability. However, as a rational function for which we allow for backward errors in the denominator, (9) does imply backward stability.

For the forward error, we can adapt the analysis of [14, Proposition 2.4.3]. Assume that the computed coefficients $\widehat{\alpha},\widehat{\beta}$ are obtained through a backward stable process,

[TABLE]

where $\kappa_{\alpha}$ and $\kappa_{\beta}$ are condition numbers associated with the matrices used to determine $\widehat{\alpha}$ and $\widehat{\beta}$. Then, if $x$ (the evaluation point) and $\{t_{k}\}$ are considered to be floating point numbers, we have

Lemma 1**.**

The relative forward error for the computed value $\widehat{r}(x)$ of (4) satisfies

[TABLE]

Proof.

This follows from [14, Prop. 2.4.3]. ∎

If the functions $|D(x)|$ and $|N(x)|$ appearing in the denominators of the right-hand side of (10) do not become too small over $[a,b]$, then we can expect the evaluation of $\widehat{r}$ to be accurate. Note that $|D(x)|$ is precisely the quantity examined in Section 2.2, where we argued that it takes values $O(1)$ or larger across the interval. Further, since $r(x)\approx f(x)$ implies $|N(x)|\approx|D(x)f(x)|$, we see that $|N(x)|$ is not too small unless $|f(x)|$ is small. Put together, we expect the barycentric evaluation phase to be stable unless $|f(x)|$ (and hence $|r(x)|$) is small. Note that since (10) measures the relative error, we usually cannot expect it to be $O(u)$ when $|r(x)|\approx|f(x)|\ll 1$.

3 The rational Remez algorithm

Initially developed by Werner [65, 64] and Maehly [38], the rational Remez algorithm extends the ideas of computing best polynomial approximations due to Remez [54, 53]. It can be summarized as follows:

• Step 1

Set $k=1$ and choose $m+n+2$ distinct reference points

[TABLE]

• Step 2

Determine the levelled error $\lambda_{k}\in\mathbb{R}$ (positive or negative) and $r_{k}\in\mathcal{R}_{m,n}$ such that $r_{k}$ has no pole on $[a,b]$ and

[TABLE]

• Step 3

Choose as the next reference $m+n+2$ local maxima $\{x_{\ell}^{(k+1)}\}$ of $\left|f-r_{k}\right|$ such that

[TABLE]

with $s\in\left\{\pm 1\right\}$ and such that for at least one $\ell\in\left\{0,\ldots,m+n+1\right\}$, the left-hand side of (12) equals $\left\|f-r_{k}\right\|_{\infty}$. If $r_{k}$ has converged to within a given threshold $\varepsilon_{t}>0$ (i.e., $(\left\|f-r_{k}\right\|_{\infty}-\lambda_{k})/\left\|f-r_{k}\right\|_{\infty}\leq\varepsilon_{t}$ [50, eq. (10.8)]) return $r_{k}$, else go to Step 2 with $k\leftarrow k+1$.

If Step 2 is always successful, then convergence to the best approximation is assured [63, Theorem 9.14]. It might happen that Step 2 fails, namely when all rational solutions satisfying the equations (11) have poles in $[a,b]$. If the best approximation is non-degenerate and the initial reference set is already sufficiently close to optimal, then the algorithm will converge [11, §V.6.B]. To our knowledge, there is no effective way in general to determine when degeneracy is the cause of failure.

We note that the rational Remez algorithm can also be adapted to work in the case of weighted best rational approximation. An early account of this is given in [22]. Given a positive weight function $w\in\mathcal{C}([a,b])$, the goal is to find $r^{*}\in\mathcal{R}_{m,n}$ such that the weighted error $\|f-r^{*}\|_{w,\infty}=\max_{x\in[a,b]}|w(x)(f(x)-r^{*}(x))|$ is minimal. Equations (11) and (12) get modified to

[TABLE]

and

[TABLE]

while the norm computations in Step 3 are taken with respect to $w$. Notice that the ability to work with the weighted error immediately allows us to compute the best approximation in the relative sense, by taking $w(x)=1/|f(x)|$, assuming that $f$ is nonzero over $[a,b]$.

We discuss each step of the rational Remez algorithm in the following sections. We first address Step 2, as this is the core part where the barycentric representation is used. We then discuss initialization (Step 1) in Section 5, and finding the next reference set (Step 3) in Section 6. Our focus is on the unweighted setting, but we comment on how our ideas can be extended to the weighted case as well.

4 Computing the trial approximation

For notational simplicity, in this section we drop the index $k$ referring to the iteration number, the analysis being valid for any iteration of the rational Remez algorithm. We begin with the case $m=n$.

4.1 Linear algebra in a polynomial basis

We first derive the Remez algorithm in an (arbitrary) polynomial basis. At each iteration, we search for $r=p/q\in\mathcal{R}_{n,n},p,q\in\mathbb{R}_{n}[x]$ such that

[TABLE]

and assume that we represent $p$ and $q$ using a basis of polynomials $\varphi_{0},\ldots,\varphi_{n}$ such that $\textnormal{span}_{\mathbb{R}}\left(\varphi_{i}\right)_{0\leq i\leq n}=\mathbb{R}_{n}[x]$:

[TABLE]

The linearized version of (13) is then given by

[TABLE]

which, in matrix form, becomes

[TABLE]

where $\Phi_{x}\in\mathbb{R}^{(2n+2)\times(n+1)}$ is the basis matrix $\left(\Phi_{x}\right)_{\ell,k}=\varphi_{k}(x_{\ell}),0\leq\ell\leq 2n+1,0\leq k\leq n$, and $c_{p}=[c_{p,0},c_{p,1},\ldots,c_{p,n}]^{T}$ and $c_{q}=[c_{q,0},c_{q,1},\ldots,c_{q,n}]^{T}$ are the coefficient vectors of $p$ and $q$. Note that in this paper, vector and matrix indices always start at zero. Up to multiplying both sides on the left by a nonsingular diagonal matrix $D=\mathop{\operator@font diag}\nolimits\left(d_{0},\ldots,d_{2n+1}\right)$, (14) can also be written as a generalized eigenvalue problem

[TABLE]

with $F=\mathop{\operator@font diag}\nolimits\left(f(x_{0}),\ldots,f(x_{2n+1})\right)$ and $S=\mathop{\operator@font diag}\nolimits\left((-1)^{k+1}\right)$.

As described in Powell [50, Ch. 10.2], solving (15) is usually done by eliminating $c_{p}$. His presentation considers the monomial basis, but the approach is valid for any basis of $\mathbb{R}_{n}[x]$. By taking the full QR decomposition of $D\Phi_{x}$, we get

[TABLE]

Since $D\Phi_{x}$ is of full rank, we have $Q_{1},Q_{2}\in\mathbb{R}^{(2n+2)\times(n+1)}$ and $Q_{2}^{T}Q_{1}=0$. By multiplying (15) on the left by $Q^{T}=\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}^{T}$, we obtain a block triangular eigenvalue problem with lower-right $(n+1)\times(n+1)$ block

[TABLE]

(The top-left $(n+1)\times(n+1)$ block has all eigenvalues at infinity, and is thus irrelevant.) In terms of polynomials, $(Q_{1})_{\ell,k}=d_{\ell}\psi_{k}(x_{\ell})$, $0\leq k\leq n,0\leq\ell\leq 2n+1$, where $(\psi_{k})_{0\leq k\leq n}$ is a family of orthonormal polynomials with respect to the discrete inner product $\left\langle f,g\right\rangle_{x}=\sum_{k=0}^{2n+1}d_{k}^{2}f(x_{k})\overline{g(x_{k})}$. Moreover, if $(\varphi_{k})_{0\leq k\leq n}$ is a degree-graded basis with $\deg\varphi_{k}=k$, then we have $\deg\psi_{k}=k,0\leq k\leq n$.

Let $\omega_{x}$ be the node polynomial associated with the reference nodes $x_{0},\ldots,x_{2n+1}$, and $\Omega_{x}=\mathop{\operator@font diag}\nolimits\left(1/\omega_{x}^{\prime}(x_{0}),\ldots,1/\omega_{x}^{\prime}(x_{2n+1})\right)$. We have [50, p. 114]

[TABLE]

where $V_{x}\in\mathbb{R}^{(2n+2)\times(n+1)}$ is the Vandermonde matrix associated with $x_{0},\ldots,x_{2n+1}$, that is, $(V_{x})_{i,j}=x_{i}^{j}$. Indeed,

[TABLE]

the divided differences of order $2n+1$ of the function $x^{i+j}$ at the $\left\{x_{\ell}\right\}$ nodes, hence [math] if $i+j\leq 2n$.

By using the appropriate change of basis matrix in (17), we have

[TABLE]

Now, by multiplying (15) on the left by $\Phi_{x}^{T}\Omega_{x}D^{-1}$ and using (18), we can eliminate the $c_{p}$ term to obtain

[TABLE]

Equation (19) is the extension of [50, Eq. (10.13)] from the monomial basis to $\varphi_{0},\ldots,\varphi_{n}$. Moreover, we have:

Lemma 2**.**

The matrix $\Phi_{x}^{T}\Omega_{x}S\Phi_{x}$ is symmetric positive definite.

Proof.

Since $\Omega_{x}S=\left|\Omega_{x}\right|$, it means that $\Omega_{x}S$ is symmetric positive definite, and the conclusion follows. See also [50, Theorem 10.2]. ∎

Since $\Phi_{x}^{T}\Omega_{x}F\Phi_{x}$ is also symmetric, it follows that all eigenvalues of (19) are real and at most one eigenvector $c_{q}$ corresponds to a pole-free solution $r$ (i.e., $q$ has no root on $[a,b]$). To see this, suppose to the contrary that there exists another pole-free solution $r^{\prime}$. Then, from (13), it follows that either $r(x_{k})-r^{\prime}(x_{k})$ are all zero or they alternate in sign at least $2n+1$ times. In both cases, $r-r^{\prime}\in\mathcal{R}_{2n,2n}$ has at least $2n+1$ zeros inside $[a,b]$, leading to $r=r^{\prime}$.

We can in fact transform (16) into a symmetric eigenvalue problem (an observation which seems to date to [49]) by considering the choice $D=\left|\Omega_{x}\right|^{1/2}$, which leads to $Q_{2}=SQ_{1}$ in view of (18). The system becomes $Q_{1}^{T}SFQ_{1}Rc_{q}=\lambda Q_{1}^{T}S^{2}Q_{1}Rc_{q},$ which, by the change of variables $y=Rc_{q}$, gives

[TABLE]

To get $c_{p}$, from (14), we have $\left|\Omega_{x}\right|^{1/2}\Phi_{x}c_{p}=(F-\lambda S)\left|\Omega_{x}\right|^{1/2}\Phi_{x}c_{q},$ or equivalently (by multiplication on the left by $Q_{1}^{T}$),

[TABLE]

The vectors $Rc_{p}$ and $Rc_{q}$ can be seen as vectors of coefficients of the numerator and denominator of $r$ in the orthogonal basis $\psi_{0},\ldots,\psi_{n}$. The (scaled) values of the denominator at each $x_{k}$ corresponding to an eigenvector $y$ can be recovered by computing

[TABLE]

From this we can confirm the uniqueness of the pole-free solution: since the eigenvectors are orthogonal, there is at most one generating a vector of denominator values of the same sign, making it the only pole-free solution candidate.

4.2 Linear algebra in a barycentric basis

An equivalent analysis is valid if we take $r$ in the barycentric form (4). Namely, (13) becomes

[TABLE]

where $C$ is now a $(2n+2)\times(n+1)$ Cauchy matrix with entries $C_{\ell,k}=1/(x_{\ell}-t_{k})$ (we assume for the moment $\{x_{\ell}\}\cap\{t_{k}\}=\varnothing$) and $\alpha=[\alpha_{0},\alpha_{1},\ldots,\alpha_{n}]^{T}$ and $\beta=[\beta_{0},\beta_{1},\ldots,\beta_{n}]^{T}$ are the column vectors of coefficients $\left\{\alpha_{k}\right\}$ and $\left\{\beta_{k}\right\}$. Again, this can be transformed into a generalized eigenvalue problem

[TABLE]

To reduce (23) to a symmetric eigenvalue problem as in (20), we form a link between the monomial and barycentric representations in terms of the basis matrices $V_{x}$ and $C$. We have:

Lemma 3**.**

Let $V_{x}$, $\omega_{t}$ be as defined above, and $V_{t}\in\mathbb{R}^{(n+1)\times(n+1)}$ be the Vandermonde matrix corresponding to the support points, i.e., $(V_{t})_{i,j}=t_{i}^{j}$. Then

[TABLE]

Proof.

If we look at an arbitrary element of the right-hand side matrix, we have

[TABLE]

where the second equality is a consequence of the Lagrange interpolation formula. ∎

In place of $\Omega_{x}$ we will use the following matrix $\Delta$:

Lemma 4**.**

If $\Delta=\mathop{\operator@font diag}\nolimits\left(\omega_{t}(x_{0})^{2},\ldots,\omega_{t}(x_{2n+1})^{2}\right)\Omega_{x}$, then $C^{T}\Delta C=0$.

Proof.

We apply Lemma 3 and use the fact that $V_{x}^{T}\Omega_{x}V_{x}=0$. Namely, $C^{T}\Delta C=\mathop{\operator@font diag}\nolimits\left(\omega_{t}^{\prime}(t_{0}),\ldots,\omega_{t}^{\prime}(t_{n})\right)V_{t}^{-T}V_{x}^{T}\Omega_{x}V_{x}V_{t}^{-1}\mathop{\operator@font diag}\nolimits\left(\omega_{t}^{\prime}(t_{0}),\ldots,\omega_{t}^{\prime}(t_{n})\right)=0$. ∎

We now take the full QR decomposition of $\left|\Delta\right|^{1/2}C=(S\Delta)^{1/2}C$. We have

[TABLE]

Based on Lemma 4, we can again take $Q_{2}=SQ_{1}$. From (23) we get

[TABLE]

Multiplying this expression on the left by $\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}^{T}$ gives a block triangular matrix pencil, whose $(n+1)\times(n+1)$ lower-right corner is the barycentric analogue of (16): $Q_{2}^{T}FQ_{1}R\beta=\lambda Q_{2}^{T}SQ_{1}R\beta.$ After substituting $Q_{2}^{T}=Q_{1}^{T}S$, we get

[TABLE]

which, by the change of variable $y=R\beta$, becomes a standard symmetric eigenvalue problem in $\lambda$ with eigenvector $y$ (recall that $S,F$ are diagonal):

[TABLE]

Hence, computing its eigenvalues is a well-conditioned operation. The values of the denominator of the rational interpolant corresponding to each eigenvector $y$ can be recovered by computing

[TABLE]

As in the polynomial case, there is at most one solution such that $q(x)=D(x)\omega_{t}(x)$ has no root in $[a,b]$; indeed, (21) and (26) represent the values of $q(x_{\ell})$ for $r=p/q$ and $x_{\ell}$ satisfying equation (13). We use this sign test involving (26) to determine the levelled error $\lambda$ that gives a pole-free $r$ in Step 2 of our rational Remez algorithm. The appropriate $\beta$ is then taken by solving $R\beta=y$. From (22), we have

[TABLE]

or equivalently (by multiplication on the left by $Q_{1}^{T}$)

[TABLE]

which allows us to recover $\alpha$ (and thus $r$).

Most of the derivations in this section can be carried over to the weighted approximation setting as well. In particular, the reader can check that the weighted versions of Equations (23) and (25) correspond to

[TABLE]

and

[TABLE]

where $W=\mathop{\operator@font diag}\nolimits\left(w(x_{0}),\ldots,w(x_{2n+1})\right)$ and all the other quantities are the same as before. While not leading to a symmetric eigenvalue problem, the symmetric and symmetric positive definite matrices appearing in the second pencil seem to suggest that the eigenproblem computations will again correspond to well-conditioned operations. Our experiments support this statement and we leave it as future work to make this rigorous. To recover $\alpha$, (27) becomes $R\alpha=Q_{1}^{T}(F-\lambda SW^{-1})Q_{1}y$.

4.3 Conditioning of the QR factorization

Since the above discussion makes heavy use of the matrix $Q_{1}$, it is desirable that computing the (thin) QR factorization $\left|\Delta\right|^{1/2}C=Q_{1}R$ is a well-conditioned operation.

Here we examine the conditioning of $Q_{1}$, the orthogonal factor in the QR factorization of $|\Delta|^{1/2}C$, as this is the key matrix for constructing (24). We use the fact that the standard Householder QR algorithm is invariant under column scaling, that is, it computes the same $Q_{1}$ for both $\left|\Delta\right|^{1/2}C$ and $\left|\Delta\right|^{1/2}C\Gamma$ for diagonal $\Gamma$ [33, Ch. 19]. We thus consider

[TABLE]

where $\mathcal{D}_{n+1}$ is the set of $(n+1)\times(n+1)$ diagonal matrices. We have

Theorem 5**.**

Let $t_{k}\in(x_{2k},x_{2k+1})$ for $k=0,\ldots,n$ and $s_{k}\in(x_{2k+1},x_{2k+2})$ for $k=0,\ldots,n-1$, $s_{n}\in(x_{2n+1},\infty)$, and define $\omega_{s}(x)=\prod_{k=0}^{n}(x-s_{k})$. Then

[TABLE]

Proof.

Let $\{y_{j}\}$ be a $(2n+2)$-element set such that $y_{j}\in(x_{j},x_{j+1}),j=0,\ldots,2n$, $y_{2n+1}>x_{2n+1}$ and let $C_{x,y}\in\mathbb{R}^{(2n+2)\times(2n+2)}$ be the Cauchy matrix with elements $(C_{x,y})_{j,k}=1/(x_{j}-y_{k})$. If we consider $D_{1}=\mbox{diag}(\sqrt{\left|\omega_{y}(x_{j})/\omega_{x}^{\prime}(x_{j})\right|})$ and $D_{2}=\mbox{diag}(\sqrt{\left|\omega_{x}(y_{j})/\omega_{y}^{\prime}(y_{j})\right|})$, then the matrix $D_{1}C_{x,y}D_{2}$ is orthogonal. This follows, for instance, if we examine the elements of its associated Gram matrix $G$ and use divided differences. Indeed, for an arbitrary element $(G)_{j,k}$ with $j\neq k$, we have

[TABLE]

Similarly, since $\prod_{j\neq k}(x-y_{j})=q(x)(x-y_{k})+\omega_{y}^{\prime}(y_{k})$, with $q\in\mathbb{R}_{2n}[x]$, we have,

[TABLE]

Now, if we take $t_{k}=y_{2k},s_{k}=y_{2k+1},$ for $k=0,\ldots,n$, there exist $D\in\mathcal{D}_{2n+2}$ and $\Gamma\in\mathcal{D}_{n+1}$ such that $\left|\Delta\right|^{1/2}C\Gamma=DD_{1}C_{x,y}D_{2}I_{t}$, where $D=\mbox{diag}(\sqrt{\left|\omega_{t}(x_{j})/\omega_{s}(x_{j})\right|})$ and $I_{t}$ is obtained by removing every second column from $I_{2n+2}$. In particular, $\Gamma=I_{t}^{T}D_{2}I_{t}$. It follows that

[TABLE]

∎

Let $\Gamma=I_{t}^{T}D_{2}I_{t}$ be as in the proof of Theorem 5. It turns out that for the choice $t_{k}=x_{2k+1}-\varepsilon,s_{k}=x_{2k+1}+\varepsilon$, for $k=0,\ldots,n$, as $\varepsilon\rightarrow 0$, the matrix $\left|\Delta\right|^{1/2}C$ has a finite limit $\widetilde{C}$ of full column rank, and similarly $\Gamma$ tends to some diagonal matrix $\widetilde{\Gamma}$ with positive diagonal entries. From Theorem 5 and its proof we know that $\widetilde{C}\widetilde{\Gamma}$ has condition number 1, and, more precisely, orthonormal columns. We thus obtain an explicit thin QR decomposition of $\widetilde{C}$ (by direct calculation):

Corollary 6**.**

In the limit $t_{k}\nearrow x_{2k+1}$, for $k=0,\ldots,n$, the matrix $\left|\Delta\right|^{1/2}C$ converges to $\widetilde{C}$, with entries

[TABLE]

and explicit thin QR decomposition $\widetilde{C}=Q_{1}R$, where

[TABLE]

and $R=\sqrt{2}\ \textnormal{diag}\left(\frac{|w^{\prime}_{t}(t_{0})|}{\sqrt{|w_{x}^{\prime}(t_{0})|}},\ldots,\frac{|w^{\prime}_{t}(t_{n})|}{\sqrt{|w_{x}^{\prime}(t_{n})|}}\right).$

Corollary 6 suggests the choice

[TABLE]

This takes us back to the interpolatory mode of barycentric representations (5), in which we take $\alpha_{k}=\beta_{k}(f(t_{k})-\lambda)$ for all $k$, instead of solving the system (27). This interpolatory mode formulation is used in [35, Sec. 3.2.3]. Our derivation provides a theoretical justification by showing that it is optimal with respect to the conditioning of $\left|\Delta\right|^{1/2}C\Gamma$. Moreover, since $\min_{\Gamma\in\mathcal{D}_{n+1}}\kappa_{2}(\widetilde{C}\Gamma)=1$ in (28), forming the QR factorization of $\left|\Delta\right|^{1/2}C$ via a standard algorithm (e.g. Householder QR) to obtain $Q_{1}$ is actually unnecessary, as the explicit form of $Q_{1}$ is given in Corollary 6. In addition, we reduce the problem to a symmetric eigenvalue problem (25), resulting in well-conditioned eigenvalues, with $\beta$ being obtained by solving the diagonal system $R\beta=y$ with $y$ as in (25). Compared to (13), where we want $q$ to have the same sign over $\{x_{\ell}\}$, we similarly require that $\beta$ and thus $y$ have components alternating in sign, which uniquely fixes the norm 1 eigenvector $y$ in (25). Our approach also allows for nondiagonal types, as we describe next.

4.4 The nondiagonal case $\boldmath{m\neq n}$

As pointed out in Section 2.1, when searching for a best approximant with $m\neq n$, we need to force the coefficient vector $\alpha$ or $\beta$ to lie in a certain subspace. This results in modified versions of (23). Namely,

[TABLE]

for $\widehat{\beta}\in\mathbb{C}^{n+1}$, and we take $\beta=P_{n}\widehat{\beta}$. Similarly,

[TABLE]

for $\widehat{\alpha}\in\mathbb{C}^{m+1}$, and we take $\alpha=P_{m}\widehat{\alpha}$.

Below we describe the reduction of the generalized eigenvalue problems (31) and (32) to standard symmetric eigenvalue problems.

Case $m>n$

In this case, $C\in\mathbb{R}^{(m+n+2)\times(m+1)}$. Since $\det|\Delta|^{1/2}\neq 0$, (31) is equivalent to the generalized eigenvalue problem

[TABLE]

Consider the (thin) QR decomposition of $\left|\Delta\right|^{1/2}C\begin{bmatrix}P_{n}&P_{n}^{\perp}\end{bmatrix}=(S\Delta)^{1/2}C\begin{bmatrix}P_{n}&P_{n}^{\perp}\end{bmatrix}$:

[TABLE]

Then we have the identity $\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}^{T}(SQ_{1})=0$, as can be verified analogously to (17) using divided differences. This implies $(SQ_{1})^{T}\left|\Delta\right|^{1/2}C=0$, so by left-multiplying (33) by $\begin{bmatrix}(SQ_{1})^{\perp}&SQ_{1}\end{bmatrix}^{T}$ we obtain a block upper-triangular eigenvalue problem with lower-right $(n+1)\times(n+1)$ block

[TABLE]

which again reduces to the standard symmetric eigenvalue problem (setting $y=R_{1}\widehat{\beta}$)

[TABLE]

From (33), we have $\left|\Delta\right|^{1/2}C\alpha=(F-\lambda S)\left|\Delta\right|^{1/2}CP_{n}\widehat{\beta}$. Left-multiplying by $\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}^{T}$ and using $\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}^{T}S\left|\Delta\right|^{1/2}CP_{n}=0$, we obtain

[TABLE]

Therefore

[TABLE]

which is obtained by computing the vector $\widehat{y}=\begin{bmatrix}Q_{1}&Q_{2}\end{bmatrix}^{T}FQ_{1}y$, then solving $R\widetilde{y}=\widehat{y}$ for $\widetilde{y}$, then $\alpha=\begin{bmatrix}P_{n}&P_{n}^{\perp}\end{bmatrix}\widetilde{y}$.

Case $m<n$

This case is analogous to the previous one; we highlight the differences. $C$ is a $(m+n+2)\times(n+1)$ matrix. Equation (32) is equivalent to

[TABLE]

Consider the (thin) QR decompositions

[TABLE]

Here $Q_{1}\in\mathbb{R}^{(m+n+2)\times(n+1)},\widehat{Q}_{1}\in\mathbb{R}^{(m+n+2)\times(m+1)}$. We have $\widehat{Q}_{1}^{T}(SQ_{1})=0$, which again can be established using divided differences. This implies $(SQ_{1})^{T}\left|\Delta\right|^{1/2}CP_{m}=0$, so left-multiplying equation (35) by $\begin{bmatrix}(SQ_{1})^{\perp}&SQ_{1}\end{bmatrix}^{T}$ results in a block upper-triangular eigenvalue problem with lower-right block

[TABLE]

which also reduces to the standard symmetric eigenvalue problem (setting $y=R\beta$)

[TABLE]

From (35), we have $\left|\Delta\right|^{1/2}CP_{m}\widehat{\alpha}=(F-\lambda S)\left|\Delta\right|^{1/2}C\beta$. Left-multiplying by $\widehat{Q}_{1}^{T}$ and using $\widehat{Q}_{1}^{T}S\left|\Delta\right|^{1/2}C=0$, we obtain

[TABLE]

Therefore

[TABLE]

obtained via $\widehat{y}=\widehat{Q}_{1}^{T}FQ_{1}y$, then solving the linear system $\widehat{R}\widehat{\alpha}=\widehat{y}$.

Analogously to our comments at the end of Section 4.2, the analysis for nondiagonal approximation presented here carries over to the weighted setting. In both the $m>n$ and $m<n$ scenarios, the standard symmetric eigenproblems (34) and (36) become

[TABLE]

where $y=R_{1}\widehat{\beta}$ when $m>n$ and $y=R\beta$ when $m<n$. Recovering the set of barycentric coefficients in the numerator corresponds to solving the systems

[TABLE]

and

[TABLE]

Stability and conditioning

We have just shown that the matrices arising in our rational Remez algorithm have explicit expressions, and the eigenvalue problem reduces to a standard symmetric problem. Indeed, our experiments corroborate that we have greatly improved the stability and conditioning of the rational Remez algorithm using the barycentric representation. However, the algorithm is still not guaranteed to compute $r^{*}$ to machine precision. Let us summarize the situation for the unweighted case. As shown in Corollary 6, the computation of $Q_{1}$ can be done explicitly, and the linear system $y=R\beta$ is diagonal, hence can be solved with high relative accuracy. The main source of numerical errors is therefore in the symmetric eigenvalue problem (25), (34) or (36). As is well known, by Weyl’s bound [57, Cor. IV.4.9], eigenvalues of symmetric matrices are well conditioned with condition number 1; thus $\lambda$ is computed with $O(u)$ accuracy, assuming for simplicity that $\|f\|_{\infty}=1$ (without loss of generality). The eigenvector, on the other hand, has conditioning $O(1/\mbox{gap})$ [57, Ch. V], where gap is the distance between the desired $\lambda$ and the rest of the eigenvalues. These eigenvalues are equal to those of the nonzero eigenvalues of the generalized eigenproblem (15), and are inherent in the Remez algorithm, i.e., they cannot be changed e.g. by a change of bases. For a fixed $f$, gap tends to decrease as $m,n$ increase, and we typically have $\mbox{gap}=O(|\lambda|)$. Hence the computed eigenvector tends to have accuracy $O(u/|\lambda|)$, and if the eigenvector $y$ has small elements, the componentwise relative accuracy may be worse. The computation therefore breaks down (perhaps as expected) when $|\lambda|=O(u)$, that is, when the error curve has amplitude of size machine precision.

4.5 Adaptive choice of the support points

Theorem 5 gives an optimal choice of support points $t_{k}=x_{2k+1}$ in terms of optimizing $\min_{\Gamma\in\mathcal{D}_{n+1}}\kappa_{2}(\left|\Delta\right|^{1/2}C\Gamma)$. In Section 2.2 we discussed another desideratum for the support points $\{t_{k}\}$: the resulting $|D(x_{\ell})|=|q(x_{\ell})\prod_{k=0}^{n}(x_{\ell}-t_{k})|$ should take uniformly large values for all $\ell$. Fortunately, this requirement is also met with this choice, as was illustrated in Figure 1.

When $m\neq n$, (30) does not determine enough support points. We take the remaining $|m-n|$ support points from the rest of the reference points in Leja style, i.e., to maximize the product of the differences (see for instance [52, p. 334]). This is a heuristic strategy, and the optimal choice is a subject of future work: indeed, in this case $\min_{\Gamma\in\mathcal{D}_{n+1}}\kappa_{2}(\left|\Delta\right|^{1/2}CP_{m,n}\Gamma)>1$.

5 Initialization

An indispensable component of a successful Remez algorithm implementation is a method for finding a good set of initial reference points $\{x_{\ell}\}$. A key element of our approach is the AAA-Lawson algorithm, which can efficiently find an approximate solution to the minimax problem (2) (to low accuracy).

5.1 Carathéodory-Fejér (CF) approximation

We first attempt to compute the CF approximant [59, 61] to $f$, and use it to find the initial reference points (as explained in Section 6). The dominant computation is an SVD of a Hankel matrix of Chebyshev coefficients, which usually does not cause a computational bottleneck. This method was also used in the previous Chebfun remez code. When $f$ is smooth, the result produced by CF approximation is often indistinguishable from the best approximation, but nonsmooth cases may be very different.

5.2 AAA-Lawson approximation

This approach is based on the AAA algorithm [43] followed by an adaptation of the Lawson algorithm. The resulting algorithm is also based crucially on the barycentric representation. To keep the focus on Remez, we defer the details to Section 8.

The output of the AAA-Lawson iteration typically has a nearly equioscillatory error curve $e=f-r$, from which we find the initial set of reference points as the extrema of $e$. For the prototypical example $f=|x|$, AAA-Lawson initialization lets our barycentric minimax code converge for type up to $(40,40)$. The entire process relies on a moderate number of SVDs (say $\max(m,n)+10$).

5.3 Using lower degree approximations

We resort to this strategy if CF and AAA-Lawson fail to produce a sufficiently good initial guess. For functions $f$ with singularities in $[a,b]$, the reference sets $\{x_{\ell}\}$ corresponding to best approximations in (3) tend to cluster near these singularities as $m$ and $n$ increase.

It is sensible to expect that first computing a type $(m^{\prime},n^{\prime})$ best approximation to $f$ with $m^{\prime}\ll m$ and $n^{\prime}\ll n$ is easier (with convergence achieved if necessary with the help of CF or AAA-Lawson). We then proceed by progressively increasing the values of $m^{\prime}$ and $n^{\prime}$ by small increments $j$, typically $j\in\left\{1,2,4\right\}$. The steps taken follow a diagonal path, as explained in Figure 2. Note that in addition to improving the robustness of the Remez algorithm, this strategy can help detect degeneracy; recall the discussion after (3). It proves useful for many examples, including some of those shown in Section 7: type $(n,n)$ approximations to $f(x)=|x|,x\in[-1,1]$ for $n>40$ and the $f_{1},f_{2}$ and $f_{4}$ specifications in Table 1.

6 Searching for the new reference

We now turn to the updating strategy for the reference points $x_{0}\ldots,x_{m+n+1}$ during the Remez iterations. These are a subset of the local extrema of the error function $e(x)=f(x)-r(x)$. To find them, we decompose the domain $[a,b]$ into subintervals of the form $[\tilde{x}_{\ell},\tilde{x}_{\ell+1}]$ (and $[a,\tilde{x}_{0}]$ and $[\tilde{x}_{m+n+1},b]$, if non-degenerate; here $\{\tilde{x}_{\ell}\}$ are the old reference points) and then compute Chebyshev interpolants $p_{e}(x)$ of $e(x)$ on each subinterval. In addition, if $f$ has singularities (identified by Chebfun’s splitting on functionality [46]), then we further divide the subintervals at those points. Since $e(x)$ is then smooth and each subinterval is small, typically a low degree suffices for $p_{e}=\sum_{i=0}^{k}c_{i}T_{i}(x)$: we start with $2^{3}+1$ points (degree $k=8$), and resample if necessary (determined by examining the decay of the Chebyshev coefficients). We then find the roots of $p_{e}^{\prime}(x)=\sum_{i=1}^{k}ic_{i}U_{i-1}(x)$ (using the formula $T_{n}^{\prime}(x)=nU_{n-1}(x)$) via the eigenvalues of the colleague matrix for Chebyshev polynomials of the second kind [28]. Typically, one local extremum per subinterval is found, resulting in $m+n+2$ points, including the endpoints. If more extrema are found, we evaluate the values of $|e(x)|$ at those points and select those with the largest values that satisfy (12).

7 Numerical results

All computations in this section were done using Chebfun’s new minimax command in standard IEEE double precision arithmetic.

Let us start with our core example of approximating $|x|$ on $[-1,1]$, a problem discussed in detail in [58, Ch. 25]. For more than a century, this problem has attracted interest. The work of Bernstein and others in the 1910s led to the theorem that degree $n\geq 0$ polynomial approximations of this function can achieve at most $O(n^{-1})$ accuracy, whereas Newman in 1964 showed that rational approximations can achieve root-exponential accuracy [45]. The convergence rate for best type $(n,n)$ approximations was later shown by Stahl [56] to be $E_{n,n}(|x|,[-1,1])\sim 8e^{-\pi\sqrt{n}}$.

This result had in fact been conjectured by Varga, Ruttan and Carpenter [62] based on a specialized multiple precision (200 decimal digits) implementation of the Remez algorithm. Their computations were performed on the square root function, using the fact that $E_{2n,2n}(|x|,[-1,1])=E_{n,n}(\sqrt{x},[0,1])$, as follows from symmetry. They went up to $n=40$. In both settings, the equioscillation points cluster exponentially around $x=0$ (see second plot of Figure 4), making it extremely difficult to compute best approximations. Our barycentric Remez algorithm in double precision arithmetic is able to match their performance, in the sense that we obtain the type $(80,80)$ best approximation to $|x|$ in less than 15 seconds on a desktop machine. The results are showcased in Figure 4, where our levelled error computation for the type $(80,80)$ approximation (value $4.39\ldots\times 10^{-12}$) matches the corresponding error of [62, Table 1] to two significant digits, even though the floating point precision is no better than $10^{-16}$.

Running the other non-barycentric codes (Maple’s numapprox[minimax], Mathematica’s MiniMaxApproximation (which requires $f$ to be analytic on $[a,b]$), and Chebfun’s previous remez) on the same example resulted in failures at very small values of $n$ (all for $n\leq 8$).

The robustness of our algorithm is also illustrated by the examples of Table 1 and Figure 5, which is a highlight of the paper. Computing these five approximations takes in total less than 50 seconds with minimax. Example $f_{4}$ is taken from [60, §5], while $f_{5}$ is inspired by [51]. The difficulty of approximating $f_{5}$ is even more pronounced than for $|x|$, since best type $(n,n)$ approximations to $f_{5}$ offer at most $O(n^{-1})$ accuracy (a stark contrast to the root-exponential behavior of $E_{n,n}(|x|,[-1,1])$) and the reference points cluster even more strongly, quickly falling below machine precision.

In Figures 6 and 7, we further illustrate minimax and its weighted variant, by revisiting some classical problems in rational approximation: the Zolotarev problems [2, Ch. 9]. Among other questions, Zolotarev asked what are the best rational approximants to the sign function (on the union of intervals $[-b,-a]\cup[a,b]$ for scalars $0<a<b$) and the $\sqrt{x}$ function (in the relative sense, i.e., minimizing $\|1-r/\sqrt{x}\|_{\infty}$) on $[1/b^{2},1/a^{2}]$. Zolotarev proved these problems are mathematically equivalent through the identity $\mbox{sign}(x)=x\sqrt{1/x^{2}}$: if $r$ is the type $(m,m)$ best approximant to $\sqrt{x}$ on $[1/b^{2},1/a^{2}]$, then $\mbox{sign}(x)-xr(1/x^{2})$ is found to equioscillate at $4m+4$ points on $[-b,-a]\cup[a,b]$, so $xr(1/x^{2})$ is the best approximant to $\mbox{sign}(x)$ of type $(2m+1,2m)$ on $[-b,-a]\cup[a,b]$. Furthermore, Zolotarev gave explicit solutions involving Jacobi’s elliptic functions. These rational functions have the remarkable property of preserving optimality under appropriate composition [42]. In Figure 6 we compute the best relative error approximant of type $(m,m)$ to $\sqrt{x}$ using the weighted variant of our rational Remez algorithm. We then compute $xr(1/x^{2})$, the type $(2m+1,2m)$ best approximant to the sign function. The error function is shown in Figure 7, confirming Zolotarev’s results.

We emphasize that the examples presented in this section are extraordinarily challenging, far beyond the capabilities of most codes for minimax approximation. Chebfun minimax not only solves them but does so quickly. For smoother functions such as analytic functions (with singularities, if any, lying far from the interval), we find that minimax usually easily computes $r^{*}$ so long as $\|f-r^{*}\|_{\infty}$ is a digit or two larger than $u\|f\|_{\infty}$.

8 AAA-Lawson algorithm

Here we describe a new algorithm for rational approximation that we call the AAA-Lawson algorithm; in practice we recommend this for computing an initial guess for the Remez iteration. It applies on a finite, discrete set rather than the continuous interval $[a,b]$ as in (2). Specifically, we consider the problem

[TABLE]

where $Z=\{z_{1},\ldots,z_{M}\}$ is a set of distinct points (sample points) in $[a,b]$. The number $M$ is usually large, $\textnormal{e.g.~{}}10^{5}$, and in particular much bigger than $m$ and $n$. The idea is that the solution for the discrete problem (37) should converge to the continuous one (2) if we discretize the interval densely enough.

AAA-Lawson proceeds as follows:

Use the AAA algorithm to find an approximant (5), in particular the support points $\left\{t_{k}\right\}$ for a rational approximation $r$ to $f$. This step is not tied to a particular norm. 2. 2.

Use a variant of Lawson’s algorithm to obtain a refined (near-best) rational approximant in the $\ell_{\infty}$ norm.

Below we first review the AAA algorithm, introduced in [43], then the Lawson algorithm, and then we present the AAA-Lawson combination.

8.1 The AAA algorithm

Given a function $f$ and sample points $Z\in\mathbb{C}^{M}$, the AAA algorithm finds a rational approximant of type $(n,n)$ represented as in (5) by r(z)=\widetilde{N}(z)/\widetilde{D}(z):=\sum_{k=0}^{n}f(t_{k})\beta_{k}(z-t_{k})^{-1}\big{/}\sum_{k=0}^{n}\beta_{k}(z-t_{k})^{-1}. Here, the support points $\{t_{k}\}$ are a subset of $Z$ chosen in an adaptive, greedy manner so as to improve the approximation as we increase $n$, exploiting the interpolatory property $\widetilde{N}(t_{k})/\widetilde{D}(t_{k})=f(t_{k})$ for all $k$ (unless $\beta_{k}=0$). AAA takes only $\beta_{k}$ as the unknowns, which are found by solving a linearized least-squares problem of the form $\operatorname*{minimize}_{\|\beta\|_{2}=1}\|f\widetilde{D}-\widetilde{N}\|_{\widetilde{Z}}$, where the subscript $\widetilde{Z}$ denotes the discrete $2$-norm at points $\widetilde{Z}:=Z\setminus\left\{t_{0},\ldots,t_{n}\right\}$. For details, see [43].

Noninterpolatory AAA

As we discussed in Section 2, the representation $\widetilde{N}(z)/\widetilde{D}(z)$ is unsuitable when the goal is to represent $r^{*}$: it is necessary to use the representation r(z)=N(z)/D(z)=\sum_{k=0}^{n}\alpha_{k}(z-t_{k})^{-1}\big{/}\sum_{k=0}^{n}\beta_{k}(z-t_{k})^{-1} as in (4). This leads to a noninterpolatory variant of AAA, discussed briefly in [43, Section 10]. The resulting least-squares problem $\operatorname*{minimize}_{\|\alpha\|^{2}_{2}+\|\beta\|^{2}_{2}=1}\|fD-N\|_{\widetilde{Z}}$ has unknowns $\alpha$ and $\beta$. Written in matrix form, it takes the form

[TABLE]

where $F=\mbox{diag}(f(\widetilde{Z}))$, and $C_{\ell,k}=1/(z_{\ell}-t_{k})$ is the Cauchy (basis) matrix as in (23), but with rows corresponding to $z_{\ell}\in\{t_{0},\ldots,t_{n}\}$ removed. We take the same support points $\{t_{k}\}$ as in AAA. We solve (38) by computing the SVD of the matrix $\begin{bmatrix}C&-FC\end{bmatrix}$ and finding the right singular vector v=\big{[}\begin{smallmatrix}\alpha\\ \beta\end{smallmatrix}\big{]}\in\mathbb{R}^{2n+2} corresponding to the smallest singular value. As in Section 4.4, the case $m\neq n$ also uses the projection matrices $P_{m},P_{n}$.

8.2 Lawson’s algorithm

Lawson’s algorithm [37] computes the best polynomial (linear) approximation based on an iteratively reweighted least-squares process. During the iteration, a set of weights is updated according to the residual of the previous solution.

Specifically, suppose that $f$ is to be approximated on $Z=\{z_{1},\ldots,z_{M}\}$ in a linear subspace $\mbox{span}(g_{i})_{i=0}^{n}$. With an initial set of weights $\left\{w_{j}\right\}_{j=1}^{M}$ such that $w_{j}\geq 0$ and $\sum_{j=1}^{M}w_{j}=1$, one solves (using a standard solver) the weighted least-squares problem

[TABLE]

and computes the residual $r_{j}=f(Z_{j})-\sum_{i=0}^{n}c_{i}g_{i}(Z_{j})$. The weights are then updated by $w_{j}:=w_{j}|r_{j}|$, followed by the re-normalization $w_{j}:=w_{j}/\sum_{i=1}^{M}w_{i}$. Iterating this process is known to converge linearly to the best polynomial approximant (with nontrivial convergence analysis [17]), and an acceleration technique is presented in [26].

8.3 AAA-Lawson

We now propose a rational variant of Lawson’s algorithm. (A similar attempt was made in [20, § 6.5], though the formulation there is not the same: most notably, adjusting the exponent $\gamma$ as done below appears to improve robustness significantly.) The idea is to incorporate Lawson’s approach into noninterpolatory AAA, replacing (39) with a weighted version of (38), and updating the weights as in Lawson.

Specifically, given an initial set of weights $w\in\mathbb{R}^{M-(\max(m,n)+1)}$, usually all ones, and initializing the Lawson exponent $\gamma=1$, we proceed as follows:

Solve the weighted linear least-squares problem

[TABLE]

via the SVD of the matrix $\mbox{diag}(\sqrt{w})\begin{bmatrix}C&-FC\end{bmatrix}$ (recall (38)). If the resulting $\left\|f(Z)-N(Z)/D(Z)\right\|_{\infty}$ is not smaller than before, then set $\gamma:=\gamma/2$. 2. 2.

Update $w$ by

[TABLE]

and return to step 1.

Note the exponent $\gamma$ in (41). In the linear case, this is $\gamma=1$. In the rational (nonlinear) case, for which experiments suggest convergence is a delicate issue, we have found that taking $\gamma$ to be smaller makes the algorithm much more robust. We repeat the steps until $w$ undergoes small changes, e.g. $10^{-3}$, or a maximum number of iterations (e.g. 30) is reached.

We refer to this algorithm as AAA-Lawson. Each iteration is computed by an SVD of an $(M-\max(m,n)-1)\times(m+n+2)$ matrix, so the cost for $k$ iterations is $O(kM(m+n)^{2})$. Convergence analysis appears to be highly nontrivial and is out of our scope. We simply note here that if equioscillation of $f-N/D$ is achieved at $m+n+2$ points in $Z_{*}\subset Z$, then by defining $w^{*}$ as $w_{j}^{*}=1/\sqrt{|D(Z_{j})|}$ for $j\in Z_{*}$ and [math] otherwise, we see that $w^{*}/\sum w^{*}$ (together with $N^{*}/D^{*}=r^{*}$, the solution of (2)) is a fixed point of the iteration.

8.4 Experiments with AAA-Lawson

Figure 8 compares AAA and AAA-Lawson (run for ten Lawson steps) for type (10,10) and (20,20) approximation of $f(x)=|x|$. The sample points are $10^{4}$ equispaced points on $[-1,1]$. Observe that the Lawson update significantly reduces the error and brings the error curve close to equioscillation.

AAA-Lawson is a new algorithm for rational minimax approximation. However, we do not recommend it as a practical means to obtain $r^{*}$ over the classical Remez or differential correction algorithms. The reason is that its convergence is far from understood, and even when it does converge, the rate is slow (linear at best). We illustrate this in Figure 9. In our Remez algorithm context, we take a small number (say 10) of AAA-Lawson steps to obtain a set of initial reference points, thereby taking advantage of the initial stage of the AAA-Lawson convergence.

We note that other approaches for rational approximation are available, which can be used for initializing Remez. These include the Loewner approach presented in [39] and RKFIT [6]. In particular, the Loewner approach is well suited when approximating smooth functions (and sometimes non-smooth functions like $f_{4}$ [36]), often achieving an error of the same order of magnitude as the best approximation. Our experiments suggest that AAA-Lawson is at least as efficient and robust as these alternatives.

8.5 Adaptive choice of support points

At an early stage of the AAA-Lawson iteration, we usually do not have the correct number ($m+n+2$) of reference (oscillation) points in the error curve. Therefore, choosing the support points $\{t_{k}\}$ as in (30) is not an option. Instead, we use the same support points chosen by the AAA algorithm, which is typically a good set. Once convergence sets in and the error curve of the AAA-Lawson iterates has at least $m+n+2$ alternation points, we can switch to the adaptive choice (30) as in Remez. We note, however, that adaptively changing the support points may further complicate the convergence, since it changes the linear least-squares problem (40).

8.6 Adaptive choice of the sample points

For solving the continuous problem (2), we take the sample point set $Z$ to be $M$ points uniformly distributed on $[a,b]$ ($M\lesssim 10^{5}$, chosen to keep the run time under control). Generally, it is necessary to sample more densely near a singularity if there is one; this is important e.g. for $f(x)=|x|$. We incorporate this need as follows: use AAA to find the support points $\{t_{k}\}$ (assume they are sorted), and take $M/n$ points between $[t_{k},t_{k+1}]$.

9 A barycentric version of the differential correction algorithm

The DC algorithm, due to Cheney and Loeb [16], has the great advantage of guaranteed global convergence in theory [3, 25], which applies whether the approximation domain $X$ is an interval $[a,b]$ or a finite set. It can also be extended to multivariate approximation problems [32]. In practice, however, it may suffer greatly from rounding errors, and its speed is often disappointing on larger problems. As we shall now describe, we have found that the first of these difficulties can be largely eliminated by the use of barycentric representations with adaptively chosen support points. The second problem of speed, however, remains, which is why ultimately we prefer the Remez algorithm for most problems.

9.1 The barycentric formulation

For an effective implementation, $X$ needs to be a finite set (e.g. obtained by discretizing $[a,b]$) to reduce each iteration to a linear programming (LP) problem. Considering the diagonal case $m=n$, a barycentric version of the DC algorithm can be defined recursively as follows. (We assume the support points are fixed to the values $t_{0},\ldots,t_{n}$, which do not belong to $X$.) Given $r_{k}=N_{k}/D_{k}\in\mathcal{R}_{n,n}(X)$, choose the partial fraction decompositions $N$ and $D$ of (4) that minimize the expression

[TABLE]

subject to

[TABLE]

and

[TABLE]

where $\delta_{k}=\max_{x\in X}\left|f(x)-r_{k}(x)\right|$. If $r=N/D$ is not good enough, continue with $r_{k+1}=r$. By imposing (44), we can establish convergence using an argument analogous to [3, Theorem 2]. In the polynomial basis setting, we know that the rate of convergence will ultimately be at least quadratic if the best approximation is non-degenerate [3, Theorem 3]. Non-diagonal approximations can be computed by adding the appropriate null space constraints as described in Section 4.4.

9.2 Choice of support points

Compared to the case of the barycentric Remez algorithm, changing the support points at each iteration of the DC algorithm makes it hard to impose a normalization condition similar to (44) or do a convergence analysis of the method. We therefore fix $\{t_{k}\}$ throughout the execution. The strategy we have adopted is based on Section 5.3: recursively construct type $(\ell,\ell)$ approximations with $\ell\leq n$. We take the set of support points of the $(\ell,\ell)$ problem based on a piecewise linear fit of the final reference points of the $(\ell-1,\ell-1)$ problem (similar to what is shown in Figure 2).

9.3 Experiments

We have implemented222The prototype code used is available at https://github.com/sfilip/barycentricDC. the barycentric DC algorithm in MATLAB using CVX [29] to specify the LP problems corresponding to (42)–(44), which are then solved using MOSEK’s [41] state-of-the-art LP optimizers. The four examples in Table 2 and Figure 10, for instance, demonstrate the effectiveness of the algorithm. For comparison, the sensitivity to the initial reference set prevented the convergence of our barycentric Remez implementation on all four of these examples. Function $f_{1}$ is particularly interesting since it is a version of Weierstrass’s classic example of a continuous but nowhere differentiable function.

Using a monomial or Chebyshev basis representation for the LP formulations quickly failed due to numerical errors, illustrating that the barycentric representation is crucial for the DC algorithm just as for the Remez algorithm.

We nevertheless echo the statement in the beginning of the section of the downsides of using the DC approach:

• •

Its overall cost. Producing the approximations in Figure 10 took several minutes in MATLAB on a desktop machine for each example.

• •

Numerical optimization tools for solving the corresponding LP problems break down at lower values of $m$ and $n$ than the ones we achieved with the barycentric Remez algorithm. We were usually able to go up to about type $(20,20)$.

10 Minimax approximation in Chebfun

We have presented many algorithmic details that have enabled the design of a fast and robust Remez implementation. In closing we remind readers that all this is available in Chebfun and readily explored in a few lines of code. Download Chebfun version 5.7.0 or later from GitHub or www.chebfun.org, put it in your MATLAB path, and then try for example

  [p,q,r] = minimax(@(x) abs(x),60,60);
  fplot(@(x) abs(x)-r(x),[-1 1])

In a few seconds a beautiful curve with 123 exponentially clustered equioscillation points will appear. Figure 11 summarizes our algorithm in a flowchart.