Dominique Gay

Vincent Lemaire

Should we Reload Time Series Classification Performance Evaluation ? (a position paper)

Since the introduction and the public availability of the ucr time series benchmark data sets, numerous Time Series Classification (TSC) methods has been designed, evaluated and compared to each others. We suggest a critical view of TSC performance evaluation protocols put in place in recent TSC literature. The main goal of this "position" paper is to stimulate discussion and reflexion about performance evaluation in TSC literature.

Introduction

The need for time series data mining benchmarks [START_REF] Keogh | On the need for time series data mining benchmarks: A survey and empirical demonstration[END_REF] has been fulfilled: with, firstly, the ucr Time Series Classification Archive [START_REF] Chen | The UCR time series classification archive[END_REF] and then the uea & ucr Time Series Classification Repository [START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF], the research community now have more than 85 data sets from various application domains to evaluate newly introduced TSC methods and to compare them to existing ones. The public availability and the wide di↵usion of the benchmark data had a strong and positive impact on the research community which has been prolific in terms of publications of TSC methods; as an example, the recent experimental evaluation in [START_REF] Bagnall | The great time series classification bake o↵: a review and experimental evaluation of recent algorithmic advances[END_REF] involves 18 recently proposed algorithms while in the same year, in 2017, two contenders, BoPF [START_REF] Li | Linear time complexity time series classification with bag-of-patternfeatures[END_REF] and weasel [START_REF] Schäfer | Fast and accurate time series classification with WEASEL[END_REF] have been presented in top data mining conferences.

In most of the research papers, the experimental evaluation section starts with "Each UCR dataset provides a train and test split set which we use unchanged to make our results comparable the prior publications" 3 . And this is where the reflexion we suggest begins. In the following, we will consider accuracy as the measure of predictive performance since it is the most widely used in TSC literature.

Discussion & reflexion

About train/test experiment. While in transactional data classification, resampling strategies (mainly bootstrap and stratified k-fold cross-validation) are the norm to estimate the expected performance of classifiers and to compare them [START_REF] Salzberg | On comparing classifiers: Pitfalls to avoid and a recommended approach[END_REF][START_REF] Delgado | Do we need hundreds of classifiers to solve real world classification problems[END_REF], there is a singularity concerning predictive performance evaluation in TSC literature: the vast majority of research work restrains predictive performance evaluation to a single train/test split experiment [START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF], also called hold-out method or test sample estimation [START_REF] Kohavi | A study of cross-validation and bootstrap for accuracy estimation and model selection[END_REF].

Unless disposing of a large and representative data set of the application domain, a single train/test experiment is generally ine↵ective in providing predictive performance estimation and valuable comparison between TS classifiers without random subsamplings (i.e., repeated hold out experiments). Indeed, different samplings may lead to results with strong variations [START_REF] Kohavi | A study of cross-validation and bootstrap for accuracy estimation and model selection[END_REF]. Thus, a classifier C A may show better predictive performance than another classifier C B just because of this particular split. And, in such cases, subsequent statistical tests based on single train/test accuracy results (such as now commonly used posthoc Nemenyi test [START_REF] Demsar | Statistical comparisons of classifiers over multiple data sets[END_REF]) do not help more in comparing classifiers performance over several data sets since other train/test splits could have led to di↵erent accuracy and mean rank results and thus potentially di↵erent conclusions.

Moreover, unless disposing of tens of thousands instances, it is common to keep more instances in the training set than in the test hold out, generally around class-stratified ratio of 2/3 for training against the rest i.e., 1/3 for testing. The TSC repository [START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF] provides various predefined train/test split ratios, going from 1.6% to 81% of the whole data set for training set (notice that 35 of the 85 TSC data sets show a train/test split ratio below 34%). In addition to small -if not very small-training sets, the train/test splits are not always classstratified; it results in very (too) few representative instances of some class labels (especially in multi-class problems). We are aware that, in some application domains, "labelled data is expensive to collect" [START_REF] Bagnall | The great time series classification bake o↵: a review and experimental evaluation of recent algorithmic advances[END_REF], however, considering the whole available and already-labelled data, splitting for such small non-stratified training sets is also questionable. Here are some singularities that arise from some data sets drawn from the TSC archive:

- Generally speaking, in addition to the weakness of single train/test experiments for performance evaluation, choosing to split such way, without classstratification in small training sets, results in "unnecessary di cult" TSC problems: indeed, the obtained training set is not always representative of the whole available data set and as explained above, it could lead to class distribution change between training and test sets (also known as prior probability shift [START_REF] Moreno-Torres | A unifying view on dataset shift in classification[END_REF]).

In such environments, intuitively, simple 1-Nearest Neighbor lazy learners [START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF] and ensemble methods that embed several weasel and ST [START_REF] Hills | Classification of time series by shapelet transformation[END_REF][START_REF] Bostrom | Binary shapelet transform for multiclass time series classification[END_REF] score the highest mean ranks when considering data sets with training set size greater than 500 while CoTE takes advantage when adding data with smaller training size. We also observe that the mean rank of ST increases as we consider more and more data sets with smaller training size. Notice that, it has to be balanced against the fact that as the number of data sets decreases the critical di↵erence (CD) value increases, making more di cult to state significant di↵erences of performance between contenders. Another illustrative example is about snb [START_REF] Boullé | Towards automatic feature construction for supervised classification[END_REF] (an improved Naïve Bayes classifier benefiting from multiple representations) which is competitive with weasel, CoTE, boss and st when not considering too small (< 200) training set size -confirming the importance of the size of the training set on the predictive performance of some classifiers. Aside from [START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF], we may regret the lack of experimental studies about the learning curves [START_REF] Perlich | Tree induction vs. logistic regression: A learning-curve analysis[END_REF][START_REF] Ting | Defying the gravity of learning curve: a characteristic of nearest neighbour anomaly detectors[END_REF] of TSC algorithms.

About resampling strategies.

As far as we know, very few attempts of resampling strategies for TSC performance evaluation have been led, e.g., :

-Grabocka et al. [START_REF] Grabocka | Invariant time-series classification[END_REF] provides some results by 5-folds cross validation on 35 ucr data sets. It allows to identify easy data sets in the TSC archive [START_REF] Chen | The UCR time series classification archive[END_REF]. Indeed, a default SVM classifier with polynomial kernel scores above 95% accuracy (often near perfect) on 18 data sets.

-Wang et al. [START_REF] Wang | Experimental comparison of representation methods and distance measures for time series data[END_REF], focusing on 1-NN with various distance measures, provides k-folds stratified cross-validation results. However, the k varies from 2 to 30 depending on the benchmark data set and the cross-validation method at use is unconventional: when splitting the data set T into k folds, the model is learnt on fold T k and tested on T \ T k -while conventional k-folds cross validation being the opposite: learning the model on T \T k and testing on T k . This unconventional cross-validation leads to the same problems explained above (with 1/30, i.e., 3% of the whole data set used for training).

In the same paper, the authors also noticed the importance of the e↵ect of the training data set size on 1-NN classifier accuracy : DTW-1-NN is better than ED-1-NN with small training data set, but providing that we have enough training instances (a few hundred/thousand depending on the simulated data set), the two classifiers show similar predictive performance.

-Bagnall et al. [START_REF] Bagnall | The great time series classification bake o↵: a review and experimental evaluation of recent algorithmic advances[END_REF] performs 100 resampling experiments on each of the 85 TSC data sets -followed by Nemenyi's statistical post-hoc test. However, the multiple resamplings fit the original size of train/test split provided by [START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF] -which leads to the same problem raised above about training size and class distribution change. This resampling strategy gives, all the same, a better idea of the performance of recent classifiers on data with various sizes of training sets, although the train/test split is still questionable. Indeed, if instances from original test set are authorized to be in training set due to resampling, why not use 10-CV or 10⇥10 CV, as in transactional data classification literature5 ?

The cross-validation (CV) method is not unknown to the TSC community; indeed, some algorithms, like e.g., weasel or DTW-CV used cross-validation on training set to set hyper-parameters (even if the training set is very small), then a single train/test split is performed to evaluate the performance of the "best hyper-parametered" model on a single hold out test set. Again, why not use CV method to evaluate and compare TSC algorithms ?

About (repeated) 10-folds CV, statistical tests and beyond.

While 10-folds CV with subsequent statistical tests [START_REF] Demsar | Statistical comparisons of classifiers over multiple data sets[END_REF] is now the gold standard for predictive performance evaluation and comparisons between classifiers on transactional benchmark data sets, recently, Vanwinckelen & Blockeel [START_REF] Vanwinckelen | On estimating model accuracy with repeated crossvalidation[END_REF][START_REF] Vanwinckelen | Look before you leap: Some insights into learner evaluation with cross-validation[END_REF] warns the Machine Learning community about pitfalls hidden in such comparisons. The take away messages are:

-"Our experiments show that when using cross-validation for choosing between two models, the best performing model is not always chosen". -"This discussion leads us to question the usefulness of statistical testing in the context of evaluating predictive models with cross-validation".

On the other hand, after almost a decade of the "10-CV + statistical tests" combination [START_REF] Demsar | Statistical comparisons of classifiers over multiple data sets[END_REF] to evaluate learner's predictive performance, J. Demsar et al. [START_REF] Benavoli | Time for a change: a tutorial for comparing multiple classifiers through bayesian analysis[END_REF][START_REF] Corani | Statistical comparison of classifiers through bayesian hierarchical modelling[END_REF] "discourage the use of frequentist null hypothesis significance tests (NHST) in machine learning and, in particular, for comparison of the performance of classifiers" and encourages the community to embrace Bayesian analysis using 10⇥10-CV for comparing classifiers. This is perhaps a change point for performance evaluation habits in the Machine Learning community. Notice that Bayesian analysis of performance is more conservative than NHST; that is, it is "more di cult" for a classifier C A to be better than a classifier C B , considering Bayesian analysis.

About the evaluation measure.

As recalled in the introduction, the vast majority of recently proposed TSC algorithms are evaluated and compared with regards to the accuracy measure, i.e., the number of correctly classified time series. Accuracy measure is suitable for roughly balanced 2-class data sets. However, for unbalanced and/or multiclass data sets, accuracy measure is inappropriate for evaluation since high accuracy results due to a bias towards the majority class could hide severely bad predictive performance on the minor class or on other classes in multiclass settings. Often, ROC or Precision/Recall curve analysis or lift curve and cumulative gain charts are prefered for unbalanced settings. The TSC archive [START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF] contains some 2-class unbalanced problems, e.g., Earthquakes, ToeSegmentation2 and Wafer with respectively 20.2%, 25.3% and 10.6% unbalanced ratio (i.e., the proportion of the minor class). The archive also contains many multiclass problems with severe unbalanced ratios, e.g., ECG5000 and Worms with respectively 0.5% and 9.7% unbalanced ratio. Learning in the presence of class imbalance or in multiclass settings is still an ongoing machine learning research topic [START_REF]Proceedings of the IJCAI 2017 Workshop on Learning in the Presence of Class Imbalance and Concept Drift (LPCICD'17)[END_REF]. Again, the presence of such data sets in the benchmark repository, when averaging ranks of classifiers based on accuracy results, could lead to flawed conclusions on the performance evaluation.

Conclusion

The still ongoing public release of benchmark TSC data sets to the data mining community through the ucr & uea repository has certainly been the best catalyst for the development of new TSC algorithms and methods. In this discussion paper, we briefly review the habitual protocols at use in TSC algorithms performance evaluation, discuss the pros and cons of experimental protocols and try to warn the community about the pitfalls and hidden problems of actual performance evaluation protocols in TSC literature. We agree that the core of this discussion paper needs more in-depth development and experimental arguments but we believe that interesting discussions on this important topic deserve to be launched and continued during the 3 rd ECML/PKDD Workshop on Advanced Analytics and Learning on Temporal Data.

  DiatomSizeReduction: a 4-class problem, with a train/test split ratio about 5% (16 instances for training, 306 for testing) and a class distribution (1,6,5,4) for training against (33,92,94,87) for testing. The class distribution is not respected from train to test set. There are relatively 73% more instances of class c 1 in the test set than in the training set. -SonyAIRobotSurface1: a 2-class problem, with a train/test split ratio about 3.22% (20 instances for training, 601 for testing) and a class distribution (6,14) for training against (343,258) for testing. The class distribution is not respected from train to test set. There are relatively 90% more instances of class c 1 in the test set than in the training set; and 38% less instances of class c 2 in the test set than in the training set. -. . . several additional similar cases of data sets can be found in commonly used benchmark [5] : either the number of training examples is very small w.r.t the whole available data, either there is class distribution change between training set and test set (sometimes both singularities arise).

Table 1 .

 1 1-NN classifiers (with various distance measures or on various data representations)[START_REF] Bagnall | Transformation based ensembles for time series classification[END_REF][START_REF] Bagnall | Time-series classification with COTE: the collective of transformation-based ensembles[END_REF][START_REF] Lines | Time series classification with ensembles of elastic distance measures[END_REF] still obtain "good" accuracy results since it is still possible to find a nearest neighbor of a test instance in very few training instances of a minor class. However, eager classifiers based on empirically observed per-class frequencies (such as Naïve Bayes or Decision Trees) often fail in characterizing the minor classes with very few representative instances.If the TSC archive[START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF] o↵ers a wide variety of data sets, the original train/test splits also involves hidden di cult well-known problems in the Machine Learning community: e.g., learning from few examples or in class distribution change environment. Unfortunately, averaging ranks over all data sets to produce critical di↵erence diagram presents a risk of hiding the reasons why a particular classifier shows better performance than another. As an example4 , in Table1, we report performance comparisons of 11 recent classifiers like in recent literature (single train/test split following by significance testing). We provide several statistical tests by integrating step by step data sets which involves smaller size of training set. Average ranks of 11 recent classifiers depending on the data sets taken into account from[START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF] (i.e., depending on the training set size). Accuracy results are taken from Schäfer & Leser paper on WEASEL[START_REF] Schäfer | Fast and accurate time series classification with WEASEL[END_REF]. Post-hoc Nemenyi's statistical test considering training set with size 1000 (only 10 data sets), then > 500 (only 22 data sets), . . . , until considering All 85 data sets from[START_REF] Bagnall | The UEA & UCR time series classification repository[END_REF]. Underlined rank is the best per line and bold results on the same line indicates that there is no statistically significant di↵erence of performance between bold results according to Nemenyi's test, considering the current benchmark datasets.

	Min. size #DB	CD	WEASEL DTWCV DTW BOSS LS TSBF ST	EE CoTE SNB BoPF
	1000	10 4.77403	3.250	8.300 8.600 5.700 6.050 7.700 3.050 6.250 4.100 3.400 9.600
	> 500	22 3.21865	2.932	7.659 9.023 5.705 6.909 6.841 4.068 6.273 3.636 4.000 8.955
	> 300	42 2.32949	3.560	7.929 8.893 5.857 6.845 6.298 4.500 6.060 3.369 4.821 7.869
	> 200	48 2.17903	3.594	7.917 8.885 5.781 6.885 6.271 4.625 5.865 3.219 5.010 7.948
	> 100	57 1.99962	3.737	7.868 8.860 5.474 6.623 6.237 4.544 5.930 3.298 5.535 7.895
	All	85 1.63748	3.847	7.806 8.647 5.382 6.135 6.388 4.847 5.871 3.412 6.429 7.235

quoted from weasel paper[START_REF] Schäfer | Fast and accurate time series classification with WEASEL[END_REF].

We did not integrate recently introduced hive-cote[START_REF] Lines | HIVE-COTE: the hierarchical vote collective of transformation-based ensembles for time series classification[END_REF] accuracy results since they are available under 100 resamples protocol.

However, even with k-folds CV, a particular attention must be given to the setting of k for the 17 small data sets, with less than 200 instances, from the TSC archive.

Acknowledgments

We would like to thank the anonymous reviewers for joining the discussion with valuable comments and arguments.


Franck Chauvel
email: franck.chauvel@sintef.no

Sébastien Mosser
email: sebastien.mosser@sintef.no

Arnor Solberg
email: arnor.solberg@sintef.no

Reconsidering

Reconsidering Quality of Service Analysis in Dynamic and Open Systems

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

Introduction

With the emergence of middleware platforms supporting run-time adaptations (e.g., WS, OSGi, JEE), software systems are becoming more and more open and dynamic. These new abilities, which improve software reuse, also directly challenge our ability to validate these systems. This is further complicated when extra-functional requirements like Quality of Service (QoS) must be met. Overlooking QoS undoubtedly leads to overpriced rollbacks in the development, or even jeopardise the delivery of final software artifacts. To mitigate QoS issues, modern approaches [START_REF] Balsamo | Model-based performance prediction in software development: a survey[END_REF] validate end-to-end QoS requirements against predictive QoS models representing the system and its environment. Precise quantitative descriptions (e.g., probability distributions, memory capacities, CPU speeds) are critical to ensure high accuracy of QoS predictions. However, the more dynamic a system is, the more difficult it becomes to precisely characterise its environment, and, in turn, the less significant existing techniques are. The same is true for end-to-end requirements, whose numerical estimation becomes more and more irrelevant as the system bounds become uncertain. This inherent incompatibility between dynamicity and accuracy challenges the worthiness of QoS analysis: the effort spent at design-time validating QoS requirements must be put in perspective using the significance of the validation results at run-time.

We are convinced that our ability to tame dynamicity and openness is bounded to the way we address the uncertainty it implies. We believe that postponing validation until run-time to minimise uncertainty while using proper mathematical foundations shall help reconcile dynamicity and QoS engineering. In this setting, the contribution of this vision paper is to examine to which extent recent models at run-time techniques [START_REF] Morin | Models@ Run.time to Support Dynamic Adaptation[END_REF] can be combined with advanced modelling/reasoning techniques addressing uncertainty in the first place [START_REF] Bandemer | Mathematics of Uncertainty[END_REF]. To this end, we first present our approach as a development process focusing on QoS requirements. We then reviewed the key properties that an ideal QoS model must exhibit to enable the validation of QoS requirements under openness and dynamicity. We then discuss the construction, the validation and the injection of such QoS requirements into running systems. P. Collet, P. Merle (eds.); Conférence en IngénieriE du Logiciel (CIEL), Juin 2012, pp. 1-6 1

Motivations

Let us consider an application server, as a practical example of above-mentioned QoS issues. In a general sense, a server processes requests for resources (e.g., files, HTML pages, service invocations) and users are often expecting its response time to be bounded, (we consider here the server-side response time, i.e., the time between the reception of a request and the emission of the related response). A typical QoS requirement describing the server is "the server responsetime must be less than 3 s for 90 % of requests in a single day".

Various models and techniques can be used to validate that this requirement will hold at run-time, such as Markovian models for instance. In a general setting, servers' response time depends on two main factors, namely, the arrival rate (i.e., the speed at which requests arrive, in req.s -1 ), and the service rate (i.e., the speed at which the server processes them, in req.s -1 ). The task of QoS expert is therefore to anticipate, at design time, the environments the system may be facing. In Markovian models, this is achieved by means of probability distributions quantifying the likelihood of possible arrival rates, and service rates. For the sake of the example, let us assume the two following distributions: f , which describes the arrival rate, is a normal distribution centred around 50 req.s -1 with a deviation of 25 req.s -1 , whereas g, which describes the service rate, is also a normal distribution centred around 100 req.s -1 with a deviation of 25 req.s -1 . The Markovian model assesses that our requirement holds. However, in real life, when the actual arrival rates will exceed the service rate, our requirement will actually not hold anymore, as the server will either start to store or to reject pending request, depending of the overload management policy.

So what did we get wrong? Actually, the original QoS requirement and the one that we validated, are not the same. The one we validated is: "assuming the distributions f and g for the arrival and service rates, respectively, the response time must be less than 3 s for 90 % of the requests in a single day" and is much more restrictive.We could obviously relax our expectations on the expected execution environment by using uniform distributions, but we would not be able to validate our requirement: there will always be a non null probability that the arrival rate exceeds the service rate, in which case the response time will tends towards infinity.

The previous example shows the intrinsic difficulty of QoS requirement modelling and validation, even on simple examples. Existing QoS requirements reflect our expectations from a user perspective, whereas our validation techniques reflects the system's ability from the "pragmatic" engineer standing point. Open and dynamic systems acerbate this conflict and a paradigm shift is needed to validate QoS in such a context.

Overview of the Approach

The paradigm shift we envision is the continuous anticipation of the inherent conflict between the user expectations and the system abilities, so as to trigger preventive adaptations. The proposed approach combines advanced mathematical modelling of uncertainty with recent advances in models@runtime. Fig. 1 illustrates the development process making concrete continuous QoS validation in open and dynamic systems. This process reflects our will to turn the traditional once-and-for-all QoS validation into a continuous validation process at runtime. It is structured as follows:

• QoS Requirements Elicitation is jointly carried out by the user and the QoS expert.

The objective here is to obtain detailed QoS requirements that can be later turned to operational artifacts, stating the expectations from a user point of view. describing the dynamics of the system's QoS, describing the system's abilities. Recent advances in QoS models shows that automated extraction is needed to cope with the difficulty of quantifying the various parameters needed to obtain accurate QoS predictions [START_REF] Balsamo | Model-based performance prediction in software development: a survey[END_REF][START_REF] Krogmann | Using Genetic Search for Reverse Engineering of Parametric Behavior Models for Performance Prediction[END_REF].

• Initial QoS validation is also carried out by the QoS expert, but is an initial comparison between the system abilities and the user expectations. When this first validation fails, the QoS expert has to determine whether they are due to ill-formed requirements, or to error in the running system. Once issues have been mitigated, the process restarts with QoS model extractions or with alternative requirements.

• Probes and Requirements Injection is the combination of the running system (extended with additional monitoring capabilities) and the user requirements: the mismatch between requirements and systems abilities can be dynamically checked while the system is running, taking into account unforeseen environments discovered at runtime.

• Continuous QoS Validation, supported by probes injection, aims at anticipating when QoS requirements might actually be violated. This is further complicated by the potential softness of QoS requirements, whose violation may deal with softness regarding amplitude, duration or frequency.

• Corrective Maintenance is needed to moderate QoS issues anticipated at run-time. This might be done automatically using self-adaptation techniques [START_REF] Chauvel | Using QoS-Contracts to Drive Architecture-Centric Self-adaptation[END_REF] or semi automatically [START_REF] Morin | Taming Dynamically Adaptive Systems using Models and Aspects[END_REF].

Towards Assume/Guarantee-based QoS Requirements

The first obstacle to the process sketched on Fig. 1 is the lack of a model for QoS requirements. By contrast with traditional QoS requirements, we believe that key point is not to quantify end-to-end QoS, but to constraint its possible values by capturing the "envelope" in which the QoS is expected to evolve, depending on the context. Fig. 2 illustrates what we envision as a good QoS requirement, regarding our server example. In Fig. 2, we express that when the service rate exceeds the arrival rate, the response time shall be bounded by the time spent processing all the requests that are already pending, and no request shall be rejected. By contrast, if the arrival rate exceeds the service rate, new requests shall be rejected at a rate which is the difference between arrival rate and service rate. This example highlights several key properties of what we think are good QoS requirements:

• Completeness. Such requirements forms a complete partition regarding the different situations that might occur, namely nominal behaviour, or overloads in our example. Such a partition discards superfluous details but retains the critical distinctions. Building such partitions is critical to tame open and dynamic systems.

• Uncertainty. Such QoS requirements are inherently uncertain. First, they are imprecise as they capture relative orders of magnitude, instead of accurate prediction. For instance, we only assume that the difference between the arrival rate and the service rate is a good approximation of what must be the rejection rate. In addition, QoS requirements are inherently permissive or soft, as they make explicit what is acceptable in which contexts. Imprecision and Softness are the counter part of completeness, to tame open and dynamic systems.

• Clear Intention & Operational Extension. They have a clear intention: "bound the response time" in our example. This critical from the user perspective. Besides, assuming that the values at play can be observed at run-time, such requirements must be turned into operational checkers. This is critical from the QoS Expert perspective: Considering the effort to elicit and track one single QoS requirement through the development process, we directly question the worthiness of non-operational requirements.

• Qualitative. The envelope that is presented here is qualitative: it does not refer to crisp values (except for 0, which can be considered as a qualitative value). Such qualitative description, uses merely basic arithmetic operators, as non-linearity emerges from the identification of distinct cases (e.g., nominal versus overloads).

How one can come up with such requirements? We believe that such requirements can be elicited using the two main abilities of human reasoning: abstraction and decomposition. Decomposition is our ability to break down complex problems into smaller ones that can be separately solved, whereas abstraction is our ability to get rid of irrelevant details and to retain only critical distinctions. In the previous example, abstraction selects the properties of interest: service rate, arrival rate, rejection rate, and maximum number of waiting requests. By contrast, 4 Actes de CIEL 2012 decomposition identifies our two cases: one characterising the nominal behaviour of the system, the other one characterising overloads.

How to manage uncertainty? Aside of the well known Probability theory, there have been several attempts during the past 50 years to provide sound mathematical foundations for reasoning under uncertainty [START_REF] Bandemer | Mathematics of Uncertainty[END_REF]: Fuzzy Sets, Evidence Theory, Rough Sets, Grey Theory and numerous hybrid theories. Each of these theories addresses specific aspects of uncertainty, e.g., eventuality, cognitive imprecision, approximation, incompleteness. We are firmly convinced that, while abstraction and decomposition will help designers elicit QoS requirements, such formal foundations will permit to effectively exploit their inherent imprecision and softness. Modern Control Theory relies for instance on Fuzzy Set Theory to tackle complex non-linear and multivariate problems.

QoS Model Extraction and Validation

QoS Models, at least performance models [START_REF] Balsamo | Model-based performance prediction in software development: a survey[END_REF], have always suffered a very low acceptance due to the difficulty to quantify the numerous parameters that they require. However, recent advances in parameter estimation now permit to extract such parameters directly from running artifacts, typically traces or source code [START_REF] Krogmann | Using Genetic Search for Reverse Engineering of Parametric Behavior Models for Performance Prediction[END_REF]. Alternatively, it is also possible to extract QoS models using approximations techniques combined with performance testing: The model then mimic the behaviour observed during testing. QoS Models extraction is a great step towards a general and effective QoS engineering. In our process, it is critical to continuously validate the QoS of the system as often as needed, by changes in the context. Coming back to our server example, and assuming an existing QoS model, it becomes possible to validate that the response time is bounded, in nominal conditions as well as during overloads.

Models extraction is however still not sufficient as it results in accurate predictive models whose complexity is far too high to be embedded at runtime. There is hence a need for models approximations and efficient parameter identification to reflect context changes observed at runtime, and enable in turn continuous validation.

Probe Injection and Continuous Validation

The last obstacle to the process introduced in Fig. 1 is the ability to deploy probes that observe the various parameters involved in QoS requirements, and the ability to decide whether a given requirements will be violated or not. Measuring values at play in QoS requirements implies, in the general case, the development of specific probes. The design and the integration of such probes undoubtedly impacts the system itself, so it exhibits relevant information. While probes represent an additional development, their integration in the rest of the system can be lightened by advanced software engineering techniques, e.g., aspect-oriented programming or dynamic deployment.

Deciding whether a given QoS requirement will be violated or not depends on two aspects: (i) our ability to calculate, at run-time, the expected envelope and trajectory of the QoS properties of interest, and (ii) our ability to handle the imprecision and softness inherent to QoS requirements. Regarding run-time evaluation, we envision two main approaches, the dynamic interpretation of QoS requirements at run-time, or if their semantic is not computationally efficient enough, the injection of some precomputed numeric approximations. Regarding imprecision, it is critical to properly handle all acceptable relative deviations, both in term of amplitude and frequency. Efficient QoS fault detection is critical to ensure the good usability and stability of the system, especially is fault detection triggers self-adaptive mechanisms. In this perspective, we aim at reusing existing classification techniques (e.g., clustering, pattern recognition), already effective in other engineering fields.

Conclusions & Road Map

In this paper, we envision a paradigm shift to support the modelling of QoS requirements in adaptive systems that are open and dynamic systems. This approach is complementary to usual QoS methods, as it aims to address systems that are typically not designed to support usual QoS models. We are especially interested in looking for convergence with the development various services architectures, including Internet of Things and Internet of Services. Regarding IoT, we are especially exploring within the ThingML project [START_REF] Fleurey | MDE to Manage Communications with and between Resource-Constrained Systems[END_REF], which defines a model-driven framework enhancing the development of embedded systems.

On the short term, we are targeting at first the formal definition of a QoS requirements semantics that would handle uncertainty while preserving computational efficiency. Another task is the definition of QoS fault detection mechanisms that would report gradual degradations, supporting deviations in frequency or amplitude and enable preventive adaptations. On the mid-term perspective, our focus will be on the development of QoS requirements framework, that would combine traditional verification with the imprecision and softness of QoS requirements. Finally, in the long term perspective, we will focus on the evaluation of the feasibility, of the effectiveness and of the performance in securing QoS in dynamic and open systems.

Figure 1 :

 1 Figure 1: Overview of the process, designed to tackle QoS validation for open/dynamic systems

Figure 2 :

 2 Figure 2: Alternative QoS Requirements bounding the response-time of a server application

Actes de CIEL 2012

Acknowledgements. This work is funded by the SINTEF strategic projects SiSaS and MOD-ERATES.


Mauro Zucchelli

Samuel Deslauriers-Gauthier

Rachid Deriche

A Closed-Form Solution of Rotation Invariant Spherical Harmonic Features in Diffusion MRI

Keywords: Diffusion MRI, Rotation Invariant, Gaunt Coefficients, Spherical Harmonics

Rotation invariant features are an indispensable tool for characterizing diffusion Magnetic Resonance Imaging (MRI) and in particular for brain tissue microstructure estimation. In this work, we propose a new mathematical framework for efficiently calculating a complete set of such invariants from any spherical function. Specifically, our method is based on the spherical harmonics series expansion of a given function of any order and can be applied directly to the resulting coefficients by performing a simple integral operation analytically. This enable us to derive a general closed-form equation for the invariants. We test our invariants on the diffusion MRI fiber orientation distribution function obtained from the diffusion signal both in-vivo and in synthetic data. Results show how it is possible to use these invariants for characterizing the white matter using a small but complete set of features.

Introduction

Brain tissue microstructure influences the diffusion Magnetic Resonance Imaging (MRI) signal in-vivo in multiple ways. Tissue microstructure in gray and white matter can be subdivided into two main categories: tissue composition and tissue orientation. Tissue composition corresponds to the heterogeneity of the neurites in the voxel, namely the presence and relative fraction of a certain neurite population in the voxel. Axonal diameter and intra-axonal signal fraction are two of the most studied tissue composition microstructural features in diffusion MRI. Tissue orientation, on the other hand, corresponds to the geometrical arrangement of the neurites in a given voxel. White matter axonal bundles (also known as "fiber bundles"), because of their high anisotropy, represent the key factor which influences the diffusion MRI signal via tissue orientation. It is common practice in diffusion MRI to model the tissue orientation information in each voxel using the so-called fiber Orientation Distribution Function (fODF) [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion mri: non-negativity constrained super-resolved spherical deconvolution[END_REF]. The fODF represents the probability distribution function of the fibers, namely the probability of having a fiber oriented in a given direction. Such a distribution is commonly obtained from the diffusion signal by deconvolving it with a single fiber response which, in theory, should correspond to the tissue composition in the same location [START_REF] Dell'acqua | A modified damped richardson-lucy algorithm to reduce isotropic background effects in spherical deconvolution[END_REF] [START_REF] Tournier | Robust determination of the fibre orientation distribution in diffusion mri: non-negativity constrained super-resolved spherical deconvolution[END_REF] [START_REF] Zucchelli | A generalized smt-based framework for diffusion mri microstructural model estimation[END_REF]. In diffusion MRI, fODF is modeled using its real symmetric Spherical Harmonics (SH) series expansion, and therefore, each fODF is completely characterized by its SH coefficients. Being a well defined numerical quantity, the idea of using the SH coefficients directly as a biomarker for tissue orientation is appealing but it suffers from the fact that even a small rotation of the same fODF will produce a different set of coefficients, which makes their use unsuitable for classification and clustering tasks. An intuitive solution to this problem is to use the SH coefficients to calculate Rotation Invariant Features (RIF) which, as the name implies, have values which are independent of the direction of the fODF. The first and most well-known fODF RIF which has been developed for diffusion MRI is the Generalized Fractional Anisotropy (GFA) [START_REF] Tuch | Q-ball imaging[END_REF]. GFA represents the fODF evolution of the diffusion tensor-derived Fractional Anisotropy (FA) with the advantage of being able to better represent crossing fibers. However, reducing the rich directional information present in the fODF into a single numerical value severely limits its descriptive power. Bloy and Verma [START_REF] Bloy | Demons registration of high angular resolution diffusion images[END_REF] proposed to use the power spectrum of the SH coefficients to register diffusion MR images. The power spectrum provides a richer description with respect to GFA but is still incomplete with respect to the theoretical number of independent invariants obtainable. Ghosh et al. [START_REF] Ghosh | Generalized invariants of a 4th order tensor: Building blocks for new biomarkers in dmri[END_REF] [13] proposed a method for recovering all 12 invariants obtainable from a 4 th order tensor while Caruyer and Verma [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF] derived a framework for obtaining all the algebraically independent invariants for any given order of SH. Such a framework initially derives all the polynomials invariant to a rotation of 1 degree with respect to the x and z axes, which are proven to be invariant to any additional rotation. After collecting all the possible invariants, a pruning algorithm is applied to the polynomials to reduce them to the minimal set of algebraically independent invariants. The main limitations of this technique is that the invariants cannot be expressed in closed-form and that it produces a huge number of polynomials to be tested for independence (more than 30,000 for a SH order of 6).

The goal of this work is to propose a new framework for generating a set of invariants in closed-form (Section 2) with an elegant mathematical framework. Our solution is able to produce fewer polynomials to be tested for independence compared to previous published related work. In Section 3 we will present the synthetic and in-vivo fODF data which we will use to test our RIF. Results show the potential for these invariants to disambiguate complex fiber geometries using only a small set of invariants (Section 4).

Theory

A new set of invariants for SH

Given a spherical function f (u) defined on the three dimensional unit vector u we can represent it as a linear combination of SH as

f (u) = ∞ l=0 l m=-l c lm Y m l (u) (1)
where Y m l are the real SH functions [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][7], and c lm ∈ R are the SH coefficients. The coefficients c lm are calculated as the inner product f and Y m l

c lm = S 2 f (u)Y m l (u)du. (2)
A rotation R applied to the SH functions can be written as

RY m l (u) = Y m l (R -1 u) = l m =-l D l m ,m (R)Y m l (u) (3)
with D l m ,m (R) defined as

D l m ,m (R) = e -im α d l m ,m (β)e -imγ (4)
where α, β, and γ are the Euler angles and d m ,m is the Wigner d function.

Therefore, we can define h(u) = Rf (u) the function f rotated by R and whose SH expansion is given by

h(u) = Rf (u) = ∞ l=0 l m=-l c lm Y m l (R -1 u) = ∞ l=0 l m =-l g lm Y m l (u). (5)
The rotated SH coefficients for a given l are given by

g lm = l m=-l c lm D l m ,m (R). (6)
A RIF of a Spherical Function f can be expressed as a function I for which

I[f (u)] = I[Rf (u)], ∀R (7)
If f is expressed in terms of its SH expansion, we have that

I[c] = I[g] (8)
where c and g are the SH coefficients of f (u) and h(u), respectively. It is straightforward to see that the integral of f (u) represents a RIF, in fact:

I[f ] = S 2 f (u)du = S 2 f (R -1 u)du (9)
this equivalence follows the natural law that the area of a function does not change while rotating the function. In particular in this case

I[f ] = S 2 f (u)du = c 00 √ 4π
, and therefore, normalizing by the constant factor √ 4π, the polynomial P (c) = c 00 represents a first RIF for f . Similarly, we can prove that if we raise the function f to the power d, the integral of such a function is still a RIF. In fact

I d [f ] = S 2 [f (u)] d du = S 2 f (R -1 u) d du ( 10
)
Although it is an interesting RIF, I d [f ] represents only a partial representation of f . In order to overcome this limitation we can define a new set of RIF from

I d [f ].
In particular, we observe that we can write I d [f ] as

I d [f ] = ∞ l1=0 • • • ∞ l d =0 S 2 d i=1 li mi=-li c limi Y mi li (u) du = ∞ l1=0 • • • ∞ l d =0 S 2 l1 m1=-l1 • • • l d m d =-l d c l1m1 • • • c l d m d Y m1 l1 (u) • • • Y m d l d (u)du = ∞ l1=0 • • • ∞ l d =0 l1 m1=-l1 • • • l d m d =-l d c l1m1 • • • c l d m d S 2 Y m1 l1 (u) • • • Y m d l d (u)du = ∞ l1=0 • • • ∞ l d =0 I d l1...l d [f ] = l I d l [f ] (11)
with l = [l 1 . . . l d ]. I d l [f ] can be evaluated as

I d l1...l d [f ] = l1 m1=-l1 • • • l d m d =-l d c l1m1 • • • c l d m d G(l 1 , m 1 | • • • |l d , m d ) ( 12
)
where

G(l 1 , m 1 | • • • |l d , m d )
represents the generalized Gaunt coefficient calculated as the integral of the product of d SH of the same argument [START_REF] Homeier | Some properties of the coupling coefficients of real spherical harmonics and their relation to gaunt coefficients[END_REF]. We can also view I d l as a polynomial of degree d of the SH coefficients c lm . Appendix A provides a proof that Equation [START_REF] Novikov | Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion mri[END_REF] generates rotation invariant features of the spherical function f (u) given the set of SH indices l. The function I d l shares some of the properties of the Gaunt coefficients, in particular it is zero for most combination of l and it is invariant to the permutations of the indices. In this work, we will refer to the set of all the non-zero invariants of maximal SH order l max and degree d as I d lmax .

Example 1: degree 2 invariant-polynomials

In the case of polynomials of degree d = 2 and considering a maximal SH order l max = 2, I d l corresponds to 7) 12 (11) 18 (12*) lmax = 6 1 (1) 7 (4) 13 ( 13) 28 (25*) 49 (25*) Table 1. Number of RIF for symmetric (even) SH basis given the maximum value of l and the degree. In bold we report the number of algebraically independent invariants, with an asterisk if the number corresponds to the maximum theoretical number of invariants.

I 2 l1,l2 [f ] = l1 m1=-l1 l2 m2=-l2 c l1m1 c l2m2 G(l 1 , m 1 |l 2 , m 2 ) ( 13
) d = 1 d = 2 d = 3 d = 4 d = 5 lmax = 2 1 (1) 2 (2) 3 (3*) 4 (3*) 5 (3*) lmax = 4 1 (1) 5 (3) 7 (
with l = [(0, 0), (0, 1), • • • (2, 2)] considering only the sorted couple in ascending order of (l 1 , l 2 ) because of the symmetry property of the Gaunt coefficients. G(l 1 , m 1 |l 2 , m 2 ) represents the integral of two SH which is equal to δ l2,m2 l1,m1 . This means that the invariant will be non-zero only when l 1 = l 2 and m 1 = m 2 . We can thus rewrite the degree two invariants as

I 2 l [f ] = l m=-l [c l,m ] 2 (14)
which corresponds to the power spectrum of the SH basis. The power spectrum is one of the most widely used invariants in pattern recognition application and has been recently introduced in diffusion MRI for tissue microstructure estimation [START_REF] Novikov | Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion mri[END_REF].

Example 2: degree 3 invariant-polynomials

The invariants of degree 3 can be calculated as

I 3 l1,l2,l3 [f ] = l1 m1=-l1 l2 m2=-l2 l3 m3=-l3 c l1m1 c l2m2 c l3m3 G(l 1 , m 1 |l 2 , m 2 |l 3 , m 3 ) (15)
It is possible to show that I 3 l1,l2,l3 [f ] are actually equivalent to the bispectrum invariants proposed in [START_REF] Kondor | A novel set of rotationally and translationally invariant features for images based on the non-commutative bispectrum[END_REF], [START_REF] Kakarala | A theory of phase-sensitive rotation invariance with spherical harmonic and moment-based representations[END_REF], and [START_REF] Negrinho | Shape representation via elementary symmetric polynomials: A complete invariant inspired by the bispectrum[END_REF] (results not shown due to space limitation). The main difference with respect to the bispectrum is that in this case we have the Gaunt coefficients obtained from the integral of three SH instead of the Clebsch-Gordan coefficients.

It is possible to show that both I 2 l1,l2 and I 3 l1,l2 are not only invariants but they also form a set of algebraically independent invariants. The next subsection will provide the formal definition of algebraically independence and how to apply it in order to verify the independence of any set of I d lmax . 6) 20 ( 20) 57 (33*) lmax = 6 1 (1) 7 ( 7) 30 (30) 94 (46*) Table 2. Number of RIF given the maximum value of l and the degree. In bold we report the number of algebraically independent invariants, with an asterisk if the number corresponds to the maximum theoretical number of invariants.

Algebraically independence

∀Q ∈ R[Y 1 , . . . , Y m ], Q(P 1 , . . . , P m ) = 0 ⇐⇒ Q = 0
Example Given P 1 = 2x 2 -4xy + 2y 2 and P 2 = x -y, the set [P 1 , P 2 ] is not algebraically independent. In fact, considering Q = 2P 2 2 -P 1 it leads to:

2(x -y) 2 -(2x 2 -4xy + 2y 2 ) = 0.
An efficient way for verifying whether the invariants I d lmax are algebraically independent is the Jacobian theorem proposed by [START_REF] Ehrenborg | Apolarity and canonical forms for homogeneous polynomials[END_REF]:

Theorem 1 The polynomials P 1 , . . . , P m ∈ R[X 1 , . . . , X n ] with m ≤ n are algebraically independent iff the Jacobian matrix (∂P i /∂X j ) 1≤i≤m,1≤j≤n has full rank.

In our case, the Jacobian matrix corresponds to the partial derivatives of our invariants I d l with respect to the SH coefficients c lm . The testing for full-rankness is performed taking advantage of the Schwartz-Zippel polynomial identity testing lemma [START_REF] Schwartz | Fast probabilistic algorithms for verification of polynomial identities[END_REF] [START_REF] Zippel | Probabilistic algorithms for sparse polynomials[END_REF] which states that given a random vector of coefficients c, the Jacobian rank corresponds to the rank of (∂P i /∂X j ) X=c with a probability that is proportional to the probability of extracting c randomly. In this work, we consider 100 random instances of c lm in order to assess the independence of I d l with sufficient accuracy. This approach has already been used for the same purpose by [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF], where the authors also theorized that the maximal number of independent coefficients is equal to the maximal number of SH coefficients used, n c , minus three, the degrees of freedom of the rotation. Table 1 show the number of invariants calculated for the symmetric SH basis [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF]. In bold we highlight the number of algebraically independent invariants and we also add an asterisk if this number corresponds to n c -3. For the complete SH basis we observed that using d = 4 we are able to obtain all the invariants up to l max = 6, while for the symmetric SH basis our technique require d = 5 in order to obtain the 12 th invariant for l max = 4. Caruyer and Verma [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF] were able to obtain all 12 invariants for l max = 4 at d = 4 using a different technique, and we are also able to obtain the complete set of invariants at d = 4 if we include also the odd terms of the basis (see Table 2). It is possible that with our technique the 12 th invariant can be expressed as a combination of two mixed even-odd invariants and appears as a function of purely even coefficients only at d = 5.

Note that with our technique all the invariants at a certain degree include also all the invariants of lower degree (i.e. I 2 2,2 represents the same invariant as I 4 0,0,2,2 ). Therefore, in order to calculate all the algebraically independent invariants it is sufficient to consider only the degree d that provides the maximal number of invariants, without having to compute all the invariants at lower degree. For example, taking Table 1 first row, we can observe that we have two invariants at d = 2 of which one is a new invariant and the other is the same invariant that appears at d = 1. In practice, in order to obtain the number of new invariants appearing only at degree d we have to subtract the numbers of invariants of degree d -1.

Materials and methods

Synthetic fODF dataset: We created a synthetic fODF dataset consisting of crossings of two delta functions on the sphere. We can obtain the SH coefficients of two delta functions with equal volume fraction centered around v 1 and v 2 , respectively, as

c crossing l,m = 1 2 Y m l (v 1 ) + 1 2 Y m l (v 2 ) (16)
In this dataset, we will use the cosine of the crossing angle, v 1 • v 2 , as the reference value for the invariants. Since the maximum absolute value of our RIF is obtained when the considered fODF is a delta function (OD = 0 or v 1 • v 2 = ±1) we can use this value as a normalization factor for our invariants.

In this work we will call the normalized fODF invariants Îd l [ρ] to distinguish them from the unnormalized version.

Human Connectome Project (HCP) dataset: In order to test our RIF in-vivo we considered one subject of the HCP diffusion MRI dataset [START_REF] Sotiropoulos | Advances in diffusion MRI acquisition and processing in the human connectome project[END_REF]. The HCP acquisition scheme presents 18 b-value 0 s/mm 2 volumes, and 90 volumes at bvalue 1000, 2000, and 3000 s/mm 2 , respectively. All the three shell were acquired using an unique gradient direction. HCP pulse separation time was set at ∆ = 43.1ms and pulse width at δ = 10.6ms. Each MRI volume consists of 145 × 174 × 145 voxels with a resolution of 1.25 × 1.25 × 1.25 mm 3 .

Results

Figure 1 shows the proposed algebraically independent RIF given a maximum SH order l max = 4 for the synthetic fODF dataset, before (left) and after (right) the normalization by the maximum theoretical value. This dataset was generated by taking 1000 fODF directions v 1 and v 2 randomly, and the relative invariants was ordered according to the value of v 1 • v 2 . Because of the randomness of the orientation of the crossings, the smoothness of the curves itself can be view as an empirical proof of the rotation invariance property of our RIF. In fact, In order to test the invariants in-vivo we consider one subject of the HCP dataset. The main advantage of this dataset is the high spatial resolution which enables us to observe the fine changes of our RIF in the white matter. From the diffusion signal, a voxel-dependent response function was obtained using the multi-compartment model presented in [START_REF] Kaden | Multicompartment microscopic diffusion imaging[END_REF]. Each response function has been used for calculating the fODF with constrained spherical deconvolution (l max = 4) as in [START_REF] Zucchelli | A generalized smt-based framework for diffusion mri microstructural model estimation[END_REF]. Figure 2 shows the 12 algebraically independent invariants after normalization for a coronal slice of the HCP subject. As for the synthetic data, the first invariant l = (0) is constant because the first SH coefficient c 00 is always 1/ √ 4π in order to have the integral of the fODF equal to 1 [START_REF] Novikov | Rotationally-invariant mapping of scalar and orientational metrics of neuronal microstructure with diffusion mri[END_REF]. The average values of the invariants generally decrease with the increase of each l i and with the degree of the polynomials. For example, l = (2, 2) has an average of 0.12 while l = (2, 2, 2, 2, 4) has an average of only 0.0009 in white matter. Conversely, the contrast between white matter and gray matter generally increases with the degree. A similar trend can be observed in the contrast between single bundles and crossings in white matter. Some of our RIF present negative values in certain areas of the brain, in particular in fiber crossing regions. In our simulations we observed that negative values for certain invariants are associated with a characteristic crossing angle (e.g 90 degrees crossing for l = (2, 2, 2)). However, our synthetic data is formed only of 2-fiber crossings, whereas crossings of three or even more fibers are common in the human brain. Therefore, a direct mapping of the value of the RIF to the fiber geometry is not trivial and will require more advanced mathematical tools.

Discussion and Conclusion

In this work, we proposed a new framework for generating algebraically independent rotation invariant features for the diffusion MRI signal. One of the main results of our technique is that, conversely with respect to other works, we were able to provide a general closed-form for our invariants. The core mathematical tool of our formulation are the Gaunt coefficients originating from the integral of multiple SH. Gaunt coefficients also appear in the invariants for diffusion MRI proposed in [START_REF] Schwab | Rotation invariant features for hardi[END_REF], although using a completely different approach and without testing the resulting invariants for completeness or independence.

Rotation invariant features are not only a subject of interest in diffusion MRI but also in several other fields including computer vision and pattern recognition [START_REF] Kakarala | A theory of phase-sensitive rotation invariance with spherical harmonic and moment-based representations[END_REF] [START_REF] Kondor | A novel set of rotationally and translationally invariant features for images based on the non-commutative bispectrum[END_REF]. We believe that our framework will be advantageous for non-diffusion applications given the simplicity of the mathematical formulation of our RIF. Our mathematical formulation makes them very fast to compute even for large SH order compared to other techniques.

Recent works [12][14] show the importance of rotation invariant features such as the power spectrum for estimating tissue microstructural properties. We will focus our future works in trying to establish whether including the complete set of invariants will result in a more precise estimation of brain microstructure.

A Relation between signal SH coefficients and fODF SH coefficients

Definition 1 .

 1 A set of m polynomials P 1 , . . . , P m ∈ R[X 1 , . . . , X n ] are called algebraically independent over the field R if d = 1 d = 2 d = 3

d c l1m1 • • • c l d m d S 2 Y

 2 I d l [f ] is rotation invariant iff I d l [f ] = I d l [Rf ] = I d l [h]. Given the rotated SH expansion h(u) we can calculate I d l [h] as I d l1...l d [h] = d =-l d c l1m1 • • • c l d m d Y m1 l1 (R -1 u) • • • Y m d l d (R -1 u)du m1 l1 (R -1 u) • • • Y m d l d (R -1 u)du d =-l d c l1m1 • • • c l d m d G(l 1 , m 1 | • • • |l d , m d ) = I d l1...l d [f ](17)which proves the invariance propertyI d l [f ] = I d l [Rf ].

Acknowledgment This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (ERC Advanced Grant agreement No 694665 : CoBCoM -Computational Brain Connectivity Mapping). Data were provided by the Human Connectome Project, WU-Minn Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience Research; and by the McDonnell Center for Systems Neuroscience at Washington University.

l= [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF] l= [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= [START_REF] Caruyer | On facilitating the use of hardi in population studies by creating rotation-invariant markers[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= [START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF][START_REF] Descoteaux | Regularized, fast, and robust analytical q-ball imaging[END_REF] l= ( 0.54/0.026 0.19/0.012 0.02/0.009 0.28/0.018 l = (2, 2, 2, 4) 0.69/0.019 0.38/0.013 -0.02/0.003 -0.03/0.002 l = (2, 2, 4, 4) 0.60/0.022 0.28/0.012 0.05/0.004 0.06/0.003 l = (2, 4, 4, 4) 0.51/0.025 0.17/0.010 -0.00/0.002 0.06/0.003 l = (4, 4, 4, 4) 0.45/0.028 0.15/0.011 0.05/0.005 0.21/0.014 l = (2, 2, 2, 2, 4) 0.66/0.022 0.34/0.014 0.03/0.003 0.02/0.001 Table 3. Mean and standard deviation of the 12 algebraically independent invariants for lmax = 4 with respect to different crossing angles.

two adjacent points in the graph had been generated by crossings oriented, up to a rotation, in two completely different directions. If any of the I d l was not invariant its curve would have appear as a broken line in the graph. We also test the robustness of the invariants with respect to noise. We simulate diffusion signal arising from crossing fibers at 0, 30, 60, and 90 degrees. We add Rician noise with a signal to noise ratio equal to 30, generating 100 noise realization per crossing angle. The fODF have been calculated using constrained spherical deconvolution algorithm with SH order equal to 4. Table 3 reports the mean and standard deviation of the twelve algebraically independent invariants. In general, our invariants result to be stable with values of standard deviation fifteen times lower than the mean in the worst case. Invariants including higher harmonics terms seem to present higher standard deviation (e.g l = (2, 2) is more stable than l = (4, 4)). This is probably due to the fact that higher harmonics tend to model signal noise more with respect to the lower ones that are too smooth to capture rapid signal changes. 0.0000 0.0106 0.0212 0.0318 0.0424 -0.0025 0.0047 0.0118 0.0190 0.0261 0.0000 0.0072 0.0145 0.0217 0.0290 -0.00004 0.00819 0.01642 0.02466 0.03289 Fig. 2. All the 12 algebraically independent invariants calculated on the fODF of a coronal slice of one of the HCP subjects.