Alexandre Ravet

Simon Lacroix

Gautier Hattenberger

Augmenting Bayes filters with the Relevance Vector

come

Augmenting Bayes filters with the Relevance Vector Machine for time-varying context-dependent observation distribution

Alexandre Ravet 1,2 , Simon Lacroix 1,3 and Gautier Hattenberger 4 Abstract-Bayesian filtering often relies on a reduced system state relating to robot internal variables only. The exogenous variables and their effects on the measurement process are then encompassed within a global observation noise model. Even if Bayes filters proved to be robust to such approximations, special care has to be taken to handle some of these exogenous effects, usually by introducing complex observation distributions or rejection rules. No matter how complex these models are, they often fail in dealing with contextual incidence which can hardly be explicitly encoded. This article shows how contextual information can be introduced within the Bayesian filtering framework by coupling a filter with classification and regression probabilistic models. The classification model provides an efficient context-dependent measurement selection mechanism and is specifically trained with respect to the filter estimation performance. This first component is enhanced by the introduction of context-dependent observation noise provided by the regression model. The performance of this is approach is evaluated and compared with other methods in the context of altitude estimation for a UAV.

I. INTRODUCTION

The issue of understanding and exploiting measurements provided by different sensors is of major importance in state estimation. It is commonly accepted that the more accurate the observation model, the better. However, the physical measurement process is often too complex to be fully understood and exactly modeled. In practice, one has to deal with a whole range of performance alterations going from the unavoidable average observation noise to completely unreliable measure values. Instead of trying to explain the whole measurement set with an accurate observation model, a common way to cope with unmodelled observation alterations is to select the measurements which are known to be positively informative for the estimation process -in other terms to reject outliers. This yields robustness of the estimation process with respect to unmodeled effects on the observations, but at the risk of under-exploiting information. It is desirable to define an estimation scheme that goes beyond outlier rejection, introducing context information for a time varying modeling of measurement contribution.

a) Related work: Most approaches for solving the estimation problem rely on the framework of Bayesian filtering which proved to be robust and highly reliable. Any implementation of a Bayes filter is based on two key components, the state transition model, and the observation model. The optimal observation model should describe how sensor measurements are generated in the current system state. In the Bayesian filtering framework, this model is described as a conditional distribution p(z t |x t ) where z t is the set of measurements at time t and x t is the state of the system. As the underlying physics explaining the process of measurement formation is often too complex to be perfectly described by the observation model, contextual information (i.e. information about exogenous effects required to explain the formation of measurements) is usually not included in the state x t . Such information should also account for effects stemming from the unpredictability of dynamic environments, resulting for example in unexpected occlusions. All the unmodeled effects are usually represented through the observation noise, by making the assumption that they are similar to random effects. This is a strong assumption that often has to be compensated, for example through the introduction of outlier rejection methods [START_REF] Lynen | A robust and modular multi-sensor fusion approach applied to mav navigation[END_REF], or through the derivation of complex observation pseudodensities modeling sensor specific behaviors, such as the well known beam model described in [START_REF] Thrun | Probabilistic Robotics (Intelligent Robotics and Autonomous Agents)[END_REF]. These approaches suffer from over-simplification as they rely on a subjective understanding of the sensor characteristics, and also remain self-contained systems that implicitly accommodate all the context influence in a single static model. A direct consequence of this implicitness is that these systems are quite likely to diverge [START_REF] Ravet | Learning to combine multi-sensor information for context dependent state estimation[END_REF].

As deriving an accurate model by hand is too complex, learning algorithms have been applied in order to build, at first, parametric models [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation[END_REF], and more recently nonparametric models [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF]. By learning an observation model without introducing any prior knowledge over the sensor behavior, the latter approach proved to be very efficient, and has been extended to a fully state-dependent observation model through the introduction of heteroscedastic noise (meaning the observation noise is state-dependent), whereas usual models make the assumption of homoscedastic (constant) noise. The resulting model is of the form z t = f (x t ) + ǫ(x t ) where f and ǫ are represented by Gaussian processes (GP). Still, all these attempts to learn a precise model rely on a system oriented reduced state, and are consequently unable to model contextual influence over the measurement process. Moreover the global performance of the filter can strongly depend on the algorithm that utilizes the model [START_REF] Thrun | Probabilistic Robotics (Intelligent Robotics and Autonomous Agents)[END_REF].

b) Approach: This paper aims at showing that the joint set of sensor measurement values z t , potentially extended with any relevant contextual information i t (i.e. other sensor measurements, or other internal data), provides a succinct yet rich representation of the perception context that can be used as the input of a context-dependent observation model. Intuitively, we can see this joint set of measurement values as the minimum information relating the influence of the current context over sensor performance. As shown in [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF], [START_REF] Pfaff | Efficiently learning high-dimensional observation models for monte-carlo localization using gaussian mixtures[END_REF], one intuitive way to learn an observation model is to apply a direct regression method. However, we rely on a different approach for the two following reasons:

• In our case the contextual information is not part of the system state. It is an additional set of variables c t = {z t , i }. Therefore learning a direct mapping from the context information to the measurement is irrelevant, as the output of the system is contained in the input. • It is legitimately more appealing to optimize the observation model with respect to the ultimate system output performance (the state estimate accuracy) rather than learning a direct observation model which can be seen as an intermediate goal (explaining at best the mapping between state and context to measurements), but which does not provide any guarantee in term of global performance.

The core component of our model is used to encode a high-level knowledge over measurement validity which, conversely to other approaches, is learned with respect to the only goal of helping the system to provide the best estimate. Learning such knowledge implies to take into account the whole prediction-update process in the optimization and requires to examine each filter algorithm separately. A simple way to achieve this with classical learning methods, and without any consideration regarding the type of filter, is to view the observation function as a selector over measurement values. The function is then context-dependent in the sense that the actual subset of measurement values used for the estimation process change over time. Practically, this is done by running at first a group of different filters, all of these filters using a distinct yet combinatorially exhaustive subset of measurement values from a training set D = X, C , where X is a vector of ground truth states and C a vector containing measurements c t = {z t , i t } made at state x t . We then apply the Relevance Vector Machine (RVM) classification technique in order to learn a discrete mapping from the context space (c t ) to the most accurate filter at time t. This provides a context-dependent observation function based on a classifier optimized with respect to the global performance of the system (measured w.r.t X). This work is thus closely related to the bank of filter approach [START_REF] Bar-Shalom | Estimation with Applications to Tracking and Navigation[END_REF] [3], if we consider a bank of filters differing only in the subset of measurement values they are using. However, we enhance this first classification model by learning a context-dependent (heteroscedastic) observation noise model for each sensor with RVM regression. The classification and regression models are then combined to obtain what can be seen as a fully context-dependent observation model. c) Outline: The next section provide background details of RVM. We then present the specific implementation for learning context-dependent observation function (Section III), and observation noise (Section IV). Finally, we provide experimental results and comparison to existing approaches.

II. BACKGROUND ON SPARSE BAYESIAN LEARNING

In the past decade, machine learning has seen the emergence of sparse Bayesian methods, a specialisation of which being known as the Relevance Vector Machine (RVM) [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF]. In this work, RVM is chosen for the regression and classification task for two main reasons. First, as a nonparametric model, RVM introduces no assumption over the functional form underlying the mapping from the context to the appropriate measurement subset, or from the context to the current observation noise. It then provides a generic approach for different robot configurations. Secondly, RVM training naturally provides sparsity thanks to the automatic relevance determination mechanism. It is a considerable feature in the context of real-time Bayesian filtering, since inference over the model is fast. This is in contrast with GP, widely used for enhancing Bayes filters and for the task of heteroscedastic regression [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF], [START_REF] Turner | State-space inference and learning with Gaussian processes[END_REF], [START_REF] Lzaro-Gredilla | Variational heteroscedastic gaussian process regression[END_REF], [START_REF] Kersting | Most likely heteroscedastic gaussian process regression[END_REF], which require to turn to more complex sparse GP techniques when the system is intended to be used in real time.

A. RVM principle

Based on the probabilistic Bayesian framework, RVM provides a powerful solution to the problem of supervised learning, i.e. learning a mapping model between a set of input vectors {x n } N n=1 and corresponding target values {t n } N n=1 considered as noisy outputs of an underlying noise free function y(x n ). RVM can be seen as the equivalent probabilistic treatment of the Support Vector Machine (SVM), whilst avoiding most of its limitations. The shared principal characteristic of these methods is that they both base their predictions upon a function y(x) defined as a weighted sum of basis functions given by kernels K, with one kernel defined for each input of the training set:

y(x, w) = N i=1 w i K(x, x i ) + w 0 (1)
where w 0 is a bias parameter and w = {w 1 , . . . , w N } are the weights. The most appealing aspect of the RVM is that the obtained prediction model is sparse (most of w i parameters are close to zero), while still providing excellent generalisation performance for new input values as long as they remain in the implicit boundaries defined by the exploited training set. A widely used kernel function is the Gaussian kernel defined as

K(x m , x n ) = exp(- 1 r 2 x m -x n 2 ) ( 2
)
where r is called the width parameter and can be easily tuned during the learning process. By experience this function provides results for our application. However the use of a different kernel function may improve drastically the performance of the learned model, for example if we prefer to use a different width parameter along the different input dimensions.

Training and inference over the model are done within the Bayesian framework. By lack of space we will not derive equations for training and predicting, but all details can be found in [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF]. The RVM classification method originally introduced by Tipping is however designed for binary classification, and requires adaptations for the multiclass case. This is why we prefer to rely on the multiclass Relevance Vector Machine (mRVM) introduced in [START_REF] Psorakis | Multiclass relevance vector machines: sparsity and accuracy[END_REF] which naturally handles the multiclass setting and achieve better sparsification.

III. LEARNING CONTEXT-DEPENDENT OBSERVATION

FUNCTION

In many cases deriving the noise free component f of the observation model z t = f (x t )+ǫ is not the most problematic task, as it can be obtained directly through deterministic physics rules. In next sections we assume that these functions are known, and then alleviate the rest of the system which can focus on capturing the remaining unmodeled aspects of the observation model. We now detail the procedure for learning the context-dependent selection scheme and give an illustration of its application on simulated data.

A. Goal-oriented learning

The complete observation vector z t is composed of N measurements {z n t }. Based on pre-defined noise free component of the observation model for each of these measurements {z n t }, we are able to derive the combinatorially exhaustive set of possible observation functions. As these functions correspond to different contexts, we note them f c where c ∈ {1, ..., C} with C the total amount of distinct functions. We then train a RVM model such that, at each filter iteration, the model selects the most appropriate subset of measurements {z n t } n∈c in the complete measurement vector z t according to the current context, i.e. we use the corresponding observation function f c (x t ):

{z n t } n∈c ⊆ z t = f c (x t )
1) Defining the training data: The definition of the training data is a crucial step in this approach. The idea is to form an alternative training set based on the original data D = X, C in order to explicitly introduce the requirements in terms of global system performance. Thus we run C Bayes filter, each of which based on a different function f c , and save the corresponding state estimates and uncertainty. Then, each sample of D is associated with an activation vector A te.g. A t = {0, 0, 1, 0, ..., 0} T , which means that the best estimate output at time t was provided by the filter based on the observation function f c=3 . It should be noted that the definition of the 'best' output requires the evaluation of a performance metric which represents the expected requirements to the system behavior. In this paper, we evaluate the output distribution of each filter at the corresponding ground truth value x t (given by N (x t |x c t , σ c ) where xc t and σ c are the output mean and variance of the filter based the observation function f c ). Choosing the highest value among all filters is then equivalent to selecting the filter with the lowest estimation error and uncertainty. However, other requirements may be easily introduced through the evaluation of different metrics or the addition of penalty terms.

A new training set D ′ = C, A is formed, and is readily used for training a mRVM classifier with C class labels. This classifier is then used at runtime to select the most appropriate subset of measurements according to the context c t = {z t , i t }. This process is illustrated in Fig. 1.

The training method of the selection model can be seen as a case of discriminative training, analogous to the one suggested in [START_REF] Abbeel | Discriminative training of kalman filters[END_REF]. This is due to the optimization of the RVM model done with respect to the ultimate filter performance, while the observation model is usually optimized through generative training, i.e in order to maximize the observation likelihood p(z t |x t ) over the training set. A strong benefit of this specific training approach is the subsequent model implicit capability in compensating for mis-modeled aspects of the real system. Discriminative training has then been chosen here as it serves the common purpose of compensating for Bayes filter model inaccuracies. Note that, alternatively, the RVM model could however be trained generatively.

B. Illustration

We demonstrate here the capabilities of the selection approach on a simulated dataset. For this test case the original training set D contains 6000 samples of the simulated altitude ground truth of a UAV and the corresponding altitude measurements provided by three sensors (Fig. 2) presenting some realistic characteristics, such as maximum range threshold and outliers occurrences -section V deals with real data. From these 3 sensors, seven observation functions f c are built in order to cover the exhaustive measurement combinations: Altitude measurement provided by three sensors. Sensor 1 reproduces typical ultrasonic measures, low observation noise, outliers occurrences and maximum range threshold. Sensor 2 permanently provides measures with high observation noise. Sensor 3 does not provide relevant measures below 2 meters. on the training set then provides useful insights on the relevant sensor combinations: some functions may appear to be rarely used and therefore can be removed from the initial set before re-training a refined classifier. This process can simplify the classification problem and directly reduce the computational cost of the model at runtime. Once the classifier is trained, it is evaluated on a validation data set.

{[z 1 ], [z 2 ], [z 3 ], [z 1 z 2 ], [z 1 z 3 ], [z 2 z 3 ], [z 1 z 2 z 3 ]}.
The final estimate value is shown Fig. 3. Clearly, the classifier recognizes all different contexts which for sensor measurements could alter the estimation, while combining the maximum of reliable data so as to provide a low estimation uncertainty. In practice, this means that the classifier behavior consists in activating as often as possible the observation functions corresponding to combined measurements. The effect of measurement rejection over the final estimate uncertainty are well illustrated for the highest altitudes where sensor 1 does not provide any useful information and hence is not used. As a consequence the estimation uncertainty is slightly increased.

IV. LEARNING CONTEXT-DEPENDENT OBSERVATION NOISE MODEL

In practice, the core component described in the previous section is a powerful mechanism which is also capable of compensating further inaccuracies made in the estimation model. This is a simple consequence of the fact that the selection function has been trained to choose the appropriate measurement combination such as obtaining the best estimate output, no matter how well tuned the observation and transition models are. As explained in section III, the parameters of the observation function f are not the most complex to obtain, but the observation noise can be difficult to estimate. Moreover, introducing an accurate and realistic noise model in the filter helps in tempering the rejection scheme operated by the classifier. Indeed, and as shown in [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF], the use of heteroscedastic observation noise models can greatly improve the estimation precision while providing a more realistic estimation uncertainty. When combined with the selection scheme described in the previous section, the use of such noise models results in a better exploitation of each filter by extending their domain of 'reliability'.

An elegant approach to learn input dependent noise models in the context of regression was introduced in [START_REF] Kersting | Most likely heteroscedastic gaussian process regression[END_REF]. A key idea of this method is to introduce a second regression model dedicated to modeling the empirical 'observation' standard deviation for each sample. In a similar approach we form again N new dataset D ′′i = C, S i where

S i t = log[ (z i t -f (x t )) 2 ]
. This value provides an empirical estimation of (the logarithm of) the observation noise level for the measurement z i at time t. N RVM regression models are then trained on these different training sets.

Used along with the different functions f c these models notably prevents the different filters from diverging during the performance evaluation step preceding the classifier training. Note that the RVM model provides Gaussian prediction distribution over the noise level. In a fully Bayesian approach, evaluation of the final observation distribution p(z t |x t ) would require to marginalize over (integrate) the noise distribution. However, this integral is analytically intractable and then requires to turn to approximation methods. For computational efficiency, the integral is then replace by the most likely approximation, meaning that the noise level is approximated by a point estimate corresponding to the maximum of the prediction distribution (its mean). At runtime, computation time can be saved by only predicting noise values for the subset of measurements required by the active function f c . These functions associated with the classification and regression models then define a complete context-dependent observation model whose performance is illustrated in the next section.

V. EXPERIMENTS

Experiments on real data are performed in the context of altitude estimation for a quadrotor UAV. The paparazzi platform [START_REF] Brisset | The paparazzi solution[END_REF] was used to collect sensor measurements provided by an ultrasonic sensor and a barometer, along with other UAV internal data. The altitude ground truth is provided by a motion capture system and its data is synchronized at 50Hz with the UAV telemetry. All data is collected by manually flying the UAV such as covering the practicable space with different dynamic behaviors (smooth and more aggressive motions). The different datasets all contain around 6000 samples.

We aim at evaluating the benefits of introducing contextdependency over the estimation process. The proposed approach is then compared with the GP-PF algorithm defined in [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF]. This state-of-the-art algorithm is based on a GP modeling the whole observation model p(z t |x t ). Both filters instantiations are trained on a common dataset similar to the one shown in Fig. 4. For all datasets, the ultrasonic sensor measurements show frequent occurrences of strong outliers that we intuitively know to be caused by perturbations originating from the actuation of the 4 motors. We however never made deeper investigations on the nature of the perturbations and thus train the RVM-based filter by learning observation noise and measurement selection for a total of 3 measurement combinations: one filter based on the ultrasonic measurement, an other one based on the barometer, and a last one based on both ultrasonic and barometer sensors. The context data c t contains both sensor measurements and the motor thrust command which, according to our intuition about ultrasonic sensor perturbation, may provide useful information for the selection task. The GP-PF is trained on the raw dataset containing the ground truth and the measurements. A Kalman filter with 3 sigma rejection scheme on all measurements is also evaluated on the common validation set. All tested filters share a common constant velocity prediction model with identical parameters.

The final estimation error for all filters is shown in Fig. 5. As can be seen, the GP-PF diverged in many situations. Analysis of the trained observation GP shows that it failed to model the occurrence of outliers in the ultrasonic measurements. For some cases outliers are ignored, while for frequent occurrences phases the observation model tends to converge to the outlier values. Less importantly, one can notice that fluctuations of the barometer measurements when the motors are turned on (noticeable on the thrust command from Fig. 4) are also not generalizing well on the validation set. The Kalman filter provides a lower estimation error thanks to the rejection scheme, but still shows some cases of strong divergences, especially for the high dynamic parts of the flight. Finally the RVM-based filter shows a noticeable improvement in both terms of error amplitude and frequency.

Clearly, both Kalman filter and GP-PF reveals the need for a context-dependent knowledge during the estimation process. While the Kalman filter manages to get rid of most outliers, the rejection mechanism can not handle some specific configurations. For example this is the case for the part of the dataset contained between samples 3000 and 4000 corresponding to the highest observed vertical speed variations. Here, ultrasonic sensor outliers and barometer measurement latency become coherent, leading to divergence. Meanwhile, the GP observation model is trained to fit at best the training set within which strong measurement variations in the mea- Intuitively we understand that some information is missing when trying to fit a precise model of measurement generation (note however that in optimal conditions, this is not the case). This is were the benefits of the RVM-based approach come into play: the contextual knowledge compensates for the missing information and allows recognition of situations where the Kalman filter was unable to apply the correct rejection. The estimation output for the RVM-based filter on the validation set is shown in Fig. 6. This figure also illustrates the role of the heteroscedastic observation noise in the reliability of the final estimate uncertainty. As can be seen, the noise models provide realistic and consistent evaluation of the observation noise such that the altitude truth remains inside the final uncertainty boundaries. Conversely, as noticed in [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF], and as observed in these experiments, filters using homoscedastic observation noise such as the GP-PF and the Kalman filter provide lower and less variable uncertainties, leading to less consistent estimate output distribution.

VI. CONCLUSION

A. Summary

This article presented a simple approach for introducing context-dependent knowledge within the observation model of a Bayes filter. As such, the approach is generic as it RVM-based filter estimation output and associated 3 sigma uncertainty boundaries can be directly applied regardless the version of the Bayes filter. Based on contextual information defined by the current measurements (and any additional information if required), a first classification model is trained. Its role consists in incorporating measurements in the estimation process when they are known to help in improving the final estimate. Conjointly, heteroscedastic (in our case context-dependent) observation noise models are also learned and help the filter in keeping a realistic estimation uncertainty level. The approach is evaluated and compared to two standard and state-of-the-art methods. Evaluation shows that contextdependency allows the approach to outperform other filters in both terms of estimation and uncertainty consistency. This work can be seen as an extension of the approach proposed in [START_REF] Ravet | Learning to combine multi-sensor information for context dependent state estimation[END_REF], but still offers room for improvement, especially through the introduction of a dynamic model which could potentially greatly improve context identification. However, the current implementation has the advantage of relying on classical RVM training and inference methods.

B. Discussion

In this work, RVM models have been chosen instead of GP for computational efficiency reasons. There is a actually a strong relationship between these two models, as RVM can be seen as a special case of GP [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF]. Both models provide probabilistic outputs, i.e. each prediction is provided with an associated uncertainty. This specificity is a key feature and could be used to protect the system when it enters part of the input space where too few or no data has been used for training. Analysis of the prediction uncertainty can then allow the designer of the system to implement a failsafe behavior, either by asking the robot to stop, or by switching to a default standard estimation filter. It should however be noticed that in this situation, the RVM model is known to decrease drastically its uncertainty when entering unexplored parts of the input space [START_REF] Rasmussen | Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning)[END_REF], while GP behave more logically by increasing the uncertainty. This specific aspect requires further investigation about the applicability of such safety systems.

GP-PF and the RVM-based filter are both based on nonparametric models capable of representing dense and complex information. However, the principal difference between both approaches is that the GP-PF is based on a statedependent model, where the RVM-based filter only exploits the contextual information. Introduction of the state value as an additional contextual information has been tested and proved to give lower performance for our approach. This is a direct consequence of the fact that relying on perceptual information is a much more robust foundation than using the state value. Clearly, a state-dependent model trained with the ground truth is more likely to produce erroneous predictions, as errors on the system state are unavoidable at runtime (hence differing from values seen in the training set). Conversely, using the measurement inputs ensures the RVMbased filter to be more reliable as its prediction capabilities do not depend on the own filter performance. Especially, one has to notice that the probabilistic models applied here do not take into account the possible existence of noise over the input variables. Investigating the exploitation of such models [START_REF] Mchutchon | Gaussian process training with input noise[END_REF] might be an interesting direction for future work.

Fig. 1 .

 1 Fig. 1. Principle of mRVM classification based measurement selection for Bayesian filtering (contextual information it is not represented for clarity)

  Fig. 2.Altitude measurement provided by three sensors. Sensor 1 reproduces typical ultrasonic measures, low observation noise, outliers occurrences and maximum range threshold. Sensor 2 permanently provides measures with high observation noise. Sensor 3 does not provide relevant measures below 2 meters.

Fig. 3 .

 3 Fig. 3. Final estimation and associated uncertainty on validation set.

Fig. 5 .

 5 Fig. 4. Validation set

  Fig. 6. RVM-based filter estimation output and associated 3 sigma uncertainty boundaries


Alexandre Ravet

Simon Lacroix

Heterogeneous Bayes Filters with Sparse Bayesian Models: Application to state estimation in robotics

This study introduces a new augmented Bayes filter model for time-varying, context-dependent emission noise. The envisaged application, robust state estimation for a robot, motivates the use of the Relevance Vector Machine to model the emission noise, as it provides sparsity and fast inference capabilities. Besides the introduction of this new model, this work also aims at comparing the final filter performance when discriminative training is used instead of the prevalent generative training. The theoretical foundations for training and running inference over the model are proposed.

Introduction

Bayes filters (BF) have been widely applied to many areas. They are notably a workhorse of robotics, where recursive filtering [START_REF] Rauch | Maximum likelihood estimates of linear dynamic systems[END_REF], a fast and simple inference procedure, has provided the most common and reliable method for real-time state estimation over the past decades. BF used for state estimation usually rely on some optimistic assumptions for two main reasons: the real physical system is actually too complex to be perfectly described through a tractable model, and physical exactitude is neglected for the sake of computational efficiency. As a consequence, filter models designed for state estimation usually rely on a minimum system state representing the robot variables required for achieving a specific task. Any unmodeled aspect of the system, among which the effects of the current environment over the filter performance, are then encompassed within additional noise terms.

In practice, when used for robot state estimation, Bayes filters require a substantial tuning phase to provide acceptable performances. This is because capturing all unspecified aspects of the system through the sole introduction of noise often results in a trade-off between output optimality (accuracy of the state estimate) and robustness (to the different unmodeled aspects). An illustration of this problem, and the core motivation of this work, is the case of an autonomous robot navigating through different environments, in which one has to deal with a whole range of alterations in sensor readings, going from average (optimal) observation noise to complete failure (unreliable data). This is especially true when sensor performances are strongly affected by the different environment characteristics, such as luminosity, or texture and materials of the surroundings objects, ground and obstacles.

The classical state-space model is defined by:

x t = f (x t-1 ) + γ y t = g(x t ) + ν (1)
where x t is the latent state at time t with the associated observation y t , f and g the transition and emission functions respectively, γ ∼ N (0, Σ γ ) the system noise and ν ∼ N (0, Σ ν ) the observation noise. f and g can either be linear functions (linear dynamical system) or nonlinear (nonlinear dynamical system). Most often, no fixed noise values Σ ν that would yield an optimal output for a large variety of operating conditions (or environments) can be determined. Other unmodeled effects producing the strongest measurement alterations are often compensated with a rejection scheme, usually relying on a model self-consistency check. In other words, designing a Bayes filter emission model consists in finding the best distribution modeling the emission process for the nominal cases, and rejecting all the data that does not fit this distribution [START_REF] Sivia | Data Analysis: A Bayesian Tutorial[END_REF]. One strong consequence of this approach is that the system might converge to erroneous but model-consistent state values [START_REF] Thrun | Probabilistic Robotics (Intelligent Robotics and Autonomous Agents)[END_REF][START_REF] Ravet | Learning to combine multi-sensor information for context dependent state estimation[END_REF]. When such divergences are observed, additional parameter tuning is required to improve global robustness, then detrimentally affecting the state estimation performance. This classical state-space model, whose graphical representation is shown in Fig. 1, is known as (time) homogeneous. It can be enhanced by making the model parameters vary in time: the trade-off usually required when tuning the parameters is then no longer needed, and the resulting model can handle the whole spectrum of alterations over measurements. As the research context motivating this study concerns robust and adaptive perception for autonomous robots, this work focuses on the emission distribution -even though the proposed approach can be straightforwardly applied to the prediction distribution of model [START_REF] Abbeel | Discriminative training of kalman filters[END_REF].

This paper aims at developing a model able to compensate for the assumptions made by describing a system through the simplified emission distribution with simple Gaussian noise Σ ν -while maintaining high computational efficiency. It results in an enhanced Kalman filter capable of dealing with both moderate alterations and outliers, without requiring the implementation of rejection rules. This is done by training an additional model for the emission noise, which relies on contextual information input. One particularly appealing consequence of this approach lies in the introduction of a second order knowledge over measurements reliability, where basic rejection schemes rely on some knowledge about the data properties. The proposed approach is consequently less likely to diverge.

Discriminative training being known to help in compensating some of the mis-modeled aspects of a system [START_REF] Abbeel | Discriminative training of kalman filters[END_REF], we also aim at analyzing the impact of a discriminative learning method versus a generative one over the performance of the resulting filter. The next section depicts the model basic principle. Sections 3 and 4 respectively describe generative and discriminative training of the model. Inference methods are then provided in section 5, and a discussion concludes the paper.

2 Heterogeneous BF with sparse Bayesian models

Background

Unless the state of simple models such as (1) encompasses all the exogenous phenomena likely to alter the system behavior, BF are by nature unable to model time varying emission and prediction processes. Recently, extensions have been proposed in order to compensate this unability. Nonparametric models such as Gaussian Processes (GP) have been integrated for modeling transition and emission distributions [START_REF] Deisenroth | Robust filtering and smoothing with gaussian processes[END_REF], and extended to fully state-dependent models in [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF], through the introduction of a heteroscedastic observation noise. If nonparametric models improve the filter robustness compared to parametric functions, they are usually designed as state-dependent models, and as such are unable to handle contextual influence over the measurement process. For the heteroscedastic observation model proposed in [START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF], the presence of outliers in the training set is then critical: based on the sole state information, the system is unable to discern the contribution of the noise free model (g(x t ) in ( 1)) from the noise model (which is then written ν(x t )) within measurements. An other disadvantage of these approaches is that one has to turn to more complex sparse GP techniques when using a large training set if the system is intended to be used in real time.

From a different perspective, optimizing the parameters of a BF has always been mostly considered as a tuning task, achieved with the intuitive goal of providing the most accurate state estimation. This allows to take into consideration some mis-modeled aspects of the system, even if they are never explicitly described, neither understood. Surprisingly, methods for learning the parameters of a BF appeared quite recently in the literature, and mostly rely on maximum likelihood. Since BF are generative models, this implies that the parameters are not determined with respect to the system ultimate performance, but so as to get the best model for the underlying prediction and emission processes. Conversely, discriminative training is similar to manual tuning, in the sense that the parameters are optimized with respect to filter performance. But this latter training approach remains uncommon, although it proved to outperform manual tuning and maximum likelihood as well [START_REF] Abbeel | Discriminative training of kalman filters[END_REF]. To our knowledge, combining an augmented model with discriminative training remains untreated, while both approaches serve a similar purpose, i.e compensating for unmodeled aspects of the real system.

Heterogeneous Bayes filters with sparse Bayesian model

To overcome BF unability to deal with context influence, an augmented model is introduced, whose particularity is to explicit the context repercussion over the measurement process. It relies on an additional observation variable c t relating to the perception context. As suggested in previous work [START_REF] Ravet | Learning to combine multi-sensor information for context dependent state estimation[END_REF][START_REF] Ravet | Context-dependent bayesian filtering: an approach based on measurement selection and supervised learning[END_REF], this additional observation can consist in the joint set -or subset -of sensor measurement values y t possibly extended with any relevant contextual information i t (any other sensor measurement, robot internal data, or any information that might influence the measurement emission process). It is assumed that this joint set of measurements defines a proper representation space for the contextual influence over measurement noise, i.e there exists a mapping from the context input space to the observation noise level.

To further avoid ambiguities in the contribution of two distinct models g and ν in the measurement process, we assume that the noise free component g of the emission model is known and homogeneous in time, since it can generally be obtained directly through physical considerations about the nominal measurement generation process. This results in an emission model of the form:

y t = g(x t ) + ν(c t )
Reminding the goal of this model is to enhance a Kalman filter, ν(c t ) is then a zero mean Gaussian noise distribution with context-dependent variance:

ν(c t ) ∼ N 0, r(c t )
To avoid making any assumption over the functional form for the variance model, and keep computational efficiency, r(c t ) is modeled with the Relevance Vector Machine (RVM) framework [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF], naturally providing sparsity thanks to the automatic relevance determination mechanism. For a given training set of T observations {c i } T i=1 , we define r(c t ) = exp(z t ) to ensure variance positivity where z t is given by

z t = T i=1 w i K(c t , c i ) + w 0 + ǫ (2)
with w 0 a bias parameter, K the chosen kernel function, and ǫ ∼ N (0, σ 2 ǫ ). To foster sparsification, a zero-mean Gaussian prior is placed over the weight vector w = (w 0 , w 1 , .., w N ) T :

p(w|α) = N i=0 N (w i |0, α -1 i )
Fig. 2. Bayes Filter with heterogeneous emission noise.

where we define uniform hyperparameters priors over α = (α 0 , α 1 , ..., α N ) T . So far, the model has been depicted for a one-dimensional observation space. Real applications however require to consider the multi-dimensional case. For an observation variable y t ∈ R D , D distinct RVM models are then used to model each component of the noise covariance matrix ν(c t ). For clarity, the next sections only consider the one-dimensional case, the extension to higher dimensions being straightforward. The resulting graphical model of this augmented model is depicted in Fig. 2.

Generative training

Bayes filter parameter optimization is usually done by minimizing the likelihood of the training set [START_REF] Ghahramani | Parameter estimation for linear dynamical systems[END_REF][START_REF] Bishop | Pattern recognition and machine learning[END_REF], considering the latent state variable remains unobserved. In this work, we assume that the training set also contains ground truth data, i.e accurate values of the state variables x = {x 1 , ..., x T }. The issue of training the emission model then turns to be analogous to the regression task, and more specifically to the heteroscedastic regression task with nonparametric models [START_REF] Goldberg | Regression with input-dependent noise: A gaussian process treatment[END_REF][START_REF] Kersting | Most likely heteroscedastic gaussian process regression[END_REF][START_REF] Lzaro-Gredilla | Variational heteroscedastic gaussian process regression[END_REF][START_REF] Khashabi | Heteroscedastic relevance vector machine[END_REF][START_REF] Ko | Gp-bayesfilters: Bayesian filtering using gaussian process prediction and observation models[END_REF]. Note however that the training task is here a bit simpler since the observation function g is fixed and only the model ν(c t ) has to be determined. If sampling and variational approximation can be also used, the chosen approach relies on hard-assignment Expectation Maximization (EM) as suggested in [START_REF] Kersting | Most likely heteroscedastic gaussian process regression[END_REF]. By using hard-assignment EM we iteratively estimate the RVM parameters and predicted log noise level z at original inputs C = {c i } T i=1 . Thanks to this approximation, we are able to make direct use of classical RVM optimization and prediction equations, providing in this context the fastest solution for a real time application.

Following Kersting et al. approach [START_REF] Kersting | Most likely heteroscedastic gaussian process regression[END_REF], the hard E-step consists in empirical estimation of the noise variance. Based on real observations y = {y 1 , ..., y T } and samples y k t provided by the current observation model (using the parameters α and σ ǫ found after last EM iteration), the set of values y t and {y k t } K k=1 are seen as independent noisy observations of g(x t ). Empirical estimation of the noise variance at x t is then provided by the mean

var t = 1 2.K K k=1 (y t -y k t ) 2
In the subsequent M-step the RVM model is trained with the new training set D = {c t , log(var t )} T t=1 with a classical optimization procedure [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF]. In other words, the optimization process considers the noise variance as the hidden variable of the model, and iteratively optimize the parameters of the RVM model based on a hard assignment of the estimated noise. This method requires to use a substantial number of samples to empirically estimate the noise variance and, as any hard-assignment EM, is prone to oscillating, requiring to monitor the likelihood of the model over the training set after each algorithm iteration. It however brings an important advantage, since the optimized model, in association with the last noise variance estimation, can be readily used for prediction using classical RVM equations.

Discriminative training

The previous learning approach aims at minimizing a loss function corresponding to the emission likelihood. In other words the optimization step finds model parameters explaining at best the measurement generation process. As suggested in [START_REF] Abbeel | Discriminative training of kalman filters[END_REF], it is however better to optimize the parameters with respect to the ultimate system performance, i.e the accuracy of filter estimates. Training the model then consists in finding α max and σ ǫmax such that:

α max , σ ǫmax = arg max α,σǫ T t=1 log(p(x t |y 1:t ))
where p(x t |y 1:t ) is provided by Kalman equations. Considering f and g are linear functions corresponding to the matrices F and G respectively (classical Kalman filter case), we have:

p(x t |y 1:t ) = N (x t |µ t , σ t ) (3)
with

µ t = F µ t-1 + K t (y t -GF µ t-1 ) σ t = (I -K t G)P t-1 P t-1 = F σ t-1 F t + Σ γ K t = P t-1 G t (GP t-1 G t + Σ ν(ct) )
Since this distribution requires the evaluation of two latent variables: the log noise level and the RVM weight parameter, a procedure similar to Type II maximum likelihood is employed. z t and w are then marginalized out and we now seek for parameters α max and σ ǫmax maximizing:

L discr = T t=1 log p(x t |y 1:t , z t )p(z t |c t , w, σ ǫ )p(w|α) dw dz t
Since p(z t |c t , w, σ ǫ ) and p(w|α) are both Gaussian, the integral with respect to w is readily evaluated to give:

L discr = T t=1 log p(x t |y 1:t , z t )N (z t |0, D t )dz t
where D t = σ ǫ + K T α -1 K, with K the vector of kernel functions such that K i = K(c t , c i ) as defined in (2), and A = diag(α). The last equation is analytically intractable and is then approximated with integration by substitution and Gauss-Hermite quadrature. α max and σ ǫmax are subsequently found by conjugate gradient ascent over L discr . Note that classical optimization of the RVM model requires the computation of a design matrix containing all kernel elements evaluated at all original locations {c i } T i=1 . Here, the optimization is done separately for each kernel vector evaluated at c i , through their influence over ultimate filter performance. This different form of training (by comparison to the one in [START_REF] Kersting | Most likely heteroscedastic gaussian process regression[END_REF]) then requires to foster sparsity by thresholding of the α i values during the optimization process.

Inference

Generative training case

For the classical model of Fig. 1, filtering consists in evaluating normalized marginal distributions α(x t ) = p(x t |y 1 , .., y t ) with the recursion equation of the form:

η t α(x t ) = p(y t |x t ) α(x t-1 )p(x t |x t-1 )dx t-1 (4)
where η t = p(y t |y 1 , ..., y t-1 ) the scaling factor and p(y t |x t ) and p(x t |x t-1 ) the emission and prediction distribution considered as Gaussian for kalman filtering.

Since in the new model the noise level is considered as an additional latent variable, the emission distribution required for the evaluation of ( 4) is now given by:

p(y t |x t , c t ) = N y t |g(x t ), exp(z t ) p(z t |c t , C, z, α, σ ǫ )dz t ( 5
)
where z is the predicted log noise level at original inputs {c i } T i=1 , and p(z t |c t , C, z, α, σ ǫ ) the predictive distribution given by: p(z t |c t , C, z, α, σ ǫ ) = p(z t |c t , w, σ ǫ )p(w|C, z, α, σ ǫ )dw [START_REF] Rauch | Maximum likelihood estimates of linear dynamic systems[END_REF] This familiar predictive distribution [START_REF] Tipping | Sparse bayesian learning and the relevance vector machine[END_REF] is also Gaussian. The evaluation of the integral ( 5) is hence analytically intractable, and requires approximation. The fastest approach is the most likely approximation, where we replace the integral by N y t |g(x t ), exp(z * t ) with z * t = arg max zt p(z t |c t , C, z, α, σ ǫ ). Note that this approximation allows p(y t |x t , c t ) to be a Gaussian distribution, a necessary condition for using Kalman recursive equations.

Discriminative training case

Since discriminative training did not involve evaluation of the posterior distribution of the noise level z, we can not make straightforward use of the RVM prediction equation. Analogously to equation [START_REF] Goldberg | Regression with input-dependent noise: A gaussian process treatment[END_REF], evaluation of the emission distribution now requires to marginalize out the distribution over z t and z:

p(y t |x t , c t ) = N y t |g(x t ), exp(z t ) p(z t , z|c t , C, χ, α, σ ǫ )dzdz t (7)
with χ = {x, y}. We can decompose the last term of the right-hand side of this equation and write:

p(z t , z|c t , C, χ, α, σ ǫ ) = p(z t |c t , C, z, α, σ ǫ )p(z|χ, C, α, σ ǫ )
and, as done previously for the generative case, replace the integral in [START_REF] Kersting | Most likely heteroscedastic gaussian process regression[END_REF] with a point estimate of z t and z. Note that this approximation is mainly motivated by its low computational cost, required for online filtering.

We turn to Gibbs sampling in order to evaluate at first a point estimate for z, reminding we have:

p(z|x, y, C, α, σ ǫ ) ∝ p(x|y, z)p(z|C, α, σ ǫ ) (8)
As for standard Gibbs Sampling procedure, we then sample by cycling through the different components of z. Noting z i the i th component of z, and z \i the remaining components, we have:

p(z i |z \i , x, y, C, α, σ ǫ ) ∝ p(x t=i |y t={1:i} , z i )p(z i |z \i , C, α, σ ǫ ) (9)
where the first term of the right-hand side of ( 9) is given by the kalman equations (3) and the last term is given by the RVM prediction equation [START_REF] Rauch | Maximum likelihood estimates of linear dynamic systems[END_REF]. At first, samples are drawn from p(z i |z \i , C, α, σ ǫ ) and rejection is done according to the term p(x t=i |y t={1:i} , z i ). Intuitively, this procedure corresponds to maximizing the discriminative likelihood of the noise level z i . Once a point estimate for z is found, it can be used straightforwardly within regular RVM prediction equation for z t . The remaining of the inference is then similar to the procedure depicted for the generative training case.

The aim of this work is to define a new Bayes filter model able to encompass a variety of contexts, and to analyse different training approaches and their expected consequences over the system performance. Its specificities rely in the introduction of an additional observation used for context identification, and in the use of a sparse RVM model for context-dependent observation noise prediction. The theoretical foundations have been presented and system performances are currently being investigated. Note that the idea of introducing an additional context variable has been tested in our previous work [START_REF] Ravet | Learning to combine multi-sensor information for context dependent state estimation[END_REF][START_REF] Ravet | Context-dependent bayesian filtering: an approach based on measurement selection and supervised learning[END_REF] and proved to be relevant for the simple task of context-dependent sensor selection (equivalent to rejection). However experiments conducted so far concerned the simple task of altitude estimation for an UAV, and the approach still has to be tested on more complex scenarios involving numerous sensors and a broad range of contexts.

While augmenting Bayes filters with time-varying noise model plays a central role in trying to compensate the optimistic assumptions usually made by the classical model, the training method might also have major consequences over the system performance. As such, discriminative training seems promising in that it requires to run the filter during optimization while generative training focuses on the underlying emission and prediction processes. Discriminative training nevertheless brings some particular issues, since at first, it does not allow to use classical training (and sparsification) method for the RVM model, but also because the optimization process of L discr is much more complex. Indeed, each term of the discriminative loss function is strongly related to the preceding one as a direct consequence of the recursive equations used for state estimation. Classical RVM models already require the optimization of a nonconvex function, and we still need to study the consequences of the additional complexity introduced along with this specific loss function. Besides this particular issue, the use of discriminative training also ends up with some additional approximations during inference. Experiments will provide a better insight on how expected benefits and inherent drawbacks of the discriminative method impact the system performance, by comparison to the much simpler generative approach.

Fig. 1 .

 1 Fig. 1. Graphical model of a Bayes filter with homogeneous emission distributions.


Farouk Mansouri

Sylvain Huet

Dominique Houzet

SignalPU: A programming model for DSP applications on parallel and heterogeneous clusters

The biomedical imagery, the numeric communications, the acoustic signal processing and many others digital signal processing (DSP) applications are present more and more in the numeric world. They process growing data volume which is represented with more and more accuracy, and use complex algorithms with time constraints to satisfying. Consequently, a high requirement of computing power characterize them.

To satisfy this need, it's inevitable today to use parallel and heterogeneous architectures in order to speedup the processing, where the best examples are today's supercomputers like "Tianhe-2" and "Titan" of Top500 ranking. These architectures with their multi-core nodes supported by many-core accelerators offer a good response to this problem. However, they are still hard to program to make performance because of many reasons: Parallelism expression, task synchronization, memory management, hardware specifications handling, load balancing . . . In the present work, we are characterizing DSP applications and propose a programming model based on their distinctiveness in order to implement them easily and efficiently on heterogeneous clusters.

I. INTRODUCTION

The digital signal processing applications like computer vision, medical imagery or acoustic signal processing take a real place in many domains. Also, cause of sensors technology evolution, they process increased data volume (data length) like the high definition imagery which reach ten or so of Go. Also the data unit is represented more and more precisely (data floating point encoding) from single precision to 32 digits to quadruple precision to 128 digits. From the other side, the DSP applications use complex algorithms (time complexity) in linear, quadratic or exponential time. And are usually constrained with latency or throughput. Because of all this reason, we can sayy that DSP applications need a big computing power.

To satisfy this need, it's recommended today to focus on parallel and heterogeneous hardware architectures which include generalist multi-core processors (Intel Xeon or AMD Opteron), supported by many-core accelerators (GPU, Xeon phi, Cell, FPGA, DSP) and structured in the form of a cluster of connected nodes with a high bandwidth network. The best examples are the supercomputers like "Tianhe-2" or "Titan" related in the Top500 list.

Certainly, these architectures can get high performance computing to satisfy DSP applications requirements. However, they present some difficulties of use. In fact, to make performance with that last, the programmer has to deal with heterogeneous computation units using different languages or API, he has to manage synchronization, memory allocation, data transfers and the load balance between the processors. Thus, programming models are necessary to hide all this hardware specifications, and produce easily and efficiently the desired performance.

In the present work, we present first (Section II) the state of the art of models of programming and their classification. Afterwards, we describe in the section III the digital signal processing applications and highlight their characteristics. Then, we give the features to efficiently implement them in the next section (Section IV). In the section V we propose a programming model and explain its conception parts. Finely in the section VI we present experimentations and results of applying our model of programming (MOP) on a real world application.

II. STATE OF THE ART OF PROGRAMMING MODELS

A programming model is a programming concept providing some abstractions to mask (hide) hardware specificities of parallel and heterogeneous architectures in order to easily and efficiently express high level applications. In other words, it's an interface between the hardware complexity and the software features which represent a base for the high level program. Thus, they increase productivity, portability and performance. They exist in the form of languages, language extension or library. Below, we present their classification under three axis:

A. Classification according to abstraction level

One of the important purposes of MOP is to make easier the implementation of high level applications on a complex architecture. That is possible by bridging between the hardware functionalities like the memory management, the tasks scheduling, the data copy and the basis stand of the high level algorithm. Therefore, a good programming model has to abstract the architecture details for the programmer. Thus, we distinguish two abstraction levels for today's MoP :

Low level abstraction: In this kind, the programmer has to handle manually a lot the hardware functionalities like allocation and freeing of memory, creation, submission and synchronization of tasks, data transfer, load balance. . . Thus, this proximity of architecture gives him more control to tune the execution but reduce productivity.

High level abstraction: By contrast, in this kind, the programmer is freed of several of these constraints because the programming model does this for him.

B. Classification according to the communication model

In this classification we focus on the communication sorts used by models of programming in order to synchronize or to copy data between the tasks executed on several computing units composing the cluster. Here, we distinguish also two kinds :

Shared memory communication: In this case, the programming model assumes the data copy and the synchronization of tasks on parallel architectures which are based on shared memory like multi-core or many-core computing units.

Message passing communication (Distributed memory):

This kind of MOP assumes the data copy and the task synchronization onto the distributed architectures by exchange messages between the tasks through the network connecting nodes of the cluster.

C. Classification according to the parallelism and heterogeneity expression

To deal with the parallel and heterogeneous architecture, the programmer has to express the part of the code which will be executed in parallel or on the accelerator by programming it inside the application. However, thanks to programming models, he can do that with more easier manners. For example, by function calls in the API, by directive insertion for language extension or by compilator's options added for some languages. We propose to classify this kind of expressions under two sorts:

Explicit expression: In this case, the programmer has to explicitly express the part of the program to be executed in parallel by adding at the portion of code inside. For example, insert some directives or call functions from an API.

Implicit expression: Unlike the previous case, here, the expression of the parallelism or heterogeneity is implicitly made by the programmer. For example, by giving the task graph of the application or by adding some compilation options.

In the figure Fig. 1, we present some programming languages, extension of language and API, representing software aspects of programming models and their characteristics classified according to two axis:the communication model and the abstraction level.

1) Languages:

CUDA [START_REF] Sanders | CUDA by Example: An Introduction to General-Purpose GPU Programming[END_REF](Compute Unified Device Architecture): Is a proprietary programming model invented by NVIDIA to implement programs in their graphical processing unit (GPU). It's proposed in the form of a language inspired from C that the programmers can insert on their C/C++ code to offload in generally the SIMD portion. It's a low level abstraction MOP, and support only the shared memory architecture. OpenCL [START_REF] Munshi | OpenCL Programming Guide, ser. OpenGL[END_REF](Open Computing Language): Unlike CUDA, it's an open standard language based on C99 with some limitations and additions. It allows for the programmers to address heterogeneous platforms consisting of CPUs, GPUs, DSPs, FPGA, and others processors. Also, this MOP offers a low level abstraction and is designed for shared memory architectures but allows the both task and data parallelism.

PGAS [START_REF] Chen | Optimizing partitioned global address space programs for cluster architectures[END_REF](partitioned global address space): Is a parallel programming model which assumes a global memory address space for distributed architectures that is logically partitioned. Each portion of it is locally affected to each process or thread. It is a base for many programming languages like Chapel [START_REF] Chamberlain | Parallel programmability and the chapel language[END_REF] developed by Cray, X10 by IBM [START_REF] Saraswat | X10 language specification[END_REF], Fortress [START_REF] Allen | The Fortress Language Specification[END_REF] by Sun or Unified parallel C [START_REF]Upc language specifications, v1.2[END_REF] by UPC consortium.

StreamIt [START_REF] Gordon | Compiler techniques for scalable performance of stream programs on multicore architectures[END_REF]: Is a programming language and a compilation infrastructure, specifically engineered for streaming systems. It is designed to facilitate the programming of large streaming applications, as well as their efficient mapping to a wide variety of target architectures like CPUs, GPUs, DSPs and FPGA.

2) Extension of language:

OpenMP [START_REF] Chandra | Parallel Programming in OpenMP[END_REF](Open Multi-Processing): Is a high level abstraction MOP in the form of a set of directives and environment variables that allow for programmers to model his algorithms for data parallelism or task parallelism on a shared memory architecture. The programmer has to annotate his code by inserting some extensions (directives) which are identified at compilation time to generate supplementary code and executed on both hardware, CPUs and accelerators.

OpenAcc [START_REF] Kirk | Programming Massively Parallel Processors: A Hands-on Approach[END_REF]: Like OpenMP, it's a programming standard for parallel computing based on PRAGMA directive but more oriented to data parallelism on shared memory architectures containing accelerators.

OmpSS [START_REF] Bueno | Productive cluster programming with ompss[END_REF]: Is another variant of OpenMP extended to support asynchrony, heterogeneity and data movement for task parallelism. As OpenMP, it is based on decorating an existing serial version with compiler directives which are translated into calls to a runtime system in order to manage the parallelism extraction and the coherence and the movement of data.

3) Application interface programming (API):

MPI [START_REF] Pacheco | Parallel Programming with MPI[END_REF] (Message Passing Interface): Is the most used MOP for distributed memory architecture. It consist of a specific set of routines (i.e., an API) directly callable from supported languages (C/C++/Fortran . . . ). Using MPI, the programmers can distribute execution of their code on many nodes connected through a TCP/IP network.

TBB [START_REF] Reinders | Intel threading building blocks -outfitting C++ for multicore processor parallelism[END_REF] (Threading Building Blocks): Is a C++ template library developed by Intel for writing software programs under task parallelism model and executing it on multi-core processor according to the dependency graph. It also schedules the processing to balance the load.

StarPU [START_REF] Augonnet | StarPU: a Runtime System Scheduling Tasks over Accelerator-Based Multicore Machines[END_REF]: In the same manner, StarPU is a task programming library but for hybrid architectures CPU-GPU, and allows also for programmers to express the data parallelism. Using some routines and data structures, the programmer constructs a graph of tasks which are optimally scheduled and executed on heterogeneous and distributed memory cluster.

III. DIGITAL SIGNAL PROCESSING APPLICATIONS

Nowadays the digital signal processing applications know a great boom and are present in a lot of domains. They are a repetitive (iterative) processing of data set of input digital signal for producing an output signal or a result, as shown in the figure Fig. 2. In the algorithmic aspect, it represents a main loop which iteratively process all input data units.

Fig. 2: DSP application's form illustration

The processing part represents the firing of each data unit composing the input signal using some operators (kernels), where each operator represents an independent function with some input and output arguments. In the algorithm, that's the function calls inside the main loop. We can model this part with a data flow graph, where the vertex of the graph represents the kernels and the edges represent the data trading between operators in the form of flow. The figure (Fig. 3) illustrates this model. In this work, we limit ourselves to the treatment of synchronous data flow graph [START_REF] Lee | Static scheduling of synchronous data flow programs for digital signal processing[END_REF], [START_REF] Bhattacharya | Parameterized dataflow modeling for dsp systems[END_REF].

Fig. 3: Data flow graph (DFG) model of DSP application

In majority of DSP applications, multiple input data units are processed with the same set of instructions. For example, in video processing, each input image is processed in the same manner using the same algorithm. So, it's possible to run multiple data flow graph (DFG) in the same time on a parallel architecture. Also, each unit of data goes through the operators as input arguments and comes out as output arguments. Thus, depending on the DFG organization, some operators (kernels) must be executed in order, but some of them may be executed separately and in parallel. On the other hand, according to the data kind and the kernel algorithm, it could be interesting in some cases to offload execution of certain of kernels on a massively parallel computation unit (accelerators) like the GPU, Xeon Phi or Cell in order to enhance performances.

Taking into account these characteristics, it's possible to deduce some rules to apply in implementation and execution of DSP applications on parallel and heterogeneous architectures:

Task parallelism: Using the DFG model, highlight the dependencies between the tasks and detect the tasks able to be executed in parallel.

Data parallelism: Identify the tasks according to their Flynn taxonomy [START_REF] Flynn | Some computer organizations and their effectiveness[END_REF]. The MISD (Memory bounds) tasks are oriented to generalist processors and the SIMD tasks (Compute bounds) towards the accelerators.

Graph parallelism: In order to optimize the occupancy of computing units composing the cluster, deal with several graphs where each of them processes one input data unit.

In the next section (Section IV), based on extracted distinctiveness and the rules cited above, we focus on the implementation side of the DSP applications, and we discuss on which programming model is more suitable for programmers to easily and efficiently porting these applications on parallel and heterogeneous architectures.

IV. IMPLEMENTATION OF DSP APPLICATIONS

As presented in the preceding section (Section III), the DSP applications have some characteristics, that the programmers must exploit in order to take advantage of targeted heterogeneous and parallel architectures. First: To highlight the kernels able to be executed in parallel (task parallelism), they have to express their algorithm in the form of a set of tasks using threads or process technologies. They have to manage these threads for communicating between them or to be synchronized according to the application's dependencies on the both shared and distributed memory architectures. In addition, to profit from the accelerator's capacity to speedup the SIMD processing (data parallelism), the user has to offload a part of their task towards these compute units. To do this, he has to deal with memory allocation on accelerators, copy-in the input data, lunch the execution, copy-out the results and finally freedom the used memory zone. Also, the DSP applications are mostly iterative, so it's a good idea to unroll the main loop of the application and therefore process a number of data units in the same time (graph parallelism) in order to increase the occupancy of computing units. To do that, the code writers must duplicate the process (thread) in charge of executing the main loop taking care to guarantee the data coherence by restricting some variables or sharing others. All this implementation features are necessary for porting DSP applications on heterogeneous clusters but not enough to optimize productivity of the hardware. In fact, the programmer must cope with others difficulties like communication cost which must be masked by overlap it with the computation, or the load balancing between the computational units which must be assumed by a good scheduling of tasks.

Applying all these implementation rules is very hard. The programmer has to combine the handling of some API, language or extension of language which are low level for certain or restricted to specific hardware for others. For example, the programmer has to use Pthread , TBB or OpenMP to generate threads and express task parallelism on each node of the cluster (shared memory architecture), but also the MPI or PGAS model to manage them by creating processes onto many nodes (distributed memory architectures). He has to use CUDA or OpenCL to address accelerators like the GPU, Cell or Xeon-Phi and offload a part of a SIMD work on it. In the other case, the higher level tools like OpenACC, OmpSS or StarPU which are based on the low level tools, must be the solution. They offer more abstraction of the hardware and can target the complete cluster. But some of them are restricted to a particular model of programming, for example OpenAcc express only the data parallelism. Others of them like OmpSS are rather oriented to decorating an existing sequential code by inserting some PRAGMA directive and transforming it at compilation time into a parallel code. The rest, based on API like StarPU is, in our opinion, the most adapted programming models to implement high level applications on heterogeneous cluster. It offers an interface based on a large routines and structures which the programmers can use to design their applications, and in addition proposes a runtime which manages the tasks, their dependencies and dynamically schedule their executions on the architecture. However, it's not adapted (specified) to DSP applications as characterized in the section III with their iterative and repetitive form, and also it's still complicated to handle because of the number of routines and data structures proposed to the user as interface to implement their applications. Because of these reasons, we propose in the next section, a programming model based on a data flow graph model to make easier the application modeling and automatize the generation of the directional acyclic graph of tasks (DAG) in order to adapt StarPU to the implementing of DSP applications on heterogeneous cluster.

V. OUR MODEL OF PROGRAMMING

In this section we propose a programming model as an extension of StarPU's application programming interface (API) [START_REF] Augonnet | StarPU: a Runtime System Scheduling Tasks over Accelerator-Based Multicore Machines[END_REF], which allows for programmers to express easily and efficiently DSP applications on parallel and heterogeneous cluster. In comparison with the state of the art classification described in section II, our MOP is a high level abstraction concept. The programmers don't have to worry about several architecture specificities, like memory management, task creation, computing units profiling etc. . . They can design their algorithms for tasks and data parallelisms. Also, because it's based on StarPU, our MOP take in charge the both shared and distributed memory architectures, and deal with many-node cluster using the messages passing interface (MPI). Finally, it's implicit parallelism and heterogeneity expression. In fact, thanks to the DFG application design, the programmer is free to specify in the code which tasks will be executed in parallel and on accelerators. 4). First level: the DFG-XML interface with which the programmer can describe his application in the form of DFG, he can specify the operators composing the processing (function, arguments, architecture kind) and the data flow information (size, type, link). Thus, he is saved to manipulate the StarPU's API for creating tasks, for managing buffers between each couple of task, or for submitting jobs onto the corresponding computation unit. Second level: the application design. In this step, the DFG is transformed into a directional acyclic graph (DAG) of task using graphs unfolding techniques [START_REF] Parhi | Static rate-optimal scheduling of iterative data-flow programs via optimum unfolding[END_REF], where tasks are linked according to data dependencies. Also here, the user doesn't have to deal with the API to unroll the main loop. Third level: the StarPU runtime is used to physically manage the set of tasks and execute them on the cluster according to an optimal dynamic scheduling.

Below, we describe with more details and through a synthetic application example the 3 steps of our programming model :

A. Level 1: SignalPU DFG-XML interface

Let's take the synthetic application described in the algorithm 1

Algorithm 1 Synthetic DSP application

Input: Number of iterations (N br). Input data set (Data in ). Output: Output data set (Data out ).

1: for each Data unit in Data in do 2:

V ar 1 ← P roducer(Data unit )

3:

V ar 2 ← kernel 1 (V ar 1 )

4:

V ar 31 ← Kernel 2 (V ar 2 )

5:

V ar 32 ← Kernel 3 (V ar 2 )

6:

V ar 4 ← Kernel 4 (V ar 31 , V ar 32 ) 7:

V ar 5 ← Kernel 5 (V ar 4 , V ar ′ 4 )

8:

V ar ′ 4 ← V ar 4

9:

Consumer(V ar 5 ) 10: end for /// Data unit may be an image or a sample /// In this step, the programmer has to express his application using the DFG-XML interface description. First, he has to describe each kernel (functions) in the algorithm in the form of graph node (vertex) using the given XML structure. He has to put the number of input and output arguments, the architecture kind corresponding to the kernel (CPU, GPU, Cell, Xeon Phi . . . ). Second, he has to describe, using an XML structure, the data flows between the kernels in the form of graph edges with a weight corresponding to the type and size of the data which is traded between kernels. As result, a data flow graph of application is produced as shown in the figure Fig. 5. T set ← Link task ()

7:
end for 8:

DAG ← Submit tasks (T set ) 9: end for

C. Level 3: StarPu runtime

In this step, we use StarPU runtime to manage the DAG of tasks generated in the previous level and to dynamically schedule them taking into account the dependencies between the tasks in order to guarantee the data coherence. The figure below (Fig. 7) show how StarPU runtime does that:

At runtime, all the submitted tasks are stored on a stack, each task has an affected status ( ready to be executed, waiting some things, in firing, terminated . . . ). StarPU runtime piles up the "ready to execute" tasks on another task list which represents the input set of tasks for scheduler, which dynamically affects them to be executed on the appropriate computation unit using some heuristic algorithms in order to determine the best time for each task. One of these algorithms is the Work Stealing (WS) [START_REF] Augonnet | StarPU: a Runtime System Scheduling Tasks over Accelerator-Based Multicore Machines[END_REF] is based on estimated transfer times Fig. 6: The DAG of tasks of synthetic DSP application of task's data to predict which compute unit will be the best for each task. Another policy like the Heterogeneous Earliest Finish Time (HEFT) [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] uses in addition, the estimation of execution times of each task on each device in order to produce more efficient scheduling.

With these three steps of our proposed MOP, the programmer can easily implement his DSP application and applies the optimizations: Tasks parallelism (TP) by extracting the tasks in the DFG which be able to be executed in the same time. Data parallelism (DP) by off-loading some tasks on (SIMD) accelerators. Graph parallelism (GP) by overlapping the processing of some graphs. And the load balance tanks to dynamic scheduling using StarPU runtime. In this section we present you the real world experimentations in order to validate our approach and demonstrate the interest of its usage. We use the saliency application to process a set of images on the heterogeneous CPU-GPU architecture. First, we describe the saliency application and give its algorithm. Then, we explain its implementation using our programming model. And finally, we give the results and discuss their impacts.

A. The saliency application

Based on the primate's retina, the visual saliency model is used to locate regions of interest, i.e. the capability of human vision to focus on particular places in a visual scene. The implementation that we use is the one proposed by [START_REF] Itti | A model of saliency-based visual attention for rapid scene analysis[END_REF] as shown in figure Fig. 8. His algorithm (Algorithm 3) is: First, the input image (r-im) is filtered by a Hanning function to reduce intensity at the edges. In the frequency domain,(cf -f im) is processed with a 2-D Gabor filter bank using six orientations and four frequency bands. The 24 partial maps (cf -maps[i; j]) are moved in the spatial domain (c-maps[i; j]). Short interactions inhibit or excite the pixels, depending on the orientation and frequency band of partial maps. The resulting values are normalized between a dynamic range before applying Itti's method for normalization, and suppressing values lower than a certain threshold. Finally, all the partial maps are accumulated into a single map that is the saliency map of the static pathway. for j ← 1 to frequencies do 5:

cf maps[i, j] ← GaborF ilter(cf f im, i, j)

6: c maps[i, j] ← IF F T (cf maps[i, j]) 7:
r maps[i, j] ← Interactions(c maps[i, j])

8:
r normaps[i, j] ← N ormalizations(r maps[i, j])

9:
end for 10: end for 11: saliency map ← Summation(r normaps[i, j])

B. The SignalPU implementation

To implement the application with our programming model, the first step is to model its algorithm (Algorithm 3) given before in the form of DFG-XML using the Sig-nalPU interface. For this, we represent each of all functions (Hanningf ilter(), F F T (), GaborF ilter(), IF F T (), Interactions(), N ormalizations(), Summation()) with a node in the graph including the characteristics of each of them (architecture kind, input arguments, output arguments). Then, we represent the data flow between each twice kernels with an edge in the graph including its characteristics (data type, data size). In the figure (Fig. 9) we present the DFG resulting of this step. At runtime, the XML-DFG description of the saliency application is analyzed and a DAG of independent tasks is iteratively generated, where each task represents the execution of each kernel's code (function's code) for each image on the corresponding computation unit (CPU,GPU). Thus, we haven't to use the StarPU's API for describing the application's tasks. Also, we have not to manage the buffer's allocating and freeing. We haven't written the main loop which processes the set of input images, and don't have to unroll it. The StarPU's API is almost entirely masked.

C. The results

In this subsection, we present results of two experimentations where we show the performance provided by the implementation with our programming model using the optimizations which we propose (graph unfolding + dynamic scheduling), compared to the same implementation with a static scheduling, and without graph unfolding. The aim is to give you an indication of the performance gained using our model of execution compared to another classical implementation. The architecture used for experimentations is a heterogeneous CPU-GPU node composed of a 4 cores CPU (intel-i7 core) and 3 GPU (1 NVIDIA GeForce GT 610, 1 NVIDIA Quadro 4000, 1 NVIDIA GeForce GTX TITAN). Note that the GPU devices are ordered from the lowest to the highest powerful. For all experimentations, we decide to put 2 tasks per device as a static scheduling. That's not the optimally placement, but the goal is to compare a dynamic scheduling to a naive user placement without profiling device's power. For the first experimentation, we present first in the figure Fig. 10 the evolution of the performance (Time of execution) as a function of number of processed images (iterations) in a static scheduling context. In the green curve (marked with 'x'), we show the evolution of execution time (in second) necessary to process each set of images using our programming model with graph unfolding which allow to process several graphs in the same times. In the blue curve (marked with '*'), we show the same result of implementation but without graph unfolding techniques. The comparison between them is represented on the dotted black curve as speedup obtained between both results which reach and is stabilized on 1.04 x. In the figure Fig. 11, we present the same experimentation (evolution of time of execution) as a function of number of processed images (iterations)) but here, we combine the dynamic scheduling and the graph unfolding optimizations. So, here the speedup obtained between both results reach 3.5 x.

We explain this variance of performance gains in both parts of this first experimentation by the increasing of the number of tasks executed in parallel. In fact, the present application is completely sequential as shown in the figure Fig. 10, but thanks for the graph unfolding, it's possible to process a number of data flow graphs and several images are executed at the same time. However, in the first part of this experiment (Fig. 11), the tasks are statically putted on the architecture and each kernel of the DFG model of the application is associated with a computing unit. Therefore, each compute unit processes all the time the same size of the task, consequently, one of them (which processes the biggest tasks) creates a bottleneck, and it delays the processing. In contrast, with using the dynamic scheduling, it's possible to more take advantage of vacancy of all computation units and enhance the processing.

In the second experiment, we process a fixed number of images (100 images) but we vary its size (height and width). The figure Fig. 12 shows the evolution of execution time in tow implementations: with dynamic scheduling, and with static scheduling. Also here, the green curve (marked with 'x') represents the evolution of execution time of the set of images processed by the application with dynamic scheduling but without graph unfolding, as function of image size. And in the blue curve (marked with '*'), we present the same result obtained by applying static scheduling. The difference between both results is shown in the black dotted line as a speedup that we get by using our programming model which reach only 1.19 x. The figure Fig. 13 represent the second part of experimentation, shows the same result where we add the graph unfolding optimization to dynamic scheduling in the the green curve (marked with 'x'), and keep the second implementation (marked with '*') in static scheduling and without graph unfolding. Also here, the difference between both results is shown in the black dotted line as a speedup which grows-up to 2 x in this case.

The explanation of the performance increasing (from 1.19x to 2x) we can give is: In the first part of the experimentation (Fig. 12), the use of dynamic scheduling without graph unfolding restrict the performance because there is not enough tasks to process inside each iteration of the main loop, in order to optimally occupy all computing units. In fact, by using the graph unfolding in the second part of experimentation (Fig. 13), it's possible to overlap the processing of several iterations (graphs), so the runtime continually distribute works, and optimally takes advantage of hardware availability. In addition, we note that the dynamic scheduling is more efficient when data is bigger. That's due to the reduction of the overhead time compared to execution and communication times of processed tasks.

From this result, we can say that it's it's profitable to use graph unfolding techniques combined with a dynamic scheduling to process the DSP applications as we propose in our programming model in order to better take advantage of computation unit availability.

VII. CONCLUSION

This paper presents our proposed programming model for parallel and heterogeneous architecture, which allows a high level abstraction from the hardware specificities, and at the Fig. 12: The performance as a function of size of processed images without graph unfolding Fig. 13: The performance as a function of size of processed images with graph unfolding same time increases the performance of the implemented application. First, we presented the state of the art of MOP and classified them according to different axis. We described the DSP applications and specified their characteristics in order to implement their with optimal manner. For that, we proposed an Xml interface to easily describe DSP applications in the form of a DFG model. Also, we proposed a strategy based on unfolding graph techniques to construct a directional acyclic graph (DAG) of tasks which we process using StarPU on a heterogeneous and parallel architecture with a dynamic scheduling. Finally, we experimented our MOP on the realworld saliency application and shown that's easier to use our programming model to design it, but at the same time, it's possible to optimally take advantage of architecture's power to speed up the execution using the optimizations: Task parallelism, data parallelism, graph parallelism, and dynamic scheduling.

Fig. 1 :

 1 Fig. 1: Architecture and models of programming

Fig. 4 :

 4 Fig. 4: SignalPU design : Three levels of processing

Fig. 5 : 2 : 5 :

 525 Fig. 5: The DFG-XML of the synthetic DSP application

Fig. 7 :

 7 Fig. 7: StarPU runtime's levels

Fig. 8 :

 8 Fig. 8: The static pathway of the visual saliency application

Algorithm 3

 3 Static pathway of visual model Input: An image r im of size w • l Output: The saliency map 1: r f im ← Hanningf ilter(r im) 2: cf f im ← F F T (r f im) 3: for i ← 1 to orientations do 4:

Fig. 9 :

 9 Fig. 9: The DFG-XML model of the visual saliency application

Fig. 10 :Fig. 11 :

 1011 Fig. 10: The performance as a function of number of processed images with static scheduling