Silviu-Ioan Filip

Yuji Nakatsukasa
email: nakatsukasa@nii.ac.jp

Lloyd N Trefethen
email: tre-fethen@maths.ox.ac.uk

AND Bernhard Beckermann
email: bernhard.beckermann@univ-lille1.fr

RATIONAL MINIMAX APPROXIMATION VIA ADAPTIVE BARYCENTRIC REPRESENTATIONS

Keywords: barycentric formula, rational minimax approximation, Remez algorithm, differential correction algorithm, AAA algorithm, Lawson algorithm AMS subject classifications. 41A20, 65D15

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.

1. Introduction. The problem we are interested in is that of approximating functions f ∈ C([a, b]) using type (m, n) rational approximations with real coefficients, in the L ∞ setting. The set of feasible approximations is

R m,n = p q : p ∈ R m [x], q ∈ R n [x] . (1.1)
Given f and prescribed nonnegative integers m, n, the goal is to compute

min r∈Rm,n f -r ∞ , (1.2)
where • ∞ denotes the infinity norm over [a, b], i.e., f -r ∞ = max x∈[a,b] |f (x)r(x)|. The minimizer of (1.2) is known to exist and to be unique [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Ch. 24]. Let the minimax (or best) approximation be written r * = p/q ∈ R m,n , where p and q have no common factors. The number d = min {m -deg p, m -deg q} is called the defect of r * . It is known that there exists a so-called alternant (or reference) set consisting of ordered nodes a x 0 < x 1 < • • • < x m+n+1-d b, where f -r * takes its global extremum over [a, b] with alternating signs. In other words, we have the beautiful equioscillation property [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]Theorem 24.1] f (x ) -r * (x ) = (-1) +1 λ, = 0, . . . , m

+ n + 1 -d, (1.3)
where |λ| = f -r * ∞ . 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 (1.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), . . . , (m, n) with k = min {m, n} [11, p. 161].

In the approximation theory literature [START_REF] Braess | Nonlinear Approximation Theory[END_REF][START_REF] Cheney | Introduction to Approximation Theory[END_REF][START_REF] Meinardus | Approximation of Functions: Theory and Numerical Methods[END_REF][START_REF] Powell | Approximation Theory and Methods[END_REF][START_REF] Watson | Approximation Theory and Numerical Methods[END_REF], two algorithms are usually considered for the numerical solution of (1.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 [START_REF] Pachón | Barycentric-Remez algorithms for best polynomial approximation in the Chebfun system[END_REF]. 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 [START_REF] Driscoll | Chebfun Guide[END_REF] remez code. The Boost C++ libraries [START_REF]Boost C++ Libraries[END_REF] 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 [START_REF] Brophy | Synthesis of spectrum shaping digital filters of recursive design[END_REF][START_REF] Deczky | Equiripple and minimax (Chebyshev) approximations for recursive digital filters[END_REF] and the representation of special functions [START_REF] Cody | The FUNPACK package of special function subroutines[END_REF][START_REF] Cody | Algorithm 715: SPECFUN-a Portable FORTRAN Package of Special Function Routines and Test Drivers[END_REF]. 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 k, converge geometrically at the rate O(9.28903 • • • -k ) [START_REF] Nakatsukasa | Rational approximation of x n[END_REF].

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 [START_REF] Varga | Numerical results on best uniform rational approximation of |x| on [-1, +1][END_REF] using 200-digit arithmetic 1 .

A similar Remez iteration using the barycentric representation was described by Ionit , ȃ [START_REF] Ionit | Lagrange Rational Interpolation and its Applications to Approximation of Large-Scale Dynamical Systems[END_REF]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 Ionit , ȃ's treatment is that we reduce the core computational task to a symmetric eigenvalue problem, rather than a generalized eigenproblem as in [START_REF] Ionit | Lagrange Rational Interpolation and its Applications to Approximation of Large-Scale Dynamical Systems[END_REF]. The bigger difference is that Ionit , ȃ 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 [START_REF] Nakatsukasa | The AAA algorithm for rational approximation[END_REF] 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 ∞ 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 [START_REF] Curtis | The construction of minimax rational approximations to functions[END_REF], 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

r(z) = N (z) D(z) = n k=0 α k z -t k n k=0 β k z -t k , (2.1)
where n ∈ N, α 0 , . . . , α n and β 0 , . . . , β n are sets of real coefficients and t 0 , . . . , t n is a set of distinct real support points. The names N and D stand for "numerator" and "denominator". If we denote by ω t the node polynomial associated with t 0 , . . . , t n ,

ω t (z) = n k=0 (z -t k ), then p(z) = ω t (z)N (z) and q(z) = ω t (z)D(z) are both polynomials in 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 (µ, ν) with µ < n and/or ν < n.) At each point t k with nonzero α k or β k , formula (2.1) is undefined, but this is a removable singularity with lim z→t k r(z) = α k /β k (or a simple pole in the case α k = 0, β k = 0), meaning r is a rational interpolant to the values {α k /β k } at the support points {t k }.

Much of the literature on barycentric representations exploits this interpolatory property [START_REF] Berrut | Rational functions for guaranteed and experimentally well-conditioned global interpolation[END_REF][START_REF] Berrut | Recent developments in barycentric rational interpolation[END_REF][START_REF] Berrut | Barycentric Lagrange interpolation[END_REF][START_REF] Brezinski | Padé-type rational and barycentric interpolation[END_REF][START_REF] Floater | Barycentric rational interpolation with no poles and high rates of approximation[END_REF][START_REF] Schneider | Some new aspects of rational interpolation[END_REF] by taking α k = f (t k )β k , so that r is an interpolant to some given function values f (t 0 ), . . . , f (t n ) at the support points. In this case

r(z) = n k=0 f (t k )β k z -t k n k=0 β k z -t k , (2.2)
with the coefficients {β k } commonly known as barycentric weights; we have r(t k ) = f (t k ) as long as β k = 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 (2.1), 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 {α k } and {β k } [43, Theorem 2.1]. We remark that variants of (2.1) also form the basis for the popular vector fitting [START_REF] Gustavsen | Improving the pole relocating properties of vector fitting[END_REF][START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF] 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 (2.1) 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 = 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 [START_REF] Berrut | Matrices for the direct determination of the barycentric weights of rational interpolation[END_REF], 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 = ∞. 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

r(z) = ω t (z)N (z) ω t (z)D(z) = n k=0 (z -t k ) n k=0 α k z -t k n k=0 (z -t k ) n k=0 β k z -t k =: p(z) q(z) . (2.3)
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 α = [α 0 , . . . , α n ] T lies in a subspace spanned by the null space of the (transposed) Vandermonde matrix

V m =      1 1 • • • 1 t 0 t 1 • • • t n . . . . . . . . . t n-1-m 0 t n-1-m 1 • • • t n-1-m n      . (2.4)
That is, to enforce r ∈ R m,n with m < n, we require α ∈ span(P m ), where P m ∈ R (n+1)×(m+1) has orthonormal columns, obtained by taking the full QR factorization

V T m = P ⊥ m P m R m 0 , where P ⊥ m ∈ R (n+1)×(n-m) , R m ∈ R (n-m)×(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 (2.1), that is,

r(z) = m k=0 α k (z -t k ) -1 m k=0 β k (z -t k ) -1
, and force β ∈ span(P n ), where span(P n ) is the null space of the matrix

V n =      1 1 • • • 1 t 0 t 1 • • • t m . . . . . . . . . t m-1-n 0 t m-1-n 1 • • • t m-1-n m      , (2.5)
obtained by the QR factorization

V T n = P ⊥ n P n R n 0 , where P ⊥ n ∈ R (m+1)×(m-n) , R n ∈ R (m-n)×(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 (2.4)-(2.5) and real-valued Vandermonde matrices are usually highly ill-conditioned [START_REF] Beckermann | The condition number of real Vandermonde, Krylov and positive definite Hankel matrices[END_REF][START_REF] Beckermann | On the singular values of matrices with displacement structure[END_REF][START_REF] Pan | How bad are Vandermonde matrices?[END_REF], 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 [START_REF] Berrut | Matrices for the direct determination of the barycentric weights of rational interpolation[END_REF] 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 T m and V T n , we propose the following simpler approach, based on an Arnoldi-style orthogonalization:

1. Let Q = [1, . . . , 1] T when m > n, and Q = [f (t 0 ), . . . , f (t n )] T when m < n,
and normalize to have Euclidean norm 1. 2. Let q be the last column of Q. Take the projection of diag(t 0 , . . . , 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 ∈ C (max(m,n)+1)×(|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. Take the orthogonal complement Q ⊥ of Q via computing the QR factorization of Q. Q ⊥ is the desired matrix, P m or P n .

Note that the matrix Q in the final step is well conditioned (κ 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 {t k } 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 {t k } 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 2.1.

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 (2.1)) shows that evaluating r(x) is backward stable in the sense that the computed value r(x) satisfies

r(x) = n k=0 α k (1 + α k ) x -t k n k=0 β k (1 + β k ) x -t k , (2.6)
where α k , β k denote quantities of size O(u), or more precisely, bounded by (1 + u) 3n+4 . In other words, r(x) is an exact evaluation of (2.1) for slightly perturbed {α k }, {β k }. Note that when r represents a polynomial (as assumed in [START_REF] Higham | The numerical stability of barycentric Lagrange interpolation[END_REF]), (2.6) does not imply backward stability. However, as a rational function for which we allow for backward errors in the denominator, (2.6) does imply backward stability.

For the forward error, we can adapt the analysis of [START_REF] Celis | Practical Rational Interpolation of Exact and Inexact Data[END_REF]Proposition 2.4.3]. Assume that the computed coefficients α, β are obtained through a backward stable process,

α k = α k (1 + δ α k ), δ α k = O(κ α u), β k = β k (1 + δ β k ), δ β k = O(κ β u), k = 0, . . . , n,
where κ α and κ β are condition numbers associated with the matrices used to determine α and β. Then, if x (the evaluation point) and {t k } are considered to be floating point numbers, we have Lemma 2.1. The relative forward error for the computed value r(x) of (2.1) satisfies 3. The rational Remez algorithm. Initially developed by Werner [START_REF] Werner | Die konstruktive Ermittlung der Tschebyscheff-Approximierenden im Bereich der rationalen Funktionen[END_REF][START_REF] Werner | Tschebyscheff-Approximation im Bereich der rationalen Funktionen bei Vorliegen einer guten Ausgangsnäherung[END_REF] and Maehly [START_REF] Maehly | Methods for fitting rational approximations, Parts II and III[END_REF], the rational Remez algorithm extends the ideas of computing best polynomial approximations due to Remez [START_REF] Remes | Sur le calcul effectif des polynômes d'approximation de Tchebichef[END_REF][START_REF] Remes | Sur un procédé convergent d'approximations successives pour déterminer les polynômes d'approximation[END_REF]. It can be summarized as follows:

r(x) -r(x) r(x) ≤ u(n+3+O(κ α )) n k=0 α k x-t k n k=0 α k x-t k +u(n+2+O(κ β )) n k=0 β k x-t k n k=0 β k x-t k +O(u 2 ).
Step 1 Set k = 1 and choose m + n + 2 distinct reference points a ≤ x (k) 0 < • • • < x (k) m+n+1 ≤ b.
Step 2 Determine the levelled error λ k ∈ R (positive or negative) and r k ∈ R m,n such that r k has no pole on [a, b] and

f (x (k) ) -r k (x (k) ) = (-1) +1 λ k , = 0, . . . , m + n + 1. (3.1)
Step 3 Choose as the next reference m + n + 2 local maxima {x

(k+1) } of |f -r k | such that s(-1) f (x (k+1) ) -r k (x (k+1) ) ≥ |λ k | , = 0, . . . , m + n + 1, (3.2)
with s ∈ {±1} and such that for at least one ∈ {0, . . . , m + n + 1}, the lefthand side of (3.2) equals f -r k ∞ . If r k has converged to within a given threshold

ε t > 0 (i.e., ( f -r k ∞ -λ k )/ f -r k ∞ ≤ ε t [50, eq. (10.8)]) return r k , else go to Step 2 with k ← k + 1. If
Step 2 is always successful, then convergence to the best approximation is assured [START_REF] Watson | Approximation Theory and Numerical Methods[END_REF]Theorem 9.14]. It might happen that Step 2 fails, namely when all rational solutions satisfying the equations (3.1) 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 [START_REF] Curtis | Theory and calculation of best rational approximations[END_REF]. Given a positive weight function

w ∈ C([a, b]), the goal is to find r * ∈ R m,n such that the weighted error f -r * w,∞ = max x∈[a,b] |w(x)(f (x) -r * (x))| is minimal. Equations (3.1) and (3.2) get modified to w(x (k) ) f (x (k) ) -r k (x (k) ) = (-1) +1 λ k , = 0, . . . , m + n + 1 and s(-1) w(x (k+1) ) f (x (k+1) ) -r k (x (k+1) ) ≥ |λ k | , = 0, . . . , m + n + 1,
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 ∈ R n,n , p, q ∈ R n [x] such that f (x ) -r(x ) = (-1) +1 λ, = 0, . . . , 2n + 1 (4.1)
and assume that we represent p and q using a basis of polynomials ϕ 0 , . . . , ϕ n such that span

R (ϕ i ) 0≤i≤n = R n [x]: p(x) = n k=0 c p,k ϕ k (x), q(x) = n k=0 c q,k ϕ k (x).
The linearized version of (4.1) is then given by

p(x ) = q(x ) f (x ) -(-1) +1 λ ,
which, in matrix form, becomes

Φ x c p =           f (x 0 ) f (x 1 ) . . . f (x 2n+1 )      -λ      -1 1 -1 . . .           Φ x c q , (4.2)
where andc p = [c p,0 , c p,1 , . . . , c p,n ] T and c q = [c q,0 , c q,1 , . . . , 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 = diag (d 0 , . . . , d 2n+1 ), (4.2) can also be written as a generalized eigenvalue problem

Φ x ∈ R (2n+2)×(n+1) is the basis matrix (Φ x ) ,k = ϕ k (x ), 0 ≤ ≤ 2n + 1, 0 ≤ k ≤ n,
DΦ x -F DΦ x c p c q = λ 0 -SDΦ x c p c q , (4.3)
with F = diag (f (x 0 ), . . . , f (x 2n+1 )) and S = diag (-1) k+1 .

As described in Powell [50, Ch. 10.2], solving (4.3) is usually done by eliminating c p . His presentation considers the monomial basis, but the approach is valid for any basis of R n [x]. By taking the full QR decomposition of DΦ x , we get

DΦ x = Q 1 Q 2 R 0 = Q 1 R. Since DΦ x is of full rank, we have Q 1 , Q 2 ∈ R (2n+2)×(n+1) and Q T 2 Q 1 = 0. By multi- plying (4.3) on the left by Q T = Q 1 Q 2 T
, we obtain a block triangular eigenvalue problem with lower-right (n + 1) × (n + 1) block

Q T 2 F Q 1 Rc q = λQ T 2 SQ 1 Rc q . ( 4.4)
(The top-left (n+1)×(n+1) block has all eigenvalues at infinity, and is thus irrelevant.)

In terms of polynomials, (Q

1 ) ,k = d ψ k (x ), 0 ≤ k ≤ n, 0 ≤ ≤ 2n + 1
, where (ψ k ) 0≤k≤n is a family of orthonormal polynomials with respect to the discrete inner product f,

g x = 2n+1 k=0 d 2 k f (x k )g(x k ). Moreover, if (ϕ k ) 0≤k≤n is a degree-graded basis with deg ϕ k = k, then we have deg ψ k = k, 0 ≤ k ≤ n.
Let ω x be the node polynomial associated with the reference nodes x 0 , . . . , x 2n+1 , and Ω x = diag (1/ω x (x 0 ), . . . , 1/ω x (x 2n+1 )). We have [50, p. 114]

V T x Ω x V x = 0, (4.5)
where ) is the Vandermonde matrix associated with x 0 , . . . , x 2n+1 , that is, (V x ) i,j = x j i . Indeed,

V x ∈ R (2n+2)×(n+1
V T x Ω x V x i,j = 2n+1 =0
x i+j 1 ω x (x ) = (x i+j )[x 0 , . . . , x 2n+1 ] = 0, i, j ∈ {0, . . . , n} , the divided differences of order 2n + 1 of the function x i+j at the {x } nodes, hence 0 if i + j ≤ 2n. By using the appropriate change of basis matrix in (4.5), we have

Φ T x Ω x Φ x = 0. (4.6)
Now, by multiplying (4.3) on the left by Φ T

x Ω x D -1 and using (4.6), we can eliminate the c p term to obtain

Φ T x Ω x F Φ x c q = λΦ T x Ω x SΦ x c q . (4.7)
Equation (4.7) is the extension of [50, Eq. (10.13)] from the monomial basis to ϕ 0 , . . . , ϕ n . Moreover, we have: Lemma 4.1. The matrix Φ T x Ω x SΦ x is symmetric positive definite. Proof. Since Ω x S = |Ω x |, it means that Ω x S is symmetric positive definite, and the conclusion follows. See also [START_REF] Powell | Approximation Theory and Methods[END_REF]Theorem 10.2].

Since Φ T

x Ω x F Φ x is also symmetric, it follows that all eigenvalues of (4.7) 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 . Then, from (4.1), it follows that either r(x k ) -r (x k ) are all zero or they alternate in sign at least 2n + 1 times. In both cases, r -r ∈ R 2n,2n has at least 2n + 1 zeros inside [a, b], leading to r = r .

We can in fact transform (4.4) into a symmetric eigenvalue problem (an observation which seems to date to [START_REF] Pelios | Rational function approximation as a well-conditioned matrix eigenvalue problem[END_REF]) by considering the choice

D = |Ω x | 1/2 , which leads to Q 2 = SQ 1 in view of (4.6). The system becomes Q T 1 SF Q 1 Rc q = λQ T 1 S 2 Q 1 Rc q , which, by the change of variables y = Rc q , gives Q T 1 SF Q 1 y = λy. (4.8)
To get c p , from (4.2), we have

|Ω x | 1/2 Φ x c p = (F -λS) |Ω x | 1/2 Φ x c q ,
or equivalently (by multiplication on the left by Q T 1 ),

Rc p = Q T 1 (F -λS) |Ω x | 1/2 Φ x c q = Q T 1 F Q 1 y.
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 ψ 0 , . . . , ψ n . The (scaled) values of the denominator at each x k corresponding to an eigenvector y can be recovered by computing

|Ω x | 1/2 Φ x c q = Q 1 y. (4.9)
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 (2.1). Namely, (4.1) becomes

Cα =           f (x 0 ) f (x 1 ) . . . f (x 2n+1 )      -λ      -1 1 -1 . . .           Cβ, (4.10)
where C is now a (2n + 2) × (n + 1) Cauchy matrix with entries C ,k = 1/(x -t k ) (we assume for the moment {x } ∩ {t k } = ∅) and α = [α 0 , α 1 , . . . , α n ] T and β = [β 0 , β 1 , . . . , β n ] T are the column vectors of coefficients {α k } and {β k }. Again, this can be transformed into a generalized eigenvalue problem

C -F C α β = λ 0 -SC α β . (4.11)
To reduce (4.11) to a symmetric eigenvalue problem as in (4.8), we form a link between the monomial and barycentric representations in terms of the basis matrices V x and C. We have: Lemma 4.2. Let V x , ω t be as defined above, and V t ∈ R (n+1)×(n+1) be the Vandermonde matrix corresponding to the support points, i.e., (V t ) i,j = t j i . Then

diag 1 ω t (x 0 ) , . . . , 1 ω t (x 2n+1 ) V x = C diag 1 ω t (t 0 ) , . . . , 1 ω t (t n ) V t .
Proof. If we look at an arbitrary element of the right-hand side matrix, we have

C diag 1 ω t (t 0 ) , . . . , 1 ω t (t n ) V t j, = n k=0 1 (x j -t k )ω t (t k ) t k = x j ω t (x j ) ,
where the second equality is a consequence of the Lagrange interpolation formula.

In place of Ω x we will use the following matrix ∆:

Lemma 4.3. If ∆ = diag ω t (x 0 ) 2 , . . . , ω t (x 2n+1 ) 2 Ω x , then C T ∆C = 0.
Proof. We apply Lemma 4.2 and use the fact that V T

x Ω x V x = 0. Namely,

C T ∆C = diag (ω t (t 0 ), . . . , ω t (t n )) V -T t V T x Ω x V x V -1
t diag (ω t (t 0 ), . . . , ω t (t n )) = 0. We now take the full QR decomposition of |∆| 1/2 C = (S∆) 1/2 C. We have

|∆| 1/2 C = Q 1 Q 2 R 0 = Q 1 R.
Based on Lemma 4.3, we can again take Q 2 = SQ 1 . From (4.11) we get

|∆| 1/2 C -F |∆| 1/2 C α β = λ 0 -S |∆| 1/2 C α β .
Multiplying this expression on the left by Q 1 Q 2 T gives a block triangular matrix pencil, whose (n + 1) × (n + 1) lower-right corner is the barycentric analogue of (4.4):

Q T 2 F Q 1 Rβ = λQ T 2 SQ 1 Rβ. After substituting Q T 2 = Q T 1 S, we get Q T 1 (SF )Q 1 Rβ = λQ T 1 S 2 Q 1 Rβ, (4.12)
which, by the change of variable y = Rβ, becomes a standard symmetric eigenvalue problem in λ with eigenvector y (recall that S, F are diagonal):

Q T 1 (SF )Q 1 y = λy. (4.13)
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 diag (ω t (x 0 ), . . . , ω t (x 2n+1 ))

Cβ = diag (ω t (x 0 ), . . . , ω t (x 2n+1 )) |∆| -1/2 Q 1 y. (4.14)
As in the polynomial case, there is at most one solution such that q(x) = D(x)ω t (x) has no root in [a, b]; indeed, (4.9) and (4.14) represent the values of q(x ) for r = p/q and x satisfying equation (4.1). We use this sign test involving (4.14) to determine the levelled error λ that gives a pole-free r in Step 2 of our rational Remez algorithm. The appropriate β is then taken by solving Rβ = y. From (4.10), we have

|∆| 1/2 Cα = (F -λS) |∆| 1/2 Cβ,
or equivalently (by multiplication on the left by Q T 1 )

Rα = Q T 1 (F -λS) |∆| 1/2 Cβ = Q T 1 (F -λS)Q 1 y = Q T 1 F Q 1 y, (4.15)
which allows us to recover α (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 (4.11) and (4.13) correspond to

C -F C α β = λ 0 -SW -1 C α β and Q T 1 (SF )Q 1 y = λQ T 1 W -1 Q 1 y
, where W = diag (w(x 0 ), . . . , w(x 2n+1 )) 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 α, (4.15) becomes Rα = Q T 1 (F -λ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

|∆| 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 |∆| 1/2 C, as this is the key matrix for constructing (4.12). 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 |∆| where D n+1 is the set of (n + 1) × (n + 1) diagonal matrices. We have Theorem 4.4. Let t k ∈ (x 2k , x 2k+1 ) for k = 0, . . . , n and

s k ∈ (x 2k+1 , x 2k+2 ) for k = 0, . . . , n -1, s n ∈ (x 2n+1 , ∞), and define ω s (x) = n k=0 (x -s k ). Then min Γ∈Dn+1 κ 2 (|∆| 1/2 CΓ) ≤ max ω s (x ) ω t (x ) • max k ω t (x k ) ω s (x k ) . (4.17)
Proof. Let {y j } be a (2n + 2)-element set such that y j ∈ (x j , x j+1 ), j = 0, . . . , 2n, y 2n+1 > x 2n+1 and let C x,y ∈ R (2n+2)×(2n+2) be the Cauchy matrix with elements (C x,y ) j,k = 1/(x j -y k ). If we consider D 1 = diag( |ω y (x j )/ω x (x j )|) and D 2 = diag( ω x (y j )/ω y (y j ) ), 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 = k, we have

-(G) j,k = ω x (y j )ω x (y k ) ω y (y j )ω y (y k ) 2n+1 =0 ω y (x ) (x -y j )(x -y k )ω x (x ) = ω x (y j )ω x (y k ) ω y (y j )ω y (y k ) ω y (x) (x -y j )(x -y k ) [x 0 , . . . , x 2n+1 ] = 0.
Similarly, since j =k (x -

y j ) = q(x)(x -y k ) + ω y (y k ), with q ∈ R 2n [x], we have, -(G) k,k = ω x (y k ) ω y (y k ) 2n+1 =0 ω y (x ) (x -y k ) 2 ω x (x ) = ω x (y k ) ω y (y k ) j =k (x -y j ) x -y k [x 0 , . . . , x 2n+1 ] = ω x (y k ) ω y (y k ) q(x) + ω y (y k ) x -y k [x 0 , . . . , x 2n+1 ] = ω x (y k ) 1 x -y k [x 0 , . . . , x 2n+1 ] = ω x (y k ) -1 ω x (y k ) = -1.
Now, if we take t k = y 2k , s k = y 2k+1 , for k = 0, . . . , n, there exist D ∈ D 2n+2 and Γ ∈ D n+1 such that |∆| 1/2 CΓ = DD 1 C x,y D 2 I t , where D = diag( |ω t (x j )/ω s (x j )|) and I t is obtained by removing every second column from I 2n+2 . In particular, Γ = I T t D 2 I t . It follows that

κ 2 (|∆| 1/2 CΓ) ≤ κ 2 (D) = max ω s (x ) ω t (x ) • max k ω t (x k ) ω s (x k ) .
Let Γ = I T t D 2 I t be as in the proof of Theorem 4.4. It turns out that for the choice t k = x 2k+1 -ε, s k = x 2k+1 + ε, for k = 0, . . . , n, as ε → 0, the matrix |∆| 1/2 C has a finite limit C of full column rank, and similarly Γ tends to some diagonal matrix Γ with positive diagonal entries. From Theorem 4.4 and its proof we know that C Γ has condition number 1, and, more precisely, orthonormal columns. We thus obtain an explicit thin QR decomposition of C (by direct calculation):

Corollary 4.5. In the limit t k x 2k+1 , for k = 0, . . . , n, the matrix |∆| 1/2 C converges to C, with entries

( C) j,k =        |w t (t k )| √ |w x (t k )| if j = 2k + 1, 0 if j = 2 + 1, = k, |wt(xj )| √ |w x (xj )| /(x j -t k ) if j = 2 ,
and explicit thin QR decomposition C = Q 1 R, where 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 4.5. In addition, we reduce the problem to a symmetric eigenvalue problem (4.13), resulting in well-conditioned eigenvalues, with β being obtained by solving the diagonal system Rβ = y with y as in (4.13). Compared to (4.1), where we want q to have the same sign over {x }, we similarly require that β and thus y have components alternating in sign, which uniquely fixes the norm 1 eigenvector y in (4.13). Our approach also allows for nondiagonal types, as we describe next.

(Q 1 ) j,k =        1/ √ 2 if j = 2k + 1, 0 if j = 2 + 1, = k, wt(xj ) w t (t k ) w x (t k ) 2w x (xj ) /(x j -t k ) if j =

4.4.

The nondiagonal case m = n. As pointed out in Section 2.1, when searching for a best approximant with m = n, we need to force the coefficient vector α or β to lie in a certain subspace. This results in modified versions of (4.11). Namely,

C -F CP n α β = λ 0 -SCP n α β , when m > n, (4.19)
for β ∈ C n+1 , and we take β = P n β. Similarly,

CP m -F C α β = λ 0 -SC α β , when m < n, (4.20)
for α ∈ C m+1 , and we take α = P m α.

Below we describe the reduction of the generalized eigenvalue problems (4. [START_REF] Cody | Algorithm 715: SPECFUN-a Portable FORTRAN Package of Special Function Routines and Test Drivers[END_REF]) and (4.20) to standard symmetric eigenvalue problems.

Case m > n. In this case, C ∈ R (m+n+2)×(m+1) . Since det |∆| 1/2 = 0, (4. [START_REF] Cody | Algorithm 715: SPECFUN-a Portable FORTRAN Package of Special Function Routines and Test Drivers[END_REF]) is equivalent to the generalized eigenvalue problem

|∆| 1/2 C -F |∆| 1/2 CP n α β = λ 0 -S |∆| 1/2 CP n α β . (4.21) Consider the (thin) QR decomposition of |∆| 1/2 C P n P ⊥ n = (S∆) 1/2 C P n P ⊥ n : |∆| 1/2 C P n P ⊥ n = Q 1 Q 2 R = Q 1 Q 2 R 1 R 12 0 R 2 .
Then we have the identity Q 1 Q 2 T (SQ 1 ) = 0, as can be verified analogously to (4.5) using divided differences. This implies (SQ 1 ) T |∆| 1/2 C = 0, so by left-multiplying (4.21) by (SQ 1 ) ⊥ SQ 1 T we obtain a block upper-triangular eigenvalue problem with lower-right (n + 1) × (n + 1) block

(SQ 1 ) T F Q 1 R 1 β = λ(SQ 1 ) T SQ 1 R 1 β,
which again reduces to the standard symmetric eigenvalue problem (setting

y = R 1 β) Q T 1 (SF )Q 1 y = λy. (4.22)
From (4.21), we have |∆|

1/2 Cα = (F -λS) |∆| 1/2 CP n β. Left-multiplying by Q 1 Q 2 T and using Q 1 Q 2 T S |∆| 1/2 CP n = 0, we obtain R P n P ⊥ n T α = Q 1 Q 2 T F |∆| 1/2 CP n β = Q 1 Q 2 T F Q 1 R 1 β = Q 1 Q 2 T F Q 1 y. Therefore α = P n P ⊥ n R -1 Q 1 Q 2 T F Q 1 y,
which is obtained by computing the vector y = Q 1 Q 2 T F Q 1 y, then solving R y = y for y, then α = P n P ⊥ n y.

Case m < n. This case is analogous to the previous one; we highlight the differences. C is a (m + n + 2) × (n + 1) matrix. Equation (4.20) is equivalent to

|∆| 1/2 CP m -F |∆| 1/2 C α β = λ 0 -S |∆| 1/2 C α β . (4.23)
Consider the (thin) QR decompositions

|∆| 1/2 C = (S∆) 1/2 C = Q 1 R, |∆| 1/2 CP m = Q 1 R.
Here m+1) . We have Q T 1 (SQ 1 ) = 0, which again can be established using divided differences. This implies (SQ 1 ) T |∆| 1/2 CP m = 0, so left-multiplying equation (4.23) by (SQ 1 ) ⊥ SQ 1 T results in a block uppertriangular eigenvalue problem with lower-right block

Q 1 ∈ R (m+n+2)×(n+1) , Q 1 ∈ R (m+n+2)×(
(SQ 1 ) T F Q 1 Rβ = λ(SQ 1 ) T SQ 1 Rβ,
which also reduces to the standard symmetric eigenvalue problem (setting y = Rβ)

Q T 1 (SF )Q 1 y = λy. (4.24)
From (4.23), we have

|∆| 1/2 CP m α = (F -λS) |∆| 1/2 Cβ. Left-multiplying by Q T 1 and using Q T 1 S |∆| 1/2 C = 0, we obtain R α = Q T 1 F |∆| 1/2 Cβ = Q T 1 F Q 1 Rβ = Q T 1 F Q 1 y. Therefore α = R -1 Q T 1 F Q 1 y, obtained via y = Q T 1 F Q 1 y
, then solving the linear system R α = 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 (4.22) and (4.24) become

Q T 1 (SF )Q 1 y = λQ T 1 W -1 Q 1 y,
where y = R 1 β when m > n and y = Rβ when m < n. Recovering the set of barycentric coefficients in the numerator corresponds to solving the systems

α = P n P ⊥ n R -1 Q 1 Q 2 T (F -λSW -1 )Q 1 y, m > n and α = R -1 Q T 1 (F -λSW -1 )Q 1 y, m < n.
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 4.5, the computation of Q 1 can be done explicitly, and the linear system y = Rβ is diagonal, hence can be solved with high relative accuracy. The main source of numerical errors is therefore in the symmetric eigenvalue problem (4.13), (4.22) or (4.24). 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 λ is computed with O(u) accuracy, assuming for simplicity that f ∞ = 1 (without loss of generality). The eigenvector, on the other hand, has conditioning O(1/gap) [57, Ch. V], where gap is the distance between the desired λ and the rest of the eigenvalues. These eigenvalues are equal to those of the nonzero eigenvalues of the generalized eigenproblem (4.3), 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 gap = O(|λ|). Hence the computed eigenvector tends to have accuracy O(u/|λ|), and if the eigenvector y has small elements, the componentwise relative accuracy may be worse. The computation therefore breaks down (perhaps as expected) when |λ| = O(u), that is, when the error curve has amplitude of size machine precision. When m = n, (4.18) 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 Γ∈Dn+1 κ 2 (|∆| 1/2 CP m,n Γ) > 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 }. A key element of our approach is the AAA-Lawson algorithm, which can efficiently find an approximate solution to the minimax problem (1.2) (to low accuracy).

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

We first attempt to compute the CF approximant [START_REF] Trefethen | The Carathéodory-Fejér method for real rational approximation[END_REF][START_REF] Van Deun | A robust implementation of the Carathéodory-Fejér method for rational approximation[END_REF] 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.

AAA-Lawson approximation.

This approach is based on the AAA algorithm [START_REF] Nakatsukasa | The AAA algorithm for rational approximation[END_REF] 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 [START_REF] Meinardus | Approximation of Functions: Theory and Numerical Methods[END_REF][START_REF] Meinardus | Approximation of Functions: Theory and Numerical Methods[END_REF]. The entire process relies on a moderate number of SVDs (say max(m, n) + 10). 7.1, with a singularity at x = 1/ √ 2. The y components of the red crosses correspond to the final references {x } for the (m , n ) = (10, 10) best approximation, while the y components of the black circles are the initial guess {x } for the (m , n ) = [START_REF] Braess | Nonlinear Approximation Theory[END_REF][START_REF] Braess | Nonlinear Approximation Theory[END_REF] problem, taken based on the piecewise linear fit at {x }. Note how the y components of both sets of points cluster near the singularity.

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 } corresponding to best approximations in (1.3) tend to cluster near these singularities as m and n increase.

It is sensible to expect that first computing a type (m , n ) best approximation to f with m m and n 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 and n by small increments j, typically j ∈ {1, 2, 4}. The steps taken follow a diagonal path, as explained in Figure 5.1. Note that in addition to improving the robustness of the Remez algorithm, this strategy can help detect degeneracy; recall the discussion after (1.3). It proves useful for many examples, including some of those shown in Section 7: type (n, n) approximations to f (x) = |x|, x ∈ [-1, 1] for n > 40 and the f 1 , f 2 and f 4 specifications in Table 7.1.

6.

Searching for the new reference. We now turn to the updating strategy for the reference points x 0 . . . , 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 [x , x +1 ] (and [a, x0 ] and [x m+n+1 , b], if non-degenerate; here {x } 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 [START_REF] Pachón | Piecewise-smooth chebfuns[END_REF]), 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 = k i=0 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 (x) = k i=1 ic i U i-1 (x) (using the formula T n (x) = nU n-1 (x)) via the eigenvalues of the colleague matrix for Chebyshev polynomials of the second kind [START_REF] Good | The colleague matrix, a Chebyshev analogue of the companion matrix[END_REF]. 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 We split this interval into subintervals separated by the previous reference points (red circles), and approximate e on each subinterval by a low-degree polynomial. We then find the roots of its derivative.

and select those with the largest values that satisfy (3.2).

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 [START_REF] Trefethen | Approximation Theory and Approximation Practice[END_REF]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 ≥ 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 [START_REF] Newman | Rational approximation to[END_REF]. The convergence rate for best type (n, n) approximations was later shown by Stahl [START_REF] Stahl | Best uniform approximation of |x| on [-1, 1[END_REF] to be E n,n (|x|, [-1, 1]) ∼ 8e -π √ n . This result had in fact been conjectured by Varga, Ruttan and Carpenter [START_REF] Varga | Numerical results on best uniform rational approximation of |x| on [-1, +1][END_REF] 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 ( √ 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 7.1), 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 7.1, where our levelled error computation for the type (80, 80) approximation (value 4.39 . . . × 10 -12 ) matches the corresponding error of [START_REF] Varga | Numerical results on best uniform rational approximation of |x| on [-1, +1][END_REF]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 ≤ 8).

The robustness of our algorithm is also illustrated by the examples of Table 7.1 and Figure 7.2, 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 [START_REF] Pushnitski | Best rational approximation of functions with logarithmic singularities[END_REF]. 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 ) accu- Table 7.1: Best approximation to five difficult functions by the barycentric rational Remez algorithm. f 1 is discontinuous at x = 1/ √ 2, f 2 is discontinuous at x = 0, f 3 is unbounded as x → 0, f4 has two sharp peaks at x = ±0.6, and f5 has a logarithmic singularity at x = 0. . Furthermore, Zolotarev gave explicit solutions involving Jacobi's elliptic functions. These rational functions have the remarkable property of preserving optimality under appropriate composition [START_REF] Nakatsukasa | Computing fundamental matrix decompositions accurately via the matrix sign function in two iterations: The power of Zolotarev's functions[END_REF]. In Figure 7.3 we compute the best relative error approximant of type (m, m) to √ 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.4, confirming Zolotarev's results.

i fi [a, b] (m, n) f -r * ∞ 1 x 2 , x < 1 √ 2 -x 2 + 2 √ 2x -1, 1 √ 2 ≤ x [0,
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 * ∞ is a digit or two larger than u f ∞ . x and a type [START_REF] Cline | Rate of convergence of Lawson's algorithm[END_REF][START_REF] Cline | Rate of convergence of Lawson's algorithm[END_REF] rational approximation. We plot the absolute error curve on the left, while the relative error (right), matching our choice of w, gives an expected equioscillating curve. This is Zolotarev's third problem. [20, § 6.5], though the formulation there is not the same: most notably, adjusting the exponent γ as done below appears to improve robustness significantly.) The idea is to incorporate Lawson's approach into noninterpolatory AAA, replacing (8.3) with a weighted version of (8.2), and updating the weights as in Lawson.

Specifically, given an initial set of weights w ∈ R M -(max(m,n)+1) , usually all ones, and initializing the Lawson exponent γ = 1, we proceed as follows:

1. Solve the weighted linear least-squares problem minimize

α 2 2 + β 2 2 =1 f ( Z)D( Z) -N ( Z) w , (8.4)
via the SVD of the matrix diag( √ w) C -F C (recall (8.2)). If the resulting f (Z) -N (Z)/D(Z) ∞ is not smaller than then set γ := γ/2.

Update w by

w j ← w j f (Z j ) - N (Z j ) D(Z j ) γ , ∀j, then w j := w j i w i (8.5)
and return to step 1.

Note the exponent γ in (8.5). In the linear case, this is γ = 1. In the rational (nonlinear) case, for which experiments suggest convergence is a delicate issue, we have found that taking γ 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 (Mmax(m, n) -1) × (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 * ⊂ Z, then by defining w * as w * j = 1/ |D(Z j )| for j ∈ Z * and 0 otherwise, we see that w * / w * (together with N * /D * = r * , the solution of (1.2)) is a fixed point of the iteration. 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 8.2. 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 [START_REF] Mayo | A framework for the solution of the generalized realization problem[END_REF] and RKFIT [START_REF] Berljafa | The RKFIT algorithm for nonlinear rational approximation[END_REF]. In particular, the Loewner approach is well suited when approximating smooth functions (and sometimes non-smooth functions like f 4 [START_REF] Karachalios | Hyperbolic function and the Loewner framework[END_REF]), 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 (4.18) 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 (4.18) 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 (8.4). 8.6. Adaptive choice of the sample points. For solving the continuous problem (1.2), we take the sample point set Z to be M points uniformly distributed on [a, b] (M 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 [START_REF] Cheney | Two new algorithms for rational approximation[END_REF], has the great advantage of guaranteed global convergence in theory [START_REF] Barrodale | The differential correction algorithm for rational ∞-approximation[END_REF][START_REF] Dua | Further remarks on the differential correction algorithm[END_REF], 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 [START_REF] Hettich | An algorithm for general restricted rational Chebyshev approxima-tion[END_REF]. 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. 3), we can establish convergence using an argument analogous to [START_REF] Barrodale | The differential correction algorithm for rational ∞-approximation[END_REF]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 nondegenerate [START_REF] Barrodale | The differential correction algorithm for rational ∞-approximation[END_REF]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 (9.3) 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 ( , ) approximations with ≤ n. We take the set of support points of the ( , ) problem based on a piecewise linear fit of the final reference points of the ( -1, -1) problem (similar to what is shown in Figure 5.1). 9.3. Experiments. We have implemented2 the barycentric DC algorithm in MATLAB using CVX [START_REF] Grant | CVX: Matlab software for disciplined convex programming[END_REF] to specify the LP problems corresponding to (9.1)-(9.3),

Fig. 2 . 1 :

 21 Fig.2.1: Linear (top) and semilogy (bottom) plots of q and |D| in r * = p/q = N/D, the best rational approximation for |x| of type (m, n) =[START_REF] Cooper | Rational Approximation of Discrete Data with Asymptomatic Behaviour[END_REF][START_REF] Cooper | Rational Approximation of Discrete Data with Asymptomatic Behaviour[END_REF]. Here p, q are the polynomials in the classical quotient representation (1.1), and D is the denominator in the barycentric representation (2.1). The dots are the equioscillation points {x }, while the set of support points {t k } consists of every other point in {x }.

(2. 7 )

 7 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 (2.7) do not become too small over [a, b], then we can expect the evaluation of 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) ≈ f (x) implies |N (x)| ≈ |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 (2.7) measures the relative error, we usually cannot expect it to be O(u) when |r(x)| ≈ |f (x)| 1.

4. 5 .

 5 Adaptive choice of the support points. Theorem 4.4 gives an optimal choice of support points t k = x 2k+1 in terms of optimizing min Γ∈Dn+1 κ 2 (|∆| 1/2 CΓ). In Section 2.2 we discussed another desideratum for the support points {t k }: the resulting |D(x )| = |q(x ) n k=0 (x -t k )| should take uniformly large values for all . Fortunately, this requirement is also met with this choice, as was illustrated in Figure2.1.

(m 1 , n 1 )Fig. 5 . 1 :

 1151 Fig. 5.1: Initialization with lower degree approximations. The left plot shows the three possible paths for updating the degrees (assuming the increment is j = 1): m < n (red), m = n (black) and m > n (blue). The right plot shows how initialization is done at an intermediate step. The function is f1 from Table7.1, with a singularity at x = 1/ √ 2. The y components of the red crosses correspond to the final references {x } for the (m , n ) = (10, 10) best approximation, while the y components of the black circles are the initial guess {x } for the (m , n ) =[START_REF] Braess | Nonlinear Approximation Theory[END_REF][START_REF] Braess | Nonlinear Approximation Theory[END_REF] problem, taken based on the piecewise linear fit at {x }. Note how the y components of both sets of points cluster near the singularity.

Fig. 6 . 1 :

 61 Fig. 6.1: Illustration of how a new set of reference points (black stars) is found from the current error function e = fr (blue curve). Shown here is the error curve after three Remez iterations in finding the best type (10, 10) approximation to f (x) = |x| on [-1, 1].We split this interval into subintervals separated by the previous reference points (red circles), and approximate e on each subinterval by a low-degree polynomial. We then find the roots of its derivative.

12 Fig. 7 . 1 :

 1271 Fig.7.1: In the first plot, the upper dots show the best approximation errors for the degree 2n best polynomial approximations of |x| on [-1, 1], while the lower ones correspond to the best type (n, n) rational approximations, superimposed on the asymptotic formula from[START_REF] Stahl | Best uniform approximation of |x| on [-1, 1[END_REF]. The bottom plot shows the minimax error curve for the type (80, 80) best approximation to |x|. Note that the horizontal axis has a log scale: the alternant ranges over 11 orders of magnitude. The positive part of the domain [-1, 1] is shown (by symmetry the other half is essentially the same).

Fig. 7 . 2 :

 72 Fig. 7.2: Error curves for the best rational approximations of Table 7.1.

Fig. 7 . 3 :

 73 Fig. 7.3: Result of the weighted version of our barycentric Remez algorithm for the function f(x) = √ x, x ∈ [10 -8 , 1] with w(x) = 1/ √x and a type[START_REF] Cline | Rate of convergence of Lawson's algorithm[END_REF][START_REF] Cline | Rate of convergence of Lawson's algorithm[END_REF] rational approximation. We plot the absolute error curve on the left, while the relative error (right), matching our choice of w, gives an expected equioscillating curve. This is Zolotarev's third problem.

Fig. 7 . 4 :

 74 Fig. 7.4: The error in type (35, 34) best approximation to the sign function on [-10 4 , -1] ∪ [1, 10 4 ], computed via xr(1/x 2 ), where r(x) ≈ √ x as obtained in Figure 7.3. This is Zolotarev's fourth problem.

8. 4 .

 4 Experiments with AAA-Lawson.

Figure 8 .

 8 1 compares AAA and AAA-Lawson (run for ten Lawson steps) for type[START_REF] Berrut | Barycentric Lagrange interpolation[END_REF][START_REF] Berrut | Barycentric Lagrange interpolation[END_REF] and[START_REF] Cooper | Rational Approximation of Discrete Data with Asymptomatic Behaviour[END_REF][START_REF] Cooper | Rational Approximation of Discrete Data with Asymptomatic Behaviour[END_REF] 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.

Fig. 8 . 1 :

 81 Fig. 8.1: Error of rational approximants to f (x) = |x| by the AAA and AAA-Lawson algorithms. The black dots are the support points. They are also interpolation points for AAA, but not for AAA-Lawson.

Fig. 8 . 2 :

 82 Fig. 8.2: Convergence of AAA-Lawson alone and AAA-Lawson followed by Remez, for f (x) = |x|, m = n = 10. The error is measured by r * -r k ∞, where r k is the kth iterate. AAA-Lawson converges linearly, whereas Remez converges quadratically.

9. 1 .

 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 , . . . , t n , which do not belong to X.) Givenr k = N k /D k ∈ R n,n(X), choose the partial fraction decompositions N and D of (2.1) that minimize the expressionmax x∈X |f (x)D(x) -N (x)| -δ k |D(x)| |D k (x)| , t (x)D(x)) = sign(ω t (y)D(y)), ∀x, y ∈ X, x = y, δ k = max x∈X |f (x) -r k (x)|. If r = N/D is not good enough, continue with r k+1 = r. By imposing(9.

Fig. 9 . 1 :

 91 Fig. 9.1: The functions of Table9.1 with error curves for best rational approximations computed by the barycentric DC algorithm.

  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 (1.2). = {z 1 , . . . , z M } is a set of distinct points (sample points) in [a, b]. The number M is usually large, e.g. 10 5 , and in particular much bigger than m and n. The idea is that the solution for the discrete problem (8.1) should converge to the continuous one (1.2) if we discretize the interval densely enough. 8.3. AAA-Lawson. We now propose a rational variant of Lawson's algorithm. (A similar attempt was made in

			Specifically, we
	consider the problem
	minimize r∈Rm,n	f (Z) -r(Z) ∞ ,	(8.1)
	where Z

Table 9 .

 9 1: Best type[START_REF] Cheney | Two new algorithms for rational approximation[END_REF][START_REF] Cheney | Two new algorithms for rational approximation[END_REF] approximations to four functions using the barycentric DC algorithm. X consists of 20000 equispaced points inside [-1, 1].

Table 9

 9 .1 with error curves for best rational approximations computed by the barycentric DC algorithm.which are then solved using MOSEK's[START_REF] Mosek Aps | The MOSEK optimization toolbox for MATLAB manual[END_REF] state-of-the-art LP optimizers. The four examples in Table9.1 and Figure9.1, 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

Chebfun's previous remez command (until version 5.6.0 in December

2016) could only go up to type[START_REF] Berrut | Recent developments in barycentric rational interpolation[END_REF][START_REF] Berrut | Recent developments in barycentric rational interpolation[END_REF].

The prototype code used is available at https://github.com/sfilip/barycentricDC.

SF and LNT were supported by the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement 291068. The views expressed in this article are not those of the ERC or the European Commission, and the European Union is not liable for any use that may be made of the information contained here. YN was supported by Japan Society for the Promotion of Science as an Overseas Research Fellow.

AAA-Lawson proceeds as follows:

1. Use the AAA algorithm to find an approximant (2.2), in particular the support points {t k } for a rational approximation r to f . This step is not tied to a particular norm. 2. Use a variant of Lawson's algorithm to obtain a refined (near-best) rational approximant in the ∞ norm. Below we first review the AAA algorithm, introduced in [START_REF] Nakatsukasa | The AAA algorithm for rational approximation[END_REF], 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 ∈ C M , the AAA algorithm finds a rational approximant of type (n, n) represented as in (2.2) by

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 N (t k )/ D(t k ) = f (t k ) for all k (unless β k = 0). AAA takes only β k as the unknowns, which are found by solving a linearized least-squares problem of the form minimize β 2=1 f D -N Z , where the subscript Z denotes the discrete 2-norm at points Z := Z \ {t 0 , . . . , t n }. For details, see [START_REF] Nakatsukasa | The AAA algorithm for rational approximation[END_REF].

Noninterpolatory AAA. As we discussed in Section 2, the representation N (z)/ D(z) is unsuitable when the goal is to represent r * : it is necessary to use the representation

where F = diag(f ( Z)), and C ,k = 1/(z -t k ) is the Cauchy (basis) matrix as in (4.11), but with rows corresponding to z ∈ {t 0 , . . . , t n } removed. We take the same support points {t k } as in AAA. We solve (8.2) by computing the SVD of the matrix C -F C and finding the right singular vector v = α β ∈ R 2n+2 corresponding to the smallest singular value. As in Section 4.4, the case m = n also uses the projection matrices P m , P n . 8.2. Lawson's algorithm. Lawson's algorithm [START_REF] Lawson | Contributions to the Theory of Linear Least Maximum Approximations[END_REF] 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 , . . . , z M } in a linear subspace span(g i ) n i=0 . With an initial set of weights {w j M j=1 such that w j ≥ 0 and M j=1 w j = 1, one solves (using a standard solver) the weighted least-squares problem minimize c0,...,cn

and computes the residual r j = f (Z j ) -n i=0 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 / M i=1 w i . Iterating this process is known to converge linearly to the best polynomial approximant (with nontrivial convergence analysis [START_REF] Cline | Rate of convergence of Lawson's algorithm[END_REF]), and an acceleration technique is presented in [START_REF] Ellacott | Linear Chebyshev approximation in the complex plane using Lawson's algorithm[END_REF]. 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 9.1 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 [START_REF] Cooper | Rational Approximation of Discrete Data with Asymptomatic Behaviour[END_REF][START_REF] Cooper | Rational Approximation of Discrete Data with Asymptomatic Behaviour[END_REF]. Step 1 consists of picking the initial reference set. This is done by applying in succession (if needed) the strategies discussed in Sections 5.1, 5.2 and 5.3. Next up in Step 2 is computing the current approximant r k and alternation error λ k . We do this by solving a symmetric eigenvalue problem (4.13), (4.22) or (4.24), depending on m = n, m > n or m < n. We then pick, if possible, the eigenpair leading to a rational approximant with no poles in [a, b] (see discussion around equation (4.14)). The next reference set is determined in Step 3 as explained in Section 6. If convergence is successful, the routine outputs a numerical approximant of r * . 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 10.1 summarizes our algorithm in a flowchart.


Keywords: single molecule toroics, lanthanide, dysprosium, cubane, ab initio calculations

A Dy 4 Cubane: A New Member in the Single-Molecule Toroics Family

Guglielmo Fernandez Garcia, [a,b] Djamilla Guettas, [c] Vincent Montigaud, [a] Paolo Larini, [d] Roberta Sessoli, [b] Federico Totti, [b] Olivier Cador,* [a] Guillaume Pilet,* [c] and Boris Le Guennic* [a] Abstract: Molecular materials that possess a toroidal moment associated to a non-magnetic ground state are known as Single-Molecule Toroics (SMTs) and are usually planar molecules. Herein, we report a Dy4 cubane, namely [Dy4(Bppd)4(µ3-OH)4(Pa)4(H2O)4]•0.333H2O (where BppdH=1,3-Bis(pyridin-4yl)propane-1,3-dione and PaH=2-Picolinic acid) for which magnetometry measurements and state-of-art ab initio calculations highlight SMT behavior in a tridimensional structure (3D-SMT). The in-depth theoretical analysis on the resulting low-lying energy states, along with their variation in function of the magnetic exchange pathways, allows to shed further light on the description of Single-Molecule Toroics and identify the coupling scheme that better reproduce the observed data.

Single Molecule Magnet (SMM) behavior was discovered more than twenty years ago in polynuclear 3d complexes [1] and ten years later, in mononuclear 4f complexes, [2] the latter being nowadays the most popular. [3] These objects fascinate chemists and physicists because they might be components of spintronic devices. [4,5] They are conceptually based on axial magnetic anisotropy to create a potential barrier separating the two orientations of the magnetic moment. Magnetic orbitals of lanthanide are largely subject to spin-orbit coupling but slightly perturbed by the coordination because of the small radial extension of 4f orbitals. Therefore, the amplitude of the magnetic anisotropy is much stronger for lanthanide ions and this turns to be an advantage to produce mononuclear SMMs. [START_REF] Tang | Lanthanide Single-Molecule Magnets[END_REF] However, the fast relaxation observed in these systems at zero magnetic field is a long-standing problem. [START_REF] Liddle | [END_REF] To overcome it, the suppression of hyperfine couplings through isotopic substitution and the dilution in diamagnetic matrices are efficient alternatives. [8,9] Another strategy resides in having inside the molecule more than one paramagnetic center coupled through exchange interactions. The magnetic interactions between lanthanide magnetic moments in homometallic polynuclear complexes are, in general, very small (of the order of few Kelvin) [10,[START_REF] Tang | Polynuclear Lanthanide Single Molecule Magnets[END_REF] and not always able to quench the tunneling relaxation process between the two low-lying exchange multiplets, even if radical ligands strongly modify this picture. [START_REF] Rinehart | [END_REF] Indeed parallel and axial principal magnetic axes of the individual centers are needed. [13] On the other hand, polynuclear lanthanide systems can show very peculiar magnetic moments configurations at low temperature which are due to a subtle balance between anisotropy directions and exchange interactions. Probably, Single-Molecule Toroics (SMTs) are the most aesthetic examples of atypical magnetic structures. [10] SMTs have attracted the attention of the scientific community, since they offer ways to design multiferroic materials at the nanoscale. [10] This class of compounds is characterized by a "vortex-like" spatial distribution (e.g. a torus disposition) of the single ion easy-axes that leads to a zero total magnetic momentum, but a non-vanishing toroidal magnetic moment. Several examples of Dy III -based polynuclear complexes (Dy3 triangles, [10,[START_REF] Tang | Polynuclear Lanthanide Single Molecule Magnets[END_REF]14,15] Dy4 coplanar, [16,17] coupled Dy3 triangles, [18][19][20] Dy6 wheel [21] or 3d-4f metallocycles [22] ) have been reported so far with a non-magnetic Kramers doublet ground state resulting from such magnetic vortex on wheel-shaped topologies. Molecular symmetry, magnetic anisotropy direction and magnetic interactions are the key ingredients to produce SMTs. [10] In such a framework, Dy III ions are ideal candidates since it easily leads to an easy magnetic anisotropy axis.

Even if the clusters that have shown toroidal-like spin structures so far are arranged in a planar configuration, in principle, toroidal-like spin structures, and so SMTs, might emerge as well from non-planar arrangements of the paramagnetic centers. Indeed, the requirement for a toroidal ground state is only a spatial distribution of weakly coupled magnetic moments that results in a quenching of the total momentum of the ground state, but also to a residual toroidal moment. In this regard, we investigated a novel Dy4 cubane complex that shows a non-trivial toroidal magnetic moment despite its tetrahedral core. In this work, we report the structure, the magnetic characterizations along with the ab initio calculations and the spin Hamiltonian model which support the existence of uncommon toroidal magnetic moment in a nonplanar complex (3D-SMT). S1. Dy4 crystallizes in the centrosymmetric P42/n tetragonal space group. The asymmetric unit is composed by one Dy III metal center, one Bppd -and one Pa -both disordered on two positions (refined to 50:50), one µ3-OH -bridge, and one terminal coordinated water molecule (Figure S1). One disordered noncoordinated water molecule completes the unit-cell. The structure of the complex is constructed as a tetranuclear architecture where one half of the vertexes of the cubane core are occupied by Dy III centers and the other half by µ3-OH - molecules (Figure 1). Each Dy III is surrounded by eight oxygen atoms (SHAPE analysis [23] : square antiprism; Table S2), two from one Bppd -ligand, two from two Pa -ligands, three from three µ3-OH -and one from a coordinated water molecule. Pa - ligands connect two lanthanide ions in a µ1,2 mode on four of the six {Dy4O4} faces. Ln-O(H2O) bond lengths are greater (2.457 Å vs. 2.369 Å) than the other Ln-O bond lengths (Table S3). All the Dy-O(µ3-OH)-Dy bond angles are very similar ranging from 106.8° to 107.1°. Likewise, Dy … Dy distances within the cubane architecture are almost the same ranging from 3.83 to 3.85 Å, in line with the deviation from the ideal point group symmetry. The crystal packing is insured by π-π interactions involving aromatic rings of Bppd ligands forming a dense 3D network (Figure S2). The shortest Dy…Dy distance between two different Dy4 (center to center) is equal to 14.0 Å ensuring no magnetic interaction between neighbors.

Room temperature χMT value, with T the temperature in Kelvin and χM the molar magnetic susceptibility, is just slightly smaller than expected for the 6 H15/2 multiplet ground state with g=4/3 (52.1 vs. 56.6 cm 3 K mol -1 ) and collapses down to 2.5 cm 3 K mol -1 on cooling down to 500 mK (Figure S3). Meanwhile, χM passes through a rounded maximum at 2 K which clearly indicates that the ground state is non-magnetic (Figure 2). At the lowest investigated temperature, a weak increase of χM is observed which is attributed to a low temperature paramagnetic Curie tail, due to the impurities present in the crystal. The nonmagnetic ground state is confirmed by S-shaped first magnetization curve at 500 mK (Figure 2) which saturates at 19.2 Nβ (Figure S3). This value is very close to what is expected for four effective spins Seff=½ with gx=gy=0 and gz=20. This pattern corresponds to the stabilization of the two MJ=±15/2 states. All these characteristics indicate that significant interactions operate in the cubane leading to a non-magnetic ground state with magnetic excited states with different contribution within few wavenumbers. Ab initio SA-CASSF/SI-SO calculations (see details in SI) have been performed to i) determine the mutual orientation of the principal magnetic axes of the Dy III ions related to their arrangement in the cubane core and ii) elucidate the magnetic interactions which stabilize the non-magnetic ground state. The magnetic properties of each Dy III ion were calculated by a diamagnetic substitution approach, i.e. substituting in turn the other three Dy III ions with Y III ions. The Kramers ground doublet is characterized by an axial g-tensor with a gzz ~ 19.6 in the effective spin ½ model (Table S4) with the first excited state at ~100 cm -1 (Table S5). |±15/2> MJ component (> 96%) prevails in the ground state of each Dy III ion (Table S6). The easy magnetic axes are symmetry related due to the almost perfect tetrahedral structure (Table S7). They make an angle of 30.5° with c axis (Figure S4). From the c axis perspective, it is possible to evidence the circular pattern described by the easy axes. The presence of magnetic dipole-dipole interactions (Jdip, see SI) and exchange coupling, Jex, have been considered to account for the magnetic inter-ions interactions. The latter was included in the framework of the Lines model. [24] Furthermore, it was impossible to properly reproduce the experiment with enough accuracy either with Jex=0 (i.e. only dipolar interactions) or with a unique Jex (Figure S5). To clarify the magnetic exchange coupling scheme, the following spin Hamiltonian was used (the subscript « ex » is omitted):

( ) ( ) = - ⋅ + ⋅ + ⋅ + ⋅ - ⋅ + ⋅ 1 1 2 2 3 3 4 4 1 2 1 3 2 4
ˆˆˆˆˆˆˆŝ pin

H J S S S S S S S S J S S S S

with J1 and J2 referring to the exchange couplings when µ3-OH and carboxylate or only µ3-OH bridging groups are involved, respectively, Si's are the effective spin operators at Dy III sites (Figure 3). The dipolar interaction is calculated routinely in the point-dipole approximation. All the calculations were performed with the POLY_ANISO routine [25,26] and homemade fitting code, with the inclusion of the g tensors for each Dy center calculated separately (Table S4). With this model both the crystallographic

Accepted Manuscript

disorder on each metal site and the non-collinearity of the singleion g tensors are considered. In this framework, the resulting low-lying 16 states are arranged in four groups of degenerate states (Figure 3). For the sake of simplicity, the spatial disposition of the single magnetic moments (left-hand side in Figure 3) for each group of states can be labeled using the collinear "up" (↑) and "down" (↓) notation. Remarkably, the total magnetic moment of the ground state |↑↓↑↓> is close to 0 (effective g=0.46), while excited states have a non-zero total magnetic moment. Looking at the spins configuration of the ground state along the c axis (Figure 3), the toroidal arrangement of the magnetic moments emerges. To determine under which condition state |↑↓↑↓> becomes the ground state, we mapped the energies of the low-lying states for J1 and J2 ranging from -0.5 to +0.5 cm -1 . The resulting surfaces, called Exchange Energy Surfaces (EES) in analogy with the widely used Potential Energy Surfaces (PES), are plotted on Figure 4. The ground state is always non-magnetic for J1<0 and J2>0. For J1<0 and J2<0 (only antiferromagnetic couplings) there is a competition between states |↑↓↑↓> and |↑↑↓↓>. In this subset of values, states |↑↓↑↓> or |↑↑↓↓> are stabilized when J1 or J2 dominates, respectively. However, on the basis of the well-known efficiency of the super-exchange pathway given by the carboxylate ligands, J1 is therefore expected to be dominating in Dy4. [27] Such results point out the necessity to verify the presence of equivalent set of exchange parameters able to reproduce in the most accurate way as possible the experiment. Even if only seldom performed, [28] such a detailed analysis of the EES's is mandatory whenever in presence of potentially high correlated set of parameters (Figure 4). We finally identified a set of J values (J1=-0.27 cm -1 and J2=-0.09 cm -1 ) that match reasonably well the maximum on the χM vs. T curve as well as the inflexion point on the M vs. H curve (Figures 2 andS5). The slight deviation at low T may be ascribed to small amount of paramagnetic impurities and structural disorder on the ligands that are not taken into account in the calculations. The selected set of exchange parameters is not, of course, the unique set which can reproduce the experimental data, but surely it is the most sound on the basis of the ab initio results. Indeed, when one of the most important parameters is fixed (the energy gap between the ground and first excited states), all other used sets of parameters fail to reproduce the experiment (Figure S5). The energy diagram and the calculated total g-values in the effective spin 1/2 framework are provided on Figure 3: The ground state is non-magnetic, even if its global g-value is not strictly zero due to the disorder in the crystal structure that lower the symmetry. To confirm this, we ran a calculation on a symmetrized spin Hamiltonian (Td symmetry) resulting in a strict diamagnetic ground state (see SI, Table S8). Back to the main model, the first excited states |↑↑↓↓> (g=28.33) lie 2.14 cm -1 above the ground state, while the third and fourth excited states are at 3.46 cm -1 (g=39.05) and 9.51 cm -1 (g=67.63), respectively. One must notice that for all excited states, global g tensors are not collinear (Figure S6). In order to confirm that state |↑↓↑↓> is not only non-magnetic but possesses a non-zero toroidal moment, we calculated the total spin

! S tot = ! S i i ∑ and the toroidal moment ! τ ∝ ! r i × ! S i i ∑ , respect
to the center of mass of the metal core (as a reference the same quantities were calculated for the archetypal Dy3 triangle, [15] see Table S9). State |↑↓↑↓> presents a quenching of the total spin but a non-zero τ (net toroidal moment), while state |↑↑↓↓> has both non-zero τ and Stot (mixed-moment SMT). Dy4 shows frequency dependent ac susceptibility in the 2-11 K temperature in the absence of an external dc field (Figure S7), the signature of the slowing down of the magnetic moment. According to the calculated energy diagram of Dy4 (vide supra), Accepted Manuscript even at 2 K not only the non-magnetic ground-state is populated but the other states are also involved. As a result, the relaxation of the magnetic moment implies several states with different processes. Nevertheless, the experimental data can be analyzed in the frame of the extended Debye model which feature one main relaxation with a relatively narrow distribution (α<0.3) (see SI, Figure S8, Table S10). While a thermally dependent process dominates at high temperature, the relaxation time tends to become thermally independent below 4 K. The energy diagram which results from the coupling between local magnetic moments impacts, in a complex manner, the infield relaxation dynamics. At 2 K, the in-field behavior of χM' in the 1-1000 Hz frequency window is plotted on Figure 6. It reflects the in-field evolution of the relaxation time which is not monotonic. The slowing down at low field is commonly observed in mononuclear SMMs. At higher field than 1 kOe the relaxation time decreases as the surface plot moves to higher frequencies (Figure S9). The black line at high frequencies (1000 Hz) shows two maxima at 2.8 and 4 kOe. These fields correspond to the crossings between the first three states (|↑↓↑↓>, |↑↑↓↓> and |↑↓↑↑>) in agreement with Zeeman energy diagrams calculated with the magnetic applied along the principal axis of state |↑↑↑↑> (4.3 kOe, Figure S10) and along one of the principal axes of state |↑↑↓↓> (3.2 kOe, Figure S11). This allows multiple relaxation pathways to the global magnetic moment. Significant toroidal magnetic moment in a Dy4 cubane complex associated to a non-magnetic ground state has been characterized in detail by magnetometry and state-of-art ab initio computational methods. A deep analysis of the low-lying energy states resulting from the specific coupling scheme allows to reliably identify the best set of parameters to properly reproduce the magnetic data. The computed energy ladder also suggests that in such a class of systems, the magnetic nature of the ground state could be easily modulated by altering the J1/J2 ratio and, therefore, by slight changes in the ligand nature or by external perturbations (i.e. pressure). As far as we know, this is one of the first cases in which a non-magnetic toroidal ground state is found in a 3D topology (3D-SMT). This strongly differs from ferrotoroidic systems composed by two connected frustrated rings (interacting 2D-SMT). [29] These findings pave the way to the design of original three-dimensional toroics as new objects in molecular spintronic.

Figure 1 .

 1 Figure 1. Representation of the crystal structure of Dy4. BppdH are in yellow rods, PaH in green rods, while Dy, O and N are in gold, red and blue balls, respectively. Hydrogen atoms have been omitted for clarity.

  Dr. G. Fernandez Garcia,V. Montigaud, Prof. Dr. O. Cador, Dr. B. Le Guennic

  Bppd)4(µ3-OH)4(Pa)4(H2O)4]•0.333H2O (where BppdH=1,3-Bis(pyridin-4yl)propane-1,3-dione and PaH=2-Picolinic acid), hereafter abbreviated Dy4, is represented on Figure 1. The synthesis and the details of the crystal structure are given in SI and Table

Figure 2 .

 2 Figure 2. Thermal variation of the magnetic susceptibility in black dots with the calculated one (red line) according to the model developed in main text; (inset) Field variation of the magnetization at 500 mK and the calculated curve (red line).

Figure 3 .

 3 Figure 3. Coupling scheme within Dy4 (top) with the calculated energy state degeneracies (bottom) along two directions (x axis view for (a); c axis view for (b)) with gzz-values in the effective spin 1/2 framework and calculated toroidal moment τ (in β.Å).

Figure 4 .

 4 Figure 4. EES for states |↑↓↑↓>, |↑↑↓↓> and |↑↑↑↑>. State |↑↓↑↑> has been omitted for clarity. The vertical black segment corresponds to the best agreement with the experiment (J1=-0.27 cm -1 and J2=-0.09 cm -1 ).

Figure 6 .

 6 Figure 6. Variation of the in-phase component (χM') at 2 K with the frequency of the oscillating field (ν) and the external dc field (H). The thick black line corresponds to χM' values at 1000 Hz and the two dots to the observed maxima at 2.8 and 4 kOe.

Acknowledgements D. G. is grateful for the PhD fellowship from University Claude Bernard Lyon 1. Additional support was provided by the Région Rhône-Alpes (CMIRA2015). G. F. G. gratefully acknowledges the European Commission through the ERC-AdG MolNano-Mas (project no. 267746) and the ANR (ANR-13-BS07-0022-01) for financial support. V. M. is thankful to ERC (project no. 725184) for fundings. B. L. G. thanks the French GENCI/IDRIS-CINES center for high-performance computing resources. COST Action Molspin (CA 15128) is also acknowledged.