Nyoman Juniarta
email: nyoman.juniarta@loria.fr

Victor Codocedo
email: victor.codocedo@inria.cl

Miguel Couceiro
email: miguel.couceiro@loria.fr

Amedeo Napoli
email: amedeo.napoli@loria.fr

Biclustering Based on FCA and Partition Pattern Structures for Recommendation Systems

This paper focuses on item recommendation for visitors in a museum within the framework of European Project CrossCult about cultural heritage. We present a theoretical research work about recommendation using biclustering. Our approach is based on biclustering using FCA and partition pattern structures. We investigate the possibility of incorporating the order information using this approach. Then, given the dataset of visitor trajectories, the result of our biclustering can be used to build a collaborative recommendation system.

Introduction

CrossCult (http://www.crosscult.eu) is a European project whose idea is to support the emergence of a European cultural heritage by allowing visitors in different cultural sites (e.g. museum, historic city, archaeological site) to improve the quality of their visit by using adapted computer-based devices. Such improvement can be accomplished by studying, among others, the possibility to build a dynamic recommendation system.

Here, our objective is to study a dynamic recommendation system for visitors in a museum (Hecht Museum for the case study). Given a new visitor V n , the task is to suggest a museum item that may be interesting for him/her. Based on how a suggestion is made for a V n , a recommendation system can be classified into one of the three categories [START_REF] Adomavicius | Toward the next generation of recommender systems: A survey of the state-of-the-art and possible extensions[END_REF]:

-Content-based recommendations: The system makes a suggestion based on the previous visited items of V n . -Collaborative recommendations: The system looks for previous users who have similar interest with V n , and makes a suggestion based on their visited items. -Hybrid approaches: The combination of content-based and collaborative approaches. Our method belongs to the second category (collaborative recommendation). First we group all previous users based on their visit trajectories using biclustering. When V n arrives, we try to find a G s , i.e. a group who shares a similar interest to V n . Then, based on the behavior of the visitors in G s , we can suggest one item that may be interesting for V n .

An approach using partition pattern structures was proposed to obtain biclusters with constant (or similar) values on the columns [START_REF] Codocedo | Lattice-based biclustering using partition pattern structures[END_REF]. Our main contribution 1 1 4 3 5 1 1 2 5 1 3 3 4 2 1 TAB. 1 -A bicluster (shaded) with constant value. here is to extend this approach to mine another type of biclusters: those with coherent evolution on the columns (CEC biclusters). We then use them to build a recommendation system. We focus on this bicluster type since we are dealing with the dataset of trajectories where each trajectory corresponds to an ordered list of items. Hence we are interested to take into account this ordering information, and CEC biclustering is one appropriate way to do that.

Biclustering

We consider a dataset composed of a set of objects, each of which has values over a set of attributes. This dataset can be represented as a numerical matrix C, where each entry c ij indicates the value of object i w.r.t attribute j. One may be interested in finding which subset of objects possesses the same values w.r.t. a subset of attributes. Regarding the matrix representation, this is equivalent to the problem of finding a submatrix c that has a constant value over all of its elements (example in Table 1). This task is called biclustering with constant values, which is a simultaneous clustering of the rows and columns of a matrix.

Other than constant values, the bicluster approach also focused on finding other types of submatrices, as shown in Table 2. A bicluster with constant columns (rows) is a submatrix where each column (row) has the same value, as illustrated in Table 2a (Table 2b, resp.). In a bicluster with additive coherent values, the value of each cell c ij follows the equation γ + α i + β j , where γ is a constant, α i is a constant value for row i, and β j is a constant value for column j. For example, if γ = 1, (α 1 , α 2 , α 3 , α 4 ) = (3, 2, 4, 6), and (β 1 , β 2 , β 3 , β 4 ) = (0, -2, 1, -1), then we can obtain the bicluster in Table 2c. Another interesting type is the CEC bicluster, also known as order-preserving submatrix [START_REF] Ben-Dor | Discovering local structure in gene expression data: the order-preserving submatrix problem[END_REF]. In this type of bicluster, each row induces the same linear order across all columns. For example, in the bicluster in Table 2d, each row follows column1 < column2 < column4 < column3. Moreover, a bicluster with coherent evolution on the rows can be defined similarly, as shown in Table 2e.

Those different types of biclusters are useful when we are interested to identify a group of people who behave similarly according to a set of attributes. This group identification is necessary in the task of collaborative recommendation, because in the process of making a suggestion to a person, we first identify the people who are similar to him/her. N. Juniarta et al. TAB. 3 -Order of visit of 7 items, observed from certain visitors.

i 1 i 2 i 3 i 4 i 5 i 6 i 7 v 1 1 2 3 4
# Visitors Items (in order) B1 v 1 i 1 , i 2 , i 3 , i 4 , i 5 , i 6 , i 7 B2 v 2 i 6 , i 1 , i 4 , i 2 , i 3 , i 7 , i 5 B3 v 1 , v 3 i 3 , i 4 , i 5 , i 7 B4 v 1 , v 3 i 1 , i 4 , i 5 , i 7 B5 v 1 , v 4 i 3 , i 5 , i 6 B6 v 2 , v 3 i 6 , i 1 , i 4 , i 5 B7 v 1 , v 2 , v 3 , v 4 i 2 , i 7
TAB. 4 -Some CEC biclusters in Table 3.

Recommendation

In the context of CrossCult, we are working on a visitor dataset that comprises several trajectories in a museum. Within this project, our main objective is to build a dynamic recommendation system for new visitors. This system should be able to suggest a museum item to them based on their trajectory, by looking at the dataset of previous visitors. Also, it should be able to update the suggestion as they move inside the museum. Since we take into account the visit ordering of items when studying the trajectories, CEC biclustering seems to be a suitable approach. Here we will define some strategies for recommendation systems based on this type of biclustering.

Consider the dataset about 4 visitors in a museum with 7 items, given by Table 3. The numbers in row x indicate the path of visitor v x . For example, the path of

v 2 is i 6 → i 1 → i 4 → i 2 → i 3 → i 7 → i 5 .
Now we have a new visitor, v a , who recently arrives in the museum. She visits i 3 , followed by i 4 . Our task is to recommend a new item to her, by studying the CEC biclusters over the first four visitors. Some of those biclusters are listed in Table 4. From B7, we can see that all four visitors visit i 7 after (not before) i 2 . Meanwhile in B3, regarding the items {i 3 , i 4 , i 5 , i 7 }, v 1 and v 3 follow the same order: i 3 → i 4 → i 5 → i 7 . These two visitors also agree in the order of items {i 1 , i 4 , i 5 , i 7 }, as seen in B4.

Those CEC biclusters can be studied to give a recommendation for v a , by focusing on those who are similar to him/her. In the following subsections, we will present some possible strategies for performing this recommendation. They start from the same step: finding a CEC bicluster with a similar path to v a .

This strategy focus on following sequential patterns in the dataset. The idea is if many people have the path i a → i b → i c , then we should recommend the item i c to a visitor who has done i a → i b .

Since v a has path i 3 → i 4 , we focus on the CEC biclusters that have those two items, i.e. B1, B2, and B3. One of those biclusters (B2) has a different ordering (i 3 after i 4 ), and thus we filter it out. Then, in B3 for example, the path is i 3 → i 4 → i 5 → i 7 . Therefore, we can recommend i 5 to v a .

In this strategy, we assume that the order of visit reflects the order of interest. In this case, for v 2 , i 6 is the most interesting (since (s)he visited it first) while i 5 is the least interesting (since (s)he visited it last). Based on this assumption, if many people have the path i a → i b → i c that reflects their interest, then we should recommend the item i a to a visitor who has done i b → i c , assuming that (s)he missed this item earlier in his/her visit.

Similar to the previous strategy, here we first identify CEC biclusters with similar path to v a . According to the visitors in B3, i 1 and i 3 are the most interesting items for v 1 and v 3 respectively. Therefore, we can suggest those two items to v a . Furthermore, we know that v a has already visited i 3 , so we recommend her to visit i 1 .

In the framework of the CrossCult project, we are working on a specific dataset about the trajectories of 254 visitors in Hecht Museum in Haifa, Israel [START_REF] Lanir | The influence of a location-aware mobile guide on museum visitors' behavior[END_REF]. In this dataset, there are 52 items over 8 rooms (A-G). Therefore, the visitor-item matrix for biclustering will have 254 rows and 52 columns.

Conclusion

In this work, we explore an approach to build collaborative recommendation strategy for visitors in a museum. This strategy takes into account the order of interest or the order of visit for each visitor, and we showed how to use CEC biclustering to obtain a set of similar visitors. Further comparisons of CEC biclustering and sequential pattern mining [START_REF] Yan | CloSpan: Mining: Closed sequential patterns in large datasets[END_REF]; [START_REF] Yu | A dynamic recurrent model for next basket recommendation[END_REF] should be investigated, in particular, regarding their complexities and their results.

  Examples of some bicluster types. (a) Constant columns, (b) constant rows, (c) additive coherent values, (d) coherent evolution on the columns, and (e) coherent evolution on the rows.


Joris Van Der Hoeven
email: vdhoeven@lix.polytechnique.fr

On high precision integration of sti di erential equations

Keywords: sti di erential equation, numerical method, reliable computation, rigorous integration, Taylor models, multiple precision computations A.M 65L04, 65L20, 65G20, 37-04

In this paper, we present various new algorithms for the integration of sti di erential equations that allow the step size to increase proportionally with time. We mainly focus on high precision integrators, which leads us to use high order Taylor schemes. We will also present various algorithms for certifying the numerical results, again with the property that the step size increases with time.

INTRODUCTION

Consider the di erential equation + = ( ), (1) where d d is a diagonal matrix with positive real entries 0 } 1 } } d , and where F 1] , &, F d] ] d is a polynomial forcing term. We are interested in analytic solutions :I d , where I is either an interval of the form 0,T] with T >0, or a larger subset of that contains an interval of this form. The d components of this mapping will be written 1] , &, d] : I , whereas subscripts will be reserved for the coe cients of local power series solutions. A similar notation will be used for the components of other vectors and matrices.

If the largest eigenvalue d of gets large with respect to the forcing term , then the equation ( 1) is said to be sti . Generic numerical integration schemes experience di culties with this kind of equations, whether we use Euler's method, a Runge Kutta scheme, or high order algorithms based on Taylor expansions. Roughly speaking, the problem is that all generic methods become inaccurate when the step size exceeds d 1 . One of the most di cult cases is when d is large, but each of the quotients k+1 / k remains reasonably close to one. This happens for instance when k = k 2 , a case that is naturally encountered when discretizing certain partial di erential equations.

Implicit integration schemes allow for larger step sizes than generic methods, but involve the computation and inversion of Jacobian matrices, which may be expensive in high dimension d. But even such more expensive schemes only seem to deal with the step size issue in a somewhat heuristic manner: one typically selects a scheme that is stable when applied to any equation 2 + =0 for all ~0, and then hopes that the scheme continues to behave well for the actual target equation (1). We refer to 25, Section 17.5] and 4, Section 6] for more information on classical integration schemes for sti di erential equations.

We are mainly interested in high precision computations, in which case it is natural to use high order integration schemes that are based on Taylor expansions. In this context, we are not aware of any numerical method that allows the step size to increase proportionally with time. The main aim of this paper is to present such a method, together with several variants, as well as a way to compute rigorous bounds for the error. Our fastest method is an explicit scheme, but its convergence deteriorates when the eigenvalues i are close. The slower and rigourous counterparts rely on the computation of Jacobian matrices.

In sections 3 and 4, we recall various approaches that can be used to compute truncated series solutions to initial value problems and how to derive high order integration schemes from this. More precisely, given a numerical approximation t for (t) at time t, we compute the rst n terms of a power series solution f t to the initial value problem

f t z + f t = ( f t ), f t 0 = t (2)
and return f t 0 + f t 1 + + f t n 1 n 1 as an approximation of (t+ ) for some suitable step size . The inaccuracy of such schemes for larger step sizes is due to the fact that the k-th coe cient of f t tends to change proportionally to d k /k! for slight perturbations of t . If d is large, then even small multiples of d k /k! quickly dominate f t k , which severely limits the step sizes that we can use.

This numerical instability is an unpleasant artifact of traditional methods for the integration of di erential equations. It is well know that solutions of sti equations tend to become very smooth after a certain period of time (in which case one says that the system has reached its steady state), but this does not lead to larger step sizes, as we would hope for.

In order to make this idea of smoothness after a certain period of time more precise, it is instructive to study the analytic continuation of in the complex plane. In section 2, for a xed initial condition (0)=c and a xed forcing term , we show the existence of a compact half disk H R = {t : |t| } R, Re t ~0}, R > 0, on which the solution of (1) is analytic, for any choice of . In particular, for t 0,R/2], the radius of convergence t of the exact Taylor expansion f t of at time t is at least t. Setting K =max z H R | (z)|, Cauchy's formula thus yields the bound | f t k |}K/t k for all k .

If t k d 1 , then this means that e cient integration schemes should allow for step sizes of the order t at time t, for some xed constant >0. In section 5, we will present various such schemes. They are all based on the computation of good approximations f t ;n = f t 0 + f t 1 z+ + f t n 1 z n 1 of the exact Taylor expansion f t ;n of at time t and order n. By what precedes, such an approximation cannot be obtained by solving (2). Instead, for a suitable so-called critical index {0,&,d} that depends on t, we will solve the following so-called steady-state problem:

f t z + f t = ( f t ), { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { f t 0 i] = t i] if i } f t n i] =0 if i > . ( 3
)
Intuitively speaking, the system has reached a steady state for the eigenvalues i with i> . Such eigenvalues i are large , so it is more accurate to determine coecients f t k i] as a function of later coe cients f t k+1 , &, f t n 1 rather than earlier ones. Furthermore, the coe cients f t n i] are so small that they can be approximated by zeros; this explains the steady-state conditions f t n i] =0. The remaining eigenvalues i with i} are small enough for the step size under consideration that the dependence of f t on the initial condition t i] is su ciently moderate for (3) to admit an accurate solution. We say that the system (1) is in a transient state for such eigenvalues i at time t.

A nal topic is the computation of rigorous error bounds for the numerical integration process. For di erential equations that are not sti , this is a classical theme in interval analysis [START_REF] Moore | Automatic local coordinate transformations to reduce the growth of error bounds in interval computation of solutions to ordinary di erential equation[END_REF][START_REF] Hn | Rigourously computed orbits of dynamical systems without the wrapping e ect[END_REF][START_REF] Butcher | Numerical methods for ordinary di erential equations in the 20th century[END_REF][START_REF] Lanfort | A computer assisted proof of the Feigenbaum conjectures[END_REF][START_REF] Lohner | On the ubiquity of the wrapping e ect in the computation of error bounds[END_REF][START_REF] Makino | Remainder di erential algebras and their applications[END_REF][START_REF] Van Der Hoeven | Newton's method and FFT trading[END_REF]. It has long been an open problem to develop e cient reliable integrators for sti di erential equations. We report on progress in this direction in section 6.

ANALYTIC CONTINUATION

Consider the equation (1) for a xed initial condition (0) = c d and a xed forcing term , but where we will allow ourselves to vary . Our main aim is to show that there exists a compact half disk H R = {t :|t|}R,Re t ~0} with R > 0 on which the solution to (1) exists and is analytic, independently of the choice of . We will both prove a weak and simpler version of this result and a stronger version that allows us to take larger radii R and that is also more convenient if we want to actually compute a radius R that works.

Notations

We will need a few more notations. We denote the interior of H R by H R ={t :|t|<R,Re t> 0}. Given a non-empty open subset U , we write (U) for the Banach space of analytic functions g: U with nite sup-norm g =sup z U

|g(u)|.

We will often use vector notation: given v d , g (U) d , and r, s ( ~)d , we write |v|

= (|v 1] |,&,|v d] |) g = ( g 1] ,&, g d] ) r }s r 1] }s 1] ' ' r d] }s d] .
In addition to such componentwise norms, we will also use more traditional supnorms: given v d , g (U) d , and an n n matrix M d d , we de ne

|v| = max(|v 1] |,&,|v d] |) g = max( g 1] ,&, g d] ) |M| = sup v d ,|v| =1 |Mv| .

Reminders from complex analysis

We will also need two well known results from complex analysis, namely an explicit version of Cauchy-Kovalevskaya's theorem and a variant of Montel's theorem. THEOREM 1. Let (F)= (F) F, m c =|c| , and consider the initial value problem = ( ), (0)=c. (

) Writing = i 1 ,&,i d i 1 ,&,i d (F 1] ) i 1 (F d] ) i d , let M c = max i 1 ,&,i d | i 1 ,&,i d |(2m c d) i1+ +i d . 4
Then [START_REF] Brent | Modern Computer Arithmetic[END_REF]

c X ) ) ) ) ) ) ) ) ) ) ) ) ) ) m c ,
of radius c . Since the formal power series solution of (4) satis es k i] } k i] =g k for all i {1,&,d} and k , it follows that has radius of convergence at least c . THEOREM 2. Given a non-empty subset U of and a bounded sequence g 1 ,g 2 ,& (U) d , we can extract a subsequence g k 1 , g k 2 , & that converges uniformly to a limit g (U) d on every compact subset of U.

Proof. If d=1, then this is Montel's theorem. The general case is obtained using an easy induction on d: we rst extract a subsequence g i 1 ,

g i 2 , & such that g i 1 d] , g i 2 d]
, & converges to a limit in (U) and then apply the induction hypothesis to this subsequence for the remaining d 1 components.

Analyticity on open half disks

Let J = / F denote the Jacobian matrix of . THEOREM 3. Let R> 0, B ( > ) d , and C <1 be such that for all z , we have

|z c|} B | (z)|R}B. (5) |z c|} B |J (z)| R} C. ( 6
)
Then, for any choice of , the initial value problem + = ( ), lim t 0 (t)=c [START_REF] Fischer | Fast on-line integer multiplication[END_REF] admits a unique analytic solution on H R . In addition, we have | (t) c| }B for all t H R .

Proof. Consider the operator R : (H R ) d (H R ) d with R (g)(t) = c +e t 5 0 t e u (g(u))du.

The xed points (H R ) d of this operator are precisely the solutions of [START_REF] Fischer | Fast on-line integer multiplication[END_REF]. Let us prove the existence of such a xed point. The uniqueness follows by analytic continuation from the well known uniqueness of the solution of the initial value problem [START_REF] Fischer | Fast on-line integer multiplication[END_REF] Let us next show that R is actually contracting on , R (c, B). Given g, h , R (c, B), consider the homotopy =(1 )g+ h , R (c,B) with 0,1]. From (6), we get (h) (g) = 5 0 1 ( ) d

} 5 0 1 J ( )(h g) d = 5 0 1 sup t H R |J ( (t))(h(t) g(t))| d } C R 5 0 1 sup t H R |h(t) g(t)| d = C R h g . It follows that 5 0 t e u ( (h(u)) (g(u)))du } C R h
Then, for any choice of , the initial value problem + = ( ), (0)=c (9) admits a unique analytic solution on H R . In addition, we have | (t) c| }B for all t H R . Proof. Let R : (H R ) d (H R ) d be as in the proof of Theorem 3 and notice that R still maps , R (c,B) into itself. Now consider the sequence g 1 ,g 2 ,& with g k T R k (c) , R (c,B). Applying Montel's theorem, we obtain a subsequence g k i that converges uniformly to a limit g , R (c,B) on every compact subset of H R .

We claim that g is a xed point of R . Indeed, for a su ciently small R2 > 0 with R2}R, we have |J (z)| R2} /

1 2 for all z with |z c|}B. As in the proof of Theorem 3, this means that R2 is contracting on , R2 (c, B). Consequently, the restrictions g 1 m , g 2 m , & of the functions g 1 ,g 2 ,& to H R2 with g k m = R2 k (c) tend to a xed point g m , and so does the subsequence g k i m . Now this xed point g m coincides with g on H R2 and also solves the initial value problem [START_REF] Fischer | Fast on-line integer multiplication[END_REF]. By analytic continuation, g therefore solves the same initial value problem on H R , which completes the proof of our claim.

It remains to be shown that g can be extended to an analytic function on H R that satis es | (t) c|}B for all t H R . Since g 2 = (g ) g on H R , we rst notice that |g 2 | is bounded on H R , whence g continuously extends to a function that is de ned on the whole of H R , and we clearly have | (t) c|} B for all t H R . Now let us consider the unique solution 0 | to the initial value problem 2(t) + (t) = ( (t)) with (0)= . Theorem 1 provides us with a lower bound 0| >0 for the radius of convergence of 0| , which depends continuously on . By the compactness of H R , it follows that there exists some > 0 with 0 | (t) ~ for all t H R . Now consider t on the boundary H R of H R and let t2 H R be such that |t2 t| < . Then 0 | (t2) ~ and 0 | (t2) (z) = (t2 +z) on some neighbourhood of t2. We conclude that 0| (t2) (z t2) is an analytic extension of to the open disk {u :|u t2|< } that contains t.

COMPUTING FORMAL POWER SERIES SOLUTIONS

Before we present algorithms for the numeric integration of (1), let us rst consider the question how to compute formal power series solutions. Since we will later be considering power series solutions at various times t, we will use the letter f instead of for such solutions. Notice that a vector f z]] d of d formal power series can also be regarded as a formal power series f = f 0 + f 1 z + f 2 z 2 + d z]] with coe cients in d . Let c d be the initial condition. Setting ( f ) = ( f ) f , we are thus interested in the computation of truncated power series solutions to the equation f z = ( f ), f 0 =c.

(10) Alternatively, this equation can be rewritten in integral form f = c +5 ( f ), [START_REF] Van Der Hoeven | Newton's method and FFT trading[END_REF] where the integration operator + sends g z]] d to g 0 z+ 1 2 g 1 z 2 + 1 3 g 2 z 3 + . In this section, we recall several approaches to the computation of power series solutions to such equations. The e ciency of each method can be measured in terms of the number of eld operations in that are required to compute the rst n terms of the solution. For the time being, we assume that all computations are done at a xed bit precision p, so that operations in are done at unit cost. One may also take into account the number s of scalar operations that are necessary for the evaluation of . A further re nement is to separately count the number s mul of multiplications that are needed. For instance, if d =1 and ( f )=( f f +3) f 100 f , then we have s =5 and s mul =3. We always assume that d =O(s). Iterative method. For the rst method, we systematically compute with truncated power series of order n in z] / (z n ), which are also called jets of order n. One addition in ]/( n ) reduces to O(n) additions in and similarly for scalar multiplications with an element in . A naive multiplication in ]/( n ) requires O(n 2 ) operations in , although this cost can be reduced to O(n log n) using FFT techniques. The integration operator + sends (g 0 + + g n 1 z n 1 ) mod z n to (g 1 + + g n 2 /(n 1) z n 1 ) mod z n and can be computed in linear time. Now assume that f ( z]/(z n )) d is an approximate solution to [START_REF] Van Der Hoeven | Newton's method and FFT trading[END_REF] whose rst k<n terms are correct. In other words, if f is the actual solution and f z] d is a preimage of f with f =( f mod z n ), then f f =O(z k ). Given such an approximation, one iteration f T c +5 ( f ) (12) of ( 11) yields a new approximation whose rst k + 1 terms are correct. Starting with f T(0 mod z n ), we thus obtain a solution modulo z n of (11) after at most n iterations. The total cost of this computation is bounded by O(s mul n 3 +s n 2 ), or by O(s mul n 2 log n+s n 2 ) when using FFT-techniques. Recurrence relations. Since ( f ) is a polynomial, it is built from the components of f using additions, multiplications, and scalar multiplications with constants in . For sums, scalar products, general products, and integrals of power series, we may extract their coe cients in z k using the following rules (g+h

) k = g k +h k (13) ( g) k = g k (14) (gh) k = g 0 h k +g 1 h k+1 + +g k h 0 (15) (+g) k+1 = 1 k+ 1 g k . ( 16
) Applying this recursively to the polynomial expression ( f ), the iteration (12) yields a recursion relation

f k+1 T 1 k +1 ( ( f )) k ( 17
)
that allows us to compute f k from f 0 , by induction on k. Proceeding in this way, the computation of the solution f 0 + + f n 1 z n 1 at order n requires O(s mul n 2 +sn) operations. The lazy power series approach. The above approach of computing the coefficients f 1 , f 2 ,& successively using the recursion relation ( 17) can be reformulated elegantly in the framework of lazy power series computations. The idea is to regard a power series g z]] as given by its stream g 0 , g 1 , g 2 , & of coe cients. Basic arithmetic operations on power series are implemented in a lazy manner: when computing the k-th coe cient of a sum g+h, a product g h, or an integral +g, we compute just the coe cients that are needed from g and h. The natural way to do this is precisely to use the relations (13 16).

The lazy power series approach has the important property that the k-th coe cient (g h) k of a product (say) becomes available as soon as g 0 ,&,g k and h 0 ,&,h k are given. This makes it possible to let g and h depend on the result g h, as long as the computation of g k and h k only involves previous coe cients (g h) 0 , &, (g h) k 1 of g h. As a consequence, the xed point equation ( 11) can be solved simply by evaluating the right hand side using lazy power series arithmetic. Indeed, the k-th coe cient (+ ( f )) k only depends on previous coe cients f 0 ,&, f k 1 of f .

The relaxed power series approach. The main drawback of the lazy power series approach is that the computation of a product at order n using [START_REF] Hn | Rigourously computed orbits of dynamical systems without the wrapping e ect[END_REF] requires O(n 2 ) operations. Here we recall that FFT techniques allow us to compute a product in z]/(z n ) using only O(n log n) operations.

One essential observation is that, in order to solve (11) using the lazy power series approach, we only relied on the fact that each coe cient (g h) k becomes available as soon as g 0 , &, g k and h 0 ,&,h k are given. In fact, it is possible to design faster multiplication methods that still have this property; such methods are said to be relaxed or on-line. A relaxed multiplication method that computes a product at order n in time O(n log 2 n) was presented in 8] and can be traced back to 6]. An even faster algorithm of time complexity n log ne O loglogn ( was given in 9]. Denoting by R(n) = n (log n) 1+o(1) the cost of relaxed multiplication at order n, the resolution of (11) at order n now requires only O(s mul R(n) +sn) operations. Newton's method. Yet another idea to speed up computations is to replace the iteration (12) with a Newton-style iteration with faster convergence. This idea was rst used by Brent and Kung 2] to show that [START_REF] Van Der Hoeven | Newton's method and FFT trading[END_REF] can be solved at order n in time O(n log n). However, this complexity analysis does not take into account the dependence on d and s. In particular, the dependence on d of Brent and Kung's method is exponential 8]. Faster algorithms were proposed in 26, 1], based on the simultaneous computation of the solution f and its rst variation. This allows for the computation of a solution of (10) at order n in time O((s mul + d) d n log n + s d n). Whenever d = o(log n), the computation time can be further reduced to O((s mul + d) n log n + s n) 10]. For small d, this leads to the asymptotically most e cient method for solving [START_REF] Van Der Hoeven | New algorithms for relaxed multiplication[END_REF]. For large d, the computation of the rst variation induces a (d) overhead, and the relaxed method usually becomes more e cient.

NUMERICAL INTEGRATION USING TAYLOR EXPANSIONS

Let ( ) = ( ) as before. The aim of this section is to present a naive algorithm for the numerical integration of the di erential equation = ( ), [START_REF] Lohner | On the ubiquity of the wrapping e ect in the computation of error bounds[END_REF] based on the computation of truncated power series solutions at successive times t 0 = 0 < t 1 < t 2 , &. We will use a xed expansion order n, a xed bit precision p ~24, but an adaptive step size t 1 t 0 ,t 2 t 1 ,&. At t=0, we are given an initial condition (0)=c d and our aim is to nd a good numeric approximation for (t) at time t = T > 0. The algorithm is not designed to be e cient when the equation gets sti and we will present a heuristic discussion on what goes wrong when this happens.

t+z) = f t ;n +O(z n ) = f t 0 + f t 1 z+ + f t n 1 z n 1 +O(z n ).
If the time t at which we expand is clear from the context, then we will simply drop the post x t and write f instead of f t . Conversely, for any other quantities that implicitly depend on t, we will use the post x t to make this dependence explicit.

So let f = f t be the power series expansion of at a xed time t. In view of [START_REF] Lohner | On the ubiquity of the wrapping e ect in the computation of error bounds[END_REF], this power series satis es the equation f z = ( f ).

(19) In the previous section, we have recalled various methods for computing truncated power series solutions f ;n as a function of the initial condition f 0 = = t = (t) at time t. In what follows, we will use any of these algorithms as a black box, and show how to device a numerical integration scheme from that.

Obviously, given (t) and an appropriate step size = t at time t, the idea is to compute f ;n = f t ;n and simply evaluate (t+ )H f t ;n ( )= f t 0 + f t 1 + + f t n 1 n 1 . We next continue with t + in the role of t and with a suitably adapted step size . The main question is therefore to nd a suitable step size . Now the expected order of magnitude of f t ;n ( ) is given by

M = max k<n | f k | k . ( 20
) Since we are computing with p bits of precision, we wish to keep the relative error of our numeric integration scheme below 2 p , approximately. Therefore, we need to ensure that the truncation error

E = | | | | | | | | | | | | | | |y k~n f k k | | | | | | | | | | | | | | | ( 21
)
remains bounded by M 2 p . Although we have not computed any of the coe cients f k for k ~n, a reasonable approximate upper bound for E is given by

E guess = max n }k<n | f k | k , ( 22
)
where >0 is a small positive integer, called the number of guard terms. In order to protect ourselves against the occasional vanishing of f n 1 , it is wise to take >1. Nevertheless, a small value such as = 2 or = 3 should usually provide acceptable upper estimates for E. We now simply take the step size to be maximal with E guess }M 2 p . This leads to the following algorithm for the numerical integration of (18): Algorithm 1 Input: an initial condition c d and T >0 Output: a numerical approximation for (T), where satis es [START_REF] Lohner | On the ubiquity of the wrapping e ect in the computation of error bounds[END_REF] with (0)=c T c, tT 0 while t< T do Compute the truncated solution f ;n to [START_REF] Makino | Remainder di erential algebras and their applications[END_REF] with f 0 = Let be maximal such that E guess }M 2 p , with M and E guess as in [START_REF] Makino | Suppression of the wrapping e ect by Taylor model-based validated integrators[END_REF] and [START_REF] Moore | Interval Analysis[END_REF] T min( ,T t) T f ;n ( ), t T t + return

Step size and error analysis

Let denote the exact solution of ( 18) with (0) = c and let denote the computed approximation by Algorithm 1. In order to analyze various aspects of the algorithm, it is instructive to look at the radius of convergence of at time t. Roughly speaking, with the notations from the previous subsection, the coe cients f k of the power series expansion f = f t of at time t grow as

| f k | H M ( ) k .
Ideally speaking, if we were able to compute f ;n with su cient accuracy, then the coecients f k should grow in a similar way. Setting for the step size in this ideal situation, we would then expect that E guess H EH M ( / ) n H M 2 p , whence H 2 p/n . [START_REF] Neumaier | The wrapping e ect, ellipsoid arithmetic, stability and confedence regions[END_REF] This suggests to take the expansion order n to be proportional to p, after which the step size should be proportional to the radius of convergence . Let us pursue this line of wishful thinking a little further. Theorem 4 implies that is analytic on a compact half disk H R that is independent of . In particular, we get that t ~t for t } R /2. It can also be shown that the radius of convergence of at the origin is of the order 0 H d 1 or more. For sti equations, we typically have R k d 1 . In order to integrate the equation until a time T }R/2, we thus hope for a step size that increases geometrically from d 1 to T 2 p/n . The entire integration would then require approximately 2 p/n log(T d ) steps. Using the most e cient algorithms from section 3, each step requires O (s n) oating point operations (here the at Oh notation O ( ) stands for O( (log ) O(1) )). Since additions and multiplications of p bit oating point numbers can be performed in time O (p), the overall bit cost of the entire integration would thus be bounded by O (snp2 p/n log(T d )).

But are we indeed able to compute f ;n with enough accuracy in order to ensure that the coe cients f k grow according to | f k | HM ( ) k ? Let us carefully analyze each of the sources of error for one step of our integration scheme. By construction, we ensured the truncation error to be of the order |( f f ;n )( )| HM 2 p . One of the most intrinsic sources of error comes from the initial condition f 0 : since we are computing with a precision of p bits, the mere representation of f 0 induces a relative error of the order of 2 p . Even when computing f ;n from f 0 with in nite precision, this intrinsic source of error cannot be removed.

Let us now study how the error in the initial condition a ects the errors in the other coe cients f k . For this, we need to investigate the rst variation of , which describes the sensitivity of the ow to the initial condition. More precisely, let f t| be the power solution of (19) at time t as a function of the initial condition f t | 0 = d . Then the rst variation V t | is the d d matrix with entries V t | i,j] = f t | i] / j] . Dropping t | su xes when they are clear from the context, the rst variation satis es the linear di erential equation

V z = (J ( f ) )V, V 0 =Id d .
Here we recall that J = / F stands for the Jacobian matrix of . If our equation is very sti , then J is small with respect to , which leads to the extremely crude approximation V H e z for the rst variation. Now the relative error in the initial condition = f 0 is at best 2 p , as explained above, which means that | | H M 2 p for = f 0 . In combination with the relation f f H V (

)He z ( ), this leads to the following errors for the other coe cients:

| f k f k | H d k k! M 2 p . Now if d k ( ) 1 , then the error | f n f n | H M d n 2 p /n!
instead of the desired step size H 2 p/n . The actual bit cost of the complete integration is therefore bounded by O (snpT/ d ) instead of O (snp2 p/n log(T d )).

An interesting aspect of this analysis is the fact that the step size still turns out to be n/e times larger than d 1 , whence larger orders n allow for larger step sizes. However, this comes at the expense of a larger relative error than 2 p . Indeed, we have

| f k k f k k | H ( d ) k k! M 2 p H n k k M 2 p .
This error is maximal for k max H n/e, in which case we have | f k max kmax f k max kmax | H e n/e M 2 p . This means that the relative error for one step of our integration method is 2 n/(elog2) p instead of 2 p . In other words, we have sacri ced n/(e log 2) bits of precision, so that the method admits an e ective precision of only p n/(elog 2) bits.

The last source of errors for Algorithm 1 comes from rounding errors during the computation of f ;n from f 0 . The nature of these errors depends on the particular method that we used for computing the power series f ;n . Nevertheless, for most methods, the rounding errors only contribute marginally to the total error. This is due to the fact that

| f k k | 2 p j | f k k f k k | for k > 0,
so the rounding errors are absorbed by the errors induced by the error in the initial condition f 0 .

FIGHTING STIFFNESS

Let us continue with the notations from section 4 and its subsection 4.2. In particular, we assume that the exact solution to (1) with initial condition (0)=c is analytic on the compact half disk H R . We also denote by K an upper bound for | | on H R , so that t ~t and | f t k | }Kt k for all t } R/2 and k , by Cauchy's theorem. From now on, we assume that R k d 1 , so that Algorithm 1 only allows us to use a step size of the order (24) instead of [START_REF] Neumaier | The wrapping e ect, ellipsoid arithmetic, stability and confedence regions[END_REF]. The aim of this section is to present a new way to compute truncated power series solutions of the equation

f z + f = ( f ) (25)
at time t, with the property that f ;n is a good approximation of the true truncated solution f ;n . In a similar way as in section 4, we will then use this to derive an algorithm for the numerical integration of (1). Contrary to before, the property that

| f k | H | f k | } K /t k
for k <n allows us to take step sizes of the desired order (23) for t }R/2.

Integrating sti equations using Taylor series

We stress once more that the reduced step size for Algorithm 1 is a consequence of our choice to compute the truncated solution f ;n of ( 25) in terms of the initial condition f 0 = (that can only be known approximately) and not of the choice of the particular algorithm that is used for this computation. Indeed, as explained in section 4.2, only a small change in of the order of | | H | | 2 p H M 2 p can have a dramatic e ect on the solution f ;n , since the coe cient f n can change by as much as 1 , it follows that a minor change in leads to a completely incorrect computation of the coe cients f k with k close to n.

| f n | H M 2 p d n /n!. Since | f n | } K /t n j M 2 p d n /n! for large tk d
For small t, the error | f n | generally does remain bounded by the actual value | f n | . In the theory of sti di erential equations, it is customary to say that the system is in a transient state for such t. As soon as the error | f n | H M 2 p d n /n! exceeds | f n |, we say that the system has reached its steady state for the largest eigenvalue d of . When this happens, the system may still be in a transient state for some of the other eigenvalues i < d . This motivates the de nition of the critical index {0,&,n} as being the largest index i such that we are in a transient state for i ; we take =0 if we reached the steady state for all eigenvalues i .

The concepts of transient state, steady state, and critical index are deliberately somewhat vague. As a tentative de nition, we might say that we reached the steady state for i if | f n | }M 2 p i n /n!. However, for computational purposes, it is convenient to interpret this inequality as an approximate one. The crucial property of reaching the steady state for i is that the smallness of |( f n ) i] | essentially determines the i-th component i] of the initial condition up to the last p-th bit.

The main idea behind our algorithm is to use this property as a boundary condition for the computation of f ;n . More precisely, we boldly impose the boundary conditions f n +1] = = f n d] =0 as a replacement for the initial conditions f 0 +1] = +1] ,&,f 0 d] = d] , while keeping the remaining initial conditions f 0 1] = 1] , &, f 0 ] = ] for the transient states. In other words, we wish to nd the truncated solution f ;n of the system

f z + f = ( f ), { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { f 0 i] = i] if i } f n i] =0 if i > . ( 26
)
We will call f n i] =0 the steady-state condition for i and (26) a steady-state problem.

In order to solve [START_REF] Press | Numerical recipes, the art of scienti c computing[END_REF], it is natural to adapt the iterative method from section 3 and introduce the operator : z] ;n

d z] ;n d with i] (g) = { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { i] ++( i] (g) i g i] ) if i } i 1 ( i] (g) g i] ) if i > mod z n . ( 27
)
The actual arithmetic is performed over z] / (z n ), but it will be more convenient to view as an operator from z] ;n d into itself. The computation of f ;n as a xed point of leaves us with two questions: how to chose a suitable ansatz for the iteration gT (g) and how to determine the critical index ? For the ansatz, we go back to the solution f t2 ;n from the previous step at time t2 = t 2 and simply use f t2 ;n ( 2 + z) as a rst approximation of the solution f t ;n (z) at time t. For t=0, we fall back to the traditional method from section 4. For the initial and steady-state conditions to propagate to the other end , at least n iterations are required in order to nd a xed point of , whereas 2 n iterations usually su ce. One may thus take f t ;n T 2n ( f t2 ;n ( 2 +z)).

As to the critical index , we determine it as a function of the step size through =max i }d: i } n e .

(28) This choice is meant to ensure the fastest convergence of the iteration g T (g) and we used n /e as an approximation of n! n ( . The step size itself is chosen in the same way as in section 4. Since f ;n is not yet available for the computation of , we may use the previous step size 2 as an approximation for in (28). Optionally, one may also wish to implement a few additional sanity checks in order to ensure convergence of the -iteration and decrease the step size in case of failure. One useful such check is to verify that

| f t 0 f t2 ;n ( 2)| }M 2 p/2 .
Altogether, this leads to the following algorithm: Algorithm 2 Input: an initial condition c d and T >0 Output: a numerical approximation for (T), where satis es (1) with (0)=c T c, tT 0 Compute the truncated solution f ;n to [START_REF] Makino | Remainder di erential algebras and their applications[END_REF] with f 0 = while t< T do

Let be maximal such that E guess }M 2 p , with M and E guess as in ( 20) and (22) T min( ,T t) T f ;n ( ), t T t + T max{i}d: i } n e } Replace f ;n by the approximate xed point 2n ( f ;n ( +z)) of as in (27) return

Convergence of the -iteration

Let us now study the convergence of as a mapping on the dn-dimensional vector space z] ;n d over . We de ne a norm on this space by g = max

k<n k |g k | = max 1}i}d,k<n k g k i] k = k! k = +1 4
.

For a linear map on this space, represented as a matrix M, we have the corresponding matrix norm M = sup g =1

Mg .

Recall that J and J stand for the Jacobian matrices of and .

THEOREM 5. Assume that g is a xed point of and that A>0 and >n/ are constants such that |(J (g)) k | }A k for all i, j =1,&,d and k =0,&,n 1. Then J (g) } + A ( ( ( ( ( ( ( ( ( (1

n ) ) ) ) ) ) ) ) ) ) 1 .
Proof. Consider an in nitesimal perturbation g of g with (g + g) = (g) + (g). Given i } and 0<k<n, we have k ( (g) ( ( ( ( ( ( ( ( ( ( ( ( y 0} }k 1 ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( y 0} }k 1 ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( y 0} }k 1 ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( y ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( y }k 1 A ( ( ( ( ( ( ( ( ( ( n ( ( ( ( ( ( ( ( ( A ( ( ( ( ( ( ( ( ( (1

) k i] = k k (J (g) g) k 1 i] i ( g) k 1 i] } k k ( (
|(J (g) g) k 1 | + i |( g) k 1 | ) ) ) ) ) ) ) ) ) ) ) ) ) ) } k k ( (
A g k 1 + i g k 1 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) = ( (
A k k k 1 k 1 k 1 + i k k k 1 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) g = ( (
}k 1 A (k 1)! (k 1 )!( ) + i ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) g } (
) ) ) ) ) ) ) ) ) ) + i ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) g } (
n ) ) ) ) ) ) ) ) ) ) 1 + ) ) ) ) ) ) ) ) ) ) g
Similarly, given i > and k <n, we have k ( (g) ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( ( y }k A i ( ( ( ( ( ( ( ( ( ( n ( ( ( ( ( ( ( ( A ( ( ( ( ( ( ( ( ( (1

) k i] = k i (J (g) g) k i] (k+1)( g) k+1 i] } k i ( (
y }k |(J (g) g) k | +(k+1)|( g) k+1 | ) ) ) ) ) ) ) ) ) ) ) ) ) ) } k i (
y }k A g k +(k+1) g k+1 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) = ( (
y }k A i k k + (k+1) k i k+1 ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) g = ( (
y }k A i k! (k )!( ) + i ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) g } ( (
) ) ) ) ) ) ) ) ) ) + i ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) g } ( (
n ) ) ) ) ) ) ) ) ) ) 1 + ) ) ) ) ) ) ) ) ) ) g
Putting both relations together, we obtain (g) } A 1 n 1 + g .

Assume that

+ A ( ( ( ( ( ( ( ( ( (1

n ) ) ) ) ) ) ) ) ) ) 1 <1. ( 29
)
Then the theorem implies the existence a small neighbourhood of the xed point g of on which the -iteration converges to g. Whenever this condition is met, the -iteration actually tends to converge on a rather large neighbourhood that includes our ansatz; see the next subsection for a more detailed discussion. Intuitively speaking, the condition requires the eigenvalues i to be su ciently separated with respect to the norm of the forcing term . Even when the condition does not hold, the -iteration usually still displays an initial convergence for our ansatz, but the quality of the approximate solution to [START_REF] Press | Numerical recipes, the art of scienti c computing[END_REF] ceases to improve after a while.

If the condition ( 29) is satis ed for all critical indices that we encounter when integrating from 0 until time T, then Algorithm 2 should produce an accurate result. The idealized analysis from section 4.2 then also applies, so the algorithm takes O(2 p/n log( d T)) steps. Since each step now requires O (s n 2 ) oating point operations at bit precision p, we nally obtain the bound O (sn 2 p2 p/n log( d T)) for the total running time.

Quality of the computed solution

Let us now investigate in more detail why fixed points f ;n of the -iteration indeed approximate the true solution f ;n quite well. For this, we will determine a more precise small neighbourhood of the xed point g on which the -iteration converges and show that this neighbourhood in particular contains f ;n . We start with two lemmas. Proof. Given v d with |v| =1, we have

|M k v| = 1 2 i 9 |z|= M(z) v z k+1 dz } 1 k sup |z|= |M(z)v| } A k . It follows that |M k | =sup |v| =1 |M k v| }A/ k . LEMMA 7. Let >0, >0, C <1, g z] ;n d , and A = sup |u|} ,| |} 1 C e |J (g(u)+ )| .
Then for all g z] ;n d with g g } 1 C and k , we have |J (g ) k | } A k . Proof. Setting g =g g, we rst notice that

sup |u|} | g(u)| } y k | g k | k } y k || g|| k k = y k ( ) k k! || g|| } 1 C e g .
It follows that sup |u|} |J (g (u))| }A and we conclude using the previous lemma.

THEOREM 8. Let >n/ and A> 0 be such that C= + A ( ( ( ( ( ( ( ( ( (1

n ) ) ) ) ) ) ) ) ) ) 1 <1. Let g z] ;n d and = (g) g be such that |u|} ,| |} 1 C e |J (g(u)+ )| }A.
Then the sequence g, (g), 2 (g),& tends to a unique xed point g on the set , g, 1 C = h z] ;n d : h g } 1 C .

Proof. A straightforward adaptation of the proof of Theorem 5 shows that J (g) } C

on the ball , g, 1 C , which means that (h 1 ) (h 2 ) } C h 1 h 2 on this ball.

By induction on i , it follows that i+1 (g) i (g) } C i and i+1 (g) g }

(1 + C + +C i ). We conclude that g, (g), 2 (g), & converges to a fixed point g +( (g) g) +( (g 2 ) (g))+ of in , g, 1 C .

Returning to the Taylor series expansion f ;n of the exact solution of (1) at time t, we notice that i] ( f

;n ) = { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { ( f ;n ) i] if i } ( f ;n ) i] + n i ( f n ) i] z n 1 if i > . It follows that ( f ;n ) f ;n } n n 1 K t n =K n! ( t) n H K n e n ( ( ( ( ( ( ( t
) ) ) ) ) ) ) n }K ( ( ( ( ( ( ( t ) ) ) ) ) ) ) n .
Now assuming that the aimed step size H t 2 p/n was more or less achieved at the previous step, it follows that T ( f ;n ) f ;n r K 2 p is of the desired order. If the condition (29) is indeed satis ed for =t, we thus should be able to apply Theorem 8 and conclude that the computed xed point is within distance /(1 C) for the norm. Using a similar reasoning, we also see that the ansatz at the next step will be su ciently close to the true solution for the -iteration to converge.

Solving steady-state problems through the rst variation

A more robust but costly approach to solve the steady-state problem ( 26) is to compute the initial values f 0 +1] ,&, f 0 d] with su cient precision using Newton's method. Given a tentative approximation = f 0 for the initial condition, we both compute f ;n+1 =f t | ;n+1 and the rst variation V ;n+1 =V t | ;n+1 , after which we update T + by solving the linear system

f n +V n = 0, 1] = = ] =0.
This method admits quadratic convergence, but it requires us to compute with a precision of n log 2 (e d /n) bits at least in order to be accurate. Indeed, this comes from the fact that |V n | grows roughly as d n /n!. On the upside, we may compute f ;n+1 and V ;n+1 using any of the algorithms from section 3. The total running time is therefore bounded by O (d s n (p+n log d ) 2 p/n log( d T)). Notice also that it actually su ces to compute the last d rows of V n , due to the requirement that 1] = = ] =0.

The main disadvantage of the above method is that the computation of the rst variation V ;n+1 along with f ;n+1 induces an overhead of d, which may be a problem for systems of high dimension. Let us now sketch how one might reduce this overhead by combining the variational approach with the -iteration. We intend to return to a more detailed analysis in a future work.

Instead of using a single critical index , the rst idea is to use a range , &, of indices starting at = , and such that the condition (29) holds for

= +1 4 .
The second idea is to only vary the components +1] ,&, ] of the initial condition and use the steady-state conditions for the indices +1,&,d.

More speci cally, for a tentative initial condition f 0 = , we rst solve the steady-state problem

f z + f = ( f ), { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { f 0 i] = i] if i} f n i] =0 if i>
using the -iteration technique. In a similar way, we next solve the following variational steady-state problem for the d ( ) matrix W:

W z + W =J ( f )W, { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { W 0 i,j] =Id i,j+ ] if i} W n i,j] =0 if i> .
As our ansatz, we may use W i,j] = Id i,j+ ] for all i, j. Having computed f ;n+1 and W ;n+1 at precision n +1, we nally update T +d by solving the linear system

f n i] +W n i,1] ( ) +1] + +W n i, ] ( ) ] = 0, i = +1,&,
and setting 1] = = ] = +1] = = d] = 0. We repeat this whole process until f n i] is su ciently close to zero for i= +1,&, . Remark 9. The algorithms in this section are reminiscent of implicit numerical schemes for the integration of (1). One interesting di erence is that our second optimized method only needs to compute a small part of the full Jacobian matrix. Remark 10. Instead of imposing the exact steady-state conditions f n +1] = = f n d] =0, yet another approach would be to minimize the norm | f n | under the condition that | f 0 | }| | 2 p . This approach admits the advantages that one does not need to know and that it might be applied to more general complex matrices . However, it requires the computation of the full Jacobian matrix.

CERTIFICATION

In the previous two sections, we have described numerical algorithms for integrating (1). An interesting question is how to compute a tight error bound such that the distance between the true and the computed solutions lies within this error bound. Ball arithmetic provides a suitable framework for such computations. This variant of interval arithmetic is well suited for high precision computations with complex numbers and we will brie y recall its basic principles in section 6.1. As a rst application, we show how to make Theorem 4 more e ective.

The certi ed integration of di erential equations that are not sti (i.e. the robust counterpart of section 4) is a classical topic in interval arithmetic [START_REF] Makino | Suppression of the wrapping e ect by Taylor model-based validated integrators[END_REF][START_REF] Moore | Automatic local coordinate transformations to reduce the growth of error bounds in interval computation of solutions to ordinary di erential equation[END_REF][START_REF] Hn | Rigourously computed orbits of dynamical systems without the wrapping e ect[END_REF][START_REF] Butcher | Numerical methods for ordinary di erential equations in the 20th century[END_REF][START_REF] Neumaier | Taylor forms -use and limits[END_REF][START_REF] Lanfort | A computer assisted proof of the Feigenbaum conjectures[END_REF][START_REF] Fischer | Fast on-line integer multiplication[END_REF][START_REF] Moore | Interval Analysis[END_REF][START_REF] Lohner | On the ubiquity of the wrapping e ect in the computation of error bounds[END_REF][START_REF] Van Der Hoeven | Evaluating straight-line programs over balls[END_REF][START_REF] Lohner | Einschlie ung der L sung gew hnlicher Anfangs-und Randwertaufgaben und Anwendugen[END_REF][START_REF] Neumaier | The wrapping e ect, ellipsoid arithmetic, stability and confedence regions[END_REF][START_REF] Makino | Remainder di erential algebras and their applications[END_REF]. For recent algorithms of good complexity, much in the spirit of the present paper, we refer to 11]. Most of the existing algorithms rely on Taylor series expansions as we do, while providing rigorous tail bounds for the truncation error.

The e ective counterpart of Theorem 4 provides suitable tail bounds in the case of sti di erential equations. In sections 6.3 and 6.4, we show how to use this for the certi ed resolution of steady-state problems using the -iteration from section 5.1. Unfortunately, our rst algorithms are rather naive and typically lead to heavily overestimated error bounds. The classical way to reduce this overestimation is to also compute the rst variation and apply the mean value theorem: see section 6.5 for how to do this in our context.

Ball arithmetic

Given a and r ~T {x :x~0}, we write B(a, r) ={z :|z a| }r} for the closed ball with center a and radius r. The set of such balls is denoted by ,( , ) or " . One may lift the ring operations +, , in to balls in " , by setting: ,(a,r) ,(b,s) T ,(a b,r +s) ,(a,r) ,(b,s) T ,(ab,(|a|+r)s+|b|r).

These formulas are simplest so as to satisfy the so called inclusion principle: given {+, , }, u ,(a,r) and v ,(b,s), we have u v ,(a,r) ,(b,s). This arithmetic for computing with balls is called exact ball arithmetic. It extends to other operations that might be de ned on , as long as the ball lifts of operations satisfy the inclusion principle. Any ordinary complex number z can be reinterpreted as a ball ,(z,0) " . Given a ball ,(a, r) " , we also notice that ,(a, r) = |a| + r and ,(a,r) = max(|a| r,0) provide us with reliable upper and lower bounds for |,(a,r)| in ~.

Another interesting operation on balls ,(a,r),,(b,s) " that we will need below is intersection. Assuming that the set intersection I T ,(a,r)),(b,s) is non-empty, we de ne the ball intersection J =,(a,r) ,(b,s) " to be the ball of smallest radius that contains I. In order to determine this ball intersection, we may assume without loss of generality that r~s. If ,(a,r) ,(b,s), then J =,(b,r). Otherwise, let u,v be the two (possibly identical) intersections of the circles ,(a, r) and ,(b, s). If |a (u+v)/2| |a b|, then we still have J =,(b,r). Otherwise, J =,((u+v)/2,|v u|/2). It will also be convenient to extend vector notations to balls. First of all, we identify vectors of balls (,(a 1] , r 1] ), &, ,(a d] , r d] )) ,( , ) d with ball vectors ,(a, r) ,( d , d ). Given z d and v " ( " ) d , we also write z v " if and only if z i] (v " ) i] for i =1,&,d. Similar remarks apply to ball matrices, ball series, etc. Remark 11. For e ective computations, one can only work with approximations of real and complex numbers of nite precision. The IEEE 754 norm provides a standard for real oating point arithmetic based on the concept of correct rounding . It naturally generalizes to oating point numbers with a mantissa of p bits and an exponent of e bits 3]. Let p,e and p,e the sets of real and complex oating point numbers of this type (strictly speaking, one also has { , } p,e for the reliable rounding of over owing operations). One may adapt the above arithmetic on exact balls to oating point balls in ,( p,e , p,e ) while preserving the inclusion principle: it su ces to slightly increase the radii of output balls so as to take care of rounding errors. For precise formulas and interesting variations, we refer to 13]. For simplicity, the sequel will be presented for exact ball arithmetic only, but it is not hard to adapt the results to take into account rounding errors.

Computing bounds on compact half disks

The polynomials i] in (1) can either be represented as linear combinations of monomials, or using a dag (directed acyclic graph). In both cases, the ball arithmetic from the previous subsection allows us to reliably evaluate at balls in " .

For the e ective counterpart of Theorem 4, there is a trade-o between the qualities of the radius R (larger radii being better) and the bound B ( > ) d (smaller values of B being better). For a given B ( > ) d , it is simple to compute the best corresponding R that satis es the condition [START_REF] Gambill | Logarithmic reduction of the wrapping e ect with application to ordinary di erential equations[END_REF], by taking R = min i B/ (,(z, B)) i] . Conversely, for a xed R > , let us recall a classical technique from interval analysis that can be used to compute an almost best corresponding bound B ( > ) d that satis es [START_REF] Gambill | Logarithmic reduction of the wrapping e ect with application to ordinary di erential equations[END_REF].

We simply construct a sequence B 0 }B 1 }B 2 } in ( ~)d by taking B 0 =0 and B k+1 = (,(c,B k )) R for all k. If there exists a nite B for which (8) holds, then the sequence B k usually converges to a minimal such B quite quickly. After say ten iterations, we should therefore have a reasonable approximation. One then slightly in ates the last value by taking B = B 10 + (B 10 B 9 ). If (,(c, B)) R } B, then we have succeeded in nding a suitable B. If not, then we return failed .

Using the above procedure, we may also compute a reasonably large compact half disk H R on which is analytic, together with a bound for | |: we simply perform a dichotomic search for the largest R for which the computation of a corresponding B does not fail. If we also wish the bound B to remain reasonably small, one may simply divide R by two at the end of the search and compute the corresponding B.

Certi ed integration of sti di erential equations

Let us now return to the integration of (1) and let be the exact solution for the initial condition (0)= c d . Let R > and B ( > ) d be computed as above such that (8) holds. Assume that we were able to reliably integrate (1) until a given time t }R/2 and let " = " (t) ( " ) d , so that (t) " . In order to adapt the -iteration to ball arithmetic, we rst deduce from Theorem 4 that f t n i] } B i] /t n for i = + 1, &, d. This provides us with the required steady-state conditions f n +1] = ,(0, B +1] / t n ), &, f n d] = ,(0, B d] /t n ), in addition to the initial

conditions f 0 1] = 1] ,&, f 0 ] = ] .
The ball counterpart of our steady-state problem thus becomes

f " z + f " = ( f " ), { { { { { { { { { { { { { { { { { { { { { { { { { { { { ( f 0 " ) i] =( " ) i] if i} ( f n " ) i] =,(0,B i] /t n ) if i> (30)
Let cen and cen denote the centers of " and " . Starting from a numerical approximation f H t| cen , cen ;n of f t| cen , cen ;n , we wish to compute the truncation f ;n " " z] ;n d of a solution to (33), with the property that f t| , ;n f ;n " for all " and " . In order to do this, the idea is again to use a suitable ball version of , , given by

" , " i] (g " ) = { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { ( " ) i] ++( i] (g " ) i (g " ) i] ) if i} i 1 ( i] (g " ) (g " ) i] n ( " ) i] z n 1 ) if i> mod z n ,
and to keep applying " , " on the ansatz f H t | cen , cen ;n until we are su ciently close to a xed point. Using a similar in ation technique as in section 6.2, we nally obtain a truncated series f ;n " with ( f ;n " ) f ;n " , or failed . Now for any " and " , we notice that this also yields , ( f ;n " ) f ;n " . Thinking of f ;n " as a compact set of series in z] ;n d , this means in particular that , admits a xed point g = f t | , ;n f ;n " , as desired.

Curbing the wrapping e ect

If = d, then the algorithm from subsection 6.3 specializes to a classical way to compute the enclosure (31) of a solution to (1) at time t + as a function of the enclosure " = ,( cen , rad ) = f 0 at time t. However, it is well known that this method su ers from a large overestimation of the errors, due to the so-called wrapping e ect . A well know technique to reduce this overestimation is to rst compute a certi ed solution for the initial condition ,( cen , 0) instead of " and then to bound the error due to this replacement, by investigating the rst variation. Let us now show how to extend this technique to general critical indices .

With cen in the role of " and ,(0,B/t n ) in the role of " , we start with the computation of a truncated certi ed solution f ;n acc," " z] ;n d to (33). Recall that this ball solution has the property that f t | cen , ;n f ;n acc," for all " . Since the initial conditions are exact, this solution is generally fairly accurate. Now for any 2 " , we may write f t | cen , = f t| cen , +5 0 1 f t | cen +( 2 cen ) , d .

(34) We next observe that f t | cen +( 2 cen ) , = V t | cen +( 2 cen ) , ( 2 cen ), where V t| , denotes the exact solution of the following equation for the rst steadystate variation:

V z + V =J ( f t| , )V, { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { V 0 i,j] =Id d i,j] if i } V n i,j] =0 if i > .
This equation again admits a ball analogue

V " z + V " =J ( f " )V " , { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { (V 0 " ) i,j] =Id d i,j] if i} (V n " ) i,j] =0 if i> ,
where f " is a solution to (30). We only compute a crude solution V cru," to this equation, for a crude solution f cru," to (30), both truncated to the order n. These solutions can be obtained using the techniques from section 6.3. From (34) we nally conclude that f ;n " = f ;n acc," +V ;n cru," ,(0, rad ) is a ball enclosure for f ;n . This enclosure is usually of a much better quality than f cru," . However, its computational cost is roughly d times higher, due to the fact that we need to determine the rst variation.

Handling close eigenvalues

In the case when +1 / is too small for (29) to hold, we may generalize the above strategy and incorporate some of the ideas from section 5.4. With = and from there, let us brie y outline how to do this. We proceed in four stages:

" We rst compute a truncated numeric solution f ;n H z] ;n d to [START_REF] Press | Numerical recipes, the art of scienti c computing[END_REF], using the technique from section 5.4.

" Let acc," ( " ) d be such that ( acc," ) i] =( f 0 H ) i] for i} and ( acc," ) i] =0 for i> .

Similarly, let acc," ( " ) d be such that ( acc," ) i] = 0 for i } and ( acc," ) i] = ,(0,B i] /t n ) for i > . We use f ;n H as an ansatz for its deformation into a reliable truncated solution f ;n acc," " z] ;n d to (33), but with in the role of , with acc," in the role of " , and acc," in the role of " .

" We further deform f ;n acc," into a reliable truncated solution f ;n aux," " z] ;n d to (33), this time with in the role of . We do this by computing f ;n cru," as above, together with a crude but reliable solution to the equation

W " z + W " =J ( f cru," )W " , { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { { (W 0 " ) i,j] =0 if i} (W n " ) i,j] =Id d i,j] if <i} (W n " ) i,j] =0 if i> .
We then take f ;n aux," = f ;n acc," +W ;n " ,(0, aux ), where ( aux ) i] =B i] /t n for <i} and ( aux ) i] =0 for i } and i > .

" We nally compute the desired truncated solution f ;n " " z] ;n d to (30) as f ;n " = f ;n acc," +V ;n cru," ,(0, rad ), where V ;n cru," is computed as above.

CONCLUSION

There are several directions for generalizations and further improvements of the results in this paper. For simplicity, we assumed to be a diagonal matrix with positive eigenvalues. It should not be hard to adapt the results to arbitrary matrices with positive real eigenvalues only (as long as the forcing term does not explode for the change of coordinates that puts in Jordan normal form).

A more interesting generalization would be to consider complex eigenvalues i with strictly positive real parts. In that case, the half disks H R would need to be replaced by smaller compact sectors of the form S R, = {z : |z|}R, |arg z| } } with < /2. Even more generally, one may investigate how to develop accurate Taylor series methods for integrating arbitrary di erential equations with the property that the step size is proportional to the radius of convergence of the exact solution; see also Remark 10. In this paper, we focussed on high precision computations using high order Taylor schemes. Another interesting question is whether it is possible to develop analogues of our methods in the same spirit as more traditional Runge Kutta schemes. How would such analogues compare to implicit integration schemes?

Concerning certi ed integration, the theory of Taylor models 18, 23, 19] allows for higher order expansions of the ow f t| as a function of the initial condition (that is, beyond the computation of the rst variation as in this paper). We think that it should be relatively straightforward to adapt our techniques to such higher order expansions. In a similar vein, we expect that the techniques from 7, 11] to further curb the wrapping e ect can be adapted to the sti setting.

From a more practical point of view, it would nally be interesting to have more and better machine implementations. For the moment, we only did a toy implementation of the main algorithm from section 5.1 in the MATHEMAGIX system 12]. We found this implementation to work as predicted on various simple examples. Detailed machine experiments with larger systems and the algorithms from sections 5.4 and 6 should make it possible to further improve the new techniques.

  admits a unique analytic solution on the open disk with center 0 and radius 0)= = 1] (0)=m c constitutes a majorant equation for (4) in the sense of Cauchy-Kovalevskaya. By symmetry, each of the components i] of this equation satis es the simpler equation

  m

LEMMA 6 .

 6 Let M(z) = k M k z k be a d d matrix of analytic functions de ned on the closed ball ,(0, ) with center zero and radius >0. Let A=sup |z|} |M(z)| . Then |M k | }A/ k for all k .

  dominates the actual value | f n | H M ( ) n , which yields E H M ( d t ) n 2 p /n! instead of E H M ( / ) n 2 p . When chosing our step size such that EH M 2 p , as in Algorithm 1, this yields ( d ) n H n! and

	H	n! n ( d	H	n e d	,

We correspondingly de ne the ball version of for g " z] ;n d by i] (g

i] ++( i] (g " ) i (g " ) i] ) if i} i 1 ( i] (g " ) (g " ) i] ,(0,nB d] /t n )z n 1 ) if i> mod z n . This map has the property that f ;n g " f ;n (g " ). In our ball context, we may actually iterate using the following improved version of :

(g " ) = g " (g " ). Here we understand that the intersections are taken coefficientwise. Since iterations using can only improve our enclosures, we do not need to worry much about the ansatz. We may deduce a reasonably good ansatz from the power series expansion f " t2 ;n at the previous time t2 < t, and the Cauchy bounds f " t2 k i] } B i] / (t2) k for all k ~n. But it is perfectly reasonable as well to use f ;n an," = " +,(0,B/t)z+ +,(0,B/t n 1 )z n 1 . Applying a su cient number of times to this ansatz, we obtain the desired ball enclosure f ;n " = 2n ( f ;n an," ) of the truncated power series expansion of at time t. Assuming that <t, we may then deduce the enclosure +,( ( ( ( ( ( (0,B t t ( ( ( ( (

of at time t + . Notice that the cost of the computation of certi ed solutions in this way is similar to the cost of Algorithm 1, up to a constant factor.

Certi cation of approximate solutions

The above method relies on the property that the steady-state problem (30) for balls specializes into a numerical steady-state problem that admits f as a solution, since f 0 " and f n ,(0,B/t n ). Given a numerical approximation of f ;n , as computed in section 5, a related problem is whether we can certify the existence of an actual solution in a small neighbourhood of the approximate solution.

In order to make this more precise, let us introduce a few more notations. Given 1] ,&, ] and +1] ,&, d] , consider the numeric steady-state problem

We will denote by f t | , any exact solution to this problem if such solution exists. It will be convenient to regard and as elements of d , through padding with zeros. Now consider the operator , : z] ;n d z] ;n d de ned by

Through the iterated application of , to a suitable ansatz, we may obtain a numerical approximation f H t| , ;n of f t| , ;n . Now consider balls ( " ) 1] ,&,( " ) ] " and ( " ) +1] ,&,( " ) d] " , again padded with zeros. Then we have the following ball analogue of (32):