Chebyshev model arithmetic for factorable functions

This article presents an arithmetic for the computation of Chebyshev models for factorable functions and an analysis of their convergence properties. Similar to Taylor models, Chebyshev models consist of a pair of a multivariate polynomial approximating the factorable function and an interval remainder term bounding the actual gap with this polynomial approximant. Propagation rules and local convergence bounds are established for the addition, multiplication and composition operations with Chebyshev models. The global convergence of this arithmetic as the polynomial expansion order increases is also discussed. A generic implementation of Chebyshev model arithmetic is available in the library MC++. It is shown through several numerical case studies that Chebyshev models provide tighter bounds than their Taylor model counterparts, but this comes at the price of extra computational burden.


Introduction
Many complete search methods for problems in global optimization and constraint satisfaction hinge on the ability to compute enclosures for the range of nonconvex functions as well as sets defined by multiple equalities and/or inequalities [44]. A variety of parameterizations can be used to describe or approximate compact sets, including convex sets such as intervals, ellipsoids and polytopes as well as nonconvex sets for instance in the form of the image of a multivariate polynomial. Chachuat et al. [13] have argued that the foregoing parameter-izations and many others can be cast as affine set-parameterizations with either convex or nonconvex basis sets. In particular, a number of arithmetics have been developed over the years that enable propagation of certain affine-set parameterizations through the atom operations of factorable functions in order to determine rigorous enclosures of their image sets. This includes interval arithmetic [40,41], ellipsoidal calculus [28,61], polyhedral relaxations [54,55], and Taylor model arithmetic [11,36], to name but a few. Other established techniques to relax sets or programs defined by nonconvex factorable functions involve constructing a pair of convex/concave bounds for each function, such as McCormick relaxations [38,39] and αBB/γ BB relaxations [1,3,23]. In the case of semialgebraic problems, hierarchies of linear matrix inequality (LMI) relaxations can also be constructed based on a semidefinite programming formulation of the sum of squares decomposition for multivariate polynomials [29,46].
The focus in this paper is on polynomial models of multivariable functions, namely enclosures in the form of a polynomial approximant of given order and an interval remainder bound on the approximation error. Taylor models (or Taylor forms) are a special case of polynomial models, whereby the polynomial approximation matches the multivariate Taylor expansion of the function at a given point in the variable domain [43]. The idea of such interval polynomial enclosures dates back to Moore [40], and methods for computing the remainder bounds concurrently with the polynomial were later developed in the early 1980s [17,51] and popularized from the mid 1990s by Berz et al. [7,35,36]. Similar in essence to interval arithmetic, a well-developed Taylor model arithmetic is available [8,34], encompassing rules for binary sums, binary products and outer-compositions with a library of univariate functions such as exp(·), log(·) and √ ·. For sufficiently smooth functions, the diameter of the remainder interval constructed with (qth-order) Taylor model arithmetic can be proved to be a high-order power (q + 1) of the diameter of the variables domain [11]. This gives Taylor models a clear advantage over traditional interval extensions or centered forms for sufficiently narrow domains, but conversely it may result in a large overestimation or may even be poorer than naive interval evaluation over wider domains. Nevertheless, this approach has proved successful in computing tight enclosures for the solutions of differential equations and implicit algebraic equations [24,33,42,49,50,53,60], and it has enabled complete search for a range of global optimization or constraint satisfaction problems that could not be tackled using interval techniques alone (see, e.g., [4,9,31,32,47,52]). Such higher-order inclusion techniques are indeed appealing in complete search applications based on branching or subdivision, where they can mitigate the clustering effect [15,62].
A rather natural idea involves replacing the Taylor series approximation with a better one, ideally the minimax polynomial approximation, or near-minimax approximations such as Chebyshev series. The development of affine arithmetics [19] yielding first-order approximations can be seen as a precursor of this approach. Regarding univariate functions, improved variants based on the Chebyshev polynomials, Bernstein polynomials and others were widely used during the early 1980s under the name of ultra-arithmetic [26]. A comparison of univariate Taylor forms and Chebyshev forms in [27] showed that expansions in Chebyshev series may be orders of magnitude more accurate than expansions in Taylor series. The use of Chebyshev series expansion as an arithmetic is also the philosophy behind the package Chebfun by Trefethen et al. [6,57], which has recently been extended to functions in two variables as well [56]. It should be noted that Chebfun relies on computations with functions to 15-digit accuracy, thus removing the need for propagating a remainder term, but the results are not validated per se. In contrast, Brisebarre and Joldeş [12,25] developed an arithmetic for univariate functions based on Chebyshev interpolation polynomials and a remainder term. More recently, Dzetkulič [16] extended the previous work to enable enclo-sure of multivariate functions using truncated Chebyshev series and an interval remainder term, namely Chebyshev models. This latter work makes use of Chebyshev models as a means of propagating bounds for differential equations, but does not study the convergence properties of this arithmetic.
Building upon these existing contributions, this article presents a Chebyshev model arithmetic in order to compute enclosures for an inclusive class of factorable functions (Sect. 3). Besides formalizing the rules of addition, multiplication and composition as well as the polynomial range bounders, we describe an improvement of the remainder term in composition operations with certain univariate functions. We also present a convergence analysis of Chebyshev models, which investigates the propagation of local convergence properties through the addition, multiplication and composition operations (Sect. 4), and we discuss the global convergence of this arithmetic as the polynomial expansion order increases. Finally, we present numerical case studies illustrating the convergence properties of Chebyshev models and comparing their tightness and computational performance with Taylor models (Sect. 5). The underlying implementation of Chebyshev model arithmetic is part of the library MC++, which is made freely available at: http://omega-icl.bitbucket.org/mcpp/.

Notation and background
The set of n-dimensional interval vectors is denoted by IR n . The midpoint and radius of an interval vector X := [x L , x U ] ∈ IR n are defined as mid X := 1 2 (x L + x U ) ∈ R n and rad X : The total variation of a (possibly discontinuous) function g : R → R on the interval Y := [y L , y U ] ∈ IR is the quantity where ξ 0 , . . . , ξ N is any finite partition of Y , and the supremum norm of g on Y is given by Moreover, by a slight abuse of the notation, we use V [g] and g when Y = [−1, 1].
Chebyshev polynomials of the first kind, denoted by T k , are defined by the three-term recurrence relation An explicit expression of these polynomials on [−1, 1] is T k (x) = cos(k arccos(x)).
The Chebyshev polynomials are extremal polynomials with respect to the property that, for each k > 0, T k (x) is the polynomial with the largest possible leading coefficient subject to the condition that its range is in [−1, 1] for x ∈ [−1, 1].
Any polynomial function P q of degree q can be written in terms of the Chebyshev polynomials as [37] An efficient way of evaluating such a linear combination of Chebyshev polynomials is using the Clenshaw algorithm [14], where the b k 's are obtained using the (reverse) recurrence formula Moreover, the product between two Chebyshev polynomials expands to The Chebyshev polynomials form a sequence of orthogonal polynomials with respect to the weight 1 Any Lipschitz continuous function f on [−1, 1] has a unique representation as an absolutely and uniformly convergent series [58], where the notation indicates that the first term is halved. In general, the Chebyshev coefficients a k can be approximated using a numerical quadrature, possibly in combination with the fast Fourier cosine transform to improve the computation speed [21].
For practical purposes, a qth-order polynomial approximant P q f for f can be obtained from the partial sums of (3) as The smoother the function f on [−1, 1], the faster the convergence of its approximants P -If f and its derivatives through f (q) are continuously differentiable on [−1, 1], then [18] f An alternative approach to obtaining a qth-order polynomial approximant is by constructing the polynomial interpolant P q f at Chebyshev points, The polynomial P q f yields a near-minimax approximation, and counterparts to the error estimates in (5) and (6) are, respectively, [18,58] While the truncated Chebyshev series expansion P q f is generally a better approximation than the Chebyshev interpolation polynomial P q f , the latter comes with the inherent advantage that the coefficients a k in (7) can be obtained from a few function evaluations only.
The (over-)approximation of multivariate functions using Chebyshev polynomials is at the heart of the present paper. Like in the univariate case, any qth-order polynomial function in n variables can be written in terms of the Chebyshev polynomials as where κ ∈ N n is a multi-index; the order of κ is given by |κ| := n i=1 κ i ; and T κ (x) is a shorthand notation for the expression n i=1 T κ i (x i ) for any x ∈ R n .

Chebyshev model arithmetic
Similar to Taylor models [11,43], Chebyshev models are estimators made by the sum of a multivariate polynomial in Chebyshev form and a remainder interval providing a bound on the approximation error. A formal definition follows:

Definition 1
Let the function f : Z → R be defined on Z ⊆ R n . For every interval vector Y ⊂ Z , let the multivariate polynomial P q f,Y : [−1, 1] n → R with P q f,Y (ξ ) := |κ|≤q a κ T κ (ξ ) and the symmetric interval R q f,Y ∈ IR be such that where • denotes the entrywise (or Hadamard) product. We call the pair (P Notice that, in the sense of Definition 1, the polynomial part P q f,Y of a Chebyshev model may be different from the actual truncated Chebyshev expansion of f . Notice also the need for scaling the variables in [−1, 1] for consistency with the definition and properties of Chebyshev polynomials of the first kind. The definition of schemes of Chebyshev models is in analogy with the scheme of estimators introduced by Bompadre et al. [10,11] to support convergence analysis. Similar in essence to interval (or Taylor model) arithmetic, our focus in the following subsections is on the development of an arithmetic that enables the construction of Chebyshev models for factorable functions, namely functions which can be represented by a finite number of binary sums, binary products and outer compositions with a univariate function. As well as the specification of these three "elementary" operations, developing such an arithmetic also calls for an initialization procedure for variables and constants, and the ability to bound the range of Chebyshev models.

Initialization
Given a set Z ⊆ R n , each variable z i ∈ Z i with i ∈ {1, . . . , n} can be represented by a scheme of qth-order Chebyshev models (P Likewise, any real constant r ∈ R is trivially represented by the scheme of Chebyshev models

Binary sum operations
Consider two functions f 1 , f 2 : Z ⊆ R n → R, and suppose that (P are corresponding schemes of Chebyshev models. Then, a scheme of Chebyshev models (P is simply obtained by addition/subtraction of the polynomial parts and of the remainder parts,

Binary product operations
Consider again two functions f 1 , f 2 : Z ⊆ R n → R, and suppose that (P In order to arrive at a valid scheme of qth-order Chebyshev models (P for the product function f 1 f 2 in Z , we start by expressing the product between the respective polynomial parts, A complication with polynomial multiplication in Chebyshev basis is that, according to the property (2), each product term T λ (ξ )T μ (ξ ) generates 2 N (λ,μ) terms with See also Sect. 5 for further discussions and related implementation considerations.
Introducing the map P q : N n → N n × N n given by it follows from (2) that Observe that the resulting product polynomial is of order 2q and that the index sets P q (κ) are non-empty for |κ| ≤ 2q. A simple way of defining the polynomial part P q f 1 f 2 ,Y in the resulting Chebyshev model is therefore by retaining all the terms of order no larger than q as Then, the corresponding remainder term R q f 1 f 2 ,Y can be obtained by assembling four terms, namely the multiplication of the two remainders, the multiplication of each remainder with a bound on the other polynomial range, and a bound on all the terms of order larger than q from the product of the two polynomials, The latter hinges on the ability to compute bounds on the polynomial parts of the Chebyshev model operands, denoted by Such range bounders for Chebyshev models are discussed later on in Sect. 3.5.
An important aspect about polynomial multiplication in Chebyshev basis is that the computed qth-order polynomial P q f 1 f 2 ,Y in (17) will not correspond to the qth-order Chebyshev expansion of f 1 f 2 in general, even when P q f 1 ,Y and P q f 2 ,Y are themselves qth-order Chebyshev expansion of f 1 and f 2 , respectively. The reason behind this discrepancy is that the Chebyshev coefficients in the qth-order expansion of f 1 f 2 depend on terms of order greater than q in the Chebyshev expansions of f 1 and f 2 , which are missing here because of the truncation. Dzetkulič [16] argue that the resulting Chebyshev models nevertheless provide a tighter approximation than by multiplying P q f 1 ,Y and P q f 2 ,Y in monomial basis.

Univariate outer-composition operations
Consider the functions f : We will discuss the computation of such a scheme at the end of this subsection.
In order to construct a valid scheme of qth-order Chebyshev models (P with: Observe that this recurrence consists of binary sum and product operations only, for which rules have already been established in Sects. 3.2 and 3.3. Observe also that when rad Finally, the polynomial and remainder parts of the desired Chebyshev models (P In order to use this approach, one needs the ability to construct Chebyshev models (P q F,W , R q F,W ) corresponding to a given library of univariate functions, such as , |x|, etc. -A systematic way of doing this is by computing a truncated Chebyshev expansion of F of order q as Then, if F is q + 1-times continuously differentiable on W , the truncation error can be bounded as Otherwise, if F and its derivatives through F (s−1) are absolutely continuous on W and F (s) has finite total variation V [F (s) ] W < ∞ with s < q, a bound on the truncation error is given by -Alternatively, a Chebyshev interpolating polynomial of order q may be computed as The error bound R q F,W in (26) remains valid for any q + 1-times continuously differentiable function F on W . Otherwise, the following bound may be used instead of (27) if F and its derivatives through F (s−1) are absolutely continuous on W and F (s) has finite In comparison with their Taylor expansion counterparts for sufficiently smooth functions, we note that the width of the remainder is smaller by a factor of 2 q . Moreover, the Chebyshev approach remains applicable for continuous, yet nonsmooth, functions. For instance, since the total variation of the derivative of |x| on the interval [−ρ, ρ] is finite, equal to 2, a bound for the interpolation error based on (29) In practice, the bounds R q F,W on the truncation error may be tightened further if F meets certain regularity and monotonicity conditions. Under these extra conditions, it can be established that the maximum absolute error between F and its polynomial approximant P q F,W -as given by (25) or (28) for any q ≥ 0, with P q F and P q F given by (4) and (7), respectively.
Proof See "Appendix 1". This method of calculating the remainder is particularly useful as functions satisfying the conditions of Lemma 1 include the exponential, logarithm, inverse and square-root functions on any finite interval in their domains of definition.

Range bounding
Such bounders are needed for computing binary product and univariate operations in Chebyshev model arithmetic-see Sects. 3.3 and 3.4 above, thus directly affecting the tightness of the computed Chebyshev models. Unfortunately, exact range bounding for multivariate polynomials is NPhard, so one has to resort to some sort of over-approximation in practice.
A number of methods exist for bounding the range of multivariate polynomials in monomial form [20,30,36,43], which can be readily adapted to their Chebyshev counterparts. Apart from the naive bounding of each term separately, which may produce weak bounds, this subsection presents adaptations of: (i) the methods in [33] that involves exact bounding of the first-and diagonal second-order terms; and (ii) the method of Bernstein.

Exact bounding of first-and diagonal second-order terms
We start by rewriting the polynomial part in the form where κ 1, j and κ 2, j are the multi-indices given by κ . . , n}, using the Kronecker δ notation; and H q f,Y is the multivariate polynomial containing the same terms as P q f,Y apart from those indexed by κ 1 j and κ 2 j . Then, a similar rearrangement to the one proposed in [33] can be obtained by using the property that This way, an interval enclosure is obtained as (30). This bounder may only be considered when |a q Y,κ 2, j | ≥ , where is a small positive number; otherwise, separate bounding of the linear and quadratic terms for the corresponding component j as in (30) may be used.

Bounding method of Bernstein
The Bernstein algorithm is a well-established tool for computing bounds on the range of multivariate polynomials over a box; see, e.g., [20,30] and references therein. The procedure involves rewriting the multivariate polynomial from Chebyshev form into Bernstein form: with r ≥ q, and where B r i , with i ≤ r , denotes the ith Bernstein polynomial of order r on [0, 1], given by In particular, the Bernstein coefficient b q Y,κ in (31) can be expressed in terms of the Chebyshev coefficients a q Y,κ as follows [48] ∀κ ∈ {0, . . . , This transformation turns out to be particularly well-conditioned from a numerical stability standpoint [48]. At this point, an interval enclosure of P q f,Y on [−1, 1] n is simply obtained as The main advantage of this approach is that the foregoing enclosure converges to the actual polynomial range in the Hausdorff sense as the order r increases [20]. Moreover, it is possible to bound the actual under-or over-estimation for a given order r ≥ q.

Convergence analysis
This section presents an analysis of the local convergence rate of Chebyshev models constructed from the application of the arithmetic operations described in Sect. 3. The global convergence of this arithmetic as the polynomial expansion order q increases is also discussed at the end of the section.

Definition 2
Let the function f : Z → R be defined on Z ⊆ R n . The scheme and for all interval vectors Y ⊂ Z with sufficiently small rad Y .
Notice the extra conditions (32) imposed on the coefficients of the polynomial approximant in the previous definition. While these conditions are not strictly necessary to establish local convergence, the propagation of (local) convergence order for the polynomial coefficients through binary sum, binary product and univariate outer-composition operations provides insight into the propagation of the convergence order for the remainder term in the following analysis.
As far as initialization is concerned, it is clear that a scheme of qth-order Chebyshev . . , n}, as given by (11) and (12), has local convergence order 1 if q = 0, and infinite order if q ≥ 1. In both cases, the schemes of qth-order Chebyshev models for the variables thus have convergence order no less than q + 1. Regarding the range bounders introduced in Sect. 3.5, the convergence condition (32) Next, we investigate the local convergence order of schemes of Chebyshev models as propagated through factorable functions, with the corresponding binary sum, binary product and univariate outer-composition operations derived in Sects. 3.2-3.4. With regards to binary sum operations given by (13)-(14), the sum of two schemes of qth-order Chebyshev models trivially preserves the least convergence order of the operands, be they smaller than, equal to, or larger than q + 1. In particular, the addition/subtraction of two schemes with local convergence order q + 1 is itself a scheme of order q + 1. The propagation of convergence orders through binary product and univariate outer-composition operations is somewhat more involved, and detailed in Sects. 4.1 and 4.2.

Local convergence rate of binary product operations
Adopting the notation introduced in Sect. 3.3, we assume here that the schemes of qth-order Chebyshev models (P The following lemma is instrumental to prove convergence of the Chebyshev coefficients in the product scheme: Lemma 2 For any q ≥ 0, any κ ∈ N q , and any pair (λ, μ) ∈ P q (κ), we have |λ + μ| ≥ |κ|.
By the construction of the product polynomial in (16) and by Lemma 2, we have for every κ ∈ N n with |κ| ≤ 2q. Therefore, the product rule (17) propagates the least convergence order of the operands through the polynomial coefficients, up to order |κ|.
so it follows from the remainder propagation rule (18) that Overall, the least convergence order of the operands thus propagates through the remainder, up to order q + 1. In particular, the product scheme has local convergence order q + 1 whenever both operands have convergence order q + 1. These results are summarized in the following theorem: be schemes of qth-order Chebyshev models, with local convergence orders r 1 ≥ 1 and r 2 ≥ 1, respectively, for the functions (17) and (18) has local convergence order min{r 1 , r 2 , q + 1}.

Local convergence rate of univariate outer-composition operations
Using the notation introduced in Sect. 3.4, we assume that the outer functions F : X ⊆ R → R is s-times continuously differentiable, and we consider the scheme of qth-order Chebyshev models (P q F,W , R q F,W ) W ⊂X , as constructed from (25), (26) [resp. from (26), (28)] if s ≥ q + 1; or otherwise constructed from (25), (27) [resp. from (28), (29)] if F (s) has a finite total variation on X .

Lemma 3
The scheme (P q F,W , R q F,W ) W ⊂X has local convergence order min{s, q + 1} on X.
To carry out the analysis, we also assume that the scheme of qth-order Chebyshev models (P q f,Y , R q f,Y ) Y ⊂Z for the inner functions f : Z ⊆ R n → R has local convergence order r ≥ 1. Notice that the Clenshaw summations (20)- (22) consists of binary sum and product operations only and is initialized with constants. However, by construction of has local convergence order r − 1 only, the following lemma shows that the original order r of the inner operand nonetheless propagates through the scheme (P Proof See "Appendix 3". A direct consequence of Lemmata 3 and 4 is that the remainder term R q F• f,Y , as defined in (24), has convergence order min{r, s, q + 1}. Therefore, the convergence order of the inner operand may propagate through the remainder, up to order q + 1, in a qth-order scheme of Chebyshev models. Moreover, if the outer operand is s-times continuously differentiable, the convergence order of the composition scheme may not be larger than s. In particular, the composition scheme has convergence order q +1 whenever the inner scheme has convergence order q + 1 and the outer function is q + 1-times continuously differentiable. These results are summarized in the following theorem: be a scheme of qth-order Chebyshev models for the s-times continuously-differentiable function F : X ⊆ R → R, as constructed from (25), (26) [resp. from (26), (28)] if s ≥ q + 1; or otherwise constructed from (25), (27) [resp. from (28), (29)] if F (s) has a finite total variation on X . Let also (P q f,Y , R q f,Y ) Y ⊂Z be a scheme of qth-order Chebyshev models for the function f : Z ⊆ R n → R with local convergence order r ≥ 1. Then, the composition scheme (P (23) and (24) has convergence order min{r, s, q + 1} on Z .

Global convergence of Chebyshev model arithmetic
The previous convergence analysis investigates the rate at which the remainder bound shrinks as the domain Y of the variables is progressively reduced to a point, for a given expansion order q. This section discusses a different type of convergence, namely the property of the remainder bound to shrink to zero upon increasing the expansion order q → ∞, for a fixed variable domain Y .

Definition 3
Let the function f : Z → R be defined on Z ⊆ R n , and let IR n Y ⊂ Z . The As far as initialization is concerned, it is clear that a scheme of qth-order Chebyshev models (P q y i ,Y , R q y i ,Y ) q≥0 for the variable y i ∈ Y i with i ∈ {1, . . . , n}, as given by (11) and (12), is globally convergent on Y . Regarding the range bounders introduced in Sect. 3.5, it is also worth noting that the sequence {rad [P q f,Y ]} q≥0 is bounded when (36) holds. Next, we adopt the notation introduced in Sects. 3.2 and 3.3, and we assume that the schemes (P q≥0 , as given in (14), is trivially globally convergent by the triangle inequality. Concerning the product operation in Sect. 3.3, we have from (16) and, in particular, It then follows from (18) that the product scheme (P These considerations are summarized in the following theorem: q≥0 given by (13) and (14) is globally convergent on Y ; and so is the product scheme (P q≥0 given by (17) and (18) for the range bounding methods in Sect. 3.5.
In order to formally establish global convergence of the Chebyshev model arithmetic presented in Sect. 3, one also needs to establish sufficient conditions on both the inner and outer operands of a composition operation in order for the resulting composition scheme to be globally convergent. The development of such conditions calls for an analysis of the propagation of the global convergence property through the Clenshaw recurrence, for which results are currently lacking. The numerical case studies presented in the following sections suggest that global convergence might be guaranteed under mild regularity conditions nonetheless, and this will be the topic of further research.

Numerical implementation and case studies
This section presents numerical case studies that illustrate the convergence properties of Chebyshev models, and makes comparisons with Taylor models in terms of tightness and computational effort. All the computations that led to these results are performed with the library MC++ (http://omega-icl.bitbucket.org/mcpp/), which features classes for the evaluation of factorable functions in Taylor and Chebyshev arithmetics (along with other arithmetics). Verified interval libraries, such as PROFIL (http://www.ti3.tu-harburg.de/) or FILIB++ (http://www2.math.uni-wuppertal.de/~xsc/), are used to perform the remainder interval computations, but any round-off caused by operations between the polynomial coefficients is not accounted for in the current implementation. We note that approaches to account for errors due to floating-point arithmetic are well-documented for Taylor models [36,43], and the very same approaches may be used for Chebyshev models in order to arrive at a fully verified implementation.
MC++ implements the Chebyshev model arithmetic and range bounding operations described in Sect. 3, based on a dense representation of the multivariate polynomials. If not otherwise noted, Chebyshev models for the intrinsic functions are computed based on the Chebyshev interpolation polynomial. Moreover, the exact remainder bounding approach established in Sect. 3.4 is used for the intrinsic functions exp(x), log(x), √ x, and 1 x , and the default range bounder uses exact bounding of first-and diagonal second-order terms as described in Sect. 3.5. Another improvement involves intersecting the polynomial model bounds with those derived from simple interval analysis at each operation as a means of tightening the remainder bound, avoiding division by zero, etc; this extension is similar in essence to the mixed affine-arithmetic/interval-arithmetic model by Stolfi and Figueiredo [19].
The most critical operation in an efficient implementation is multiplication of Chebyshev models, since it is used extensively in univariate outer-composition operations besides bivariate product operations. The approach used in MC++ involves constructing the index sets P q (κ) given by (15), for all κ ∈ N n with |κ| ≤ q, prior to propagating the Chebyshev models through a DAG of the factorable function. The number of floating-point operations (FLOP) required to compute the product polynomial P (17) is shown on the left plot of Fig. 1 for Chebyshev expansions of order q = 1, . . . , 4 and for functions with up to n = 19 variables. Even for low expansion orders does performing polynomial products in Chebyshev basis result in an increase by 1 or 2 orders of magnitude in the number of FLOP, compared with monomial basis which is shown on the right plot. The better approximation capability of Chebyshev models may thus come at the price of a much higher computational burden compared with Taylor models, at least in a dense implementation. This refinement calls for the development of sparse implementations for both Taylor and Chebyshev model arithmetics in order to overcome this limitation and tackle larger-scale problems, a topic for future research. 18 20 derivative-based remainder exact remainder

Example 1
We consider the univariate function f given by for the variable x ∈ X := [0. 3,2]. A number of Chebyshev model enclosures for various expansion orders q in the range 2-7 are shown on the top-left plot of Fig. 2, along with the corresponding approximation errors on the opposite plot. The scheme of Chebyshev models appears to be globally convergent as q → ∞, both with the derivative-based remainder formula and with the exact remainder formula for the univariate outer-composition operations, as shown on the bottom-right plot of Fig. 2. However, the use of the exact remainder significantly tightens the bounds and improves the rate of convergence. In contrast, classical Taylor models fail to converge to the function as the order q is increased, despite the function f being analytic on the set X . This is because the variable range of interest is partially outside the radius of convergence for this function, which also means that the width of the interval remainder increases with q. Finally, the bottom-left plot on Fig. 2 shows that local convergence of qth-order Chebyshev model is of order q + 1 around x = 1. This observation is consistent with the analysis conducted in Sect. 4 since f is analytic. 18 20 derivative-based remainder reformulation as |x| = √ x 2

Example 2
We consider the univariate function f given by for the variable x ∈ X := − π 2 , π 4 . Notice that classical Taylor models may not be used to bound the function f on X given the nonsmoothness at x = 0. In contrast, Chebyshev models can be constructed since the univariate sin(x) is analytic in R, whereas the derivative of |x| has total bounded variation on X . The top plots in Fig. 3 shows how the bounds improve with a larger expansion order q, up to q = 16. It is evident that a much higher-order expansion is necessary in comparison with the smooth function in the previous example (Sect. 5.1), and that the approximation error is the largest close to x = 0. In agreement with the theory in Sect. 4, the local convergence of any qth-order Chebyshev model at this point is of order 1, regardless of the expansion order q ≥ 0; a prediction that is confirmed in the bottom-left plot of Fig. 3. But despite the nonsmoothness at x = 0 and the linear local convergence around this point, Chebyshev models are found to be globally convergent on X as q → ∞, when the remainder formula (29) is used for bounding the remainder of the univariate term |x| (bottom-right plot on Fig. 2).
Notice that this function may also be reformulated by substitution of |x| = √ x 2 , thus providing an alternative route to creating the Chebyshev model. This reformulation makes the use of the exact remainder estimate for the square-root term possible, which explains why tighter bounds for Chebyshev models of order q = 5 or less are obtained here with the reformulation approach (bottom-right plot of Fig. 3). However, the reformulation also causes Chebyshev models to diverge on X as q → ∞, a behavior attributed to the fact that the total variation of the square-root term is unbounded on intervals containing zero.

Example 3
As a final example, we consider the multivariate function f given by for the variable x ∈ X := [−0.6, 0.6] n , with the dimension n taken as 2, 3 or 4. The Chebyshev model enclosure and corresponding approximation error in the 2-variable case are presented on the left part of Fig. 4 for an expansion of order q = 8, showing a tight approximation of the function. From the right plot of Fig. 4, it is seen that the computational burden increases with both the expansion order q and the dimension n of the problem. Moreover, computing Chebyshev models (solid lines) is more demanding than Taylor models (dashed lines), due to the extra complexity of the polynomial product in Chebyshev basis (see Fig. 1 above and related discussion). However, this extra complexity pays off in terms of the tightness of the resulting bounds. For n = 2, one must consider expansions of order greater than q = 18 in order for a Taylor model's remainder bound to be as tight as those given by a second-order Chebyshev models only. For n = 3 or n = 4 likewise, although the Taylor models appear to be globally convergent on X , an expansion order far greater than q = 20 would be necessary in order to be competitive with just second-order Chebyshev models. Overall, Chebyshev models thus outperform their Taylor counterparts. Finally, it is interesting to note that as the function f being even, odd ordered polynomial models prove to provide tighter bounds than even ones, thus explaining why there is sometimes little improvement (or even deterioration) in the bounds even though the computational time increases.

Conclusions
This paper has formalized an arithmetic for the propagation of Chebyshev models through binary sums, binary products and univariate outer-composition operations in factorable functions. A simple method to calculate the exact remainder bounds for certain univariate terms, including exp(·), log(·), √ · and (·) −1 , is established. Adaptations of existing range bounders for Taylor models in order to use them with Chebyshev models are also discussed, which are a vital part of polynomial model arithmetic. The local convergence of Chebyshev models has been analyzed and proven to be equivalent to Taylor models, although the Chebyshev model remainder bounds are often found to be orders of magnitude better than their Taylor model counterparts over larger variable domains. This behavior is supported by the result that Chebyshev models are globally convergent through addition/subtraction and multiplication operations. The global convergence property is also conjectured for composition of Chebyshev models under mild conditions, although sufficient conditions are yet to be formally established. Such convergence properties have been illustrated through several numerical examples, based on an implementation in the authors' in-house library MC++. These examples also illustrate some of the advantages of Chebyshev models over Taylor models, including being able to bound functions with points of non-smoothness or a divergent Taylor expansion. Even though Chebyshev models of equivalent order are computationally more expensive to create than their Taylor counterparts, mainly due to the much greater number of operations required for the binary product operation, they provide benefits by being significantly tighter. This may allow for a much lower order polynomial model to be used and creates a net benefit in terms of computational effort.
On the computational side, the current implementation creates a polynomial which is of the same order for all the variables and is stored using a dense representation. For many applications however, a large number of coefficients may be equal to zero, thereby calling for a sparse implementation of Chebyshev model arithmetic. Such an implementation would reduce both the computational time and the memory requirements. Another area of improvement would concern the binary product operation, whose complexity currently scales exponentially in the number of variables and the expansion order. Several studies have explored approaches to speeding up the multiplication of univariate and bivariate polynomials in Chebyshev basis (e.g., [2,5,22,45]). Ways to use some of these developments for multivariate Chebyshev polynomials could be explored as part of future work. Binary product operations are also where much of the overestimation occurs, given that the polynomial part of the Chebyshev model product may not match the Chebyshev expansion of the product function. Instead, it might be possible to directly create multivariable Chebyshev approximations, for instance by following a similar approach as in chebfun2 [56]. This idea is also similar to the recent work on multivariate McCormick relaxations [59], where instead of decomposing the factorable function down to binary sums, products and univariate compositions, it would be possible to directly bound some multivariate terms. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix 1: Proof of Lemma 1
Since F in Lipschitz continuous on [−1, 1], it has a unique representation as an absolutely and uniformly convergent Chebyshev series, and we have with the Chebyshev coefficients given by Observe that F − P q F = ∞ k=q+1 |a k | under the following scenarios: (i) the coefficients a k with k ≥ q + 1 are of the same sign, either non-negative or nonpositive, in which case |F(1) − P q F (1)| = ∞ k=q+1 |a k | by the property that T k (1) = 1, ∀k ≥ 0; (ii) the coefficients a k with k ≥ q +1 alternate sign, either non-negative for even coefficients and non-positive for odd coefficients or vice-versa, in which case |F(−1) − P q F (−1)| = ∞ k=q+1 |a k | by the property that T k (−1) = (−1) k , ∀k ≥ 0. Without loss of generality, we shall prove that case (i) holds with non-negative Chebyshev coefficients a k when F has an analytic extension on the closed unit disk and all its successive derivatives are nonnegative on [−1, 1]; the other three cases follow readily by symmetry, considering −F(ξ ), F(−ξ) or −F(−ξ).
Under the assumption that F has an analytic extension on the closed unit disk, its Taylor series at any pointξ ∈ [−1, 1] converges to F on [−1, 1]. Atξ = 0 in particular, we have Recall that any polynomial P q of order q ≥ 1, given in monomial basis by P q (ξ ) := q k=0 α k ξ k , has a unique representation in Chebyshev basis as P q (ξ ) = If follows from the nonnegativity of all the successive derivatives F (i) that a k ≥ 0 for every k ≥ 1.
Regarding Chebyshev interpolation, on the other hand, we have with the coefficients a k given by a k := 2 q + 1 q j=0 F(ζ j )T k (ζ j ) withζ j := cos π 2 Like previously, we shall prove that in the case that F has an analytic extension on the closed unit disk and all its successive derivatives are nonnegative on [−1, 1]-the other cases following by symmetry. Since a k ≥ 0 for every k ≥ 1 under these assumptions, it is sufficient to prove that a k ≥ a k for each k = 1, . . . , q. By the aliasing property of the Chebyshev polynomials [37], we have ∀k ∈ {0, . . . , q}, a k = a k − a 2q+2−k − a 2q+2+k + a 4q+4−k + a 4q+4+k − · · · (41) and, therefore, the result follows by showing that a k ≥ a k+2 for every k ≥ 1. Considering the case of odd coefficients first, (39) gives and a similar result holds for even coefficients.

Appendix 2: Proof of Lemma 3
Let the Chebyshev coefficients of the function F on the subset W ⊂ X be denoted as

Appendix 3: Proof of Lemma 4
We start by showing, using finite recursion on k = q + 2, q + 1, . . . , 1, that the intermediate where b q k,Y,κ denote the coefficients of P q β k ,Y .
-This property is trivially satisfied for the schemes (P  Y ) min{r,q+1} ), and so R k β k ,Y satisfies condition (43). Concerning the polynomial part, the effect of multiplying P q β k+1 ,Y with P q Φ,Y is a reduction of the order of all the coefficients b q k+1,Y,κ by one order. Moreover, subtracting P q β k+2 ,Y does not change the order of the coefficients and by Lemma 3 the lead coefficient ϕ q B q Y ,k has order O((rad Y ) min{k,s} ). Therefore, the coefficients b q k,Y,κ satisfy condition (42), too.
The same reasoning applies in taking the recurrence one step further through (20), thereby giving the result.