L Bouchet
email: lbouchet@irap.omp.eu

P Amestoy

A Buttari

F.-H Rouet

M Chauvin

Simultaneous analysis of large INTEGRAL/SPI 1 datasets: Optimizing the computation of the solution and its variance using sparse matrix algorithms

INTEGRAL/SPI X/γ -ray spectrometer data analysis.

• Large astronomical datasets arising from the simultaneous analysis of years of data.

• Resolution of a large sparse system of equations; the solution and its variance.

• The Multifrontal Massively Parallel Solver (MUMPS) to solve the equations.

• MUMPS A -1 feature to compute selected inverse entries (variance of the solution etc.).

a b s t r a c t Nowadays, analyzing and reducing the ever larger astronomical datasets is becoming a crucial challenge, especially for long cumulated observation times. The INTEGRAL/SPI X/γ -ray spectrometer is an instrument for which it is essential to process many exposures at the same time in order to increase the low signalto-noise ratio of the weakest sources. In this context, the conventional methods for data reduction are inefficient and sometimes not feasible at all. Processing several years of data simultaneously requires computing not only the solution of a large system of equations, but also the associated uncertainties. We aim at reducing the computation time and the memory usage. Since the SPI transfer function is sparse, we have used some popular methods for the solution of large sparse linear systems; we briefly review these methods. We use the Multifrontal Massively Parallel Solver (MUMPS) to compute the solution of the system of equations. We also need to compute the variance of the solution, which amounts to computing selected entries of the inverse of the sparse matrix corresponding to our linear system. This can be achieved through one of the latest features of the MUMPS software that has been partly motivated by this work. In this paper we provide a brief presentation of this feature and evaluate its effectiveness on astrophysical problems requiring the processing of large datasets simultaneously, such as the study of the entire emission of the Galaxy. We used these algorithms to solve the large sparse systems arising from SPI data processing and to obtain both their solutions and the associated variances. In conclusion, thanks to these newly developed tools, processing large datasets arising from SPI is now feasible with both a reasonable execution time and a low memory usage.

Introduction

Astronomy is increasingly becoming a computationally intensive field due to the ever larger datasets delivered by observational efforts to map ever larger volumes and provide ever finer details of the Universe. In consequence, conventional methods are often inadequate, requiring the development of new data reduction techniques. The SPI X/γ -ray spectrometer, aboard the INTEGRAL observatory, perfectly illustrates this trend. The telescope is dedicated to the analysis of both point-sources and diffuse emissions, with a high energy resolution (Vedrenne et al., 2003). Its imaging capabilities rely on a coded-mask aperture and a specific observation strategy based on a dithering procedure (Jensen et al., 2003). After several years of operation, it also becomes important to be able to handle simultaneously all the data, in order, for example, to get a global view of the Galaxy emission and to determine the contribution of the various emission components.

The sky imaging with SPI is not direct. The standard data analysis consists in adjusting a model of the sky and instrumental background to the data through a chi-square function minimization or a likelihood function maximization. The related system of equations is then solved for the intensities of both sources and background. The corresponding sky images are very incomplete and contain only the intensities of some selected sky sources but not the intensities in all the pixels of the image. Hence, images obtained by processing small subsets of data simultaneously cannot always be combined together (co-added) to produce a single image. Instead, in order to retrieve the low signal-to-noise ratio sources or to map the low surface brightness ''diffuse'' emissions (Bouchet et al., 2011), one has to process simultaneously several years of data and consequently to solve a system of equations of large dimension. Grouping all the data containing a signal related to a given source of the sky allows to maximize the amount of information on this specific source and to enhance the contrast between the sky and the background.

Ideally, the system of equations that connects the data to the sky model (where the unknown parameters are the pixels intensities) should be solved for both source intensities and variability timescales. This problem, along with the description and treatment of sources variability, is the subject of another paper (Bouchet et al., in press).

It is mandatory, for example when studying large-scale and weak structures in the sky, to be able to process large amounts of data simultaneously. The spatial (position) and temporal (variability) description of sources leads to the determination of several tens of thousands of parameters, if ∼6 years of SPI data are processed at the same time. Consequently, without any optimization, the systems to be solved can exceed rapidly the capacities of most conventional machines. In this paper we describe a technique for handling such large datasets.

Material and methods

The SPI spectrometer

SPI is a spectrometer provided with an imaging system sensitive both to point-sources and extended source/diffuse emission. The instrument characteristics and performances can be found in Vedrenne et al. (2003) and [START_REF] Roques | [END_REF]. Data are collected thanks to 19 high purity Ge detectors illuminated by the sky through a coded-mask. The resulting Field-of-View (FoV) is ∼30°and the energy ranges from 20 keV to 8 MeV. The instrument can locate intense sources with an accuracy of a few arc minutes [START_REF] Dubath | [END_REF].

Functioning of the ''spectro-imager'' SPI

The coded mask consists of elements which are opaque (made of tungsten) or transparent to the radiation. Photons coming from a certain direction cast a shadow of the mask onto the detectors plane. The shadowgram depends on the direction of the source (Fig. 1). The recorded counts rate in each detector of the camera is the sum of the contribution from all the sources in the FoV. The deconvolution consists of solving a system of equation which relates a sky model to the data through a transfer function. In the case of SPI, the imaging properties rely on the coded aperture, but also on a specific observing strategy: the dithering.

Dithering and sources variability

The reconstruction of all the pixels of the sky image enclosed in the FoV is not possible from a single exposure. Indeed, dividing the sky into ∼ 2°pixels (the angular resolution), we obtain, for a 30°F oV, ∼(30°/2°) 2 = 225 unknowns. However, a single exposure contains only 19 data values which are the number of observed counts in the 19 Ge detectors and does not permit us to obtain the parameters necessary to determine the model of the sky and background. The related system of equations is thus undetermined. The dithering observation technique aims to overcome this difficulty. By introducing multiple exposures for a given field that are shifted by an offset that is small compared to the size of the FoV, it is possible to increase the number of equations, by grouping exposures, until the system becomes determined and thus solvable. An appropriate dithering strategy (Jensen et al., 2003) has been used where the spacecraft continuously follows a dithering pattern throughout an observation. In general, the pointing direction varies around a target by steps of 2°within a five-by-five square or a seven-point hexagonal pattern. A pointing (exposure) lasts between 30 and 60 min. Thus, the dithering allows to construct a solvable system of equations.

However, in addition to the variable instrumental background, sources are also variable on various timescales ranging from hours (roughly the duration of an exposure) to years. This is not a major problem at high energy (E 100 keV), since there are only few emitting sources, whose intensities are rather stable in time with respect to the statistics. At lower energies (E 100 keV) and in most cases, the source intensities vary during the time spanned by the all the exposures. The chi-square, of the associated least-square problem, for this group can be relatively high, if sources intensity variations are not taken into account. In spite of this, it is possible to include a model of the source intensity variations in the formulation of the problem and to re-optimize the system of equations accordingly (Bouchet et al., in press). Nevertheless, including sources variability in the system of equations increases the number of unknowns to determine (2.2.3) since intensities, in each ''time-bin'' (a segment of time where the intensity of a given source does not change statistically), are to be determined simultaneously along with the parameters which model the instrumental background.

Cases where it is better to process large amount of data simultaneously

It is impossible from a single exposure (19 data values) to obtain the sky image in the 30°FoV; only a coarse image containing at most 19 sources can be obtained. This coarse image is undersampled and contains information on only 19 pixels (there is no information on the other pixels). Hence, images cannot always be combined together (co-added) to produce a single image. Furthermore, crowded regions like the Galactic Center contain hundreds of sources and thus a single exposure cannot provide the amount of information needed, even to build only a coarse image. The grouping of the exposures, by selecting all those that contain signal on a specific sky target, can provide the necessary information. The FoV spanned by these exposures is large (30°-60°) and contains numerous sources. Finally, as the ''diffuse'' emission study involves the entire sky, it is necessary to process large amounts of data together to enhance the contrast between the sky and the instrumental background, which implies to solve a large system of equation.

Problem formulation

The signal (counts and energies) recorded by the SPI camera on the 19 Ge detectors is composed of contributions from each source in the FoV plus the background. For n s sources located in the field of view, the data D raw dp obtained from detector d during an exposure (pointing) p, for a given energy band, can be expressed by the relation:

D raw dp = n s j=1 R dp,j I s p,j + B bg dp + ǫ dp (1)
where R dp,j is the response of the instrument for source j (function of the source direction relative to the telescope pointing axis 1)). To reduce the number of free parameters related to background, we take advantage of the stability of relative counts rate between detectors and rewrite the background term as:

B bg dp = I bg p × U d × t dp (2)
where I bg p is a normalization coefficient per pointing related to the background intensity, U d is a background count rate pattern (uniformity map) on the SPI camera for detector d, and t dp the effective observation time for pointing p and detector d. The number of parameters necessary to model the background reduces to n p if U is assumed to be known. 3 However, in some cases it can be determined while processing the data (Appendix A.4).

The two extreme cases, in terms of number of parameters to be determined, are

• First, when the sources and background intensities are assumed to be constant throughout all the observation (time spanned by the exposures), the relation between the data and the sky model can be written, omitting the detector indices, as 2 The number of functioning detectors could be n d = 16, 17, 18 or 19 for single events and up to 141, when all the multiple events are used in addition to the single events [START_REF] Roques | [END_REF].

3 Derived from ''empty-field'' observations (Bouchet et al., 2010). (3)

The aim is to compute the intensities I s (j) of the n s sources and the background relative intensity I bg . Therefore, the relation can be written in matrix form, as

     P 1 R 1,1 . . . R
            I bg I s 1 I s 2 . . . I s n s        =         D raw 1 D raw 2 . . . D raw n p -1 D raw n p         +       ǫ 1 ǫ 2 . . . ǫ n p -1 ǫ n p      
.

We can rewrite the system in a more compact form as

y = H 0 x + ǫ or y i = n s +1 j=1 h ij x j + ǫ i for i = 1, . . . , M (4)
where • Second, if the background or the sources are variable on the exposure timescale, the number of unknowns (free parameters) of the set of n p × n d equations is then (n s + 1) × n p (for the n s sources and the background intensities, namely I s and I bg ). This leads, unless the number of sources is small, to an underdetermined system of equations. 4Fortunately, in real applications, many sources vary on timescales larger than the exposure. This leads to a further reduction of the number of parameters compared to the case where all sources vary on the exposure timescale. In addition, many point sources are weak enough to be considered as having constant flux within the statistical errors, especially for higher energies (E 100 keV).

H 0 (elements h ij ) is an M × (n
Then the n p × n s parameters related to sources will reduce into N eff s parameters and, similarly, N b for the background. As these parameters have also a temporal connotation, they will hereafter be referred to as ''time-bins''.

If the source named or numbered J is variable, then the total duration covered by the n p exposures is divided into K J sub-intervals where the source intensity can be considered as stable/constant regarding the data statistics. The solution x J is expanded in K J segments, it takes the value ''time-bins'' s J k in segment k, and can be written in compact notation

x j = K J k=1 s J k I J k with I J k = 1 if t ∈ [t J k-1 , t J k [ I J k = 0 otherwise.
Actually the instants t J k correspond to the exposure acquisition time (exposure number), with t 0 = 1 and t J k = n p + 1. There is at least one and at most n p time segments for each source J

(x J = [s J 1 , . . ., s J K J ] becoming a vector of length K J ).
The matrix H 0 (Eq. ( 4)) is to be modified accordingly.

When expanding matrix H 0 , column J is expanded in K J new columns, hence the number of intensities (unknowns) increases.

H 0 =         h 11 h 12 h 13 . . . h 1N h 21 h 22 h 23 . . . h 2N h 31 h 32 h 33 . . . h 3N . . . . . . . . . . . . . . . h M1 h M2 h M3 . . . h MN         -→ H =       h 11 0 0 0 h 12 0 0 h 13 0 h 1N . . . 0 0 h 21 0 0 h 22 0 0 0 h 23 h 2N . . . 0 0 0 h 31 0 0 h 32 0 0 h 33 0 . . . h 3N . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 h M1 0 0 h M2 0 h M3 0 . . . h MN       Box I. Schematically H 0 (M × (n s + 1)) is mapped into a matrix H (M × N), N being the sum of all sources intervals (N = n s J=0 K J ),
that is the number of ''time-bins'' (the index J = 0 correspond to the background). Matrix H(1 : M, 1 : K 0 ) is related to the background while H(1 : M, K 0 + 1 : N) is related to the sources response. Parameters x(1 : K 0 ) and x(K 0 + 1 : N) are related to background and source intensity variations with the time (number of exposures). Box I illustrates schematically how the matrix H is derived from the matrix H 0 .

Finally, the relation between the data and the sky model, similarly as in Eq. ( 4), is

H x = y + ǫ. (5)
Physically, H corresponds to the transfer function or matrix, y to the data and x to the unknown intensities (sources plus background) to be determined (a vector of length N).

Taking into account the variability of sources and instrumental background increases the size of the system of equation and the number of unknowns, but also increases the sparsity of the matrix related to the system of equations, which means that the underlying matrices have very few non-zero entries. In our application, the matrix H 0 is sparse, thus matrix H is even sparser. Objective methods to construct the matrix H from H 0 are described in Bouchet et al. (in press).

To give an idea, for the dataset which corresponds to (∼6 years of data, the number of free parameters N = N eff s + N b to be determined are between N ∼ 5000 and N ∼ 90 000 depending on the energy band considered and hypotheses made on sources and background variability timescale (Section 2.3).

Material

The material is related to the analysis of data accumulated between 2003 and 2009 with the spectrometer SPI. The astrophysical application is the study of diffuse emission of the Galaxy. The details can be found in Bouchet et al. (2011). The goal is to disentangle the ''diffuse'' emission (modeled with 3 spatially extended components) from the point-sources emission and instrumental background. We need to sample these large-scale structures efficiently over the entire sky and consequently use the maximum amount of data simultaneously, since a single exposure covers only one-hundredth of the total sky area. The datasets consist of 38 699 exposures that yield M = 649 992 data points. In most cases considered here, the background intensity is considered to be quite stable on a ∼6 h timescale, which corresponds to N b ≃ 5870 unknowns. ''time-bins'' intensity are to be determined.

In some configurations, essentially used to assess the results, background intensity and/or strongest variable sources vary on the exposure timescale, and the number of unknowns could be as high as N ≃ 55 000 to N ≃ 90 000. Nevertheless, the matrices associated with these problems remain relatively structured. (c) To avoid excessively structured matrices, we generate also matrices H, with a variable number of columns, the number of segments K J for a particular source being a random number between 1 and n p . This results in a different number of parameters N.

Another astrophysical application is the study of a particular source or sky region, here the crowded central region of the Galaxy. In this case, it is possible to use a smaller number of exposures. We use 7147 exposures which cover a sky region of radius 30°around the Galactic center. We measure the intensity variations of a set of 132 sources. The number of parameters to determine N = 3578 is relatively small. Details can be found in Bouchet et al. (in press).

A second matrix, used for verification purposes, has N = 9437. It corresponds to the case where some sources are forced to vary on shorter timescales. The material consists of rectangular matrices H and symmetric square matrices A (A = H T H) related to the above physical problems (2.2.3). The characteristics of some of these matrices are described in Table 1.

The system we use in the experiments consists of an Intel i7-3517U processor with 8 GB main memory. We ran the experiments on a single core, although our algorithms are amenable to parallelism.

Methods

The mathematical problem described in Section 2.2.3 and developed in 2.4.1 requires the solution of several algebraic equations. First, if the chi-square statistic is used, a linear leastsquares problem has to be solved to estimate the parameters of the model. Second, elements (entries) of the inverse of a matrix have to be computed in order to determine the error bars (variances of these parameters). Third, in some cases, the parameters are adjusted to the data through a multi-component algorithm based on likelihood tests (Poisson statistics); this problem leads to a nonlinear system of equations (Appendix A.1).

These three problems can be reduced to solving a linear system with a square matrix: a linear least-squares problem min x Hxy can be transformed into a square system Ax = b by use of the normal equations5 (A = H T H and b = H T y). Similarly, computing entries of the inverse of a matrix amounts to solving many linear systems, as described in detail in Section 3.3.1. For the above mentioned non-linear problem, we chose a Newton-type method; this involves solving several linear systems as well. Our problems are large, but sparse (cf. Table 1), which justifies the use of sparse linear algebra techniques. In Section 3.1, we describe how we selected a method suitable for our application.

The least-square solution (LSQ)

The system is, in most cases, overdetermined (there are more equations -or measures here -than unknowns), therefore there is (generally) no exact solution, but a ''best'' solution, motivated by statistical reason, obtained by minimizing the following merit function, which is the chi-square6 :

χ 2 = M i=1      y i - N j=1 h ij x j σ i      2 (6)
y = (y 1 , . . . , y M ) is vector of length M representing the data, [Σ ] a diagonal matrix of order M whose diagonal is (σ 1 , . . . , σ M ), where σ i is the measurement error (standard deviation) corresponding to the data point y i . These quantities are assumed to be known (formally

σ i = √ y i ). H = h ij is a matrix of size M × N. The least-square solution x = (x 1 , . . . , x N ) is obtained
by solving the following normal equation:

(H T Σ -2 H)x = H T [Σ -2 ]y or as Ax = b. (7)
Once the solution has been computed, the uncertainties on the estimated solution x are needed as well. The corresponding variance can be obtained by computing the diagonal of A -1 :

Var(x i ) ∝ a -1 i,i
where a -1 i,j refers to (A -1 ) i,j .

(8)

Theory

Processing large datasets: efficient solution of large sparse systems of equations

Sparse matrices appear in numerous industrial applications (mechanics, fluid dynamics,. . . ), and the solution of sparse linear systems has been an active field of research since the 1960s. Many challenges still arise nowadays, because industrial applications involve larger and larger number of unknowns (up to a few billions nowadays), and because hardware architectures are increasingly complex (multi-core, multi-GPU, etc.).

Exploiting sparsity can significantly reduce the number of operations and the amount of memory needed to solve a linear system. Let us take the example of a partial differential equation to be solved on a 2D physical domain; the domain can be discretized on a k × k 2D grid and using, say, finite differences, the equation can be transformed into a sparse linear system with N = k × k unknowns. Without exploiting sparsity, this system would be solved in O(N 3 ) operations (using an exact method), with a memory usage in O(N 2 ). It has been shown that, for this particular case, the number of arithmetic operations can be reduced to O(N 3/2 ), and space complexity to O(N log N) by exploiting the sparsity of the matrix [START_REF] Hoffman | Complexity bounds for regular finite difference and finite element grids[END_REF].

Many methods exist for solving sparse linear systems [START_REF] Duff | Direct Methods for Sparse Matrices[END_REF][START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. Two main classes can be distinguished: direct methods, that rely on a matrix factorization (e.g., A = L U), and iterative methods, that build a sequence of iterates that hopefully converges to the solution. Direct methods are known to be numerically robust but often have large memory and computational requirements, while iterative methods can be less memory-demanding and often faster but are less robust in general. Iterative methods often need to be preconditioned, i.e., to be applied to a modified system M -1 Ax = M -1 b for which the method will converge more easily; a trade-off has to be found between the cost of computing and using the preconditioner M and how the preconditioner improves the convergence. The choice of a method is often complicated and strongly depends on the application. In our case, we choose to use a direct method for the following reasons:

• Memory usage is often a bottleneck that prevents the use of direct methods, but with the matrices arising from our application, direct and iterative methods have roughly the same memory footprint. This is explained in the next section.

• The matrices from our application are numerically challenging; we found that unpreconditioned iterative methods (we tried GMRES) have difficulties converging and that a direct method that does not implement robust numerical features is also likely to fail (we illustrate this in Section 5).

• We need to compute error bars, which amounts to solving a large number (O(N)) of linear systems with different right-hand sides but the same matrix. This is particularly suitable for direct methods; indeed, once the matrix of the system is factored (e.g., A = L U), the factors can be reused to solve for different right-hand sides at a relatively inexpensive cost. We describe this in Section 3.3.1.

In this work, we use the MUMPS (Multifrontal Massively Parallel Solver) package. MUMPS [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic scheduling[END_REF][START_REF] Amestoy | Hybrid scheduling for the parallel solution of linear systems[END_REF] aims at solving large problems on parallel architectures. It is known to be very robust numerically, by offering a large variety of numerical processing operations, and provides a large panel of features. In the following section, we briefly describe how sparse direct methods work. We introduce the basic material needed to understand the algorithm used for the computation of error bars (described in Section 3.3.1). D is a diagonal matrix if A is symmetric (which is the case in our application). Sparse direct methods depend on the non-zero pattern of the matrix and are optimized in that sense; specialized mathematical libraries for tridiagonal, banded, cyclic matrices are common. If the pattern is more complex, then the method usually consists of three phases: analysis, factorization and solution.

Sparse direct methods

Direct

Analysis

The analysis phase applies numerical and structural preprocessing to the matrix, in order to optimize the subsequent phases. One of the main preprocessing operations, called reordering, aims at reducing the fill-in, i.e., the number of non-zero elements which appear in the factors but do not exist in the initial matrix; this step consists in permuting the rows and columns of the initial matrix in such a way that less fill-in will occur in the permuted matrix. Table 2 shows the amount of fill-in for different problems coming from our astrophysical application when the matrices are permuted using the nested-dissection method. For each matrix, the number of non-zero elements in the original matrix A and in the L factor of the L D L T factorization of A are reported. Note that in our application, the fill-in is not very large: the number of nonzero elements in the factors is of the same order of magnitude as in the original matrix. As a result, the use of sparse, direct methods is likely to provide a good scalability with respect to the size of the matrix produced by the application. Moreover, this implies that, for our application, direct and iterative methods will have roughly the same memory requirements; indeed, in an unpreconditioned iterative method, the memory footprint is mainly due to the storage of the matrix A, while the major part of memory requirements of direct methods comes from the factors. Note that, while our application exhibit low amount of fill-in, this not the case in other applications; in many problems, especially those involving PDEs on 3D physical domains, the number of non-zero coefficients in the factors can be as big as one hundred times more than in the original matrix. In this case, using an iterative method can be interesting memory-wise.

An important step of the analysis phase is the symbolic factorization: this operation computes the non-zero pattern of the factors, on which the numerical factorization and the solution will rely. The symbolic factorization computes the structure of the factors by manipulating graphs, and also a structure called the elimination tree, a tree-shaped graph with N vertices. This tree represents tasks dependencies for the factorization and the solution phases. We describe in more details the elimination tree since it is a key structure to explain and understand (see Section 3.3.1) how to accelerate the solution phase since computing entries in the inverse of the matrix corresponds to incomplete traversals of the elimination tree. Fig. 2(b) shows an elimination tree and we use it to illustrate some definitions: one of the nodes is designated to be the root; in the figure, this is node 6. For our purpose, the root is the node corresponding to the variable of the linear system that is eliminated last.

An ancestor of a vertex v is a vertex on the path from v to the root. The parent (or father) of v is its first ancestor; all the nodes but the root have a parent. For example, on the figure, nodes 6 and 5 are ancestors of 4; 5 is the parent of 4. A child of a vertex v is a vertex of which v is the parent. For example, 4 and 3 are the children of 5. A vertex without children is called a leaf ; 1 and 2 are leaves. Descendants of a vertex v are all the nodes in the subtree rooted at v; for example, 1, 2, 3 and 4 are descendants of 5.

In the following subsections (factorization and solution phase), we describe briefly how a sparse direct solver uses elimination trees; we will also rely on this notion in Section 3.3.1 for the computation of error bars. Further details about the construction and the role of elimination trees in sparse solvers are given in [START_REF] Liu | The role of elimination trees in sparse factorization[END_REF].

Factorization

After the preprocessing performed during the analysis phase, the numerical factorization takes place and the matrix A is transformed into a product of factors (e.g., L U). The factorization consists in traversing the elimination tree following a postorder, that is a topological ordering (i.e. each parent is visited after its children) where the nodes in each subtree are visited consecutively. In Fig. 2(b), 1-4-2-3-5-6 is, for example, a postorder. At each node, a partial factorization of a dense matrix is performed. Note that nodes that belong to different branches can be processed independently, which is especially useful in a parallel setting.

The factorization phase tries to follow as much as possible the preparation from the analysis phase, but sometimes, because of numerical issues (typically, division by a ''bad pivot'', i.e. a very small diagonal entry that could imply round-off errors), it has to adapt dynamically: the structure of the factors and the scheduling of the tasks can be modified on the fly.

Solution phase

Once the matrix has been factored, the linear system is solved.

For example, in the case of the L U factorization, the system Ax = b becomes LUx = b and is solved in two steps (two solutions of triangular systems):

z = L -1 b ''Forward substitution'' x = U -1 z ''Backward substitution''.
The forward substitution follows a bottom-up traversal of the elimination tree as in the factorization, while the backward substitution follows a top-down traversal of the tree. At each node, one component of the solution is computed, and some updates are performed on the dependent variables (ancestor nodes for the forward phase, descendant nodes for the backward phase).

Computing error bars: partial computation of the inverse of a sparse matrix

In our astrophysical application, once the solution, either for the linear or non-linear problem, has been found, it is necessary to compute the variances of the parameters of the fitted function. In the case of multiple regressions such as least squares problems, the standard deviation of the solution can be obtained by inverting the Hessian or covariance matrix. However, since the inverse of a sparse matrix is structurally full, it is impractical to compute or store it [START_REF] Duff | Sparsity structure and Gaussian elimination[END_REF]. In our case, the whole inverse of the covariance matrix is not required: since we only want the variances of the parameters (not their covariances), we only need to compute the diagonal elements of the inverse.

Some work has been done since the 1970s in order to compute a subset of elements of the inverse of a sparse matrix. One of the first works is [START_REF] Takahashi | Formation of a sparse bus impedance matrix and its application to short circuit study[END_REF] which has been extended in [START_REF] Campbell | Computing the sparse inverse subset: an inverse multifrontal approach[END_REF]; this approach relies on a direct method (i.e. on a factorization). An iterative method has been proposed in [START_REF] Tang | A probing method for computing the diagonal of the matrix inverse[END_REF] for matrices with a decay property. Some methods have also been developed for matrices arising from specific applications; a more detailed survey is given in [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF]. Many of these methods provide sophisticated ideas and interesting performance on specific problems, but no software package is publicly available, with the exception of the approach implemented within MUMPS solver, that we describe in the next section.

MUMPS A -1 feature

The A -1 feature in MUMPS has been described in Slavova ( 2009) and was motivated by the INTEGRAL/SPI application, among other applications that require the computation of inverse entries, or, more generally, applications that involve sparse right-hand sides (as explained in this section). This feature is able to compute any set of entries of A -1 , relying on a traditional solution phase, i.e. by computing every required entry a -1 ij as:

a -1 ij = A -1 e Ď j i
.

Using the LU factors of A, this amounts to solving two triangular systems:

Lz = e Ď j a -1 ij = U -1 z i .
The first triangular system in the equation above is particular because its right-hand side e Ď j is very sparse (only one non-zero entry). Furthermore, we do not need the whole solution of the second triangular system, but only one component. This information can be exploited to reduce the traversal of the elimination tree; while a regular solution phase would consist in visiting the whole elimination tree twice (a bottom-up traversal followed by a topdown traversal), computing a -1 ij consists in two partial traversals of the tree: the first triangular system is solved by following the path from node j to the root node, and the second triangular system is solved by following the path from the root node to node i; this is referred to as pruning the elimination tree. Since each node of the tree corresponds to operations to be performed (arithmetic operations, or expensive accesses to the factors in the out-of-core case), this leads to significant improvements in computation time. Moreover, since we do not have to manipulate dense solution vectors, this also leads to significant savings in memory usage.

We illustrate this technique in Fig. 3: entry a -1 23 is required, thus the only nodes of tree that have to be visited lie on the path from node 3 to the root node (6) and on the path from the root node to node 2. Therefore, one does not have to perform operations at nodes 4 and 1.

When many entries of the inverse are requested, they cannot generally be computed all at once (mainly because of memory 23 . The traversal of the tree is reduced to the path from 3 to 6 and the path from 6 to 2; no computation is performed at nodes 1 and 4. usage), but they can be computed by blocks, which allows to take advantage of efficient dense linear algebra kernels. Work has been performed in order to find optimal ways to form the blocks in different contexts [START_REF] Amestoy | On computing inverse entries of a sparse matrix in an out-of-core environment[END_REF] and to improve the parallel efficiency.

Calculation

A substantial time is spent in computing A = H T H with a basic algorithm. The use of an appropriate algorithm to perform the operation A = H T H helps to reduce the computation time (see Section 4.1). The MUMPS solver is used to solve the system of equations as described in Section 4.2. Finally, the error bars on the solution are computed, which means the calculation of the diagonal elements of inverse matrix. The new A -1 feature of MUMPS is compared with several algorithms, in terms of computation time in Section 4.3.

Improvements of the computation of A = H T H

The computation of the normal equation A = H T H is of paramount importance in many problems, yet is a very challenging operation due to the considerable amount of symbolic operations needed to compute the sparsity structure of A. For this reason efficient algorithms have been developed in the past. To perform this operation we decided to use part of a larger code developed by [START_REF] Puglisi | QR factorization of large sparse matrix overdetermined and square matrices with the multifrontal method in a multiprocessing environment[END_REF] for computing the QR factorization of sparse matrices. The used part was originally developed to compute only the structure of the A matrix and, thus, we had to extend it in order to compute the coefficients values. This was possible thanks to the help of the original code developer.

One important feature of this code offers the possibility to update the elements of A that are changed after modification of a Based on an original package from [START_REF] Puglisi | QR factorization of large sparse matrix overdetermined and square matrices with the multifrontal method in a multiprocessing environment[END_REF] and improved as suggested by the author.

b N matrix vector product are used following the Compressed Column Storage scheme, but for each operation a dense vector of length M (with many zero element), that represents a column of H is built in place. some numerical values of the columns of H without recomputing the whole matrix (the technique used to compute simultaneously the solution and the background pattern in the algorithm is described in 2.2.3 and Appendix A.4).

Table 3 shows the time reduction achieved for both the computation of A = H T H and its update after the modifications of some columns of H. The results in the first column are obtained with the code extracted from the software package by [START_REF] Puglisi | QR factorization of large sparse matrix overdetermined and square matrices with the multifrontal method in a multiprocessing environment[END_REF] and improved as suggested by the author. The results on the second column, instead, are obtained by computing N matrix vector products where, for each product, a dense vector of length M (with many zero elements) corresponding to a column of H is built in place.

The gain over a simple basic algorithm is significant (a factor ∼300) and demonstrates the interest of using specialized libraries dedicated to sparse matrix computations.

Solving a sparse linear system

Here we briefly illustrate the interest of exploiting sparsity of the matrix when solving a linear system. In Table 4, we compare the time for solving linear systems arising from our application using a dense solver (LAPACK [START_REF] Anderson | LAPACK: a portable linear algebra library for high-performance computers[END_REF]) and a sparse solver (MUMPS). Times are in seconds and include the L D L T factorization of a symmetric matrix A of order N and the computation of the solution x of the system Ax = b (x and b are vectors of length N). In the results related to MUMPS, the time for the analysis phase is included. In the second row of the table, instead, the matrix is treated as dense, hence its full storage is used and no analysis phase is performed. For the largest two problems, the dense algorithm cannot be used as the memory requirements are roughly 23 GB and 167 GB respectively. We can extrapolate that on this system, the run time would be around 22 h for the largest problem (instead of 6.7 s using a sparse algorithm).

The results in Table 4 confirm that sparse, direct solvers achieved a good scalability on the problems of our target application whereas dense linear algebra kernels quickly exceed the limit of modern computing platforms.

Time to compute error bars

In this section we present experimental results related to the computation of error bars or, equivalently, of the diagonal entries 0.28 0.9 4.9 36.0 9.5

Note: Execution times (in seconds) for the computation of all the diagonal entries of the A -1 matrix with the left-looking, brute force and MUMPS A -1 methods. For the brute force approach results are provided for blocks of size 1 and 128.

of the inverse matrix A -1 . Our approach that relies on the pruned tree, presented in Section 3.3.1, is compared to the basic, leftlooking approach described in [START_REF] Stewart | Matrix Algorithms[END_REF]. In the case of a symmetric matrix, this approach computes the diagonal entries of the inverse matrix as

a -1 ii = N k=1 1 d kk l -1 ki 2
where we denoted with a -1 ij and l -1 ij the coefficients of A -1 and L -1 , respectively. This amounts to computing, one at a time, the columns of L -T and then summing the corresponding contribution onto the a -1 ii coefficients. In this algorithm, the sparsity of the right-hand side and of the factor matrix L is exploited but not completely, and the experimental results discussed below show that this results in a higher execution time. Furthermore, because of memory issues, this simple algorithm does not allow to simultaneously compute many diagonal entries of A -1 ; clearly this is also a limiting factor for performance. Our implementation of this method is based on the LDL package [START_REF] Davis | User guide for LDL, a concise sparse Cholesky package[END_REF]. As a second term of comparison we also provide experimental results for a brute force approach with no exploitation of sparsity of the right-side and solution vectors. For this purpose, we use directly the MUMPS package and solve several systems of equations in order to compute the inverse matrix. In addition, we analyze the influence of grouping the computation of the diagonal entries (1 right-hand-side (RHS) at a time or 128 at a time).

The experimental results for the three methods described above are reported in Table 5. For the sake of this comparison, all these methods are executed in sequential mode although the code of the brute force approach and of the MUMPS A -1 feature are parallel. The experiments were carried out on the above-mentioned system.

These results show that the brute force algorithm becomes competitive with respect to the simple algorithm when the entries are processed by blocks. The MUMPS A -1 feature described in Section 3.3.1 is significantly faster than all other approaches and the gain increases with the size of the problem. Pruning leads to clear gains over a strict traditional solution phase. The gain is even larger for the largest problem due to the good scalability of the A -1 algorithm with respect to the problem size. The simple, left-looking approach shows reasonable performance for small problems, but could not be tested on our largest matrix because numerical pivoting, not available in LDL package, is needed during factorization to obtain an accurate solution.

Results and discussion

The MUMPS solver and its A -1 functionality are the core tools to solve systems of equations related to the measurements of the sources intensity. Fig. 4 shows the application to the determination of the different components of the Galaxy spectrum. The related analysis is performed in 24 consecutive energy bands in order to extract counts spectra. The counts spectra are then converted into photon spectra. The details are given in Bouchet et al. (2011). Another application is the study of the intensity variations of a peculiar source or sky region. Fig. 5 shows the intensity in function of the time (exposure) of some of the sources located in the central crowded region of the Galaxy. For this application, the end and start of the ''time-bins'' are determined by a segmentation algorithm, which is based on the efficient L D L T factorization of symmetric matrix provided by MUMPS, details can be found in Bouchet et al. (in press).

We have demonstrated that even for the basic operation such as sparse matrix product, it is better to use dedicated algorithms or libraries (4.1).

The MUMPS solver is very effective on the sparse matrix structure arising from astrophysical problems encountered with SPI. This solver is robust and the matrix factorization is stable against rounding errors. It provides many numerical preprocessing options and implements robust pivoting strategies, which make it one of the most numerically stable solvers available. The matrices arising from the INTEGRAL/SPI application are symmetric and indefinite; they are not extremely challenging numerically, but they do require two-by-two numerical pivoting for stability (in Table 5 the LDL package could not provide an accurate solution on our largest matrix). The proposed approach not only leads to substantial time reduction but also enables the resolution of large sparse system of equations which could not be solved using basic algorithms.

Other sparse linear systems solvers exists and have been used to validate the performance reported in the experimental section; for an exhaustive list see for example [START_REF] Bai | Templates for the Solution of Eigenvalue Problems: A Practical Guide[END_REF], but they all lack a function to compute also the error bars on the solution quickly, which is mandatory in our astrophysical application.

The A -1 feature in MUMPS (computation of selected inverse entries) did not exist before the beginning of this study, the IN-TEGRAL/SPI application was actually one of the motivating applications for developing techniques for the computation of inverse entries, and for releasing a publicly available code. This functionality allows to compute easily and rapidly the error bars on the solution. The gain in time over already optimized algorithms is important.

Among other methods to solve the problem completely, solution and error bars, one should mention alternative methods such as Monte Carlo Markov Chains (Metropolis et al., 1953;[START_REF] Hastings | [END_REF][START_REF] Neal | Probabilistic inference using Markov chain Monte Carlo methods[END_REF] or Simulated Annealing (Kirkpatrick et al., 1983). Such advanced statistical tools can provide the best fit and the variances of the solution of both linear and non-linear systems of equations. In particular MCMC methods could be useful when computing error bars, in case of complex constraints on the function. However, these methods may be very prohibitive in time, especially if high precision on the parameters is required; they have in general a weak or non-guaranteed convergence and are not the best suited for our needs, given the complexity of our problem.

Conclusions

We have developed algorithms to process years of data and to enhance the production of INTEGRAL hard X/soft γ -ray survey catalogs. These methods have been successfully applied to a set of ∼6 years of data (Bouchet et al., 2011). We have shown that, for processing efficiently and accurately years of data, it is critical to use algorithms that take advantage of the sparse structure of the transfer function (matrix), such as those implemented in the MUMPS software.7 It was also demonstrated that error bars can be obtained at a relatively inexpensive cost (the same order of magnitude as a simple problem solution) thanks to a recently developed algorithmic feature that efficiently computes selected entries of the inverse of a matrix. In addition, thanks to many efforts in optimization, gains are achieved both in memory usage and in computation time. Hence, it is possible to process large datasets in a reasonable time and to compute the diagonal of the covariance matrix, and thus error bars, in a rather short time. More generally, the ideas described here can be applied to a large variety of problems. Finally, we are today able to solve sparse linear systems with millions, sometimes billions, of unknowns. Although we have not used explicitly parallel computing but instead performed many sequential computations at the same time (for each energy band, etc.), the proposed approach can also be used directly in a parallel setting on massively parallel machines.

In the near future, instruments will commonly create datasets of a few tens to a few hundreds of Terabytes for a single observation project. Many of the current tools and techniques for managing large datasets will not scale easily to meet this challenge. Surveys of the sky already require parallel computing in order to be performed. New surveys will demand orders of magnitude increases in the available data and therefore in data processing capabilities. It is also a challenge for scientists who need to extract a maximum of science from the data. Exciting scientific breakthroughs remain to be achieved as astronomers manipulate and explore massive datasets, but they will require advanced computing capabilities, infrastructure and algorithms.

The INTEGRAL/SPI project has been completed under the responsibility and leadership of CNES. We are thankful to ASI, CEA, CNES, DLR, ESA, INTA, NASA and OSTC for support. We are very grateful to Chiara Puglisi, research engineer at INPT-IRIT, for her contribution to improve the performance of the H T H product in Section 4.1.

Appendix. Schematic solution of the system of equations

A.1. Maximum likelihood estimator (MLE) of the solution

In the case of a low number of counts, it is recommended to use the MLE of the solution instead of the χ 2 solution. Following Cash (1979), we maximize the likelihood function,

L = -2 × M i=1 e i -y i ln e i (A.1)
where e i is the model of the data obtain through the relation e = Hx.

A.2. Optimization of the non-linear problem

To optimize this non-linear problem, potentially with bound constraints (such as positivity of the solution), there are at least four approaches: (a) Newton type methods (or modified Newton methods): they use the Hessian matrix to define a search direction and hence need the solution of a large linear system of equations at least at each few iterations. They are the most powerful methods available and can find the solution in a few iterations. (b) Quasi-Newton methods: they build an approximation of the Hessian at each iteration. They optimize a quadratic function in at most n iterations (n being the number of unknowns). (c) Conjugate-gradient methods: unlike the Newton-type and quasi-Newton methods, conjugate gradients methods do not require the storage of an n by n Hessian matrix and thus are suited to solve large problems. The gradient of the function (a vector of length n) is used to define the search direction. They are not usually as reliable or efficient as Newton-type methods and often a relatively large number of iterations has to be performed before obtaining an acceptable solution. (d) Simplex [START_REF] Nelder | [END_REF], simulated annealing (Kirkpatrick et al., 1983) or Monte Carlo Markov Chain (MCMC) [START_REF] Neal | Probabilistic inference using Markov chain Monte Carlo methods[END_REF] can also be considered, but they are often prohibitive in time.

Methods (a) and (b) are known as order-2 optimization methods (gradient and Hessian used), (c) as an order-1 optimization method (gradient used), while method (d) can use only the function value.

To use a Newton type method (order-2), we need to compute the gradient G and the Hessian H ess of the function . As for the LSQ case, the variance of the solution is obtained thanks to the inversion of the Hessian matrix (note that in the limit lim e →d

G j = ∂L ∂x j = 2 × M i=1
d i e 2 i = 1 d i = 1 σ 2 i
, the likelihood (H ess ) and chisquare (A) Hessian matrices are the same). A guess solution to this non-linear optimization problem is the LSQ solution.

A.3. Codes for non-linear optimization

The fitting algorithm, based on the likelihood test statistic, is a non-linear optimization problem. To optimize a non-linear problem, potentially with bound constraints, a Newton type method, known for its efficiency and reliability can be used, as we already have a solver for large sparse systems at hand. A software package for large-scale non-linear optimization such as IPOPT 8 (Interior Point OPTimizer) can be used. IPOPT uses a linear solver such as MUMPS or MA57 [START_REF] Duff | MA57-a code for the solution of indefinite sparse symmetric linear systems[END_REF] as a core algorithm. For more details on this large-scale nonlinear algorithm, see Wächter and Biegler (2006). A few similar software packages for large-scale non-linear optimization exist, among them LANCELOT [START_REF] Conn | Numerical experiments with the LANCELOT package (Release A) for large-scale nonlinear optimization[END_REF], MINOS [START_REF] Murtagh | A projected Lagrangian algorithm and its implementation for sparse nonlinear constraints[END_REF] and SNOPT [START_REF] Gill | SNOPT: an SQP algorithm for large-scale constrained optimization[END_REF].

A.4. ''Empty-field'' auto-computation

Sometimes the ''empty-field'' or ''uniformity map'' U has to be computed with the solution. In order to preserve the linearity of the problem, we have adopted the algorithm described below. We consider that if the solution x is known,

y i = K 0 j=1 h ij x j + N j=K 0 +1
h ij x j = y B i + y S i . (A.4)

In the above formula y S i ≡ D S dp is the counts due to the sources, assumed to be known. y B i ≡ U d I bg p t dp is the background contribution, I bg is assumed to be known and U is to be estimated. At this stage, using the model of the sky described by (A.4) For the MLE statistic, we do not have to preserve the linearity of the problem and hence the computation of the improved ''emptyfield'' pattern can be done during the non-linear optimization process. On another side, the algorithm is simplified if we proceed similarly as in the LSQ case. Then, we wish to maximize the following quantities for each of the working detectors, One should mention that it is possible to compute, similarly, an ''empty-field'' pattern on some restricted time interval instead of the whole dataset; the best ''empty field'' for pointing intervals p k to p k+1 is then, (A.9)

L(d) = -
                           U MLE (d, k) =

A.5. ''Empty-field'' schematic construction

A sub-optimal algorithm to obtain both the sources and the background fluxes, as well as the improved ''empty-field'' pattern is described in Algorithm 1. We start with an approximation U 0 and apply some iterative refinement. In practice, the algorithm converges in a few iterations.

Algorithm 1 Computation of the ''Empty field'', the solution and its variance 1: U = U 0 , compute the structure of the Hessian (A or H ess ) 2: for i = 1 to itermax do {Iterative computation of U and x} If χ 2 stops decreasing or the likelihood function stops increasing then go to step 8 7: end for 8: Compute H at the solution (if not already done) and the diagonal of H -1 to obtain the uncertainties on the solution

Fig. 1 .

 1 Fig. 1. SPI imaging principle. The mask consists of elements transparent or opaque to the radiation. The opaque elements (made of tungsten) are shown in black. The shadowgram of the mask casts onto the detector plane (camera) depends on the source direction. Here the counts in the different detectors of source-1 and source-2 are shown in black and red. The counts recorded by the detectors are the sum of all the contributions from all the sources in the FoV. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

  p I bg + ǫ p with P p = t dp × U d .

  (a) The highest energy bands (E 100 keV) are less problematic in terms of number of parameters to determine, as illustrated by the 200-600 keV band. The sky model contains only 29 sources which are essentially not variable in time (given the instrument sensitivity). The number of unknowns is N ≃ 5900. (b) The lowest energy bands (E 100 keV) are more problematic. We use the 25-50 keV band. The sky model contains 257 sources variable on different timescales. When the background intensity is assumed to vary on ∼6 h timescale, N ≃ 22 500

  methods are commonly based on Gaussian elimination, with the aim to factorize the sparse matrix, say A, of the linear system Ax = b into a product of ''simpler'' matrices called factors. Typically, A can be factored into A = L U where L and U are lower and upper triangular matrices respectively, or A = L D L T , where

  (a) Factors L + L T . (b) Elimination tree of A.

Fig. 2 .

 2 Fig. 2. The factors and the elimination tree of a symmetric matrix A. (a) Pattern of the L + L T factors of A with filled-in entries shown with squares, (b) the elimination tree of A where the children of a node are drawn below the node itself.

Fig. 3 .

 3 Fig. 3. Computation of a -123 . The traversal of the tree is reduced to the path from 3 to 6 and the path from 6 to 2; no computation is performed at nodes 1 and 4.

Fig. 4 .Fig. 5 .

 45 Fig.4. Different contributions to the total emission at hard X-ray and soft gammaray energies in the central radian of the galaxy. The data points shown in black (plus filled circle) correspond to the contribution due to 270 point sources. The average spectrum of these sources can be viewed at http://sigma-2.cesr.fr/integral/. The data points shown in blue correspond to the diffuse emission. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

  approximation of U by minimizing again LSQ or maximizing MLE statistics 5: Update H (The first K 0 columns of H and update the new Hessian matrix (Section 4.1)) 6:

  s + 1) matrix and M = n d × n p . The parameters to be determined, x = (I bg , I s 1 , . . . , I s n s ) is a vectors of length n s +1. The data y = (D raw 1 , D raw 2 , . . . , D raw n p ) and the associated statistical errors ǫ = (ǫ 1 , ǫ 2 , . . . , ǫ n p ) are vectors of length M.

Table 1

 1 Sparsity of matrices H and H T H. Matrix) is the so-called density of the matrix: the ratio between the number of non-zero elements in the matrix and the total number of elements in the matrix (M × N for H and N 2 for A = H T H, where M is the number of rows of H. The matrix H arising from the diffuse emission study have M = 672 495 rows. The number of non-zero elements is constant nz = 27 054 399 for the matrices with N ≥ 22 503 corresponding to the 25-50 keV band and nz = 4 677 821 for the 200-600 keV band. The matrix with N = 3578 has M = 124 921 rows and nz = 11 948 840 non-zero elements.

	N	ρ(H) (%)	ρ(A) (%)
	3 578	2.67	2.96	Central galaxy (27-36 keV)
	9 437	1.01	1.05
	5 900	0.12	0.13	Diffuse emission 200-600 keV
	22 503	0.18	0.28	Diffuse emission 25-50 keV
	55 333	0.07	0.09
	149 526	0.03	0.04	Simulation (25-50 keV)
	Notes: ρ(

Table 2

 2 Number of non-zeros in the original matrix A and in the L factor of the A = L D L T factorization for different problems of our experimental set.

	Matrix size	3578	9437	22 503	55 333	149 526
	nz(A)	378 475	932 143	1436 937	2705 492	9 379 127
	nz(L)	519 542	1380 444	2885 821	9189 447	14 432 264

Table 3

 3 Time for the computation of A = H T H.

	Matrix A = H T H	Improved algorithm a	Simple algorithm b
	N = 22 503 Full matrix	28.2	5 779
	5000 H columns modified	0.43	258
	N = 149 526 Full matrix	41.2	18 585
	5000 H columns modified	0.13	31.7
	Note: Times are in seconds. H is an N by M matrix; here M = 672 495 and H has 27 054 399 non-zero elements.

Table 4

 4 Times (in seconds) for the computation of the solution.

	Matrix size	3578	9437	22 503	55 333	149 526
	Sparse	0.2	0.7	1.6	8.0	6.7
	Dense	1.2	20.1	169.9	/	/

Table 5

 5 Time to compute the diagonal elements of the inverse of a symmetric matrix.

	Matrix size	3578	9437	22 503	55 333	149 526
	Left-looking	28.2	376.1	2567.9	489.1	/
	MUMPS (1 RHS)	3.77	38.4	204.1	1324.9	8230.5
	MUMPS (128 RHS)	1.32	7.34	45.5	245.6	2833.5
	MUMPS A -1

  , a rough estimate of the pattern is U d ≃ Expression for the detector patternFor the LSQ statistic, we wish to minimize the following quantities for each of the working detectors,χ 2 (d) = LSQ solution U LSQ (d) is U LSQ (d) =

	n p p=1	(D raw dp -D S dp ) × I bg p t dp /σ 2 dp p=1 p t dp ) 2 /σ 2 dp (I bg n p	.	(A.6)
					np p=1 (D raw dp -D S dp ) np p=1 t dp	.
				A.4.1. n p p=1	D raw dp -D S dp -U d I bg p t dp σ dp

2 for d = 1, . . . , n p . (A.5) 8 IPOPT is available at https://projects.coin-or.org/Ipopt.

The

  2 U d I bg p t dp -D raw dp ln [D S dp + U d I bg p t dp ] for d = 1, . . . , n p .

		n p
		p=1	D S dp + (A.7)
	The MLE solution U MLE (d) is
	U MLE (d) =	n p p=1	(D raw dp -D S dp ) I p t dp bg n p	.	(A.8)
			p=1

With the Compressed Sensing approach[START_REF] Bobin | [END_REF],Wiaux et al. (2009) and references therein), it is possible to find a sparse solution even if the system is underdetermined and for systems in which the matrix is sparse.

For clarity, we omit to weight the matrix H and the data by the inverse of the data standard deviation, see Section 2.4.1

The number of counts per detector is high enough to use the Gaussian statistics.

Available at http://mumps.enseeiht.fr/ or http://graal.ens-lyon.fr/MUMPS/.

Acknowledgments

funded by ESA member states (especially the PI countries: Denmark,


Jérome Bobin

Jean-Luc Starck

Jalal M Fadili

Yassir Moudden

Sparsity and Morphological Diversity in Blind Source Separation

published or not. The documents may come   L'archive ouverte pluridisciplinaire


R Aimino

S Vaienti
email: vaienti@cpt.univ-mrs.fr

R Hervé

Aimino

A NOTE ON THE LARGE DEVIATIONS FOR PIECEWISE EXPANDING MULTIDIMENSIONAL MAPS

Introduction

There are different ways to establish large deviation principles (LDP) for dynamical systems. One of them is the so-called "Laplace Method", which relies on the spectral properties of the Perron-Frobenius, or transfer, operator. This strategy has been developed in a very general and abstract setting by Hennion and Hervé in [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF]. They assume that the transfer operator acts on a Banach spaces of measurable functions, and it is quasi-compact on it, i.e. it has a spectral gap. The existence of an invariant probability measure follows immediately, and, by using perturbation theory for linear operators, they derive a few others statistical properties. This approach covers a lot of systems, for instance expanding maps of the interval [15, XII.1], Gibbs measures for subshift of finite type [15, XII.2], and expanding Young towers [START_REF] Young | Statistical properties of dynamical systems with some hyperbolicity[END_REF]. Nevertheless, this theory seems inappropriate for expanding (discontinuous) maps in higher dimension, like those treated by several authors [START_REF] Adl-Zarabi | Absolutely continuous invariant measures for piecewise expanding C 2 transformations in R n on domains with cusps on the boundaries[END_REF][START_REF] Blank | Stochastic properties of deterministic dynamical systems[END_REF][START_REF] Boyarsky | Absolutely continuous invariant measures for piecewise expanding C 2 transformations in R n[END_REF][START_REF] Buzzi | Absolutely continuous invariant probability measures for arbitrary expanding piecewise R-analytic mappings of the plane[END_REF][START_REF] Buzzi | Zeta functions and transfer operators for multidimensional piecewise affine and expanding maps[END_REF][START_REF] Cowieson | Stochastic stability for piecewise expanding maps in R d[END_REF][START_REF] Liverani | Multidimensional expanding maps with singularities: a pedestrian approach, to appear in Ergodic Theory and Dynamical Systems[END_REF][START_REF] Poppe | Inadequacy of the bounded variation technique in the ergodic theory of higher-dimensional transformations in R n[END_REF][START_REF] Saussol | Absolutely continuous invariant measures for multidimensional expanding maps[END_REF][START_REF] Thaler | Estimates of the invariant densities of endomorphisms with indifferent periodic points[END_REF][START_REF] Thomine | A spectral gap for transfer operators of piecewise expanding maps[END_REF][START_REF] Tsujii | Absolutely continuous invariant measures for expanding piecewise linear maps[END_REF]. Indeed, Hennion and Hervé assume that the Banach space on which the transfer operator acts, consists of measurable functions defined everywhere, but for the higher dimensional systems quoted above, the functional spaces usually considered (bounded variation, quasi-Hölder or Sobolev spaces) consist of classes of equivalence of functions modulo the reference measure, and hence they are only defined almost everywhere. Furthermore, in [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF], the Dirac masses must belong to the topological dual of the Banach space, so this theory cannot be applied directly for those systems.

Nevertheless it appears that one could slightly modify the proofs from [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF] in order to deal with a Banach space consisting of classes of functions. In particular we will SV warmly thanks E. Ugalde for the kind invitation to participate to the Conference in honor of Valentin Afraimovich; RA and SV express their sincere gratitude to I. Melbourne who helped them to simply the proofs.

do it for the functional space of quasi-Hölder functions mostly investigated in [START_REF] Saussol | Absolutely continuous invariant measures for multidimensional expanding maps[END_REF] and which verify an additional algebraic assumption which also plays a role in the Hennion-Hervé approach. As a consequence we will get the large deviations principle for such systems and, as far as we know, this result is not present in the literature. We also prove the central limit theorem, but for the latter one already disposes of the Gordin-Liverani theorems [START_REF] Gordin | The central limit theorem for stationary processes[END_REF][START_REF] Liverani | Central Limit Theorem for Deterministic Systems[END_REF]. Actually a weaker result for the large deviations of systems like those considered above has been recently obtained in [START_REF] Alves | From rates of mixing to recurrence times via large deviations[END_REF]. We will comment about the difference with the spectral technique presented in this note in the Remark 2 below. We anticipate here that the paper [START_REF] Alves | From rates of mixing to recurrence times via large deviations[END_REF] furnishes an upper bound for the deviation functions and whenever the correlation functions involving L 1 observables decay to zero with a summable rate. In order to check these assumptions for our systems we should further require that the density of the invariant measure is essentially bounded from below, but this assumption is not necessary in the spectral approach discusses later on.

Assumptions and statement of the results

We now give the precise assumptions under which the LDP is valid. Let (X, A, m) a probability space, and T : X → X a measurable transformation, non singular with respect to m. Under these conditions, the Perron-Frobenius operator P : L 1 (m) → L 1 (m) is well defined by P f = dm f dm , where m f (A) = T -1 A f dm is absolutely continuous with respect to m. We stress here the fact that the functions under consideration are complex valued, as required by the spectral theory we are going to use below. The transfer operator enjoys some classical properties that we resume below; see [START_REF] Boyarsky | Laws of Chaos : Invariant Measures and Dynamical Systems in One Dimension, Probability and its Applications[END_REF] or [START_REF] Lasota | Chaos, Fractals and Noise: Stochastic Aspects of Dynamics[END_REF] for more details.

(1) Linearity : P is a linear operator on L 1 (m), satisfying ||P f || 1 ≤ ||f || 1 for all f ∈ L 1 (m); (2) Positivity : For all f ∈ L 1 (m) such that f ≥ 0 m-ae, we have P f ≥ 0 m-ae;

(3) Preservation of integrals : For all f ∈ L 1 (m), we have P f dm = f dm; (4) Duality : For all f ∈ L 1 (m) and g ∈ L ∞ (m), we have f (g The last assertion means that there exists a decomposition

• T ) dm = (P f )g dm; (5) Invariant Measures : f ∈ L 1 (m)
P = s i=1 λ i Π i + Q
where λ i are complex numbers of modulus 1, Π i are finite-rank projections satisfying Π i Π j = 0 when i = j and Q is a bounded operator on B with spectral radius strictly less than 1 and satisfying QΠ i = Π i Q = 0 for all i. The spectrum of P consists then of a finite number of eigenvalues of modulus 1, with finite multiplicity, and the rest of the spectrum lies in a disc centered at 0 with radius strictly less than 1. When B is compactly injected in L1 (m), this can be deduced from a Lasota-Yorke type inequality, by means of the Ionescu-Tulcea and Marinescu theorem [START_REF] Ionescu-Tulcea | Théorie ergodique pour des classes d'opérations non complètement continues[END_REF][START_REF] Hennion | Sur un théorème spectral et son application aux noyaux Lipschitziens[END_REF]. See [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF] for precise definitions and results about quasi-compactness.

Under those conditions, the existence of an T -invariant probability µ absolutely continuous w.r.t m, such that dµ dm ∈ B is a classical result : for every f ∈ B such that f ≥ 0, f dm = 1, and so in particular for f = 1, quasi-compactness implies that the sequence 1 n n-1 k=0 P k f converges in B to a function f ⋆ such that the measure µ f with dµ f dm = f ⋆ is an acip. Furthermore, 1 is an eigenvalue of P . If we assume that 1 is a simple eigenvalue of P , then there exists an unique acip µ such that dµ dm ∈ B. From now, we will always assume that 1 is a simple eigenvalue, and that there is no other eigenvalue of modulus 1. µ will denote the unique acip, and v ∈ B its density. We then have 1 , for all n ≥ 1 and f ∈ B

P n f =< m, f > v + Q n f
As a consequence, we get exponential decay of correlation : there exists C ≥ 0 and 0 ≤ λ < 1 such that for every f ∈ B and every g ∈ L ∞ (µ) we have

f (g • T n ) dµ -f dµ g dµ ≤ Cλ n ||f || B ||g|| L ∞ µ
Let now φ : X → R a bounded observable which lie in B, with zero mean φ dµ = 0. Denote by S n the Birkhoff sums :

S n = n-1 k=0 φ • T k
We are now able to state the LDP : Theorem 1. (Large Deviation Principle) Under the above conditions, the limit σ 2 = lim n→∞ ( Sn √ n ) 2 dµ exists, and if σ 2 > 0, then there exists for some ǫ 0 > 0 a rate function c : ]ǫ 0 , +ǫ 0 [→ R, continuous, strictly convex, vanishing only at 0, such that for every 0 < ǫ < ǫ 0 and every probability measure ν with ν ≪ m and dν dm ∈ B, we have

lim n→∞ 1 n log ν(S n > nǫ) = -c(ǫ)
As an easy consequence of the techniques introduced in the next section, we also get the central limit theorem. We denote with N (0, σ 2 ) the Dirac mass δ 0 if σ 2 = 0, and the probability with density

1 σ √ 2π e -t 2 2σ 2 with respect to Lebesgue if σ 2 > 0.
Theorem 2. (Central Limit Theorem) Under the same assumptions of Theorem 1, Sn √ n converges in distribution to N (0, σ 2 ) in the probability space (X, A, ν) for every probability ν with ν ≪ m and dν dm ∈ B : for every bounded continuous function g : R → R, we have

lim n→∞ g( S n √ n )dν = g dN (0, σ 2 )
Remark 3.

(1) Theorems 1 and 2 apply in particular for ν = m and ν = µ, so the LDP and the CLT are valid for both reference and invariant measures. (2) As we anticipated in the Introduction, the paper [START_REF] Alves | From rates of mixing to recurrence times via large deviations[END_REF] gives an upper bound for the large deviation function and under related assumptions. In particular Th. E in [START_REF] Alves | From rates of mixing to recurrence times via large deviations[END_REF] states the following, with our notations. Let us suppose that T preserves an ergodic probability measure µ; then let B ⊂ L 1 (µ), φ ∈ B, and assume that there exists

ξ(n) with ∞ n=0 ξ(n) < ∞ such that for all ψ ∈ L 1 (µ) we have φ(ψ • T n ) dµ -φ dµ ψ dµ ≤ ξ(n)||f || B ||g|| L ∞ µ .
Then there exists τ = τ (φ) > 0 and, for every ǫ > 0, there exists

C = C(φ, ǫ) > 0 such that µ(S n > nǫ) ≤ Ce -τ n .
The proof of this result relies on a martingale approximation. Section C in [START_REF] Alves | From rates of mixing to recurrence times via large deviations[END_REF] provides examples of systems for which one gets exponential decay of correlations (and hence ξ(n) is summable) and against ψ ∈ L 1 (µ): this last assumption requires the density to be bounded from below. In conclusion the result of our note extends the previous one in two directions : first, we obtain a LDP, and not only an upper bound; secondly, we do not have to assume anything about the density, but the fact that it is in the L 1 norm with respect to the conformal (reference) measure.

Proofs

We begin by observing that the existence of the variance σ 2 results from a straightforward computation : since the sequence φ(φ • T n ) dµ decays exponentially fast, it is absolutely summable, and then we see, by expanding the term S 2 n , that the limit σ 2 = lim n→∞ ( Sn √ n ) 2 dµ exists and we have

σ 2 = φ 2 dµ + 2 +∞ n=1 φ(φ • T n ) dµ
We assume from now σ 2 > 0. Our proof of the LDP follows closely [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF] except for a minor modification, which will be mentioned later. The same approach had been employed in [START_REF] Rey-Bellet | Large deviations in non-uniformly hyperbolic dynamical systems[END_REF]. Let f ∈ B the density of the measure ν with respect to m. We will apply Gartner-Ellis theorem [START_REF] Dembo | Large Deviations, Techniques and Applications[END_REF][START_REF] Ellis | Entropy, Large Deviations and Statistical Mechanics[END_REF], so we are interested in the convergence of the sequence 1 n log e θSn f dm for θ ∈ R small enough. We introduce the "Laplace transform" operators P z , for z ∈ C, defined by

P z (f ) = P (e zφ f ), f ∈ B
Assuming for a moment that P z is well defined, we see immediately that we have e θSn f dm = P n θ (f ) dm. In order to prove that P z is a bounded operator on B, we just have to check that e zφ ∈ B. Since B is a Banach algebra, the sequence n k=0

(zφ) k k!
converges in B, and hence in L 1 (m). On the other hand, this sequence converges uniformly, and hence in L 1 (m), to e zφ , and so we get that

e zφ = +∞ n=0 (zφ) n n! in B.
It also proves that the map z → P z is holomorphic and we have the expansion

P z = +∞ n=0 C n n! z n
where C n (f ) = P (φ n f ).

We can now apply perturbation theory for linear operator to prove the following result. The proof relies on analytic functions of operators, see [START_REF] Dunford | Linear Operators, Part I : General Theory[END_REF], or on the implicit function theorem, see [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF]. For θ > 0, we denote D θ = {z ∈ C / |z| < θ}.

Proposition 4. There exist θ 0 > 0, C > 0, η 1 , η 2 > 0 and holomorphic functions λ(.) :

D θ 0 → C, v(.) : D θ 0 → B, ϕ(.) : D θ 0 → B ⋆ and Q(.) : D θ 0 → L(B) such that for all z ∈ D θ 0 (i) λ(0) = 1, v(0) = v, ϕ(0) = m, Q(0) = Q; (ii) P z (f ) = λ(z) < ϕ(z), f > v(z) + Q(z)f for all f ∈ B; (iii) < ϕ(z), v(z) >= 1; (iv) Q(z)v(z) = 0 and ϕ(z)Q(z) = 0; (v) |λ(z)| > 1 -η 1 ; (vi) ||Q(z) n || ≤ C(1 -η 1 -η 2 ) n .
So, for all n ≥ 1, we have

P n z (f ) = λ(z) n < ϕ(z), f > v(z) + Q(z) n f
We can say much more on eigenvalues and eigenvectors when z = θ is real. At this point, we need to show that for every positive function f ∈ B with f = 0, there exists a positive linear form ϕ ∈ B ⋆ such that < ϕ, f > > 0. In the context of [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF], since functions are defined everywhere, there exists x ∈ X such that f (x) > 0, and so the Dirac mass δ x does the job. In our context, Dirac masses are not available, but the reference measure is usable , since necessarily < m, f >> 0, otherwise, f would be 0 m-ae, and so f = 0 in B. This was not the case in [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF] because they consider functions defined everywhere, and not classes of equivalence. We can also use arguments from complex Banach lattice theory [START_REF] Meyer-Nieberg | Banach Lattices[END_REF] : a modification of the Hahn-Banach theorem shows that there exists a positive bounded linear form ϕ on B R = {f ∈ B / f (x) ∈ R m-ae}, such that < ϕ, f >= 1, and then we can extend it on all B. This argument could be employed in more abstract contexts, where the Banach space B consists of distributions-like objects and when we don't have a good knowledge of its topological dual. Proposition 5. There exists 0 < θ 1 < θ 0 such that for every θ ∈ R with |θ| < θ 1 , we have λ(θ) > 0. Furthermore, v(.) and ϕ(.) can be redefined such that v(θ) ≥ 0, ϕ(θ) ≥ 0.

Proof. As P θ is a real operator, we have P θ f = P θ f for all f ∈ B. So, we have

P θ v(θ) = P θ v(θ) = λ(θ) v(θ). Since λ(θ)
is the unique eigenvalue of P θ with maximal modulus, we get λ(θ) = λ(θ), and hence λ(θ) ∈ R. Since λ(0) = 1, by a continuity argument, we obtain λ(θ) > 0 for small θ. For z ∈ C small enough, < ϕ(z), 1 > = 0. We define ṽ(z) =< ϕ(z), 1 > v(z) and φ(z) =< ϕ(z), 1 > -1 ϕ(z). Those new eigenfunctions satisfy obviously the conclusions of the previous proposition. We have just to prove that ṽ(θ) and φ(θ) are positive for θ ∈ R small enough. By the spectral decomposition of P θ , we see that λ(θ) -n P n θ 1 goes to ṽ(θ) in B, and hence in L 1 (m). We then get ṽ(θ) ≥ 0 because P θ is a positive operator and λ(θ) is positive too. Now, let ψ(θ) ∈ B ⋆ positive such that < ψ(θ), ṽ(θ) >= 1. Then, λ(θ) -n (P ⋆ θ ) n ψ(θ) goes to < ψ(θ), v(θ) > ϕ(θ) = φ(θ), which proves that φ(θ) is a positive linear form.

We denote Λ(θ) = log λ(θ)

We then have Proposition 6. There exists 0 < θ 2 < θ 1 such that for every θ ∈ B with |θ| < θ 2 and every f ∈ B with f ≥ 0 and f dm = 1, we have

lim n→∞ 1 n log e θSn f dm = Λ(θ)
Proof. We have the identity

e θSn f dm =< m, P n θ (f ) > = λ(θ) n < ϕ(θ), f > < m, v(θ) > + < m, Q(θ) n f > = λ(θ) n (< ϕ(θ), f > < m, v(θ) > +λ(θ) -n < m, Q(θ) n f >)
All involved quantities are positive, hence we can write

1 n log e θSn f dm = log λ(θ)+ 1 n log(< ϕ(θ), f > < m, v(θ) > +λ(θ) -n < m, Q(θ) n f >) Since lim θ→0 < ϕ(θ), f > < m, v(θ) >= 1
and since the spectral radius of Q(θ) is strictly less than λ(θ), it's easy to see that for θ small enough, we have

lim n→∞ 1 n log(< ϕ(θ), f > < m, v(θ) > +λ(θ) -n < m, Q(θ) n f >) = 0
In order to apply Gartner-Ellis theorem, we just have to show that Λ is differentiable function, strictly convex in a neighborhood of 0. Since λ is real-analytic, Λ is too. Computations from perturbation theory2 show that λ ′ (0) = φ dµ = 0 and λ ′′ (0) = σ 2 , so we have Λ ′′ (0) = λ ′′ (0)λ(0)-λ ′ (0) 2 λ(0) 2 = σ 2 > 0 and we can now apply the following local version of Gartner-Ellis theorem, whose proof can be found in lemma XIII.2 in [START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF] : Proposition 7. For all n ≥ 1, denote by P n a probability measure on some measurable space (Ω, T ), by E n the corresponding expectation operator and by S n a real valued random variable. Assume that on some interval [-θ Λ , θ Λ ], θ Λ > 0, we have

lim n→∞ 1 n log E n [exp(θS n )] = Λ(θ),
where Λ is a strictly convex continuously differentiable function satisfying Λ ′ (0) = 0.

Define ǫ + = Λ(θ Λ ) θ Λ > 0, ǫ -= Λ(θ Λ ) θ Λ < 0 and c(ǫ) = sup |θ|≤θ Λ {θǫ -Λ(θ)}.
Then c is a positive function, strictly convex on [ǫ -, ǫ + ], continuous, vanishing only at 0, and, for all 0 < ǫ < ǫ 0 = ǫ + , we have

lim n→∞ 1 n log P n (S n > nǫ) = -c(ǫ)
We now prove our Theorem 2. We have

Central

e it Sn √ n f dm =< m, P n it √ n (f ) >= λ( it √ n ) n < ϕ( it √ n ), f > < m, v( it √ n ) > + < m, Q( it √ n ) n f >
We just have to prove that

lim n→∞ λ( it √ n ) n = e -t 2 σ 2 2
But the Taylor's expansion says that in a complex neighborhood of 0

λ(z) = λ(0) + λ ′ (0)z + λ ′′ (0) 2 z 2 + z 2 η(z) = 1 + σ 2 z 2 2 + z 2 η(z)
where lim z→0 η(z) = 0. Then, a standard computation concludes the proof.

Application to uniformly expanding maps

The main application of our Theorem 1 will be to multidimensional piecewise uniformly expanding maps, in particular when we equip them with the space of the quasi-Hölder functions. This space, introduced by Keller [START_REF] Keller | Generalized bounded variation and applications to piecewise monotonic transformations[END_REF], developed by Blank [START_REF] Blank | Stochastic properties of deterministic dynamical systems[END_REF] and successfully applied by Saussol [START_REF] Saussol | Absolutely continuous invariant measures for multidimensional expanding maps[END_REF] and successively by Buzzi [START_REF] Buzzi | Absolutely continuous invariant probability measures for arbitrary expanding piecewise R-analytic mappings of the plane[END_REF] (see also [START_REF] Buzzi | Zeta functions and transfer operators for multidimensional piecewise affine and expanding maps[END_REF]) and Tsujii [START_REF] Tsujii | Absolutely continuous invariant measures for expanding piecewise linear maps[END_REF], reveals to be very useful to control the oscillations of a function under the iteration of the transfer operator across the discontinuities of the map. Moreover it verifies the algebraic assumption 3 in Section 2 above; it is not straightforward to replace this condition and in order to fit with the Hennion-Hervé theory, if one uses the more conventional spaces of bounded variation functions (see for instance [START_REF] Adl-Zarabi | Absolutely continuous invariant measures for piecewise expanding C 2 transformations in R n on domains with cusps on the boundaries[END_REF][START_REF] Boyarsky | Absolutely continuous invariant measures for piecewise expanding C 2 transformations in R n[END_REF][START_REF] Cowieson | Stochastic stability for piecewise expanding maps in R d[END_REF][START_REF] Liverani | Multidimensional expanding maps with singularities: a pedestrian approach, to appear in Ergodic Theory and Dynamical Systems[END_REF][START_REF] Poppe | Inadequacy of the bounded variation technique in the ergodic theory of higher-dimensional transformations in R n[END_REF]) or the Sobolev spaces [START_REF] Thomine | A spectral gap for transfer operators of piecewise expanding maps[END_REF], and this topic deserves to be investigated in the future3 .

Let us now recall the precise definitions of our system by following closely the assumptions imposed in [START_REF] Saussol | Absolutely continuous invariant measures for multidimensional expanding maps[END_REF]. Let M ⊂ R d be a compact subset with intM = M and piecewise C 1 boundary. We denote by d the Euclidean distance and by m the Lebesgue measure on R d . We can assume without loss of generality that m(M) = 1. For A ⊂ M and ǫ > 0, we denote B ǫ (A) = {x ∈ R d / d(x, A) ≤ ǫ}. Let T : M → M a measurable application, and suppose there exists 0 < α ≤ 1 such that for some small enough ǫ 0 we have :

(1) There are finitely many disjoint open sets

U i ⊂ M with m(M \ ∪ i U i ) = 0 such that for each i, T i := T | U i → M is C 1+α and can be extended on a neighborhood V i of U i to a C 1+α map T i : V i → R d such that B ǫ 0 (T i U i ) ⊂ T i (V i ). Moreover, each T i : V i → R d is injective with C 1+α inverse; (2
) There exists c > 0 such that for any i, and any x, y ∈ T (U i ) with d(x, y) ≤ ǫ 0 we have The last condition can be greatly weakened, but the condition in [START_REF] Saussol | Absolutely continuous invariant measures for multidimensional expanding maps[END_REF] is of a very abstract nature, and it's more easy to handle with this one when the boundaries of the U i are smooth. We define then the functional space on which acts the transfer operator. Let f ∈ L 1 (R d ). If A ⊂ R d is a Borel subset, we define the oscillation of f over A by osc(f, A) = ess sup

| det DT -1 i (x) -det DT -1 i (y)| ≤ c| det DT -1 i (x)|d(x, y) α ; (3) There exists s(T ) < 1 such that sup i sup x∈T i (V i ) ||DT -1 i (x)|| < s(T ); (4
x 1 ,x 2 ∈A |f (x 1 ) -f (x 2 )|
where the essential supremum is taken with respect to the product measure m × m on A × A. We get a lower semi-continuous and hence measurable function x → osc(f, B ǫ (x)). We set

|f | α = sup 0<ǫ≤ǫ 0 1 ǫ α R d osc(f, B ǫ (x))dx
We define

V α (R d ) = {f ∈ L 1 (R d ) / |f | α < ∞} and V α (M) = {f ∈ V α (R d ) / supp f ⊂ M} both endowed with the norm ||f || α = ||f || L 1 m + |f | α .
Adapting proofs from [START_REF] Keller | Generalized bounded variation and applications to piecewise monotonic transformations[END_REF], we can show that V α (M) is Banach space, with compact injection in L 1 (M). It's proven in [START_REF] Saussol | Absolutely continuous invariant measures for multidimensional expanding maps[END_REF] that V α (M) also a Banach algebra, and it's obviously a Banach lattice. So, if we want to apply our previous results to those maps, we are left to prove that the transfer operator for T acts on V α (M) and is quasi-compact of diagonal type. But Saussol proved in this context a Lasota-Yorke inequality (lemma 4.1 in [START_REF] Saussol | Absolutely continuous invariant measures for multidimensional expanding maps[END_REF]) which implies the quasi-compactness of the Perron-Frobenius operator. Hence, assuming that the system is mixing, we get the central limit theorem and the large deviations principle for bounded real observables in V α (M).

)

  Boundaries of U i are piecewise C 1 codimension one embedded compact submanifolds and we have η 0 (T ) < 1 whereη 0 (T ) = s(T ) α + 4s(T ) 1s(T ) Y (T ) γ d-1 γ d Y (T ) = supx∈R d i ♯{smooth pieces intersecting ∂U i and containing x} and γ d = π d/2 (d/2)! is the d-volume of the d-dimensional unit ball of R d .

  is the density of a T -invariant probability if and only if f ≥ 0, f dm = 1 and P f = f . B is a Banach algebra : there exists C > 0 such that for all f, g ∈ B we havef g ∈ B with ||f g|| B ≤ C||f || B ||g|| B ;(4) B is a complex Banach lattice : for every f ∈ B, we have f , |f | ∈ B;

	(3) (5) B is stable under P : P (B) ⊂ B; (6) P is a bounded operator on B, with spectral radius equal to one; (7) P is quasi-compact of diagonal type on B.
	Let us suppose now that we have a subspace B ⊂ L 1 (m), equipped with a norm || . || B such that (1) (B, || . || B ) is a complex Banach space with continuous injection B → L 1 (m); (2) Constant functions lie in B;

  Limit Theorem. By Levy's continuity theorem, it suffices to show that for all t ∈ R

	lim n→∞	e it Sn √ n f dm = e -t 2 σ 2 2

When ϕ ∈ B ⋆ belongs to the topological dual of B, we denote < ϕ, f >= ϕ(f ). The linear form f → f dm belongs to B ⋆ , and we denote it by m.

See corollaries III.11 and III.6 in[START_REF] Hennion | Limit theorems for Markov chains and Stochastic Properties of Dynamical Systems by Quasicompactness[END_REF].

In fact, if we check the proof, we only need the fact that B is Banach algebra and φ ∈ B to prove that the operators P z are well defined and holomorphic in z. So we can suppose the weaker assumption that φ is such that P z define a holomorphic family of bounded operators on B for z in a complex neighbourhood of 0.