Tamas Elteto
email: tamas.elteto@lri.fr

Nikolaus Hansen
email: nikolaus.hansen@lri.fr

Cecile Germain-Renaud

Pascal Bondon
email: pascal.bondon@lss.supelec.fr

T Éltető

N H A N S E N A

C G E R M A I N -R E N A U D A

P B O N D O N

Scalable structural break detection

Keywords:

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The so-called data deluge highlights two issues for data analysis: firstly, modelling non-stationary behaviour; second, scalability. This paper addresses a specific case for scaling model selection for structural break detection in time series. More precisely, building on the the AutoPARM procedure [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF] for piecewise autoregressive (AR) time seriesd e fi n e db yD a v i se ta l . ,w ep r o p o s ea n alternative optimisation technique for the model selections t e p . T h ep e r f o r m a n c ei m p r o v e m e n to f our alternative solution scales the size of the data sets by one to two orders of magnitude at constant computational budget.

In the Minimum Description Length context, whether in the perspective of machine learning or model fitting, the learning of fitting goal encompasses two tasks: the first one is to define an encoding of the data using some hypotheses/models; the second one is enacting an optimisation procedure to explore the associated hypothesis/model space. The optimisation problem for structural break detection is difficult, because of high dimensionality and a complex objective function. Then, evolutionary algorithms are a method of choice. Our approachc o n s i s t st or e v i s i tt h eo p t i m i s a t i o n strategy in the light of, first general advances in evolutionary computation, second the opportunities for a separable representation, and finally, an analysis of the limits of the AutoPARM optimisation method proposed in [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF].

The major contributions of this paper are the following 1. The single-level optimisation problem is decoupled into abileveloptimisation. Theupper level is the problem of finding the optimal number and location of theb r e a kp o i n t s . T h elower level optimises the autoregressive models given the number and locations of break points. 2. At the upper level, our optimisation strategy exploits thes t a t e -o f -t h e -a r tC M A -E S( C o v a r iance Matrix Adaptation -Evolutionary Strategy) instead of ar e l a t i v e l ys t r a i g h t f o r w a r dG enetic Algorithm. The associated representation addresses an important shortcoming of [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF]: the distance in the chromosomes space better maps to the distance in the model space. 3. The representation becomes scalable. More precisely, it scales linearly with the length of the data set, independently of the cost of the objective function. [START_REF] Bai | Multiple structural change models: asimulationanalysis.InEconometric theory and practice[END_REF]. At the lower level, finding the optimal autoregressive models is solved by an exhaustive search, and a significant amount of computations can be spared using the Yule-Walker method. 5. Regarding the problem of finding the optimal number and location of the break points, an alternative optimisation algorithm, the Covariance MatrixA d a p t a t i o nE v o l u t i o nS t r a t e g y (CMA-ES), and an adapted chromosome representation is proposed. This representation expresses the complexity of a piecewise autoregressive model better in the sense that a small change in the model maps to a small change in the representation. 6. The problem of finding the optimal autoregressive models int h el o w e rl e v e lo ft h ed e c o upled optimisation problem is solved by an exhaustive search,a n das i g n i fi c a n ta m o u n to f computations can be spared using the Yule-Walker method.

The content of the paper is organised as follows. Section 2 provides some background information for the piecewise stationary models, the Minimum Description Length (MDL) principle, its application to piecewise autoregressive model fitting, the actual model fitting procedure to be improved and briefly summarises the CMA-ES optimisation algorithm. Section 3 presents the main contributions of this paper. Section 4 provides the details of the numerical experiments comparing the performances of the different model fitting procedures. Section 5 concludes the paper, and the Appendix contains detailed numerical results of the performance comparisons.

Background

Stationarity plays a fundamental role in time series analysis: meaningful statistical inference about a stochastic process requires that its joint probability distribution does not change when shifted in time [START_REF] Brockwell | Time series: Theory and Methods[END_REF].

In many cases however, real-life observations are not compatible with the stationarity assumption. Trends and periodicity are not an issue, for that they can be eliminated through simple data transformations. A considerably more challenging issue is structural non-stationarity, where the time series is represented by a parametric model in which the parameters change values at unknown times called break points. In a nutshell, structural break analysis is required when the hypothesis does not stand that, at least at some relevant scale, the underlying source that creates the observed time-series either does not change at all, or evolves along a slow drift. Some important applications of structural break analysis can be highlighted as follows.

1. The objective is to explain the structural properties of the time series. For example, whether or not the system responds to a known/suspected change in the environment. In this case, the location of the change is very important.

2. The correlation or spectral structure of the stationary segments are analysed using classical methods. Then, better quality conclusions can be drawn compared to the case when the structural breaks are neglected. Furthermore, the change of thesec o n c l u s i o n sf r o mo n es e g m e n t to the other can also be studied. 3. The fitted model is used for generating "artificial" time series and the random series can be used as an input for a system study, a simulator or any other further analysis.

Objective 1 is demonstrated by an example in [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF] with a time series analysis focusing on the location of a change. The goal of the analysis is whether or nott h ei n t r o d u c t i o no ft h es e a t -b e l t legislation has an effect on the monthly deaths and serious injuries on UK roads. Another example in [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF][START_REF] Ombao | Automatic statistical analysis of bivariate nonstationary time series[END_REF] demonstrating the significance of the search for locations of structural changes is the analysis of speech signals, which are clearly non-stationary as a whole but can be broken into approximately stationary intervals corresponding to the phonemes.

The analysis of structural breaks is the main subject of [START_REF] Yamanishi | Dynamic model selection with its applications to novelty detection[END_REF],wherethegoalistodetectmalicious use of computer resources by looking for changes in the structure of consecutive best fitting models over a nontstationary time series. [START_REF] Yamanishi | Dynamic model selection with its applications to novelty detection[END_REF] refers to its approach as Dynamic Model Selection and uses an MDL coding scheme to find the best break point locations.

Objective 1 is also important in a structural break analysis of log returns of closing stock prices for 12 companies between 1993 and 2009, that is presented in [START_REF] Aue | Break detection in the covariance structure of multivariate time series models[END_REF]. According to [START_REF] Aue | Break detection in the covariance structure of multivariate time series models[END_REF], the detected changes can be readily associated with major historical events like the Asian financial crisis, the Russian financial crisis, September 11 attacks or the beginning of the recent financial crisis.

Objective 2 is demonstrated in [START_REF] Ombao | SLEX analysis of multivariate nonstationary time series[END_REF] by the analysis of the spectral structure of electroencephalograms (EEGs) recorded during an epileptic seizure. Since thes e r i e sc a n n o tb er e g a r d e dt ob e stationary, its Fourier transform cannot reveal useful information on the power conditions in various frequencies, consequently one should look for stationary segments in the time series whose spectral properties can be properly analysed.

Structural break detection related to Objective 3 is presented in [START_REF] Elteto | Discovering piecewise linear models of grid workload[END_REF][START_REF] Elteto | Towards non-stationary grid models[END_REF] where Grid workload measurements are analysed for structural breaks. The goal istoestablisharobustworkloadgenerator model that can be used in Grid architecture studies. Besides fitting a structural break model to the workload measurements, they also present empirical results on the robustness of the estimated parameters.

The rich literature on the approach of non-stationarity by structural break models (besides the above references, other recent results can be found for example in [START_REF] Bai | Estimating and testing linear modelsw i t hm u l t i p l es t r u c t u r a lc h a n g e s[END_REF][START_REF] Bai | Multiple structural change models: asimulationanalysis.InEconometric theory and practice[END_REF][START_REF] Basseville | The Detection of abrupt changes: Theory and applications[END_REF][START_REF] Brodsky | Nonparametric methods in change-point problem[END_REF]) underlines the importance of this subject. Once such a model is established,i ts i g n i fi c a n t l yi m p r o v e st h eq u a l i t y of the further analysis. However, to carry out the parameter fitting, one should face of costly and time consuming (e.g. in the order of 1-3 hours as our examples show) calculations as the structural break models (particularly when only a few prior assumptionsa r eu s e d )u s u a l l yh a v ev e r yr i c h parameter space. The exploration of the parameter space is ofc o u r s ee s s e n t i a lf r o mt h es a k eo f practical usefulness of these models. Improvements of theirp a r a m e t e rfi t t i n gm e t h o d sc a ng r e a t l y extend the scope of the models and it is therefore an importantt a s ko ni t so w nr i g h t .

Piecewise AR models

The piecewise AR model, formalised in (1), describes a finite length discrete-time locally stationary time series as segments of stationary time series that are concatenated. The stationary time series of each segment is assumed to be an AR process. The edgeso ft h es e g m e n t sa r er e f e r r e d to as break points. An important argument for focusing on piecewise AR models is that they are dense in the class of locally stationary processes with regular time varying spectrum and because efficient algorithms exist for fitting an AR model to a stationary process.

More precisely, let n be the length of the time series, m the number of break points. For all t such that 1

≤ t ≤ n, Y t = X t,j ,τ j-1 ≤ t<τ j , X t,j = γ j + φ j,1 X t-1,j + ...+ φ j,pj X t-pj ,j + σ j t,j , (1)
where τ j , j =1 ,...,m denotes the break point between the j th and (j +1) th AR process, τ 0 =1 , τ m+1 = n,( γ j ,φ j,1 ,...,φ j,pj ,σ 2 j )i st h ep a r a m e t e rv e c t o rc o r r e s p o n d i n gt oac a u s a la u t o r e g r essive process with order p j ,and t,j are independent, identically distributed random variableswithmean 0a n dv a r i a n c e1 .

The time series is explained as a set of consecutive segments,w h e rea na u t o re g re s s i v em od e ld escribes the evolution of the time series within each segment.S e c t i o n2 . 3w i l ld e s c r i b ec o r r e s p o n d i n g model selection task in the framework of Minimum DescriptionL e n g t ho p t i m i s a t i o n .

The MDL principle

The principal motivation of inductive inference is to searchf o rr e g u l a r i t yi nt h ed a t a . F r o ma machine learning perspective, what is needed is a decision rule -connected to a set of hypothesesthat uses this regularity. Given a training data, a hypothesis is then selected such that the decision rule generalises to some test data the best with respect to a loss function. Among the many kinds of loss functions, there is one, the so called log-loss, whose" o p t i m i s a t i o n "l e a d st oad e c i s i o nr u l e that best compresses the test data.

In a nutshell, the Minimum Description Length (MDL) principle [START_REF] Grünwald | The Minimum Description Length Principle[END_REF] also seeks for regularity in the data where regularity is interpreted directly as abilitytocompress. Thehopeisthat,oncesuc h compression scenario is found by selecting a proper hypothesis in the class, the inferred regularity generalises or predicts well. Furthermore, if the selected hypothesis is happen to be the "real" data source then the MDL principle leads to consistent model fitting methods.

[14] categorizes data compression data w.r.t. a model in fourma i nt ypeso fu ni v ersa lcodes. W e defer the discussion of the relevance of these categories to the problem at hand to Section 2.4.

1. NML (Normalised Maximum-Likelihood) code The code length is described as the sum of (a) the negative log-likelihood of the maximum likelihood estimate in a model class given the observation and (b) the complexity of the model class (when finite).

2. Bayes code. The code length is the weighted sum of log-likelihood of the data using the parametric distributions of the models in the model class M,w h e r et h ew e i g h ti sap r i o r distribution over the models in the class. 3. Two-part code. Given an observation y • Choose a parametric distribution F from a model class M of distributions and encode F.

• Encode y using F.

• Select the model F that gives the shortest total description length which is thes u mo f the description length of the model and the description length of y.

4. Sequential predictive. [START_REF] Grünwald | The Minimum Description Length Principle[END_REF] names it prequential: the encoding is done sequentially. After N observation steps, the maximum-likelihood model is selected for calculating the negative log-likelihood for the N +1 st observation.

We now turn to some consideration related to Information Theory, which will be useful in the next sections for motivating our choices of both a two-part code, and CMA-ES as the optimisation procedure.

First, the relationship between code lengths on one hand, and likelihood maximization on the other hand, is rooted in the standard interpretation of complete prefix codes as probability distributions. For the sake of completeness, we recall Kraft's inequality: given an alphabet S = {s 1 ,s 2 ,...,s n } and the corresponding code lengths l 1 ,l 2 ,...,l n of a uniquely decodable binary code (that is the number of bits needed to encode s k is l k ), n k=1 2 -l k ≤ 1. When equality holds, the code is called complete, and a natural distribution over S gives probability 2 -l k to symbol s k .

Conversely, the Shannon-Fano code corresponding to this distribution has code lengths l 1 ,l 2 ,...,l n .T h e r e f o r et h ed e s c r i p t i o nl e n g t hf u n c t i o nc a nb er e g a r d e das the negative log-likelihood of this probability source given the observation for which the code length is calculated.

The second important idea derives from the well-known fact that, under certain regularity conditions, the Hessian matrix of the log-likelihood function at the optimum (that is, at the maximum likelihood estimate) is a consistent estimate of the Fisher information matrix. This way, the observed information approximates the asymptotic covariance matrix of the maximum likelihood estimate, see for instance [START_REF] Efron | Assessing the accuracy of the maximum likelihood estimator: Observed versus expected fisher information[END_REF].

The maximum likelihood estimation motivates the idea that the optimisation of the fitness functions considered in certain MDL contexts can also be usedf o ra s s e s s i n gt h er e l i a b i l i t yo ft h e model selection. In the case of two-part codes with continuous parametrisation, the parameter vector minimising the negative log-likelihood of the model is per se a Maximum-Likelihood estimate, thus its Observed Information can be used for constructing confidence intervals for the parameter estimates. As a matter of fact, Section 14.4 of [START_REF] Grünwald | The Minimum Description Length Principle[END_REF] has a relevant note, that is, it is important to examine the code length estimates in the neighbourhood of the selected model. For continuous parametrisation this is explained by the Hessian matrix at the optimum, because this provides valuable information on the quality of the model selection.

The considerations so far can be summarized as follows:

• The MDL principle formulates the model (or model class) selection task as an optimisation problem of a fitness function representing a code length.

• If properly normalised, the code length function can be considered as a log-likelihood function of a probabilistic model.

• The neighbourhood of the optimum is also important with respect to the quality of the model selection. In particular, the Hessian matrix of the fitness function at the optimum is useful and we prefer optimisation algorithms that can provide an estimate for it.

The AutoPARM representation of piecewise AR process

[9] proposes the following code length function, which will be called AutoPARM-code in the following.

CL

F (y)=CL F ( F)+CL F (ê| F), (2)
where CL F ( F)denotesthecodelengthofthefittedmodel F and CL F (ê| F)denotesthecodelength of the time series using F.I no t h e rw o r d sa c c o r d i n gt o [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF] ,CL F (ê| F)i sw h a ti sn o te x p l a i n e db y the fitted model: the residuals of the AR process when the underlying piecewise AR model is F. More precisely:

CL F ( F)=log 2 m +(m +1)log 2 n + m+1 j=1 log 2 p j + m+1 j=1 p j +2 2 log 2 n j , (3)
where log 2 m and log 2 n are the length for encoding m and a segment size, respectively, log 2 p j refers to the encoding of the AR order p j , n j = τ jτ j-1 is the length of segment j,a n d pj +2

2 log 2 n j corresponds to the encoding of the AR parameters.

For the unexplained part,

CL F (ê| F) ≈ m+1 j=1 n j 2 log(2πσ 2 j ), (4)
where σ2 j is the Yule-Walker estimate (see [START_REF] Brockwell | Time series: Theory and Methods[END_REF] for details on this estimate) of the variance σ 2 j of the j th segment. [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF] shows that, when the number of break points is known, the relative break point locations selected by the two-part AutoPARM-code converge almost surely to their true values.

Function (3) enco des the descriptions of the length of the segments by (m+1) log 2 n.T h i se n c o ding implicitly bears the assumption that the length of the segments are independent and identically distributed, because the common log-likelihood (the code length) of lengths of the segments is the sum of the individual log-likelihoods.

We note, for further use, that the co de length estimation of the time series can be decomposed into a sum of m +1 t e r m s i n ( 2 ) w i t h e a c h t e r m r e f e r r i n g t o a s e g m e n t . C o n s e q u ently, if the code length estimate (2) is viewed as an objective function with AR orders as parameters, then the minimisation problem of this function is additively separable, that is, each segment only contributes summands that are independent of the others. More details on this observation are elaborated in Section 3.1.

Discussion

As said before, model selection covers both representation and optimisation. The resulting palette of possible methods is thus large, amounting to combining two relatively independent research directions. [START_REF] Grünwald | The Minimum Description Length Principle[END_REF] shows both theoretically and in practice that the accuracy of the parameter selection method highly depends on the quality of the compression method, i.e. the calculation of the code length. This section will briefly discuss alternative representations, both in the perspective of the four universal coding methodologies presented in 2.2,andwithinAutoP ARM-code. Ho w ev er, exploring these code alternatives goes beyond the scope of the present paper, which focuses on the performance of the alternative optimisation algorithm described in Section 3.

Considering the coding methodology, AutoPARM-code is a two-part coding scheme and [9] contains a proof of consistency of this approach. We note thati fo n ei sw i l l i n gt oc h o o s ea m o n g alternative encodins (to be detailed below) the same consistency property have to be proved as well.

[13] presents Normalised Maximum-Likelihood codes for the AR and ARMA model classes, and give empirical evidence of high efficiency for model selection. However, the calculation of the complexity of the model class involves Monte-Carlo integration, which is not a good candidate for scalability. Moreover, we found it difficult to estimate the values above AR order 6, an observation already presented in [START_REF] Giurcaneanu | Estimation of ar and arma models by stochastic complexity[END_REF]. There is an alternative encoding,whic histobebasedontheC onditional Normalised Maximum-Likelihood code of [START_REF] Rissanen | Conditional nml universal models[END_REF] for generalisedG a u s s i a nr e g r e s s i o nf a m i l i e sf o rA R models. Finally, [START_REF] Rissanen | Order estimation by accumulated prediction errors[END_REF] presents an ARMA model selection method based on a predictive MDL code. The model parameters have to be fitted step-by-step, therefore this method would be too slow for our scenarios with many long AR segments.

Fo cusing on the AutoPARM-co de as describ ed by ( 3) and ( 4), some improvements, that would probably improve the rate of convergence in the asymptotics,m i g h tb ea sf o l l o w s .

First, the number of break points is encoded using log 2 m bits. However, the prefix encoding of an arbitrary integer needs more bits, for example, 2 log 2 m is a better estimation. The same applies to the encoding of the AR order p j .S e c o n d ,t h ee n c o d i n go ft h eb r e a kp o i n tl o c a t i o n sb y (m +1)log 2 n bits is not the shortest one. Assuming that the length of all segments is greater than 1, the break point locations could be encoded using log 2 n m bits. Last, but not least, according to (3), the number of bits to encode the model parameters in the j th segment is (p j +2) 1 2 log 2 n j . Theorem 10.1 in [START_REF] Grünwald | The Minimum Description Length Principle[END_REF] states that under certain conditions (e.g. on the parameter set of an exponential family with k parameters), the regret (the amount of needed extra code length compared to the shortest maximum-likelihood code) of a two-partc od eh a sap a r t i c u l a rf o r m ,w h i c hi s asymptitoctically k 2 log n 2π if n,t h el e n g t ho ft h et i m es e r i e s ,g o e st oi n fi n i t y . T h ee x t e n s i o no ft h i s Theorem to the case of AR models would theoretically justify the above outlined approach, but to the best of our knowledge, this has not yet been done in the literature.

AutoPARM-GA

AutoPARM carries out the model fitting task by a genetic algorithm, which will be named AutoPARM-GA in the rest of the paper. Here, we briefly discuss its chromosome representation and some other details of the optimisation procedure.

Definition 1. Ac h r o m o s o m ei sav e c t o ro fi n t e g e r sδ

=(δ 1 ,...,δ n ) of length n where δ t is defined as

δ t =
-1, if no break point at t, p j , if there is a break point at t and the AR order for the j th segment is p j .

AutoPARM-GA follows the island strategy [START_REF] Whitley | T h ei s l a n dm od e lg e n e t i ca l g o r i t h m : On separability, population size and convergence[END_REF]: rather than running the search in the population as a whole, it considers N islands with a fixed number of chromosomes in each island. After generating the offspring, the new and old chromosomes are compared and the best chromosomes are kept independently in each island. After a certain number of generations, migrations occur among the islands. Together with the migrations, the best chromosome (representing the model with the shortest description length among all chromosomes seen so far) is updated. The stopping criterion is that the best chromosome does not change after a certain number of consecutive migrations, or the number of migrations reaches a limit. In the following, this particular kind of genetic algorithm will be termed NI.

Afi r s tl i m i t a t i o no ft h er e p r e s e n t a t i o n so fD e fi n i t i o n1i st h ea s s o c i a t e dc o m p u t a t i o n a lc o mplexity. The crossover and mutation operations step throught h ew h o l ec h r o m o s o m e( o rap a i ro f chromosomes). Since the length of a chromosome is equal to thel e n g t ho ft h ed a t a s e tt ob es e gmented, and the crossover and mutation operations are considering the elements of the chromosome vector almost one-by-one, the number of computation steps int h e s eo p e r a t i o n sa r ep r o p o r t i o n a l to the length of the dataset. Thus, the number of computation steps per objective function evaluation during the optimisation scales linearly with the length of the dataset even when the cost of the objective function evaluation and the total number of function evaluations are disregarded. This is a limitation when the algorithm is going to be applied to a large series. The island model offers opportunity for speedup through parallelisation. Thes pe e d u pf a c t o rw o u l dber e l a t e dt ot h e number of islands or chromosomes depending on the level of parallelisation.

Am o r ef u n d a m e n t a l ,l i m i t a t i o no fA u t o P A R M -G Ai sr e l a t e dt othe chosen combination of encoding and the operators in the algorithm. Particularly for long samples, the exact location of a break point is not at all critical: a small shift of the break point location implies only a small change in the overall piecewise AR model. For instance,i nt h ee x p e r i m e n t a ls t u d y[ 1 1 ] , bootstrapped samples from segmented large time series showed some small variability in break point locations, but consistent behaviour.

The issue with the representation of Definition 1 is that the change of the location of a break point involves (at least) two operations on a chromosome. Assume, for example, that there is a break point at t δ t = p and δ t+1 = -1.

Shifting this break point from t to t +1 involves two independent operations: deleting the current break point, i.e., δ t = -1; and the new break point should be established, i.e., δ t+1 = p.I tc a nb e expected, that if the chromosome representation is modified,s u c ht h a tt h ec h a n g eo fb r e a kp o i n t location can be performed in one simple operation, then the efficiency of the optimisation should improve.

The CMA-ES optimisation algorithm and the probabilistic interpretation of its output

The goal of CMA-ES [START_REF] Hansen | The CMA evolution strategy: A comparing review[END_REF] is to minimise an objective function f : R n → R.T h eo n l ya c c e s s i b l e information on f are function values at evaluated search points. CMA-ES performs a stochastic search in the parameter space by adapting the direction and step size using inferred dependencies between the parameters. These dependencies between all parameter pairs are learned by updating ac o v a r i a n c em a t r i xo fam u l t i v a r i a t en o r m a ls e a r c hd i s t r i b ution.

After generating the next offspring of λ samples in the parameter space, the best1 µ of them are selected and are used to update the mean and the covariance matrix of the sampling distribution. The updated mean is a weighted average of the selected samples. The covariance matrix is updated by averaging the old covariance matrix estimate, the covariance matrix estimate coming from the evolution path history (rank-one update) and the covariancem a t r i xe s t i m a t ef r o mt h es e l e c t e d samples (rank-µ update).

When the CMA-ES algorithm has converged, the mean of the last update of the multivariate normal distribution is as close to the optimum as possible depending on the stopping rule, while the covariance matrix approximates, in case of a smooth objective function, the second order smoothness properties of the objective function around the optimum. If an e g a t i v el o g -l i k e l i h o o df u n c t i o no fa parametric probability distribution satisfying certain regularity requirements is optimised such that the maximum likelihood estimation is consistent and asymptotically normal, then the covariance matrix estimate of CMA-ES will also be an estimate of the observed information (that is the Hessian matrix of the negative log-likelihood function, see [START_REF] Efron | Assessing the accuracy of the maximum likelihood estimator: Observed versus expected fisher information[END_REF]) apart from a constant term. Thus, CMA-ES satisfies the 3 rd requirement on the list of considerations ath the end of Section 2.2.

Bilevel optimisation and the role of CMA-ES

This section introduces two alternative optimisation methods for fitting a piecewise AR model. The first one is the bilevel method presented in Section 3.1. Empirical evidence (see Section 4.2) shows that it has superior performance over AutoPARM-GA, which is a single-level method using the NI optimisation algorithm. The second one comes from our observation that its chromosome representation has a significant limitation in that multipleo p e r a t i o n sa r en e e d e df o re x p l o r i n g models in the parameter space that are otherwise close to eacho t h e r . S e c t i o n3 . 2i n t r o d u c e st h e application of CMA-ES and the corresponding chromosome representation that do not have this limitation. The additional performance gain will be presented in section 4.3.

The bilevel method

In AutoPARM-GA, the break point locations and the AR orders are optimised at the same time.

In the following we shall refer to this as the single-level method. The search space here is rather large, and consequently it can be costly to find the optimal solution. We saw in Section 2.1 that the optimisation problem is additively decomposable. Because of this, the problem can be formulated as a special bilevel optimisation problem (a recent overviewo nb i l e v e lo p t i m i s a t i o ni sg i v e ni n [START_REF] Colson | An overview of bilevel optimization[END_REF] ) .

We first present the decomp osition idea as follows. The ob jective function is the code length function of [START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF]. We make the parameters in CL F (y)e x p l i c i tb y

CL F (y)=f (m, κ, p),
where m is the number of break points, κ denotes the vector of relative break point locations, κ =(κ 1 ,...,κ m )w h e r eκ j = τ j /n,a n dp ∈{0, 1,...,p max } m+1 denotes the vector of AR orders of the segments.

Then, the optimisation problem can be formalised as min m,κ,p f (m, κ, p).

In the bilevel setup, we formulate this optimisation problema s min

m,κ {min p f (m, κ, p)}. (5)
The graphical interpretation of the proposed optimisation method is shown in Figure 1.

Upper level: localising the break points

As it was indicated at the beginning of the section, the upper level optimisation is the selection of the number and locations of the break points. We address this problem in two different ways. First, we used the same NI algorithm as in AutoPARM-GA, described in 2.5, in order to investigate the contribution of the bilevel setting in the original scenario. We keep the problem encoding in the upper level essentially unchanged: each locus in the chromosome δ t ,f o rt =1 ,...,n,d e t e r m i n e s whether a break point is realized, which is the case if δ t = -1. The AR order is provided by the lower level. Second, we encode each break point location as an integer. Then u m b e ro fi n t e g e r sd e fi n e st h e number of break points and also the dimension of a search problem with a fixed number of break points. Changing the number of break points means to change the search space dimension of this search problem. In this approach, we repeatedly choose a break point number and run CMA-ES, as explained in Section 3.2.2.

This problem can be solved using a generic optimisation method. In this paper, we used the same NI algorithm as in AutoPARM-GA, described in 2.5, in order to highlight the specific contribution of the bilevel setting.

In any case, the important information in a chromosome of the NI algorithm is whether the elements δ t are -1ornotfort =1,...,n,w h e r en is the length of the chromosome. The AR order, which is important in the crossover and mutation operations as it might impose constraints on the length of the segment, is not randomly chosen but it is provided by the lower level.

Lower level: determining the AR models within the segments

When a particular number and location of break points are selected, then a lower level optimisation is performed. The objective function detailed in (3) and ( 4) can be minimised segment by segment. Let

f * (κ j ,κ j+1 ,p j )=log 2 n +log 2 p j + p j +2 2 log 2 n j + n j 2 log(2πσ 2 j ), (6)
we have

f (m, κ, p)=log 2 m + m j=0 f * (κ j ,κ j+1 ,p j+1 ),
and then,

min p f (m, κ, p)=log 2 m +min p1 f * (κ 0 ,κ 1 ,p 1 )+...+m i n pm+1 f * (κ m ,κ m+1 ,p m+1 ). (7)
Fo cusing on segment j,t h ec od el e n g t he s t i m a t ef o re a c ho r d e rp j ,0≤ p j ≤ p max is calculated. The order giving the shortest code length estimate for segment j is selected and the lower level optimisation procedure steps over to segment j +1. The constraint on the length of the segment (in order to avoid high order AR parameter fitting on short segments) is taken into consideration by limiting p max .

The most expensive part of the code length estimate for segment j in terms of computational complexity is the last summand in [START_REF] Brockwell | Time series: Theory and Methods[END_REF]. In particular, the calculation of the variance estimate σ2 j is the most demanding. The variance estimates are calculatedb yt h eD u r b i n -L e v i n s o na l g o r i t h m which fit AR models of successively increasing orders 1, 2,...,p max to the data, see [START_REF] Brockwell | Time series: Theory and Methods[END_REF] Section 8.2. The complexity of this algorithm is O(p 2 max ).

Piecewise AR model fitting using CMA-ES

This Section shows how CMA-ES can be incorporated into the piecewise AR model fitting framework. There are two important issues that need further considerations: the varying number of parameters and the discrete versus continuous parameter space.

The problem formulation for CMA-ES

The following alternative chromosome representation, thatc o n t a i n st h es a m ei n f o r m a t i o na s Definition 1, is established for the optimisation procedure using CMA-ES. The chromosome representation of Definition 2 is an ordered pair whose first component is the number m of break points. The dimension of the vector in the second component of the chromosome is m.T h e r e f o r e ,t h eo b j e c t i v ef u n c t i o n ,t h a tc a l c u l a t e st h ec o de length estimate for an arbitrary chromosome, has varying number of input parameters depending on the number of break points. However, CMA-ES is designed for objective functions with real vector arguments of afi x e dd i m e n s i o n . C o n s e q u e n t l y ,t h eo p t i m i s a t i o np r oc e d u r eu s i n gC M A -E Ss h o u l dbec a r r i e do u t with a fixed number of break points. That is, several optimisation steps using different number of break points should be done.

Definition 2. Ac h r o m o s o m ei sa no r d e r e dp a i ro fa ni n t e g e ra n dav e c t o ro fo r dered real numbers in

The relative locations of the break points are discrete because of the discrete-time AR model, but CMA-ES is a continuous-space optimisation method therefore the real parameters are rounded for the evaluation in the objective function. In other words,thecon tin uous-spaceobjectiv efunction is piecewise constant.

The optimisation procedure

This Section gives a high level description of the optimisation procedure; Table 1 gives the pseudo-code. Denote by f : N→R + the following function:

f (m)= min κ∈(0,1) m min p∈{0,...,pmax} m+1 f (m, κ, p).
The value of f (m)i so b t a i n e db yfi t t i n gp i e c e w i s ea u t o r e g r e s s i v em o d e l sw i t h m break points to the time series. The model fit uses CMA-ES to find the optimal location of the break points by searching in the parameter space (κ 1 ,...,κ m ), where the constraint 0 <κ 1 <...<κ m < 1i s enforced by penalisation.

The basic idea of the search algorithm is that starting from m pointer =2 ,w ea l w a y st e s tt h r e e cases with m pointer -1, m pointer and m pointer +1 number of break points in one loop. When the estimated code length shows a decreasing trend then m pointer is increased. E.g. if

f (1) > f (2) > f (3),
then the number of analysed break points is doubled, m next pointer =2• 2=4 ,a n dt h en e wf (•)a r e tested in the new loop. (That is, we test the neighbourhood of 2, 4, 8, 16, . . . )

Oppositely, when the estimated code length values show an increasing trend then m pointer is reduced. E.g. if we found decreasing trend for 8, but increasing for 16, that is The number of break points is selected to be the one which givesthesmallestcodelengthestimation L m among all that were tested. Some technical issues should be discussed here. Firstly, ther o u n d i n gi n s i d ef (m, κ, p)m a y divert CMA-ES for short time series. However, in the range of 1k -2k or above, this effect is not apparent. Next, given f (.)i su n i m o d a l( a n dn o i s e l e s s ) ,w ec a ne a s i l yp r o v et h a tt h ea l gorithm will find the global optimum for the number of break points fast. Otherwise, only a local optimum is guaranteed to be found. In fitting piecewise AR models to empirical time series, unimodality seems to be a prerequisite for the overall goal to be reasonable: the underlying hypothesis is that there is some physical process driving the ruptures. For instance, our experiments in bootstrapping segmented time series [START_REF] Elteto | Discovering piecewise linear models of grid workload[END_REF] exhibited correlation between times series corresponding to the activity of distributed agents, corresponding to a common driving mechanism. Finally,we note, that in order to avoid models with generally too short segments, an upper limit on the possible number of break points was set to 1 100 of the length of the series.

f (7) > f (8) > f (9),
Table 1: Pseudo-code of the search for the optimal number of break points.

1. mpointer =2,m upper limit = series length 100 , m low =1,m high = m upper limit .

2.

CL-1 = f (mpointer -1), CL0 = f (mpointer), CL1 = f (mpointer +1) 3.

Store the f (•)v a l u e st o g e t h e rw i t ht h e#o fb r e a kp o i n t si nT a b l eA.

4.

If not (CL-1

≤ CL0 ≥ CL1)o r( CL-1 ≤ CL0 ≤ CL1) then 4.1 m low = mpointer. 4.2 If CL-1 >CL0 >CL1 then mpointer =2mpointer. 4.3 While not(CL-1 ≥ CL0 ≤ CL1)a n d( m low <m high ) 4.3.1 CL-1 = f (npointer -1), CL0 = f (mpointer), CL1 = f (mpointer +1) 4.3.2 Store the f (•)v a l u

Performance comparison of the single-and bilevel methods

Figures 2 to 7 show that the single-level and bilevel optimisation methods explore the parameter space with similar efficiency after an initial start up phase. Furthermore, at the b eginning of the optimisation, the bilevel method is significantly more efficient as the range of explored code length values is lower compared to its single-level counterpart. This effect is particularly apparent in Figure 2: the bilevel method reaches the 2100-2150 range of code length values at ca. 100 seconds of CPU running time, while the single-level method reaches the same range at ca. 300 seconds. After ca. 400 seconds, the two historical curves run together.

We also note, that b oth exp eriments use the NI algorithm, whose running time scales together with the data length and this is also reflected in these performance figures.

Performance comparison of the NI and CMA-ES algorithms

The straightforward way to find the number of break points minimising the code length is to test all possible values m =0 ,...,m max .H o w e v e r ,d e p e n d i n go nt h ev a l u eo fm max this might be rather time-consuming. If instead our search algorithm is used, then we will be able to find the optimal number of break points m when f (•)isunimodal. Inthiscase,theoptimaln um berofbreak points can be found within log 2 m max steps by halving the region to be explored in each test of the f (•)v a l u e s .

The number of chromosomes, λ,u s e db yC M A -E Ss h o u l dbec h o s e nc a r e f u l l y ,k e e p i n ge ffi c i e n c y in mind. According to our practical experience, a moderate (e.g. in the order of 40-60) population size is sufficient for finding the optimal number of break points, because the exact break point locations are not yet needed for this task. Instead, after then u m b e ro fb r e a kp o i n t sh a sb e e n selected, a second optimisation run can be performed with an increased population size (e.g. in the order of 500-1500) so that the break point location estimations are refined.

As a further enhancement of the model fitting method, after theo p t i m i s a t i o nh a sc o n v e r g e d ,i t is restarted with increased population sizes [START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF] until a limit on the cumulative number of function evaluations is reached.

According to our performance tests, CMA-ES needs considerably smaller number of function evaluations than the NI algorithm. Figure 8 shows the performance comparison of the bilevel method with NI and CMA-ES through different number of functione v a l u a t i o n s . O n ec a ns e ei n Figure 8, that though initially the NI algorithm performs better, after the start up phase, the code length estimates of the bilevel method with the CMA-ES algorithm can reach the range of 2000-2080 in ca. 14K function evaluations, while the bilevel method with the NI algorithm can reach the same range after ca. 520K function evaluations. Considering the slight performance gain that was experienced for the same model in Section 4.2 we can conclude that, in this test example, the performance of the bilevel method with the CMA-ES algorithm is superior by more than 1 order of magnitude over the single-level method with the NI algorithm(AutoP ARM-GA).Thisperformance gain appears to be even more appealing if we remember that the search space dimension for the NI algorithm, as proposed by [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF], is proportional to the length of the time series whereas CMA-ES is independent of this length.

Average CPU times per function evaluation

This section presents empirical comparison results concerning the average CPU time per function evaluation for all three model fitting methods and highlight the main effects that lead to our performance observations. Table 2 shows the results. The values were obtained by dividing the total CPU running times by the total number of function evaluations for each run, where a run is defined by the triple test case, algorithm, method. The most significant effect that determines the CPU time per function evaluation is, of course, the length of the series. Table 2 clearly demonstrates this, as the average CPU time for the PAR many model is about twice as much as for the other models in all the columns. The reason for this is, that for this model, the length of the series is 2048, while it is 1024 for the rest of the models.

The effect of the difference between the single-level method and the bilevel method can also be studied by comparing the first and second columns of Table 2. The bilevel method performs an exhaustive search over all possible AR orders up to p max =2 0 ,w h i l et h es i n g l e -l e v e lm e t h odd oe s not. This increases the CPU time per function evaluation for the bilevel method by at most a factor of ca. 4 as it can be seen in Table 2.

Another effect is the running time overhead of the optimisation procedure. This can be quite significant as it can be seen by comparing the second and third columns of Table 2, where the objective functions are exactly the same (in fact, the same C implementation was used in the tests), but the optimisation algorithms are different. Though a large part of the NI algorithm was implemented in C and the CMA-ES algorithm was completely implemented in Matlab, the comparison shows that the CPU time, and thus the overhead of the optimisation algorithm, is smaller for CMA-ES.

Finally, we note that the rows of Table 2 show different values for the different models. The main reason for this is that the contribution of the optimisation algorithm and the contribution of the objective function evaluation to the average CPU time depend very much on the optimisation path (e.g. the number and locations of the break points, the ARo r d e r s ,e t c . ) ,w h i c hi sd i ff e r e n tf o r each model.

Accuracy

So far we showed, that the minimal code length estimate, and thus the optimal model fit, is more efficiently approached by our proposed solution. Table 4 shows the parameter values such as the number of break points, ranges of break point locations and distribution of AR orders selected by the different runs for our illustrative scenario, the PAR many test case (Appendix B gives similar tables for the other test cases). For the sake of brevity, we show only the distribution of AR orders collected from runs with the most frequent number of break points. Also, in case of many segments, we show only the results for a few of them. The fitted values are quite similar for each of the three procedures. We note, however, that the bilevel method with the NI algorithm produces consistently slightly better solution than the single-level algorithm and the CMA-ES algorithm is slightly even more precise than the two others.

Technical details

The length of the series for all the benchmark examples were 1024, except for the first one, where it is 2048. The maximum of the AR orders to be considered was set to 20.

The dimensionality of the problem for the NI algorithm with single-level optimisation is 1000-2000, where each single parameter can take 22 different values. The dimensionality of the NI algorithm with bilevel optimisation is the same, however with effectively only two different values. The dimensionality of the problem for CMA-ES with bilevel optimisation depends on the number of break points explored on the optimisation path (it is at most [START_REF] Efron | Assessing the accuracy of the maximum likelihood estimator: Observed versus expected fisher information[END_REF][START_REF] Elteto | Towards non-stationary grid models[END_REF][START_REF] Elteto | Discovering piecewise linear models of grid workload[END_REF][START_REF] Giurcaneanu | Estimation of ar and arma models by stochastic complexity[END_REF][START_REF] Grünwald | The Minimum Description Length Principle[END_REF][START_REF] Hansen | The CMA evolution strategy: A comparing review[END_REF][START_REF] Ombao | SLEX analysis of multivariate nonstationary time series[END_REF][START_REF] Ombao | Automatic statistical analysis of bivariate nonstationary time series[END_REF][START_REF] Rissanen | Order estimation by accumulated prediction errors[END_REF][START_REF] Rissanen | Conditional nml universal models[END_REF][START_REF] Ros | A simple modification in CMA-ES achieving linear time and space complexity[END_REF].

The parameters for the NI algorithm are the parameters proposed in [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF] except that the number of islands is increased to 100 (from 40), the population size on each island is increased to 100 (from 40) and the limit on the number of migrations is set to a large value so that the optimisation always stops due to that the overall best chromosome has not changed for 10 consecutive migrations. Our reasons for changing these parameters (e.g. the approximately 4x increase in the overall size of the population) is that, according to our experience, this speeded up tho optimisation process because of the following main reasons:

• It turned out, that the larger population resulted generallysomuchgainintheobjectivefunction value that the same level could be reached using 90%-95% number of function evaluations compared to an optimisation with the original population sizes.

• It turned out that the original population sizes sometimes led the optimisation stuck into a sub-optimum. Then, neither the crossover nor the mutation operations moved the process out of this sub-optimum (at least on the intended timescale ofo u rr u n s ) .

Regarding CMA-ES, we use the default optimisation parameters except that the covariance matrix estimate is set to be diagonal [START_REF] Ros | A simple modification in CMA-ES achieving linear time and space complexity[END_REF]. The population size is set to 50 during the phase when

NI algorithm

CMA-ES algorithm single-level method bilevel method relative frequencies where different number of break points were found 4> break points 0.0000 0.0000 0.0000 4b r e a kp o i n t s 0 . 0 0 0 0 0 . 0 0 0 0 0 . 0 5 9 1 5b r e a kp o i n t s 0 . 0 2 4 6 0 . 0 0 4 9 0 . 0 2 4 6 6b r e a kp o i n t s 0 . 3 8 9 2 0 . 1 2 3 2 0 . 0 2 9 6 7b r e a kp o i n t s 0 . 5 4 6 8 0 . 7 9 3 1 0 . 8 0 3 0 8b r e a kp o i n t s 0 . 0 3 9 4 0 . 0 7 8 8 0 . 0 6 9 0 9b r e a kp o i n t s 0 . 0 0 0 0 0 . 0 0 0 0 0 . 0 1 Segment 1 order <10 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 1 0.9882 0.9882 0.9882 order 2 0.0118 0.0118 0.0118 order >20 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 2 order 0 0.0235 0.0118 0.0118 order 1 0.8588 0.9765 0.9765 order 2 0.0941 0.0118 0.0118 order 3 0.0118 0.0000 0.0000 order 4 0.0118 0.0000 0.0000 order >40 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 the number of break points are determined and it is set to 1000 during the phase when the break point locations are refined. When the CMA-ES optimisation is restarted, the population sizes are increased by a factor of 2. We randomly generate 200 time series samples following the prescribed models for each experiment. Then, we fit piecewise AR models using the single-levela n db i l e v e lo p t i m i s a t i o nm e t h o d s with the NI algorithm and using the bilevel method with the CMA-ES algorithm. The code length values of the model fittings are recorded for each experiment during the optimisation and the observations are stored in historical tables. One historical entry contains the CPU-time, the number of function evaluations and the code length value at that instant. Approximately, 20K-30K historical entries are saved. Finally, each historical curve is interpolated and sampled at certain time instants and the quantile statistics of the CPU seconds and the number of function evaluations are collected and plotted at ca. 0.003 (lower whiskers), 0.25 (bottoms of the boxes), 0.5 (middle points of the boxes), 0.75 (tops of the boxes) and at ca. 0.993 (upper whiskers).

Appendix B provides the numerical details for all experiments. Generally we observe that the bilevel method using the NI algorithm shows the same or superior performance over the single-level method and the final objective function values after convergence are always better for the bilevel optimisation method.

Extending the segmentation model

The code length function of the residual part (4) uses an approximation of the maximum likelihood estimation by using the Yule-Walker estimate. For larges e g m e n ts i z e s ,t h i sa p p r o x i m a t i o ni s satisfactory. However, the Yule-Walker estimator might be replaced by another estimator for several reasons. On the one hand, particularly for short segments, the maximum likelihood estimate might be a better choice. On the other hand, the use of the Durbin-Levinson algorithm is restricted to AR models. In the case of autoregressive and moving average (ARMA) models or nonlinear models like conditionally heteroscedastic models, effective closed-form parameter estimation methods similar to the Yule-Walker estimator do not exist in general.

In these cases, one should use numerical optimisation techniques to obtain relevant estimates. In the framework of the NI algorithm, this means that for each chromosome and for each segment one should run a second-(or third-) level optimisation and then use the obtained negative log-likelihoods for calculating the code length estimate. This can be rather time consuming and therefore the selection of an efficient optimisation algorithm is essential. In such cases, CMA-ES is a promising alternative especially for optimisation problems involving a large number of parameters and difficult optimisation functions.

Conclusions

In this paper, we consider the problem of fitting a piecewise autoregressive model to a time series using the MDL principle. We are dealing with two main aspects of the fitting procedure. The first one is the method that explores the parameter space of piecewise autoregressive models. The second one is the algorithm that actually carries out the optimisation.

Our study is based on simulations on five test examples given in [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF]a n do nas i x t he x a m p l e which is motivated by practical experience from [START_REF] Elteto | Discovering piecewise linear models of grid workload[END_REF] with a longer series and larger number of break points.

We show that the p erformance of the single-level metho d usingtheNIalgorithmcanbeimproved to a moderate extent by replacing it with a bilevel method. An empirical comparison of the total CPU-time of the single-level and bilevel method using the NI optimisation algorithm is presented. We find that the CPU-time until convergence is the same or slightly smaller in favour of the bilevel method compared to its single-level counterpart.

The same optimisation problems are used to compare the bilevel method using the NI and CMA-ES algorithms. In these cases, the number of upper level objective function evaluations are presented. Here, the upper level objective function incorporates the (lower level) optimisation of the autoregressive models. We find that CMA-ES algorithm has at least one order of magnitude faster convergence in the examples. Therefore, based on our empirical experience, we propose the bilevel method using CMA-ES for piecewise autoregressive model fitting purposes.

NI algorithm

CMA-ES algorithm single-level method bilevel method relative frequencies where different number of break points were found 1b r e a kp o i n t 0 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 2b r e a kp o i n t s 0 . 4 7 0 0 0 . 3 3 5 0 0 . 4 1 5 0 3b r e a kp o i n t s 0 . 5 2 0 0 0 . 6 4 5 0 0 . 5 8 5 0 4b r e a kp o i n t s 0 . 0 1 0 0 0 . 0 2 0 0 0 . 0 0 0 0 other 0.0000 0.0000 0.0000 relative locations of break points in cases where 2 break points were found average (standard deviation) break point 1 0.3591 (7.2e-2) 0.3574 (7.1e-2) 0.3545 (7.0e-2) break point 2 0.6801 (7.2e-2) 0.6827 (7.5e-2) 0.6817 (6.9e-2) relative frequencies of AR orders when 2 break points were found Segment 1 order <10 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 2 0.9688 0.9688 0.9792 order 3 0.0208 0.0208 0.0104 order 4 0.0104 0.0104 0.0104 order >40 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 2 order <20 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 2 0.9479 0.9583 0.9688 order 3 0.0521 0.0313 0.0208 order 4 0.0000 0.0104 0.0104 order >40 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 3 order <20 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 2 0.9583 0.9688 0.9688 order 3 0.0313 0.0208 0.0208 order 4 0.0104 0.0104 0.0104 order >40 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 NI algorithm CMA-ES algorithm single-level method bilevel method relative frequencies where different number of break points were found 1b r e a kp o i n t 0 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 2b r e a kp o i n t s 0 . 9 9 5 0 0 . 9 8 0 0 0 . 9 9 0 0 3b r e a kp o i n t s 0 . 0 0 5 0 0 . 0 2 0 0 0 . 0 1 0 0 other 0.0000 0.0000 0.0000 relative locations of break points in cases where 2 break points were found average (standard deviation) break point 1 0.4999 (6.2e-3) 0.5002 (6.3e-3) 0.5001 ( 0.0051 0.0051 0.0051 order >60 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 2 order <10 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 1 0.9796 0.9898 0.9898 order 2 0.0102 0.0102 0.0102 order 3 0.0051 0.0000 0.0000 order 4 0.0051 0.0000 0.0000 order >40 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 3 order <10 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 1 0.0051 NI algorithm CMA-ES algorithm single-level method bilevel method relative frequencies where different number of break points were found 1b r e a kp o i n t 0 . 0 4 5 0 0 . 0 2 0 0 0 . 0 3 0 0 2b r e a kp o i n t s 0 . 8 7 0 0 0 . 7 8 5 0 0 . 7 8 5 0 3b r e a kp o i n t s 0 . 0 8 5 0 0 . 1 9 5 0 0 . 1 8 0 0 4b r e a kp o i n t s 0 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 5 0 other 0.0000 0.0000 0.0000 relative locations of break points in cases where 2 break points were found average (standard deviation) break point 1 0.2400 (6.4e-2) 0.2395 (6.4e-2) 0.2438 (6.9e-2) break point 2 0.5525 (7.9e-2) 0.5571 (7.7e-2) 0.5565 (8.3e-2) relative frequencies of AR orders when 2 break points were found Segment 1 order <10 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 NI algorithm CMA-ES algorithm single-level method bilevel method relative frequencies where different number of break points were found 0b r e a kp o i n t 0 . 0 0 5 0 0 . 0 0 5 0 0 . 0 0 0 0 1b r e a kp o i n t s 0 . 9 9 5 0 0 . 9 8 5 0 1 . 0 0 0 0 2b r e a kp o i n t s 0 . 0 0 0 0 0 . 0 1 0 0 0 . 0 0 0 0 other 0.0000 0.0000 0.0000 relative locations of break points in cases where 1 break point was found average (standard deviation) break point 0.0491 (5.1e-3) 0.0486 (4.5e-3) 0.0486 (4.5e-3) relative frequencies of AR orders when 2 break points were found Segment 1 order 0 0.0102 0.0102 0.0102 order 1 0.9492 0.9543 0.9543 order 2 0.0305 0.0254 0.0254 order 3 0.0101 0.0102 0.0102 order >30 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 2 order 0 0.0000 0.0000 0.0000 order 1 0.9695 0.9949 0.9949 order 2 0.0254 0.0051 0.0051 order 3 0.0051 0.0000 0.0000 order >30 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0

Figure 1 :

 1 Figure 1: Graphical interpretation of the bilevel optimisation method.

  (0, 1): δ =(m, κ),κ ∈ (0, 1) m ,m ∈{0,...,m max }, where m max is the maximum number of break points allowed. If m =0 ,t h ec h r o m o s o m er e f e r st o an AR model without break point.

and f ( 15 )

 15 < f (16) < f (17), then we set the pointer to the middle of the pointers already visited, that is, m next pointer = 1 2 (8 + 16) = 12 and go on with the next loop. The loops are finished when f (m pointer )i sn o tl a r g e rt h a nt h en e i g h b o u r i n gf (•)v a l u e s : f (m pointer -1) ≥ f (m pointer ) ≤ f (m pointer +1).

Figure 2 :

 2 Figure2: Performance comparison of the single-level and the bilevel method in the case of the PAR many model (200 independent runs -8 AR segments were generated, 7-8 were found; ca. 30K samples of objective function values in the optimisation history per one run; the boxes represent the 0.25, 0.5, 0.75 quantiles, the whiskers represent ca. the 0.07, 0.993 quantiles of the objective function values at chosen sample points). The average code length of the random generated resies using the PAR many model evaluated at the true parameter values was 2052 (std: 42) over the 200 independent runs. The NI algorithm was used here.

  Figure 3: PAR dyad

  Figure 5: P ARMA

Figure 7 :

 7 Figure 7: Short

Figure 8 :

 8 Figure 8: Performance comparison of the NI and the CMA-ES algorithm in the case of PAR many model (200 independent runs -8 AR segments were generated, 7-8 were found; ca. 20K samples of objective function values in the optimisation history per one run; the boxes represent the 0.25, 0.5, 0.75 quantiles, the whiskers represent ca. the 0.07, 0.993 quantiles of the objective function values at chosen sample points).

Figure 9

 9 Figure 9: PAR dyad

Figure 11

 11 Figure 11: P ARMA

Figure 12

 12 Figure 12: Tvar MA

Figure 13 :

 13 Figure 13: Short

Table 2 :

 2 Comparison of CPU seconds per objective function evaluation. The values are the aggregate CPU time divided by the total number of function evaluations for each experiment.

	Model	NI alg.		CMA-ES alg.
		single-level method	bilevel method
	PAR dyad	5.7e-4	21.7e-4	8.5e-4
	SlowAR	6.0e-4	22.1e-4	9.7e-4
	P ARMA	5.3e-4	21.3e-4	8.0e-4
	Tvar MA	7.1e-4	29.2e-4	12.6e-4
	Short	5.9e-4	17.6e-4	6.5e-4
	PAR many	12.6e-4	47.2e-4	17.7e-4

Table 3 :

 3 Comparisons of the number of function evaluations until the bilevel optimisation method with NI and CMA-ES algorithms reached the final range of code length values.

	Model	NI algorithm CMA-ES algorithm	Figure
	PAR dyad	140K	2K	Figure 9
	SlowAR	140K	4K	Figure 10
	P ARMA	140K	4K	Figure 11
	Tvar MA	140K	2K	Figure 12
	Short	70K	1K	Figure 13
	PAR many	520K	14K	Figure 8

Table 4 :

 4 Comparison of some statistics of the experiments (PAR many model).

Table B . 6 :

 B6 SlowAR

  Table B.7: Comparison of some statistics of the experiments (P ARMA).

  Table B.8: Comparison of some statistics of the experiments (Tva MA).

			0.0051	0.0051
	order 2	0.1888	0.2296	0.2296
	order 3	0.4388	0.4184	0.4184
	order 4	0.2857	0.2755	0.2755
	order 5	0.0612	0.0612	0.0612
	order 6	0.0153	0.0102	0.0102
	order >60 . 0 0 5 1	0 . 0 0 0 0	0 . 0 0 0 0

  Table B.9: Comparison of some statistics of the experiments (Short).

	order 1	0.1644	0.1644	0.1644
	order 2	0.3562	0.3151	0.3014
	order 3	0.1849	0.1849	0.1849
	order 4	0.2397	0.2534	0.2740
	order 5	0.0479	0.0685	0.0685
	order >50 . 0 0 6 8	0 . 0 1 3 6	0 . 0 0 6 8
			Segment 2
	order <20 . 0 0 0 0	0 . 0 0 0 0	0 . 0 0 0 0
	order 2	0.4726	0.4247	0.4110
	order 3	0.1575	0.1644	0.1849
	order 4	0.3082	0.3493	0.3425
	order 5	0.0548	0.0548	0.0548
	order >50 . 0 0 6 8	0 . 0 0 6 8	0 . 0 0 6 8
			Segment 3
	order <20 . 0 0 0 0	0 . 0 0 0 0	0 . 0 0 0 0
	order 2	0.1507	0.1712	0.1712
	order 3	0.0685	0.0822	0.0822
	order 4	0.6712	0.6233	0.6301
	order 5	0.0685	0.0890	0.0822
	order 6	0.0274	0.0205	0.0205
	order >60 . 0 1 3 7	0 . 0 1 3 7	0 . 0 1 3 7

The samples are ordered according to the corresponding objective function values.

Acknowledgements

This work was partially supported by the Île-de-France council/DIGITEO under the Software and Complex Systems program. We thank the anonymous refereesf o rt h e i rh e l p f u lr e c o m m e n d ations.

Experimental results

The goal of the experiments is to assess, respectively, the contribution of the bilevel methodology first, and then of the introduction of a new encoding for the break points and of CMA-ES as the optimisation strategy. Thus, the comparison was split into two parts.

First, we show the empirical comparison results concerning the single-level method (AutoPARM-GA) versus the bilevel method using the NI algorithm. Though in this case the problem parametrisation was different, the optimisation algorithm is the same:t h eN Ia l g o r i t h m .T h eq u e s t i o ni nt h i s case is which one of the different exploration strategies exhibit better performance. The numbers of function evaluations are not comparable due to the fundamentally different objective functions. Therefore we chose the target performance metric to be the CPUt i m ev e r s u st h eo p t i m i s a t i o nh i story. This metric is indeed comparable for the two optimisation methods if they are run in exactly the same computer environment and the same implementation ofth eN Ia lg o rith mimp leme n ta tio n is used for both cases.

Second, we show the empirical comparison results concerningt h eb i l e v e lo p t i m i s a t i o nu s i n gt h e NI algorithm versus the CMA-ES one. Because of the different implementations of the optimisation algorithms, a fair comparison should exclude the implementation-related effects of the evolutionary algorithms. Thus, in this case we compared the number of objective function evaluations versus the history of the best objective function values.

The technical details of the experiments are given in Section4 . 6 .

The benchmarks

All experiments are conducted on the following six test cases:

• PAR many -P i e c e w i s eA Rm o d e lw i t hm a n ys e g m e n t s • Short -S h o r ts e g m e n t s PAR many is motivated by [START_REF] Elteto | Discovering piecewise linear models of grid workload[END_REF], with a relatively large number of segments. The other ones have been defined by Davis in [START_REF] Davis | Structural break estimation for nonstationary time series models[END_REF]; for the sake of completeness,t h e ya r er e c a l l e di nA p p e n d i xA . PAR many model has 8 segments defined by

Appendix A. Experiment details Appendix A.1. PAR dyad

This example follows a piecewise AR model given as

Table B.5 shows the detailed results.

Appendix A.2. SlowAR

This example follows a time-dependent AR model given as

where a t =0.8(1 -0.5cos(πt/1024)) and t ∼ iid N (0, 1). Table B.6 shows the detailed results.

Appendix A.3. P ARMA

This example follows a P ARMA model given as

Table B.9 shows the detailed results.

Appendix B. Detailed results for the fitted parameters

NI algorithm

CMA-ES algorithm single-level method bilevel method relative frequencies where different number of break points were found 1b r e a kp o i n t 0 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 2b r e a kp o i n t s 0 . 9 6 0 0 0 . 9 4 5 0 0 . 9 9 0 0 3b r e a kp o i n t s 0 . 0 4 0 0 0 . 0 5 5 0 0 . 0 1 0 0 4b r e a kp o i n t s 0 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 other 0.0000 0.0000 0.0000 relative locations of break points in cases where 2 break points were found average (standard deviation) break point 1 0.4973 (1.0e-2) 0.4985 (1.0e-2) 0.4969 (1.1e-2) break point 2 0.7512 (1.0e-2) 0.7502 (9.4e-3) 0.7514 (1.2e-2) relative frequencies of AR orders when 2 break points were found Segment 1 order 0 0.0000 0.0000 0.0000 order 1

1.0000 1.0000 1.0000 order >10 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 2 order <20 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 2 0.9572 0.9893 0.9893 order 3 0.0321 0.0053 0.0053 order 4 0.0107 0.0053 0.0053 order >40 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Segment 3 order <20 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 order 2 0.9091 0.9840 0.9893 order 3 0.0856 0.0160 0.0107 order 4 0.0053 0.0000 0.0000 order >20 . 0 0 0 0 0 . 0 0 0 0 0 . 0 0 0 0 Tamás Éltető currently is a postdoc at Paris-Sud University. His researchi n t e r e s t sa r ei n Probability Theory and Mathematical Statistics.

Nikolaus Hansen is a researcher at The French National Institute for Researchi nC o m p u t e r Science and Control (INRIA). Before he joined INRIA, he has been working in applied artificial intelligence, in genomics, in evolutionary computation andi nc o m p u t a t i o n a ls c i e n c e . H i sm a i n research interests are learning and adaptation in evolutionary computation and the development of algorithms applicable in practice.

Pascal Bondon is a researcher at The French National Center for Scientific Research (CNRS). He works in the Laboratory of Signals and Systems (L2S) at Supé l e c . H i s r e s e a r c h i n t e r e s t s a r e i n time series analysis and prediction theory.

Cécile Germain-Renaud is a full Professor at Paris-Sud University. Her research interests are in models for High Performance Computing, and in Autonomic Grid and Clouds. She is involved with the European Grid Initiative, and leads the Grid Observatory project. Her web site is at http://www.lri.fr/~cecile/communications-english.html.


Emmanuel Caruyer

Rachid Deriche

Diffusion MRI Signal Reconstruction with Continuity Constraint and Optimal Regularization

Keywords: Diffusion MRI, Laplace Regularization, Q-space imaging

In diffusion MRI, the reconstruction of the full Ensemble Average Propagator (EAP) provides new insights in the diffusion process and the underlying microstructure. The reconstruction of the signal in the whole Q-space is still extremely challenging however. It requires very long acquisition protocols, and robust reconstruction to cope with the very low SNR at large b-values. Several reconstruction methods were proposed recently, among which the Spherical Polar Fourier (SPF) expansion, a promising basis for signal reconstruction. Yet the reconstruction in SPF is still subject to noise and discontinuity of the reconstruction. In this work, we present a method for the reconstruction of the diffusion attenuation in the whole Q-space, with a special focus on continuity and optimal regularization. We derive a modified Spherical Polar Fourier (mSPF) basis, orthonormal and compatible with SPF, for the reconstruction of a signal with continuity constraint. We also derive the expression of a Laplace regularization operator in the basis, together with a method based on generalized cross validation for the optimal choice of the parameter. Our method results in a noticeable dimension reduction as compared with SPF. Tested on synthetic and real data, the reconstruction with this method is more robust to noise and better preserves fiber directions and crossings.

Introduction

In diffusion MRI, the acquisition and reconstruction of the signal attenuation on the 3D Q-space allows reconstruction of the full probability of water molecules displacement, known as the ensemble average propagator (EAP). The radial and angular information contained in the EAP opens a wide range of applications, such as the definition of new biomarkers [START_REF] Cluskey | Mechanisms of neurodegeneration 547 in amyotrophic lateral sclerosis[END_REF][START_REF] Piven | An mri study of the corpus callosum in autism[END_REF], or the characterization of axon diameters in the brain white matter [START_REF] Assaf | Ax-514 caliber: a method for measuring axon diameter distribution from 515 diffusion mri[END_REF][START_REF] Özarslan | Nmr characterization of general compartment size distributions[END_REF]. The reconstruction techniques are based on the acquisition of diffusionsensitized MR signals, with the acquisition sequence described in [START_REF] Stejskal | Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient[END_REF], in which a pair of diffusion encoding magnetic field gradient are applied before and after the 180 • pulse. There exists a

Fourier relation between the diffusion attenuation E(q)

Email address: Emmanuel.Caruyer@inria.fr (Emmanuel Caruyer) URL: http://www-sop.inria.fr/members/Emmanuel. Caruyer (Emmanuel Caruyer) and the EAP 18 P(r) = R 3 E(q)e -2ıπq•r d 3 q,

(1)

where the wave vector q is directly related to the applied 19 magnetic field gradient pulse magnitude, direction, and 20 duration.

21

The diffusion tensor [START_REF] Basser | Estimation of the ef-525 fective self-diffusion tensor from the nmr spin echo[END_REF] is the first tensor, such as high angular resolution diffusion imag-30 ing [START_REF] Tuch | Q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]Aganj et al., 31 2010), or higher order tensors ( Özarslan and Mareci, 32 2003). Beyond these approaches, it is possible to recon-33 struct the model-free diffusion propagator, through Dif-

34

fusion Spectrum Imaging (DSI) [START_REF] Wedeen | Mapping complex tissue architecture with diffusion spectrum magnetic resonance imaging[END_REF], basis [START_REF] Assemlal | Efficient and robust 517 computation of pdf features from diffusion MR signal[END_REF]. DSI relies on the sampling of the diffusion signal on a regular Cartesian grid, and reconstructs the EAP through fast Fourier transform.

The main limitation of DSI is its huge demand in acquisition time, and gradient pulse strength to fulfill the Nyquist conditions [START_REF] Callaghan | Principles of nuclear magnetic resonance mi-528 croscopy[END_REF][START_REF] Tuch | Q-ball imaging[END_REF]. DPI [START_REF] Descoteaux | Multiple q-shell diffusion propagator imaging[END_REF] is a more natural method to describe the diffusion signal by a basis of functions solution to the 3D Laplace equation by parts. Though this method enables analytical reconstruction of the diffusion propagator, it cannot represent the diffusion signal in the whole Q-space. Indeed, DPI represents the signal using the 3D Laplace equation by part [START_REF] Descoteaux | Multiple q-shell diffusion propagator imaging[END_REF])

E(q • u) = l,m c l,m q l+1 + d l,m q l Y l,m (u), (2)
where Y l,m is the real, spherical harmonic function. The basis functions in DPI diverge both for q → 0 and q → ∞.

The SPF basis functions instead have a radial profile with a Gaussian-like decay, which is similar to the commonly observed diffusion signal. Besides, it is possible to recover the EAP (Cheng et al., 2010b) and the Orientation Distribution Function (ODF) (Cheng et al., 2010a) from the coefficients of the signal reconstructed in the SPF basis. The SPF basis is thus a unique, modelfree approach for the reconstruction of the full signal E, the estimation of EAP and its derived characteristics. It has been introduced in [START_REF] Assemlal | Efficient and robust 517 computation of pdf features from diffusion MR signal[END_REF] together with a regularization method to overcome ill-condition of the estimation problem.

However, the definition of the 3D functions of the SPF basis makes use of the parameterization q ∈ R 3 = q • u, where q ∈ R + and u ∈ S 2 . Near the origin, the corresponding u is not unique, and we show in Section 2.1 that continuity problems near the origin may arise if this parameterization is not used with care. Adding to that, the regularization method introduced in [START_REF] Assemlal | Efficient and robust 517 computation of pdf features from diffusion MR signal[END_REF] is based on a pair of empirical angular and radial low-pass filters. This regularization method fully relies on the choice of the basis of functions. Besides, its implementation requires to tune two separate regularization weights, which is impractical.

In this work, we propose original and efficient solutions to solve all these important problems. First, we show that continuous functions reconstructed in the classical SPF basis lie in an affine subspace which has a significantly reduced dimension. This means that the signal diffusion could be represented in this subspace with less coefficients, leading to an estimation process with less measurements than those required when repre-

B n,l,m (q • u) = R n (q)Y l,m (u). ( 3
)
Y l,m is the real, symmetric spherical harmonic intro-112 duced in [START_REF] Descoteaux | Apparent diffusion coefficients from high angular resolution diffusion imaging: Estimation and applications[END_REF], and the radial func-

113 tion R n is reported below for the record 114 R n (q) = κ n L 1/2 n q 2 ζ exp - q 2 2ζ (4) κ n = 2 ζ 3/2 n! Γ(n + 3/2) , (5)
where L 1/2 n is the generalized Laguerre polynomial, and The SPF basis is orthonormal for the dot product

115 Γ is the Gamma function Γ(z) = ∞ 0 t z-1 e -t dt
f, g = R 3
f (q)g(q)d 3 q.

(7)

The construction of this basis was motivated by the need for a complete orthonormal basis of antipodally symmetric and real functions. Besides, the radial profiles R n have a quasi-Gaussian decay, so that even a low radial truncation order leads to an accurate reconstruction and extrapolation beyond the sampling domain of the diffusion weighted attenuation E(q). reviewFrom the reconstruction of the signal in this basis, we can estimate the EAP following Cheng et al. (2010b) and the ODF following Cheng et al. (2010a).

However, a closer look at the functions B n,l,m near the origin reveals rapid oscillations and a discontinuity.

Moreover, by definition the value of the attenuation E is equal to 1 when q = 0, but there is nothing in the SPF basis to impose this. In this work, we show that the subset of functions verifying these properties of continuity and imposed value at the origin is an affine subspace of Ω N,L . We propose mSPF, an orthonormal basis for this subspace, and we give for convenience the relation between this modified SPF (mSPF) basis and the SPF basis B n,l,m introduced in [START_REF] Assemlal | Efficient and robust 517 computation of pdf features from diffusion MR signal[END_REF].

We also derive the Laplacian regularization functional expression in the mSPF basis, for a robust reconstruction of the diffusion signal. Indeed, the dimension of the basis grows rapidly with the angular and radial orders, and diffusion weighted images have a very low SNR. For the reconstruction of a smooth function, the Laplacian operator is a commonly proposed approach for regularization [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]. We derive the calculation of the Laplacian operator in the mSPF basis. The method is analytical, which ensures a fast implementation and reconstruction.

In this section, we use indifferently a notation with three indices for the bases elements, such as B n,l,m , or a notation with a simple index i, convenient for matrix notation. The link between both indexing systems is given by the functions n(i), l(i) and m(i).

Continuity in

Ω N,L Theorem 1. A function f = n,l,m a n,l,m B n,l,m of the SPF basis is continuous if and only if ∀l > 0, ∀|m| ≤ l, n a n,l,m R n (0) = 0. ( 8
)
The proof of this theorem is detailed in Appendix A.

The linear constraint in Eq. 8 imposes that the poly-

nomial part of f l,m = n a n,l,m R n has no constant term.
This linear constraint can be imposed while estimating 176

In addition to this continuity constraint, we empha-177 size that the diffusion attenuation signal is defined as 178 E(q) = S (q)/S (0), and therefore should verify

179 f (0) = 1. ( 9
)
The set of continuous functions in Ω N,L verifying Eq. 9

180
is the solution of an inhomogeneous linear equation, and

181
therefore is an affine subspace of Ω N,L . This affine space

f (q) = exp - ||q|| 2 2ζ + n,l,m x n,l,m C n,l,m (q), (10)
where {C n,l,m } is a basis of Ω 0 N,L , the subspace of con- 195

We first construct a basis of radial functions {F n }, ex-196 pressed as F n (q) = χ n q 2 /ζP n (q 2 /ζ) exp(-q 2 /2ζ). This

F n , F p R 3 = ∞ 0 F n (q)F p (q)q 2 dq = δ n,p . (11)
The substitution u = q 2 /ζ in Eq. 11 gives

201 ∞ 0 χ n χ m ζ 3/2 2 P n (u)P p (u)u 5/2 e -u du = δ n,p . (12)
The generalized Laguerre polynomial L 5/2 n suits this or-

F n (q) = χ n q 2 ζ L 5/2 n q 2 ζ e -q 2 /2ζ , (13)
and the normalization constant

χ n = 2 ζ 3/2 n! Γ(n + 7/2) . ( 14
)
The diffusion attenuation E(q)exp(-||q|| 2 /2ζ) is reconstructed through the functions

C n,l,m (q) = F n (||q||)Y l,m q ||q|| . ( 15
)
The family of functions

{C n,l,m , n = 0 . . . N -1, l = 0 . . . L, m = -l . . . l} is the modified SPF (mSPF) basis,
an orthonormal basis of Ω 0 N,L .

The coefficients x n,l,m are estimated by minimization of the squared error criterion ||y -Hx|| 2 , where y is the

vector of observations y k = E(q k ) -exp(-||q k || 2 /2ζ)
measured at wave vectors q k . The observation matrix

has entries H k,i = C n(i),l(i),m(i) (q k ).
This new space has a substantially reduced dimen-

sion: dim(Ω N,L ) = (N + 1) • L(L + 1)/2, whereas dim(Ω 0 N,L ) = N • L(L + 1)/2
. This dimension reduction comes from the two systems of linear constraints of Eq. 8 (L(L + 1)/2 -1 equations), and Eq. 9 (1 equation). As an example, when the angular truncation order L = 4 is used, the reconstruction in Ω 0 N,L requires 15 less coefficients, to represent the same signal. This simplifies the implementation, reduces the demand in storage capacity, and improves computational efficiency.

Link with the SPF basis

In this section we give the link between SPF and mSPF bases. This relationship is useful as SPF (As- between Laguerre polynomials detailed in (Abramowitz and Stegun, 1970, p. 783), we have:

F n (q) = n i=0 3χ n 2κ i R i (q) - (n + 1)χ n κ n+1 R n+1 (q). ( 16
)
If the function f (q) = E(q)exp(-||q|| 2 /2ζ) is expressed in this basis, f (q) = x n,l,m C n,l,m (q), then the coefficients a n,l,m of E in the SPF basis are obtained by a = Mx + a 0 , where In this section, we propose to introduce a regulariza-251 tion term in the fitting procedure. We choose as a regu-

242 M i j = δ l(i),l( j) δ m(i),m( j) •                    3χ n( j) 2κ n(i) n(i) ≤ n( j) - n(i)χ n( j) κ n(i) n(i) = n( j) + 1 0 n(i) > n( j) + 1 and a 0 = [ √ 4π/κ 0 0 0 . . .] T , as exp(-||q|| 2 /2ζ) = 243 √ 4π/κ 0 B 0,0,0 (q).
252 larization functional 253 U(x) = R 3 |∆E x (q)| 2 d 3 q, ( 17
)
where

E x (q) = exp(-||q k || 2 /2ζ)+ i x i C i (q) is
U(x) = (x -x 0 ) T Λ(x -x 0 ) + U 0 . (18)
Hence the penalized least squares has a unique mini-

266 mum 267 x = x 0 + (H T H + λΛ) -1 (y -Hx 0 ). (19)
In what follows, we give explicit directions how to com-268 pute the matrix Λ and the vector x 0 .

269

When E x (q) -exp(-||q k || 2 /2ζ) is expressed in the 270 mSPF basis with coefficients x i , 271 U(x) = R 3        i x i ∆C i (q) + ∆e -||q k || 2 /2ζ        2 d 3 q (20) = i j x i x j R 3 ∆C i (q) • ∆C j (q) d 3 q +2 i x i R 3 ∆C i (q) • ∆e -||q k || 2 /2ζ d 3 q + . . . ( 21
) v i = R 3 ∆C i (q) • ∆e -||q k || 2 /2ζ d 3 q (23)
The Laplace operator ∆ can be written in spherical 276 coordinates, with the Laplace-Beltrami operator ∆ b ,

277 ∆C n,l,m (q u) = χ n 1 q 2 ∂ ∂q q 2 F ′ n (q) Y l,m (u) + F n (q) q 2 ∆ b Y l,m (u) (24)
Since the spherical harmonics are eigenfunctions of the 278 Laplace-Beltrami operator with eigenvalue -l(l + 1), we

279 have 280 ∆C n,l,m (q u) = χ n F ′′ n (q)+2 F ′ n (q) q - l(l + 1)F n (q) q 2 Y l,m (u)
(25) As the spherical harmonics form an orthonormal basis 281 for the canonical dot product on S 2 , the entries of the

282 matrix Λ are 283 Λ i, j = δ l(i),l( j) δ m(i),m( j) ∞ 0 h i (q)h j (q) dq, (26)
where

284 h i = χ n(i) qF ′′ n(i) + 2F ′ n(i) - l(i) (l(i) + 1) q F n(i) . (27)
Similarly, the vector v has entries

285 v i = δ l(i),0 δ m(i),0 ∞ 0 h i (q) • q 3 ζ 2 - 3q ζ exp - q 2 2ζ
dq.

(28) The computation of the integrals in Eq. B.1 and 28 is 286 analytical and needs no numerical integration. It is de- algorithm [START_REF] Craven | Smoothing noisy data with spline func-549 tions[END_REF] to find the regu-292 larization weight λ which guarantees the best balance 293 between the smoothness of the reconstruction, and the 294 data fit. This algorithm, as well as the L-curve method 295 [START_REF] Hansen | The l-curve and its use in the numerical treatment of inverse problems[END_REF], have already been applied successfully 296 for other applications in Q-ball diffusion MRI (Koay 297 et al., 2009;[START_REF] Descoteaux | Analytical q-ball imaging with optimal regularization[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]. The GCV 298 method has the major advantage to be generalizable to 299 the situation where there is more than one λ parameter 300 to optimize. It is the case in [START_REF] Assemlal | Efficient and robust 517 computation of pdf features from diffusion MR signal[END_REF],

301
where there are two regularization matrices N and L, 302 which act respectively as radial and angular low-pass 303 filters, with corresponding weights λ N and λ L .

304

The GCV method is based on a one-fold cross valida- For the experiments, we compare the diffusion signal, the ensemble average propagator (EAP) reconstructed from the SPF coefficients by the method in (Cheng et al., 2010b), and the orientation distribution function (ODF)

reconstructed in constant solid angle, implementing the technique in (Cheng et al., 2010a).

Exact and empirical continuity constraints

We presented in Section 2.1 a linear constraint to impose the continuity of the reconstructed signal. An alternate solution proposed in (Cheng et al., 2010b) is to artificially add P virtual data points q k , k = K + 1 . . . P close to zero, verifying E(q k ) = 1. As P goes to infinity, it is possible to show that the solution of this system tends to the exact solution (see Golub and Van Loan, 1983, pp. 410-412). We study the convergence of this empirical continuity approach. As a measure of discontinuity of the reconstructed signal Ê about 0, we define d( Ê) the difference between extremal values of the set {lim q→0 + Ê(qu), u ∈ S 2 }. We also compare the relative difference between the solution c AC of the least squares problem with analytical constraint, and the solution c EC (P) of the system with empirical constraint with P virtual measurements.

Results and discussion

Continuity constraint

We compare the solution c AC and c EC (P), for a single Gaussian distribution. To focus on the continuity constraint, we do not impose any other kind of regularization. The signal is corrupted by Rician noise, with corresponding SNR = 25. An example of signal and its reconstruction is reported on Fig. 1.

We evaluate the difference of the signal reconstructed with exact continuity constraint and with empirical constraint. We plot on Fig. 2 the relative squared difference between the coefficients estimated with a strict continuity constraint, ĉAC , and the coefficients estimated with an empirical continuity constraint, ĉEC . The convergence is pretty fast, and P = 60 virtual measurements give good results. This confirms the intuition in (Cheng et al., 2010b); however the minimum number of virtual measurements P for an acceptable accuracy heavily depends on the angular order of the SPF basis, as reported

L =4 L =6 L =8 L =10
Figure 2: Relative difference between reconstruction with a strict continuity constraint, and reconstruction with a loose continuity constraint. Results on a synthetic Gaussian diffusion signal, from K = 150 measurements on 3 Q-shells, plus P virtual measurements at q = 0, for various angular orders L of the SPF basis. Depending on the radial order, the number of additional measurements needed for an accurate reconstruction may become huge, and really impractical.

on Fig.

Strict continuity constraint

Figure 1: Diffusion signal corresponding to a single fiber oriented along the x-axis, reconstructed from 120 samples in the Q-space. The signal is shown on the (q x , q y )-plane, and the grey levels correspond to signal range from 0.0 (white) to 1.0 (black). q values are understood in mm -1 . This illustrates the discontinuity at the origin inherent to the SPF basis, and how the reconstruction in mSPF solves this problem.

1 fiber 2 fibers, 90 • 2 fibers, 60

• (λ 0 L , λ 0 N ) (4.0 • 10 -7 , 8.1 • 10 -9 ) (3.2 • 10 -7 , 1.2 • 10 -8 ) (5.1 • 10 -8 , 5.5 • 10 -8 ) GCV 0 L,N 5.7 • 10 -1 3.4 • 10 -1 4.8 • 10 -1 λ 0 Λ 1.6 • 10 -1 1.7 • 10 -1 2.4 • 10 -1 GCV 0 Λ 5.3 • 10 -1 3.1 • 10 -1 4.2 • 10 -1
Table 1: Optimal λ parameters and corresponding GCV minimum, for various synthetic diffusion models. The sampling consists in 200 diffusion weighted measurements on 3 Q-shells, with a max b-value of 3000s • mm -2 . Radial and angular orders were set to 5 and 6, respectively. 1st row: separate Laplace-Beltrami and radial low-pass filter smoothing, 2nd row: Laplace regularization.

-

Laplace regularization

Figure 4: Reconstruction and extrapolation of a diffusion signal, for a Gaussian diffusion model, from 120 measurements on 3 Q-shells. We plot the reconstructed (solid lines) and ground truth (dashed lines) radial profiles of the signal on selected lines in the Q-space. The maximum q value of the sampling scheme was set to 60mm -1 , the hatched area represents the no-sample area. We compare the reconstruction without regularization, with separate Laplace-Beltrami and radial filter, and with Laplace regularization. Laplace regularization performs better in smoothing radial profiles, and we avoid oscillations outside the sampling area.

Ground truth

Laplace-Beltrami and radial low-pass filter Laplace regularization

Mean discontinuity

Figure 3: Discontinuity, measured about the origin, of a synthetic Gaussian diffusion signal, reconstructed from K measurements on 3 Q-shells, plus P virtual measurements at q = 0. The discontinuity remains very high, even for a large number of additional, virtual measurements (P = 150).

We also compare the reconstruction with both regu-416 larization constraints on our synthetic diffusion field in The function h i can be written as

492 h i (q) = χ n(i) q ζ exp - q 2 2ζ G n(i),l(i) q 2 ζ , (B.3)
where G n,l = k g n,l k X k is a polynomial. It is hard 493 to express the coefficients g n,l k in a compact form. In-494 stead of manually deriving these coefficients, we com-495 pute them using polynomial algebra facilities, provided 496 in the SciPy library [START_REF] Jones | SciPy: Open source scientific tools for Python[END_REF] in Python TM .

497

The coefficients g n,l k are algebraically computed on de- Hence the integrand h i (q)h j (q) can be written as

502 h i (q)h j (q) = χ n(i) χ n( j) ζ exp - q 2 ζ T i, j q 2 ζ (B.4)
where T i, j (X) is the polynomial XG n(i),l(i) (X)G n( j),l( j) (X).

503

The coefficients a i, j k of T i, j are simply obtained from the 504 coefficients of G n(i),l(i) and G n( j),l( j) . Therefore, the en-

22

  model historically proposed to describe the EAP. De-23 spite its wide acceptance into the research and clini-24 cal communities, this model restricts the diffusion EAP 25 within the family of Gaussian probability density func-26 tions, and is limited for the description of complex tis-27 sue structure. Since then, several models and meth-28 ods were described to extend the results of diffusion 29

35

  Diffusion Propagator Imaging (DPI) (Descoteaux et al., 36 2011), Diffusion Order Transform ( Özarslan et al., 37 2011) or reconstruction in Spherical Polar Fourier (SPF)

170

  the coefficients by constrained least squares estimation. 171 Alternatively, we will derive a new basis of functions 172 to span the subspace of continuous functions. This ap-173 proach greatly simplifies the Laplace regularization for-174 mulation and implementation, as we show in the next 175 section.

  182 is fully characterized by an underlying linear subspace, 183 and an origin. It is underlain by Ω 0 N,L , the kernel of the 184 associated homogeneous equation f (0) = 0. As for the 185 origin of the affine subspace, we can choose any solu-186 tion of Eq. 9. For the sake of simplicity, we choose a 187 simple Gaussian as the origin. 188 To sum up, any function f ∈ Ω N,L verifying the conti-189 nuity property, together with the property f (0) = 1 can 190 be expressed as 191

  192 tinuous functions f in Ω N,L verifying f (0) = 0. In the 193 remaining of this section, we give a construction for the 194 orthogonal basis {C n,l,m }.

  197 verifies F n (0) = 0; the polynomials P n and the normal-198 ization constant χ n are to determine, provided that the 199 following orthogonality property is fulfilled 200

  semlal et al., 2009) is a now a state-of-the-art method in diffusion MRI. We can therefore reconstruct the ensemble average propagator (EAP) following Cheng et al. (2010b), the orientation distribution function (ODF) following Cheng et al. (2010a), or the apparent fiber population dispersion following Assemlal et al. (2011). The SPF basis is built on Laguerre polynomials L 1/2 n while we use L 5/2 n in this work. Using the recurrence relations

244M

  is the change-of-basis matrix from mSPF to SPF, 245 two orthonormal bases. Therefore, this matrix is orthog-246 onal: the orthogonal projection of any function in Ω N,L , 247 represented by its coefficients a in the SPF basis, onto 248 the subspace Ω 0 N,L has coefficients x = M T a. 249 2.3. Laplace regularization in the mSPF basis 250

  the recon-254 structed signal. This continuous operator is rotational 255 invariant, and independent on the choice of a specific 256 basis. Besides, the Laplace operator was already ap-257 plied successfully for several applications ranging from 258 natural image denoising (You and Kaveh, 2000; Chan 259 and Shen, 2005) to diffusion MRI analysis (Descoteaux 260 et al.[START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Koay | A signal transformational framework for breaking the noise floor and its applications in mri[END_REF][START_REF] Descoteaux | Analytical q-ball imaging with optimal regularization[END_REF].261We minimize ||y -Hx|| 2 + λU(x), where the observa-262 tions are y k = E(q k )exp(-||q k || 2 /2ζ) and H is the ob-263 servation matrix. In this section, we write the Laplace 264 penalization as a quadratic form 265

ω

  305tion: among K samples, we use K -1 samples to fit the 306 model parameters, and predict the K-th left-apart sam-307 ple. The process is repeated K times, and the mean pre-308 diction error is the value we want to minimize. Fortu-309 nately, the mean prediction error, called the GCV func-310 tion, has a simple expression 311 GCV(λ; y) = ||yŷλ || 2 K -Tr(S λ ) , (29) which makes this method very efficient. The matrix 312 S λ = H(H T H + λΛ) -1 H T is the smoother matrix, and 313 ŷλ = S λ y. With the GCV method, it is possible to 314 adapt the regularization parameters to the data. How-315 ever, there is no analytical solution for the minimization 316 of the GCV function and for computational efficiency, 317 we compute the optimal λ parameters once. This choice 318 is validated in the next section, and results show it is 319 p exp(-2πτq T D p q), (30) where P ∈ 1, 2, 3 is the number of compartments, ω p is 324 the relative compartment size and D p the corresponding 325 diffusion tensor. The diffusion weighted signal is cor-326 rupted by Rician noise, with controlled variance param-327 eter σ. Using this diffusion model locally, we created 328 a synthetic diffusion field simulating a sin-shaped and a 329 straight fiber, crossing each other at 90 • . 330 The wave vectors q k for synthesis are arranged on 3 331 shells, with the strategy recently proposed in (Caruyer 332 et al., 2011a,b). In short, this method is a generaliza-333 tion of the electrostatic repulsion, introduced in (Jones 334 et al., 1999) for single Q-shell experiment design, to the 335 multiple Q-shell case. 336 The experiments on real data were carried out on the 337 publicly available phantom (Poupon et al., 2008; Fillard 338 et al., 2011) which served as the data for a tractography contest, held at the DMFC MICCAI workshop, London (2009). The diffusion signal was sampled on 3 Q-shells, with b-values ranging from 650 to 2000 s • mm -2 , and 64 directions per shell.

402

  The regularization also impacts on the extrapolation 403 capacity of the method. Hardware limitations often re-404 strict the sampling to a bounded region in the Q-space. 405 Increasing the radial order of the mSPF basis will allow 406 better signal reconstruction within the sampled area of 407 the Q-space. It might however introduce undesirable os-408 cillations outside this area, as reported on Fig. 4, where 409 the radial truncation order was set to N = 5. Adding 410 a regularization constraint greatly improves the extrap-411 olation of the diffusion signal. Laplace regularization 412 performs slightly better in this task, though a more com-413 plete study, involving real data and outside the scope of 414 this paper, should be carried out to further validate this.

Figure 5 :

 5 Figure 5: Reconstruction of a diffusion propagator field, from 120 measurements on 3 shells (max b-value was 3000s • mm -2). We compare the diffusion EAP profile (top row) P(r 0 u), for r 0 = 15µm, and the diffusion ODF ψ(u) (bottom row). Fiber crossing are better resolved with Laplace regularization, and isotropic regions are smoother.

  417

Fig. 5 .

 5 Fig. 5. Laplace regularization performs better in cross-418

422

  Similar results are obtained on the real data experi-423 ment, depicted on Fig. 6. We have overlaid the ground 424 truth fiber orientations, as provided by Fillard et al. 425 (2011). The reconstruction results with optimal Laplace 426 regularization show slighly sharper EAP and ODF pro-427 files. We acknowledge that the reconstruction of this 428 dataset was very challenging, due to the low anisotropy 429 of the signal.4305. Conclusionsconsider 0 < l ≤ L and -l ≤ m ≤ l. By continuity of R n , ∀u ∈ S 2 , |Y l,m (u)| ||Y l,m || ∞ ≤ 1, (A.6) hence Λ i, j = δ l(i),l( j) δ m(i),m( j) χ n(i) qF ′′ n(i) + 2F ′ n(i) -l(i) (l(i) + 1) q F n(i) . (B.2)

  l(l + 1) 7(3l(l + 1)/2) 15.75(3l(l + 1

  2. This makes this empirical solution imprac-

	386
	387	tical. Besides, discontinuity is not strictly imposed: as
	388	experimented and reported on Fig. 3, the value of d( Ê)
	389	remains unacceptably high while we impose the value
	390	on P = 150 virtual measurements.
	391	4.2. Laplace regularization
	392	Laplace regularization was implemented in the mSPF
	393	basis, and we compare it with separate Laplace-
	394	Beltrami and radial low-pass filter, proposed in (Assem-
	395	lal et al., 2009). The GCV function is significantly lower
	396	for the optimal Laplace regularization (Table 1). This
	397	result suggests that Laplace regularization is more suit-
	398	able than separate Laplace-Beltrami and radial low-pass
	399	filtering. Furthermore, the optimal λ Λ parameter does
	400	not vary much from one diffusion model to another. We
	401	can therefore select a unique λ Λ parameter for the regu-
		larization of a whole volume.

the minimization. Thus we have the quadratic form of 273 Eq. 18, where 274 Λ i j = R 3 ∆C i (q) • ∆C j (q) d 3 q, (22) and x 0 = Λ -1 v, with

We have proposed a novel orthonormal basis for the 432 reconstruction of the diffusion signal in the complete 3D 465 we can write ∀ǫ ′ > 0, ∃α > 0 such that

This proves the continuity of f l,m about 0, and by linear-485 ity the continuity of f .

Appendix B. Laplace regularization matrix 487

In this appendix, we derive the general expression of 488 the Laplace regularization matrix Λ in the mSPF basis.

489

The entries of the matrix Λ are