Abderrahim Kasmi
email: a.kasmi@sherpa-eng.com

Johann Laconte

Romuald Aufrere

Dieumet Denis
email: d.denis]@sherpa-eng.com

Roland Chapuis

End-to-End Probabilistic Ego-Vehicle Localization Framework

Keywords: Autonomous vehicle, localization, lane marking, map matching, road prior, fusion framework. !

Locating the vehicle in its road is a critical part of any autonomous vehicle system and has been subject to different research topics. In most works presented in the literature, ego-localization is split into three parts: Road level-localization consisting in the road on which the vehicle travels, Lane level localization which is the lane on which the vehicle travels, and Ego lane level localization being the lateral position of the vehicle in the ego-lane. For each part, several researches have been conducted. However, the relationship between the different parts has not been taken into consideration. Through this work, an end-to-end ego-localization framework is introduced with two main novelties. The first one is the proposition of a complete solution that tackles every part of the ego-localization. The second one lies in the information-driven approach used. Indeed, we use prior about the road structure from a digital map in order to reduce the space complexity for the recognition process. Besides, several fusion framework techniques based on Bayesian Network and Hidden Markov Model are elaborated leading to an ego-localization method that is, to a large extent, robust to erroneous sensor data. The robustness of the proposed method is proven on different datasets in varying scenarios.

INTRODUCTION

O VER the past decades, the automotive industry has been growing strongly. The demand for the driver and pedestrians security coupled with the technological advances are the main factors behind this exponential growth. In that regard, Advanced Driver Assistance Systems (ADAS) have spread out. The main mission of these systems is to ensure that driver safety is constantly guaranteed. For this purpose, multiple applications have been deployed, such as lane departure warning, lane keeping assist, pedestrian detection, collision avoidance, or lane change assist system. To achieve this mission, the faultless knowledge of the localization of the ego-vehicle with regards to the surrounding environment is necessary. The ego-localization of a vehicle is a mandatory component for a safe autonomous driving. Direct applications range from the decision-making system [START_REF] Iberraken | Multi-level bayesian decision-making for safe and flexible autonomous navigation in highway environment[END_REF] regarding the decisions to make in order to keep lane or to change lane, to navigation: including path planning and vehicle control [START_REF] Claussmann | A review of motion planning for highway autonomous driving[END_REF]. Ultimately, the lane/road understanding demands in terms of precision and false alarm rate [START_REF] Hillel | Recent progress in road and lane detection: a survey[END_REF] vary from one application to another. Therefore, the ego-localization solution must suit perfectly the localization requirement for each application. The ego-localization task has been widely tackled over the years. The current literature is teeming with solutions that address this issue in a variety of manners. However, one interpretation of ego-localization consists of the knowledge of three key components:

FirstName.Lastname@uca.fr I Road level localization: the road on which the vehicle travels. II Lane level localization: the position of the host lane within the road (the lane on which the vehicle travels). III Ego-lane level localization: the position of the vehicle in the lane in terms of lateral and longitudinal position.

For the road level localization, digital maps (Google, Open-StreetMap (OSM) or Waze) are used to perform this task, GPS receivers are used to retrieve the geographic (latitude, longitude) coordinates, and map-matching procedure is performed in order to match the position of the ego-vehicle with the correct road ('link'). However, the accuracy of the localization obtained is in the order of meters. Indeed, according to the Federal Aviation Administration (FAA) GPS Performance Analysis Report [START_REF] Team | Global positioning system (gps) standard positioning service (sps) performance analysis report[END_REF], the accuracy of a standard GPS device is within 3m with a 95% confidence, which can not be sufficient for most ADAS that require a more precise localization.

For some applications like lane-keeping, knowing the road on which the vehicle is traveling is not sufficient. These systems must be informed about the position of the host lane in the road to provide the adequate maneuver instruction and maintain the vehicle safety. Further, autonomous vehicle applications need a more accurate localization, which can be translated by the knowledge of the lateral and longitudinal position of the vehicle in the ego-lane. For instance, overtaking maneuvers need a faultless knowledge of the lateral position of the ego-vehicle with respect to the ego-lane marking in order to decide whether the vehicle should overtake the obstacle or not. The task of vehicle localization is still challenging for an autonomous vehicle, a complete ego-vehicle localization must perform all the three key components described above. Thus, in this work, an end-to-end solution for egolocalization from Road level localization to Ego-lane level localization is presented.

RELATED WORK

Due to the vastity of this topic, this section will discuss the main vehicle localization techniques relevant to our presented work. To do so, the current literature is divided into three distinct topics: Localization on a map, Localization on a road and Localization on a lane.

Localization on a Map

For Autonomous Vehicles, map-matching (MM) algorithms are a critical piece on any localization system. In essence, MM algorithms integrate the geographical position from a GPS receiver with the spatial road network to identify the correct road on which the vehicle is traveling. In [START_REF] Quddus | Current map-matching algorithms for transport applications: State-of-the art and future research directions[END_REF] a comprehensive literature survey of common MM methods is presented. The basic method combines an inertial measurement unit (IMU) with a classic GNSS receiver (e.g., GPS, Galileo) into a Kalman Filter to choose the correct road [6] [7]. In the same manner, authors in [START_REF] Bierlaire | A probabilistic map matching method for smartphone gps data[END_REF] propose a probabilistic MM that takes into account both the spatial (geographical coordinates) and temporal information (speed and time) to determine the likelihood road on which the ego-vehicle is traveling. However, in more complicated cases, such as urban scenarios, multi-path GPS signals lead to a lack of accuracy. In this case, the use of topological and geographical information enhances the proposed algorithm. Towards this end, many studies proposed a Hidden Markov Model (HMM) to take into account the measurement noise [START_REF] Newson | Hidden markov map matching through noise and sparseness[END_REF]. The authors in [START_REF] Goh | Online map-matching based on hidden markov model for real-time traffic sensing applications[END_REF] introduce a HMM to robustify the proposed MM algorithm. The HMM incorporates the topological and the geographical road network information in the transition state to ensure the temporal connectivity between links. A more recent approach presented in [START_REF] Atia | A low-cost lane-determination system using gnss/imu fusion and hmm-based multistage map matching[END_REF] utilizes the HMM for the MM task. In essence, the two approaches are similar, except for the HMM modeling. Indeed, in [START_REF] Atia | A low-cost lane-determination system using gnss/imu fusion and hmm-based multistage map matching[END_REF] the transition from one 'link' or segment to another must obey some rules, these rules include the dynamics of the vehicles and the legal and logical connectivity of the road map. In a more specific way, the latest developments use the OpenStreetMap (OSM) database to perform the MM. In our previous work [START_REF] Kasmi | Map matching and lanes number estimation with openstreetmap[END_REF], a multi-criteria map-matching algorithm based on multiple probabilistic criteria has been introduced. Nevertheless, the road map topology has not been properly operated in the MM process, which consists of one of the contributions of this paper.

Localization on a road

Knowing the position of the ego-lane is still a challenging task for any localization system and subject-matter of research. Indeed, lane-level localization has been widely discussed in recent literature. According to [START_REF] Kim | End-to-end ego lane estimation based on sequential transfer learning for self-driving cars[END_REF], ego-lane localization can be performed in two different manners: model-driven and deep leaning approaches.

Model-driven approaches

The current literature abounds with ego-lane recognition techniques based on model-driven approaches. In these approaches, road level features are extracted from images via feature extraction, mainly two features (edge and color) are used. For the first feature, filters are used to extract lanes edges, then the resulted outputs are fed into Hough Transforms [START_REF] Borkar | A layered approach to robust lane detection at night[END_REF] or Ransac methods [START_REF] Jiang | Computer vision-based multiple-lane detection on straight road and in a curve[END_REF] to detect lanes in the road. For the color features, it exploits the primary color or the direction of lane components (arrow, zebra marking). Once these features extracted, they are fed into a high-level fusion framework. As an example, in [START_REF] Nieto | Robust multiple lane road modeling based on perspective analysis[END_REF] ego-lane recognition is achieved from multiple-lanes detection. Adjacent lanes are first hypothesized assuming the same curvature and lane width in the road, then tested using a video-based system. This approach shows interesting results. However the occlusion of lanes marking by obstacle vehicles can not be explained and will jeopardize the output of the egolane. More high level features systems fuse visual features from the road scene: surrounding vehicles [START_REF] Svensson | Ego lane estimation using vehicle observations and map information[END_REF], lanes marking [18] [19], road clues such as arrows marking [START_REF] Popescu | Lane identification and ego-vehicle accurate global positioning in intersections[END_REF], lane marking colors [START_REF] Lee | Accurate ego-lane recognition utilizing multiple road characteristics in a bayesian network framework[END_REF], lane marking and adjacent vehicles [START_REF] Kasmi | Probabilistic framework for ego-lane determination[END_REF]. These visual cues are fed into a fusion framework: Bayesian filter (BF) [START_REF] Svensson | Ego lane estimation using vehicle observations and map information[END_REF], Bayesian Network (BN) [20] [18] [21], Hidden Markov Model (HMM) [START_REF] Ballardini | Ego-lane estimation by modeling lanes and sensor failures[END_REF] or combination of BN and HMM [START_REF] Kasmi | Probabilistic framework for ego-lane determination[END_REF]. A worth mentioning work is the one presented in [START_REF] Li | Lane-level map-matching with integrity on high-definition maps[END_REF] in which a lane-level accurate map is used in order to match the correct localization of the ego-vehicle. The map contains center lines of every lane of the road in addition to the lanes marking. Therefore to know the lane on which the vehicle is travailing a map matching procedure is presented. The presented method deals with the ambiguities encountered in choosing the right 'link'. Nevertheless, this method is tributary to the accurate map. In contrast, we present a method for lane level localization that relies only on a camera and a coarse map (OSM).

Deep learning approaches

Convolutional Neural Networks (CNNs) have also been widely adopted for lane level localization. These techniques do not take any prior about the surrounding environment, which makes these unrelated to the type of road-scenario (e.g., highway, urban). However, a learning phase is needed to set the weights of the used network.

In [START_REF] Huval | An empirical evaluation of deep learning on highway driving[END_REF] local context of the scene is used to consider occlusion to determine the end point of the local lane segment, egolane is deduced once lane segments are extracted. Another research [START_REF] Li | Deep neural network for structural prediction and lane detection in traffic scene[END_REF] uses a combination of CNN and RNN to detect lane boundaries. Similar to this approach a CNN is trained in an end-to-end way to estimate the ego-lane [START_REF] Gurghian | Deeplanes: End-to-end lane position estimation using deep neural networksa[END_REF]. In [START_REF] Kim | End-to-end ego lane estimation based on sequential transfer learning for self-driving cars[END_REF] an end-to-end ego lane estimation is presented using a deep learning network SegNet. For all the mentioned methods before, no road prior is exploited. However, a learning stage is needed in all cases, which may require a large database to perform it.

Localization within a lane

In order to guarantee the ego-lane level localization, several data fusion methods have been investigated. The first one uses a GPS receiver with a digital map to locate the egovehicle in the lane. The lack of accuracy provided by a classic GPS that can be caused by poor satellite signals, high degree dilution of precision, or multi-path in urban scenes, is first compensated with proprioceptive sensors, such as Inertial Measurement Unit (IMU). These methods are well known as Dead Reckoning [START_REF] Milanés | Intelligent automatic overtaking system using vision for vehicle detection[END_REF]. Nevertheless, the accuracy provided by this kind of method is not sufficient for autonomous navigation, where centimetric accuracy is needed. Furthermore, dead reckoning methods suffer from integration errors that tend unbounded in time.

Leafing through the literature, it appears that most researchers use lane marking detection from a camera to provide an accurate ego-lane level localization. There are two main advantages of using the camera: currently, this sensor is the cheapest and most versatile modality for automotive applications and it provides dense information of the environment. Common approaches are based on a twostep process [START_REF] Hillel | Recent progress in road and lane detection: a survey[END_REF]. First, a road marking feature detection [START_REF] Sawano | A road extraction method by an active contour model with inertia and differential features[END_REF], then a fitting procedure [START_REF] Wang | Lane detection and tracking using b-snake[END_REF] is used. This fitting procedure can be done with a polynomial fitting model based on the road geometry design [START_REF] Hernández | Real-time lane region detection using a combination of geometrical and image features[END_REF] or non-parametric model [START_REF] Nefian | Detection of drivable corridors for off-road autonomous navigation[END_REF].

The above-mentioned approaches suffer from the camera's weaknesses: high luminosity variation and occlusion issues.

In an attempt to provide a more robust and more accurate localization, the trend for manufacturers and researchers is to equip the ego-vehicle with multiple sensors such as Velodyne lidars, multiples cameras, IMUs and high defined digital map that contains the precise location of high features such as the lane marking or landmarks [START_REF] Schreiber | Laneloc: Lane marking based localization using highly accurate maps[END_REF]. In that way, authors in [START_REF] Gruyer | Accurate lateral positioning from map data and road marking detection[END_REF] used two lateral cameras combined with an IMU and a GPS receiver to accurately determine the lateral distance of the ego-vehicle with respect to the lane marking. To do so, a vision-based ego-lane marking detector was used [START_REF] Revilloud | A new multi-agent approach for lane detection and tracking[END_REF] and fused with the IMU sensors through an Extended Kalman Filter (EKF) to locate the ego-vehicle. Furthermore, a sub-decimeter accuracy was achieved by adding lanes marking from a created map. Yet again, an EKF was introduced for the fusion stage between the lane marking extracted from the map and the marking obtained through the vision process. In order to bound the integration errors coming from dead reckoning, authors in [START_REF] Hage | Student's t information filter with adaptive degree of freedom for multi-sensor fusion[END_REF] present a localization system based on the fusion of multi-sensor and a high definition map. The integrity of the localization is insured by a Student's distribution that replaces the classical Gaussian distribution used in fusion filter like Kalman filters in an Informational Filter (IF). Close to our work, an atypical approach presented in [START_REF] Alvarez | Combining priors, appearance, and context for road detection[END_REF] relies on the use of road priors and contextual information for road detection (roadway detection). OSM map is used to create the road backbone depending on the lanes number and the width of the lane. This road shape is then projected on the image taking into account the uncertainties related to the ego-vehicle pose. The result is used as prior to roadway detection. The road shape is then delimited by estimating the vanishing point and the lane marking. Naturally, the mentioned approaches are just a small part of the ego-lane localization techniques, an interesting survey is presented in [START_REF] Hillel | Recent progress in road and lane detection: a survey[END_REF] for more detailed approaches.

To the best of the authors' knowledge, there is no end-toend localization solution proposed that tackles every aspect of the ego-localization task. Thus, in this paper, we present a complete solution for the overall localization problem. Unlike existent work, this contribution distinguishes by the information-driven approach. Indeed, we exploit the available data from the OSM map in order to extract prior about the road geometry. This strategy allows us to focalize the region of interest in the recognition procedure and ultimately reduce the computation complexity. In addition, we propose a probabilistic model that will ensure the geometric coherency between lanes marking. Finally, the main contributions can be summarized as follows:

1) An enhanced map matching method based on OSM datasets and a Hidden Markov Model 2) An ego-lane detection based on a prior OSM map and a vision based recursive recognition ego-lane marking in an information-driven fashion. 3) Once the ego lane marking is detected and the lateral distance between the ego-vehicle and the lane marking estimated, we proceed to the ego-lane determination (the postion of the ego-lane in the road) using a combination of a Bayesian Network and a Hidden Markov Model. We perform the complete localization steps in a sequential fashion. Thereby, the robustness of the proposed approach is highly improved. In addition to that, the split between the different parts of the localization allows the use of any part of the complete solution individually.

OVERALL LOCALIZATION ALGORITHM

In this work, we present an end-to-end algorithm for ego localization as highlighted on Figure 1. The algorithm is intentionally divided into three modules: Road-level localization (RLL), Ego-lane level localization (ELL) and Lane level localization (LLL). This split allows us to be in phase with the current literature. Moreover, the pipeline of the algorithm allows to have modular algorithms that can be changed in the future without changing the entire algorithm architecture. In the following, we will first introduce our map-matching module based on OSM datasets, then we will use this map as a prior for our ego-lane marking detection in order to locate the vehicle within its lane. Once the egolane level localization is performed, we will proceed to the identification of the ego-lane in the road.

Road-level localization (RLL)

In this section, we present our road-level localization algorithm using OSM datasets. As shown in Figure 2 the proposed module is an upgrade of our work presented in [START_REF] Kasmi | Map matching and lanes number estimation with openstreetmap[END_REF]. Indeed, an HMM is added to robustify the proposed Map-Matching. Therefore, the algorithm is divided into two different stages:

• Discrimination stage elimination of incongruent links based on distance, orientation difference between vehicle steering, and links heading. Additionally, a third factor was implemented, it is based on the maximum speed allowed on the road.

• Selection stage in the case where the remaining candidate links are greater than one, we proceed into a selection procedure based on multiple probabilistic criteria to eliminate the ambiguity. Fig. 2: Overall algorithm for the RLL using GPS data and OSM in order to extract the correct road and the corresponding number of lanes.

Discrimination stage

Before starting with this stage, we must present the OSM data. OSM data are composed of three keys components: Nodes, Ways, and Relations [START_REF] Ramm | OpenStreetMap: using and enhancing the free map of the world[END_REF]. Nodes are the geometrical elements that represent GPS points. Ways are an ordered list of nodes that represents roads network, every Way has a number of tags, the latter provide information about the characteristics of the Way (number of lanes, limitation speed...). In addition, each Way is composed of a set of segments. In other words, belonging to a segment is equivalent to belonging to an OSM Way. So the map matching task can be reformulated as matching a GPS point with a segment. In the remaining sections, the words segment and Way are switchable. Thus, the discrimination stage consists of the exclusion of segments which are not compatible with the current vehicle and the road network topology configuration. To do so, we introduced three discrimination criteria [START_REF] Kasmi | Map matching and lanes number estimation with openstreetmap[END_REF]:

•
The first one is based on the distance between the ego-vehicle position estimated with the GPS data and the segment.

•

The second one is based on the angle difference between the ego-vehicle steering angle and the segment heading (which expresses the traffic flow).

•

The third one is based on the speed limit. The assumption made is that the ego-vehicle is respecting the limitation of speed to some extent i.e. the ego speed can not be higher than the limitation plus

40 km h -1 .
Once the discrimination stage completed, if the number of candidates segment is higher than one, we proceed to the selection of the correct Way.

Selection of the correct Way

We define the state of ego-vehicle at time step t as x t = (x w , y w , θ w ) T . The selection procedure can be formulated as finding the highest conditional probability of belonging to a Way W i knowing the pose of the ego-vehicle x t at each time t:

arg max p(W it |x t ) (1)
The main goal of this formulation is to select the correct Way on which the vehicle is traveling. Consequently, a selection criterion must be introduced to determine the correct Way.

Based on equation (1) three probabilistic criteria [START_REF] Kasmi | Map matching and lanes number estimation with openstreetmap[END_REF] were developed:

• C e : Criterion based on Euclidean distance.

• C m : Criterion based on Mahalanobis distance.

• C p : Criterion based on the probability of belonging to a segment.

Depending on the data available, one of these criteria is calculated for each GPS data and for each candidate Way. Furthermore, to take into account the relation between two consecutive frames, a convenient and effective approach is to use a HMM. Indeed, the change in Ways over time is governed by topology constraints that can be embedded in the transition state of the HMM. Thus we enhance the proposed multi-criteria algorithm with an HMM. The HMM in this work is an upgrade of the proposed HMM by [38] [39]. We fitted the HMM to the way selection problem. However, the big difference lies in the probabilistic reasoning. Indeed, we use criteria that represent probability to model the observation probability, which is not the case in [START_REF] Kang | An enhanced map-matching algorithm for real-time position accuracy improvement with a low-cost gps receiver[END_REF] [START_REF] Luo | Enhanced map-matching algorithm with a hidden markov model for mobile phone positioning[END_REF], where arbitrary functions have been introduced leading to non-intuitive criteria.

By definition a HMM consists of two stochastic processes, the first one is a Markov Chain to model the change of a state vector over time. This change is governed by a probability that describes the transit probability over time, which is called the transition probability. In the HMM, the states of the chain are not visible but observable, for this reason they are called "hidden". The second process is called the observation space, it produces emission of the observation at each time. Although the elements of the state vector are hidden, there is a relation between the hidden elements of the state and the observations, this relation is referred as an emission probability. Figure 3

P e (W i ) = C k (W i ) j ∈ Si C k (W j ) (2)
where C k is the probabilistic criterion used with k ∈ {e, p, m}.

Transition probability

The transition probability reflects the probability of the transition matrix that a state will move from one state to another. In the map matching procedure, the topology of the network is used in order to determine the transition matrix. Indeed, the vehicle can only move on Ways that are physically connected. We will refer to P(W t i , W τ j ) as the transition probability from Way W t i to Way W τ j given the state space S t and S τ for time t and τ :

P(W t i , W τ j ) = w i j ∈ St w j (3)
Where w i is a criterion calculated for each Way, this criterion obeys some rules defined as follows:

1) w i = 0 if the way W t i and W τ j are not connected; 2) w i = s tτ ij if the Way W t i and W τ j are the same; 3) w i = c tτ ij if the Way W t i and W τ j are connected. Figure 4 shows a simplified road network to illustrate the transition matrix probability. These probabilities are shown on Table 1. It can be noticed that in this example the number of Ways is the same for two connected states, which is not always the case. Indeed, the number of Ways may vary from two consecutive frames. The transition probability s tτ ij is used whenever the two Ways are the same. Let P (x w , y w ) be the vector describing the Cartesian coordinates of the ego-vehicle in the Universal Transverse Mercator coordinate system frame (UTM), we note P the projection of the point P on the Way and ψ the Way's heading, the Way W τ i is composed of two nodes n 1 , n 2 as represented on Figure 5. The criterion s ij is computed as a Gaussian distance to the middle of the segment [n 1 , n 2 ]. Indeed, the probability to change the Way will be higher if P is in the middle and will be lower in the limits, thus it is calculated as follows:

s ij = e - 1 2   t -u 0 σ 0   2 (4)
With P = t.n 1 n 2 , u 0 = 0 and we take σ 0 = 0.25.

If the Way W i and W j are connected then we introduce a criterion based on [START_REF] Luo | Enhanced map-matching algorithm with a hidden markov model for mobile phone positioning[END_REF], this criterion depends on the Fig. 4: Example of road network to illustrate the transition probability calculation. In this example the number of Way candidates is the same between two consecutive frames

with ∆θ being the vehicle's heading change over time, ∆ψ the Way's heading change between the two Ways W i and W j and β a chosen coefficient. Once the Way selected, the number of lanes is extracted from the OSM database.

Ego-lane level localization (ELL)

Once the RLL has been performed, the ELL is initiated. We choose to start the ELL before the LLL for several reasons: the first is that the ego-lane marking is, for the most part, the easiest one to detect. The second one is due to the use of the information driven approach. Indeed, when the ELL is completed we have an estimate of the ego-vehicle's lateral position in its own lane. So, knowing the lane's width and the lanes number allows us to interpolate research zones for other lane marking.

To perform the ELL, the road is modeled in a 2D vehicle's frame (x v , y v ) as a cubic polynomial [START_REF] Abramov | A flexible modeling approach for robust multi-lane road estimation[END_REF] [41]:

x v = 1 6 c 0 y 3 v + 1 2 c 1 y 2 v + c 2 y v + c 3 (6)
The parameters c 0 is the curvature's derivative of the road, c 1 the curvature of the road, c 2 the vehicle's heading with respect to the tangent of the road, c 3 the lateral shift of the ego-vehicle with regards to the road model. Once the road has been modeled in the vehicle's frame, the next step is to project this model on the image taking into account the intrinsic and extrinsic parameters of the camera. According to [START_REF] Aufrere | A model-driven approach for real-time road recognition[END_REF] the projection of the road model (left and right lines) in the image frame (u i , v i ) is defined as follows:

u i = e u e v Z 0 (v i -e v α) 2 c 0 6 - e v Z 0 2(v i -e v α) c 1 (7)
+ v i -e v α e v Z 0 (x 0 ± L w 2 ) -x 0
where e u = f /du, e v = f /dv, f is the focal distance of the camera, du and dv are the width and height of a pixel in the image, Z 0 is the height of the camera, x 0 the lateral distance of the ego-vehicle with respect to the ego-lane marking, the ± sign indicates whether the ego-marking is right (+) or left (-), α is the camera tilt angle and L w is the road width. i = 1, ..., n Roi with n Roi the number of Region Of Interest (ROI) for each lane marking. As mentioned in [START_REF] Kasmi | Map matching and lanes number estimation with openstreetmap[END_REF], the OSM map does not provide information about the accuracy of its data. However, it is well known that OSM is a collaborative project in which volunteers provide the geospatial data. If we consider that most of the volunteers have a classic GNSS receiver, the accuracy of the OSM is thus metric. That being said, the polynomial representation of the road will be affected by this unknown accuracy. For the parameters c 0 and c 1 , they represent the shape of the road. Given that we are working mostly with highway roads, their values will not be affected by the inaccurate geospatial accuracy of the OSM. In contrast, the parameter c 2 is more sensitive to inaccurate accuracy of the OSM, and its value will be affected. Hence, to compensate for the inaccuracy of its value, we will define X v = [x 0 , c 2 ] T as the state vector of the road model in the ego-vehicle frame. We associate a covariance matrix to this vector C Xv . This matrix expresses the allowed dispersion around the average parameters previously defined by the vector X v . Concerning the parameters c 3 , the value extracted from equation ( 6) is not used. Indeed, it express the lateral shift of the vehicle to the road. Hence it is not used in equation ( 7), as we use x 0 , which is the lateral shift regarding the center of the ego-lane. It assumed that the ego-vehicle travel in the center of the ego-lane thus x 0 is around 0.

Using equation ( 7), we can express the probabilistic model (X v , C Xv ) into the image space. We will note this proba- bilistic model as u , C u , with u being the average values for the pixel in the image and C u its corresponding covariance matrix:

C u = J u C Xv J T u (8)
with:

J u =        ∂u 1 ∂x 0 ∂u 1 ∂c 2 . . . . ∂u n ∂x 0 ∂u n ∂c 2        , C Xv = σ 2 x0 0 0 σ 2 c2 (9)
The resulting projection of the equation ( 7) is shown in Figure 6. Taking into account the prior about the road geometry allows focalizing the zones of research in a Top-Down process fashion. Indeed, not all the image is used in order to perform the processing task. Consequently, the recognition process is faster, and subject to less noise considering it takes into account only the regions in the image that most likely contains a lane marking. Once the probabilistic model is defined, we proceed to the recognition phase. This stage is inspired by [START_REF] Aufrere | A model-driven approach for real-time road recognition[END_REF]. For clarity's sake, we will discuss each aspect of this stage as presented in Figure 7.

i Initialization in this stage we define the probabilistic model (u , C u ) as already presented in Figure 6. ii Selection of the best Region Of Interest (ROI) once the probabilistic model is defined, we proceed to the selection of the most informative ROI in the Image in a Top-Down fashion To do so, we present an informational criterion based on the Shannon entropy [START_REF] Delobel | A real-time map refinement method using a multi-sensor localization framework[END_REF]. For each ROI, we use an a priori selection based on a entropic criterion H sel :

H sel H - v -H * v ( 10
) = p(Dk = 1, Db = 1). 1 2 log|2π e C - Xv | -log|2π e C * Xv |
with: For the probabilities, p(D b ), p(D k ), a Bayesian Network (BN) is introduced. The presented nodes are an adaptation of the BN presented in [START_REF] Aynaud | Real-time vehicle localization by using a top-down process[END_REF]. The scheme of the network is illustrated in Figure 8. The BN is composed of the following nodes:

       X - v T he
• X k-the confidence before the detection is attempted.

• Z 0 is the chosen landmark (the white strip) is observable in the ROI. This network has two uses, the first one is to calculate the a priori probabilities p(D k ) and p(D b ), the second one is to determine the confidence p(X k+ ) obtained after proceeding to a detection. iii Detection Once the most informative ROI is chosen, we proceed to the detection. The implemented method is based on a row filter as presented in [START_REF] Aufrere | A model-driven approach for real-time road recognition[END_REF]. The idea is to compute the gradient of the image row by row with the aim to find pixels that correspond to the lane marking.

Once the patterns have been selected for the complete ROI, a Ransac method is used in order to detect the segment in the ROI. Thus, for each ROI two points p 1 (u u , v u ) and p 2 (u d , v d ) are defined. iv Update for the state vector x p = (u u , u d ) we associate a covariance matrix C p :

C p = σ 2 uu 0 0 σ 2 u u ( 11
)
With σ 2 uu being the accuracy of the segment detected. In order to update the probabilistic model, we use a Kalman filter:

u + = u -+ K u [x p -x p ] C + u = C - u -K u H u C p (12)
K u being the Kalman Gain, (u -, C - u ) is the model before detection and (u + , C + u ) is the model after detection. v End after each detection, we compute the corresponding entropic gain H gain :

H gain = p X k+ . 1 2 log|2πeC + u | -log|2πeC u | (13)
With C u the covariance from the initialization stage. The ego-lane marking has been detected if the following condition is satisfied :

H gain λH max (14)
with

H max = p X k+ . 1 2 log|2π e C + u | -log|2π e C u |
With λ a fixed coefficient (< 1). If the condition is satisfied, it means that enough detections have been attempted successfully, and the obtained probabilistic model (u + , C + u ) is sufficiently precise to consider the ego-lane recognition finished. Otherwise, the ego-lane recognition process carried out, and the second most informative ROI is selected. The end of the marking detection involves the end of the ego-lane level localization. Indeed, the update of the model u, C u leads to the update of the model X v , C Xv . This means that an estimation of the road parameters c 2 , x 0 is performed.

Lane level localization (LLL)

Once the ego-lane localization is estimated, we have to perform the lane-level localization to correctly choose the right lane on which the vehicle travels. To do so, we proposed in [START_REF] Kasmi | Probabilistic framework for ego-lane determination[END_REF] a probabilistic framework that is split into three stages. Thus, the presented LLL algorithm is an extension of our work proposed in [START_REF] Kasmi | Probabilistic framework for ego-lane determination[END_REF]. Indeed, the Hidden Markov Model's architecture is the same. The difference lies in the input of the algorithm, as we take the information about the road from OSM. In addition, the parameters of the road, such as the curvature, are also extracted from OSM, which will be useful for interpolating other lanes marking. In the first stage, we extrapolate adjacent lanes by assuming that lanes in the same road have the same width L w and the variation in the curvature c 1 is very small. The second step consists of a Bayesian Network (BN), that takes as input the results of the hypothesized adjacent lane-marking detection, whether these detections succeed or fail. Furthermore, since the proposed BN is modular, we also use an adjacent vehicle detector based on deep learning. The third step includes a filtering process, using a Hidden Markov Model (HMM).

Adjacent lanes extrapolation

The adjacent lanes are extrapolated by taking advantage of the estimated ego-lane localization. Thus, each adjacent lane is described by a probabilistic model {X l , C X l }, where X l contains the parameters (x 0 , c 2 ) T described in Section 3.2 and C X l refers to the corresponding covariance matrix. The assumption made is that the curvature and lane width stays constant for all the lanes in the same road. Thereby, the value of the vector X l remains the same as x v described in equation ( 6), only the value of x 0 will be shifted by a (±iL w ), with (-i) indicates that the edge is at the right of the ego-lane and (+i) the left. The number of adjacent lanes extrapolated is equal to the lanes number from OSM. Thus, from a perspective view, the hypothesized lanes are shown in Figure 9. As mentioned in the ELL section, the model {X l , C X l } can be transferred to the model image {u l , C u l }. As a consequence, u l represents the horizontal pixel of the edges in the image and C u l its interval confidence. In Figure 10, the adjacent lane regions of interest resulting from the extrapolation are shown.

For each adjacent lane marking extrapolated, the detection is performed and the results are fed into a Bayesian Network (BN).

Bayesian network for ego-lane determination

The proposed BN is designed to be flexible and modular for other detection results from any type of sensor, i.e. vehicle detector, guardrail detector. In order to show the flexibility of the BN, we will first use only adjacent lanes detection to determine the ego-lane. Thereafter, an adjacent vehicle detector based on Deep Learning (YOLO [START_REF] Redmon | You only look once: Unified, real-time object detection[END_REF]) is introduced. The general architecture of the BN used is illustrated in Figure 11. The described nodes are the following:

• Z ki the element i is observable, this element can represent any element of the road scene: {vehicles, lane-marking, traffic signs...},

• D ki the detection of the element of type k is successful. Later we will use D l for adjacent lanes and D v for vehicles,

• L BN the lane on which the ego-vehicle is traveling, L BN is a form of categorical distribution, L BN ∼ Cat(ξ) with L BN = {l 1 , ..., l n }, where l 1 indicates the leftmost lane and n the lanes number.

Depending on the results of the detection, we infer the probability to belong to a lane l i as follows:

P (L BN =li ) = P (L BN =li |D k1 , Z k1 , ..., D kn , Z kn ) (15)
The lane with the highest probability is chosen to be the lane on which the vehicle travels. As mentioned before, the proposed BN is aimed to be modular for other detectors.

To start, we only use the Adjacent Lanes Detection (ALD) as input into our BN (BN + ALD), where D li indicates whether the detection of adjacent lane i is successful. After, we add the information about the adjacent vehicles D vi using the Yolo detector, where D vi indicates whether the detection of the adjacent vehicle i is successful. Knowing that, this detection is performed in the regions of the image bounded by the neighboring marking-lanes.

HMM for ego-lane determination

The proposed BN in the previous section is applied on a per-frame basis. However, the dynamic relationship between two consecutive frames is not taken into account. Indeed, the BN does not take into consideration the dynamic constraints of the ego-vehicle (i.e. the ego-vehicle can only change lane to an adjacent lane). In order to take into account these constraints, we filter the output of the BN by a Dynamic Bayesian Network which is the Hidden Markov Model (HMM). We will use L t to denote the set of ego-lane state variables at time t, which depends on the lanes number n lanes and Fig. 9: Detected ego-lane in solid red and the hypothesized lanes in dash which are assumed to be observable. e t denotes the observable evidence variable. The aim of the filter algorithm is to estimate the probability P (L t+1 |e 1:t+1 ). According to [START_REF] Russell | Artificial intelligence-a modern approach[END_REF], this probability can be formulated as follows: Where η denotes a normalizing constant to make probabilities sum equal to 1. By arranging Equation ( 16): The probability P (e t+1 |L t+1 ) comes from the observation model. Hence, in this paper from the BN described previously, thus:

Z k 1 D k 1 Z k 2 D k 2 Z k 3 D k 3 . . . . . . Z k n D k n L BN
P (L t+1
P (L t+1
P (e t+1 |L t+1 ) = P (L BN ) (18)
With regard to the probability P (L t+1 |lt), it comes from the transition model. It expresses the probability of the egolane to change its current state, which is the lane change probability. Finally, the third term P (Lt, e 1:t ) expresses the current state distribution. Graphically, we can illustrate the corresponding HMM as in Figure 12.

With the recursive formulation obtained in Equation ( 17), we can estimate the current ego-lane state given the observation obtained from the BN, but before we have to compute the lane change probability.

Transition probability (Lane change probability)

To calculate the lane change probability, we model the lateral position x 0 estimated in the ego-lane level localization as a normal distribution with mean µx 0 and variance σx 0 :

x 0 ∼ N (µx 0 , σ 2 x 0 ) (19)
The value µx 0 and variance σx 0 are obtained from the estimated probabilistic model (xv, Cx v ). Therefore, we predict the lateral position x 0 at time t k+1 as shown on Fig 13.

Accordingly, the lane-change probability is calculated as Fig. 13: Lane change probability, with the blue area representing the right change probability. The coordinate system for x 0 is centered on the vehicle follows:

P (Cr) = +∞ Lw /2 1 σx 0 √ 2π e - (x -µx 0 ) 2 2σ 2 x 0 dx (20)
P (C l ) = -Lw /2 -∞ 1 σx 0 √ 2π e - (x -µx 0 ) 2 2σ 2
x 0 dx [START_REF] Lee | Accurate ego-lane recognition utilizing multiple road characteristics in a bayesian network framework[END_REF] with P (Cr) the probability to right change lane and P (C l ) the probability to left change lane. Now that we have designed the HMM, we will introduce the real-world experimental results in the following section.

REAL-WORLD EXPERIMENTAL RESULTS

The different parts of the localization algorithm have been introduced in the previous section. In order to prove the effectiveness of the presented algorithms, real-world experiments have been carried out.

To do so, we tested our algorithm on real driving datasets. The first one was collected in the region of Clermont-Ferrand in France, where we drove our acquisition vehicle on two lanes and three lanes road in a national highway. In to collect these datasets, the acquisition vehicle was equipped with a front camera, with the addition of an IMU and a classic GPS receiver. The data have been stored on a bag using using ROS. In addition to that, the ego-vehicle provides information about its own odometery. Finally, these datasets will be made available for researchers 1 . In the following, they will be referred to as "2-lanes" for two lanes road and "3-lanes" for three lanes road (1000 images frames). Furthermore, we wanted to compare the performances of our algorithm to the literature on more challenging scenarios. Naturally, we turned our attention to the KITTI datasets [START_REF] Geiger | Vision meets robotics: The kitti dataset[END_REF]. Although these dataset contain few lanechanging scenarios, we will however use some part of those in order to show the effectiveness of our RLL module. In addition, we tested our algorithm on some datasets referred to in [START_REF] Ballardini | Ego-lane estimation by modeling lanes and sensor failures[END_REF]. Unfortunately, all the data were not ready yet. We were able to test on some of them, noted as the A4-Highway Italy, for a total of 9528 frames. The collected datasets were manually annotated in a per-frame basis in order to determine the correct ego-lane classification. In the following, we will refer to these datasets as "4-lanes". Each part of the ego-localization: RLL, ELL and LLL will be tested on the presented datasets, the results will be discussed in terms of accuracy. In the end, the overall algorithm will be discussed in terms of computation time. In order to assess the presented RLL algorithm, the correct Way has been manually annotated for each GPS frame. The resulting ground truth was compared with HMM-MM. The results are summarized on Table 2. As highlighted, the Map-Matching with the addition of the HMM improves the overall accuracy precision. This can be explained by the inclusion of the topology of the road network in the HMM. Indeed, the transition states are governed by the connectivity between links. Furthermore, we test the reliability of our algorithm in more challenging scenarios. Thus, we used the well-known KITTI [START_REF] Geiger | Vision meets robotics: The kitti dataset[END_REF] database in urban scenarios. The used sequences are the following: "sequence_2001_09_26_0001" and "se-quence_2001_09_26_0001". Even if the two sequences were taking in urban areas, the correct map-matching obtained shows the effectiveness and the robustness of our proposed RLL algorithm in varying scenarios. Finally, it can be noticed that for the "4-lanes", the GPS datasets have not been mentioned. Indeed, practically the entirety of the datasets is on the same Highway. Thus, the results on MM would not be relevant. As discussed in [START_REF] Revilloud | A lane marker estimation method for improving lane detection[END_REF] [13], in most of the work presented in the literature, the aim is to find the white strips in the image without seeking for consecutive elements of lane marking. Cnversely, in our work, we are more interested in the curvature of the lane marking than the white strips. However, by doing so, the characterization of the algorithm is considerably more difficult. Some examples are presented in Figure 14 showing the robustness of the presented egolane marking detector on different lighting conditions. Once the ego-lane localization completed, the Lane level localization is performed. So as to point out the increment of each added part, we first, determined the ego-lane using solely the BN with adjacent lanes. Within the second instance, we introduced adjacent vehicle detection in the BN. In all instances, we filtered the outcome of the BN with the HMM. All the results obtained are summarized in Table 3.

Considering the results, the increment provided by each module is clearly illustrated. Indeed, altogether cases the Fig. 14: Some examples of correct ego-lane marking recognition (in Blue). As highlighted, the ego-lane marking are correctly detected in varying imaging conditions: low light level for the top images and high brightness for the bottom images. BN+ALD+AVD provides more accurate classification than the BN+ALD, which suggests that addition of another information source also will improve the obtained accuracy. After investigation on the inaccurate results, it appears that the false classifications obtained using the BN+ALD are due to two main reasons: either the lanes marking are not detected or the lanes marking are wrongly detected. For the primary case, this can be explained if the lanes marking are missing or hidden by an object. For the second case, it shows the limitation of the used lane marking detector. Moreover, even if the introduction of the adjacent vehicle detection shows excellent results, there are some cases where the vehicles detection are not relevant, for instance, if the detected vehicle is on the nearby road. To overcome this issue, we would have to determine the localization of the detected vehicle relative to the ego-vehicle, which is not the case since we use solely images as input. Finally, despite that the authors 2 in [START_REF] Ballardini | Ego-lane estimation by modeling lanes and sensor failures[END_REF] did not take into account the lane change scenarios and designed empirically the HMM. Nevertheless, we manage to outperform their results on the same datasets. Indeed, they achieved 77% correct classifications, where we were able to reach 85.35% on 9528 frames. Some ELL in different images are presented in Figure 15. The presented results have been carried using Python 3.7 under a Dell G3 3579 Core i7 8th generation equipped with an NVIDIA GeForce GTX 1050 Ti. The computation time results presented in Figure 16 shows that even if the entire algorithm was coded in Python, real-time implementation is possible. Indeed, the sum of all parts of the ego-localization algorithm is under 500 ms (446.95 ms on average). In addition, if we glance in-depth at the time consumed in the RLL, we found that on average, it takes 238.16 ms to query the local server containing the OSM data. Furthermore, from the 143.67 ms dedicated to the LLL, 65.60 ms are spent on the YOLO detector. These results lead us to assume that implementation on C/C++ will divide the calculation time by 10.

CONCLUSION

In this paper, we presented an end-to-end ego-vehicle localization. Starting from Road-level localization (RLL) using OSM datasets and a classic GPS receiver, to the Ego-lane level localization (ELL) using a recognition lane marking in an information-driven fashion and finally, a Lane-level localization (LLL) using the well-known YOLO detector in a probabilistic framework composed of a Bayesian Network and Hidden Markov Model to correctly determine the egolane in the road. The distinction from other work lies in the overall aspect of the ego-localization that has been tackled

2. The authors would like to acknowledge the authors of [START_REF] Ballardini | Ego-lane estimation by modeling lanes and sensor failures[END_REF] for their help with their datasets. in stock in the paper. In addition, this work is distinguished by the information-driven approach. For our future work, we are currently working on adding detectors from different sensors, i.e., lidar, radar to enhance the proposed framework. On the other hand, we are still working on the OSM datasets in order to update this database when it provides false information about the environment.

Fig. 1 :

 1 Fig.1: Overall algorithm for complete localization of the ego-vehicle starting from the Road-level localization (RLL): to locate the ego-vehicle on the map Using OSM and GPS datasets, to the Lane-level localization (LLL): to locate the egovehicle within its own lane using prior road model from OSM and camera images and the Ego-lane level localization (ELL): to determine the ego-lane position in the road.

Fig. 3 :

 3 Fig. 3: Modeling of the Hidden Markov model Map-Matching (HMM-MM) algorithm, with S t the state space at time t, W t the Way candidate at time t and O t the observations at time t. Note that the number of Way candidates at each time may vary.

Fig. 5 :

 5 Fig. 5: Example of calculation for transition probability s tτ ij , the vehicle is illustrated by the box with black triangle, ψ represents the Way heading and θ the vehicle's heading

Fig. 6 :Fig. 7 :

 67 Fig. 6: Projection of the initial probabilistic model (u , C u ) on the image. The mean values of the lane marking u are presented in Red, Blue illustrates the dispersion around these values. Finally, green boxes are the ROI for each lane marking, in this example the number of ROI is 9 for each lane marking.

  state vector bef ore simulation of the detection, C - Xv Covariance bef ore simulation of the detection, X * v T he state vector af ter simulation of the detection, C * Xv Covariance af ter simulation of the detection.

•

  D k a landmark is detected in the focal zone,• D b the correct landmark has been detected (which manages landmark ambiguity).

Fig. 8 :

 8 Fig. 8: Bayesian Network used for the confidence estimation. Yellow nodes are the input nodes, Purple node is the node we are seeking to estimate.

Fig. 10 :

 10 Fig. 10: Results of the adjacent lane regions of interest resulting from the extrapolation

Fig. 11 :

 11 Fig. 11: Final architecture Bayesian Network (BN) used for ego-lane determination.

  |e 1:t+1 ) = P (L t+1 |e 1:t , e t+1 )(dividing up the evidence) = ηP (e t+1 |L t+1 , e 1:t )P (L t+1 |e t+1 )(Baye's rule) = ηP (e t+1 |L t+1 )P (L t+1 |e t+1 )(Markov assumption)[START_REF] Nieto | Robust multiple lane road modeling based on perspective analysis[END_REF]

  |e 1:t+1 ) = ηP (e t+1 |L t+1 ) l t P (L t+1 |Lt, e 1:t )P (Lt, e 1:t ) = ηP (e t+1 |L t+1 ) l t P (L t+1 |Lt)P (Lt, e 1:t ) (17)

Fig. 12 :

 12 Fig. 12: Hidden Markov Model for n lanes case, with L = {l 1 , l , ..., l n being the set of hidden states and O = {e 1 , e 2 , ..., e n } being the set of observations resulting from the BN. we can estimate the current ego-lane state given the observation obtained from the BN, but before we have to compute the lane change probability.

Fig. 15 :

 15 Fig. 15: Lane Level Localization on 3 lanes and 4 lanes scenarios. The orange boxes are the output of the Yolo detector.Cyan is used for adjacent lanes markings. Red is used for the ego-lane marking. A video is available at shorturl.at/lAR28.

Fig. 16 :

 16 Fig. 16: Computation time in ms for each part of the presented algorithm: RLL in Light-blue, ELL in pink and LLL in light-green, with mean values of 244.53ms for the RLL, 58.74ms for the ELL and 143.68ms for the LLL. Since the computation time for the ELL and LLL depends on the observability conditions of the markings in the image, the dispersion in time is much bigger than the RLL.

  State space the state of the system describes the list of the Ways candidates for each observation x t . We will use S t to denote the set of candidates Ways at time t S t = {W 1 , W 2 , ..., W nt } with n t the number of candidates Ways at each time t, S t is a form of categorical distribution, S t ∼ Cat(ξ). Observation space for each candidate Way composing the state space, an emission probability is made. This emission probability is directly calculated from the multi-criteria algorithm. For each Way candidate the emission probability P

	describes the proposed HMM
	for the map-matching task.
	So, in order to model the HMM-MM, three components
	must be defined:

e (W i ) is calculated as follows:

TABLE 1 :

 1 Transition matrix for example shown on Figure4

		W 1	W 2	W 3	W 4	W 5
	W 1	s 11	0	c 13	0	c 15
	W 2	0	s 22	0	0	0
	W 3	c 31	0	s 33	c 34	c 35
	W 4	0	0	c 43	s 44	0
	W 5	c 51	0	c 53	0	s 55

TABLE 2 :

 2 Map Matching results on the entire datasets without and with the HMM, KITTI-1 and KITTI-2 refer to KITTI sequences:'_2011_09_26_0001' and '_2011_09_26_0002'

	1. shorturl.at/lAR28

TABLE 3 :

 3 Classification accuracy for ego-lane determination. BN+ALD refers to the BN feed with the adjacent lanes detection, BN+ALD+HMM indicates the HMM with the corresponding BN, BN+ALD+AVD refers to the BN feed with adjacent lanes and vehicles detection and BN+ALD+AVD+HMM refers to the HMM with the corresponding BN.

ACKNOWLEDGMENT

This work has been sponsored by Sherpa Engineering and ANRT (Conventions Industrielles de Formation par la Recherche). This work has also been sponsored by the French government research program Investissements d'avenir through the RobotEx Equipment of Excellence (ANR-10-EQPX-44) and the IMobS3 Laboratory of Excellence (ANR-10-LABX-16-01), by the European Union through the program Regional competitiveness and employment 2014-2020 (FEDER -AURA region) and by the AURA region.

*This work has been sponsored by Sherpa Engineering and ANRT (Conventions


Olivier Gouvert

Thomas Oberlin

Cédric Févotte

Ordinal Non-negative Matrix Factorization for Recommendation

We introduce a new non-negative matrix factorization (NMF) method for ordinal data, called OrdNMF. Ordinal data are categorical data which exhibit a natural ordering between the categories. In particular, they can be found in recommender systems, either with explicit data (such as ratings) or implicit data (such as quantized play counts). OrdNMF is a probabilistic latent factor model that generalizes Bernoulli-Poisson factorization (Be-PoF) and Poisson factorization (PF) applied to binarized data. Contrary to these methods, Ord-NMF circumvents binarization and can exploit a more informative representation of the data. We design an efficient variational algorithm based on a suitable model augmentation and related to variational PF. In particular, our algorithm preserves the scalability of PF and can be applied to huge sparse datasets. We report recommendation experiments on explicit and implicit datasets, and show that OrdNMF outperforms BePoF and PF applied to binarized data.

Introduction

Collaborative filtering (CF) is a popular recommendation technique based only on the feedbacks of users on items. These feedbacks can be stored into a matrix Y of size U × I, where U and I are the number of users and items respectively. Matrix factorization (MF) methods [START_REF] Hu | Collaborative filtering for implicit feedback datasets[END_REF][START_REF] Koren | Matrix factorization techniques for recommender systems[END_REF][START_REF] Ma | Probabilistic factor models for web site recommendation[END_REF] aim to approximate the feedback matrix Y by a low-rank structure WH T where W ∈ R U ×K + corresponds to user preferences and H ∈ R I×K + to item attributes.

Poisson factorization (PF) [START_REF] Canny | GaP: A factor model for discrete data[END_REF][START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF][START_REF] Gopalan | Scalable recommendation with hierarchical Poisson factorization[END_REF] is a non-negative matrix factorization (NMF) model [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF][START_REF] Lee | Algorithms for non-negative matrix factorization[END_REF][START_REF] Févotte | Algorithms for nonnegative matrix factorization with the beta-divergence[END_REF] Proceedings of the 37 th International Conference on Machine Learning, Vienna, Austria, PMLR 119, 2020. Copyright 2020 by the author(s).

2011

) which aims to predict future interactions between users and items in order to make recommendations. For this purpose, PF is often applied to a binarized version of the data, i.e., Y ∈ {0, 1} U ×I , containing only the information that a user is interacting with an item or not. A variant of PF, called Bernoulli-Poisson factorization (Be-PoF) [START_REF] Acharya | Nonparametric Bayesian factor analysis for dynamic count matrices[END_REF], has been proposed to explicitly model binary data. However, for both PF and BePoF, the binarization stage induces a loss of information, since the value associated to an interaction is removed. Although several attempts in the literature tried to directly model raw data, both for explicit [START_REF] Hernandez-Lobato | Probabilistic matrix factorization with non-random missing data[END_REF] and implicit data [START_REF] Basbug | Hierarchical compound Poisson factorization[END_REF][START_REF] Zhou | Nonparametric Bayesian negative binomial factor analysis[END_REF][START_REF] Gouvert | Recommendation from raw data with adaptive compound Poisson factorization[END_REF], this remains a challenging problem.

In an attempt to keep as much information as possible, we propose in this paper to consider ordinal rather than binary data. Ordinal data [START_REF] Stevens | On the theory of scales of measurement[END_REF] are nominal/categorical data which exhibit a natural ordering (for example: cold ≺ warm ≺ hot). This type of data is encountered in recommender systems with explicit data such as ratings. It can also be created by quantizing implicit data such as play counts. Such a pre-processing remains softer than binarization and stays closer to the raw data, as soon as the number of classes is chosen big enough. In this paper, without loss of generality, we will work with ordinal data belonging to {0, . . . , V }. Note that for this type of data, the notion of distance between the different classes is not defined. For example, this implies that the mean is not adapted to these data, unlike the median.

There are two naive ways to process ordinal data. The first one consists in applying classification methods. The scale of the ordering relation which links the different categories is then ignored. The second way considers these data as real values in order to apply regression models. By doing this, it artificially creates a distance between the different categories. These two naive methods do not fully consider the specificity of ordinal data, since they remove or add information to the data. Threshold models [START_REF] Mccullagh | Regression models for ordinal data[END_REF][START_REF] Verwaeren | Learning partial ordinal class memberships with kernel-based proportional odds models[END_REF] are popular ordinal data processing methods that alleviate this issue. They assume that the data results from the quantization of continuous latent variables with respect to (w.r.t.) an increasing sequence of thresholds. The aim of these models is then to train a predictive model on the latent variables and to learn the sequence of thresh-olds. Threshold models can thus be seen as an extension of naive regression models, where the distances between the different classes are learned through quantization thresholds. For a comprehensive and more detailed review on ordinal regression methods, we refer to [START_REF] Gutierrez | Ordinal regression methods: Survey and experimental study[END_REF].

In this paper, we develop a new probabilistic NMF framework for ordinal data, called ordinal NMF (OrdNMF). Ord-NMF is a threshold model where the latent variables have an NMF structure. In other words, this amounts to defining the approximation Y ≈ G(WH T ), where Y is the ordinal data matrix, W and H are non-negative matrices and G(•) is a link function. OrdNMF allows us to work on more informative class of data than classical PF method by circumventing binarization. Contrary to ordinal MF (OrdMF) models [START_REF] Chu | Gaussian processes for ordinal regression[END_REF][START_REF] Koren | An ordinal model for predicting personalized item rating distributions[END_REF][START_REF] Paquet | A hierarchical model for ordinal matrix factorization[END_REF][START_REF] Hernandez-Lobato | Probabilistic matrix factorization with non-random missing data[END_REF], OrdNMF imposes non-negativity constraints on both W and H. This implies a more intuitive part-based representation of the data [START_REF] Lee | Learning the parts of objects by non-negative matrix factorization[END_REF], and were shown to improve results in recommendation [START_REF] Gopalan | Scalable recommendation with hierarchical Poisson factorization[END_REF]. OrdNMF can efficiently take advantage of the sparsity of Y, scaling with the number of non-zero observations. Thereby, it can be applied to huge sparse datasets such as those commonly encountered in recommender systems. As opposed to learning-to-rank models, the aim of OrdNMF is to model ordinal data, via a generative probabilistic model, in order to predict the class of future interactions. Learning-to-rank models do not seek to predict a class but to rank items relatively to each other. For example, Bayesian personalized ranking [START_REF] Rendle | Bpr: Bayesian personalized ranking from implicit feedback[END_REF] is based on binary pairwise comparisons of the users' preferences and not on the raw matrix Y. Although such models can also be used for recommendation, they are not generative.

The contributions of this paper are the following.

• We propose a new NMF model for ordinal data based on multiplicative noise. In particular, we study an instance of this model where the noise is assumed to be drawn from an inverse-gamma (IG) distribution. We show that this instance is an extension of BePoF [START_REF] Acharya | Nonparametric Bayesian factor analysis for dynamic count matrices[END_REF] and PF [START_REF] Gopalan | Scalable recommendation with hierarchical Poisson factorization[END_REF] applied to binarized data.

• We use a model augmentation trick to design an efficient variational algorithm, both for the update rules of the latent factors W and H, and for those of the thresholds b. In particular, this variational algorithm scales with the number of non-zero values in Y.

• We report the results of OrdNMF on recommendation tasks for two datasets (with explicit and implicit feedbacks). Moreover, posterior predictive checks (PPCs) demonstrate the excellent flexibility of OrdNMF and its ability to represent various kinds of datasets.

The rest of the paper is organized as follows. In Section 2, we present important related works on cumulative link models and on BePoF. In Section 3, we present our general OrdNMF model and detail a particular instance. In Section 4, we develop an efficient VI algorithm which scales with the number of non-zero values in the data. In Section 5, we test our algorithm on recommendation tasks for explicit and implicit datasets. Finally, in Section 6, we conclude and discuss the perspectives of this work.

Related Works

2.1. Cumulative Link Models (CLMs)

CLMs were one of the first threshold models proposed for ordinal regression [START_REF] Agresti | Categorical data analysis[END_REF]. These models have been adapted to deal with the MF problem, leading to OrdMF models. They amount to finding the approximation Y ≈ G(WH T ), where Y ∈ {0, . . . , V } U ×I is an ordinal data matrix, W ∈ R U ×K and H ∈ R I×K are latent factors, and G(•) is a parametrized link function described subsequently. OrdMF has been applied mainly to explicit data in order to predict users feedbacks [START_REF] Chu | Gaussian processes for ordinal regression[END_REF][START_REF] Paquet | A hierarchical model for ordinal matrix factorization[END_REF].

The idea behind threshold models is to introduce a continuous latent variable x ui ∈ R that is mapped to the ordinal data y ui . This is done by considering an increasing sequence of thresholds b

-1 = -∞ < b 0 < • • • < b V -1 < b V = +∞,
denoted by b, which fully characterize the following quantization function, illustrated in Figure 2, by:

G b : R → {0, . . . , V } x → v such as x ∈ [b v-1 , b v ). (1)
Therefore, ordinal data result from the quantization of the variable x ui by the step function G b , i.e., y ui = G b (x ui ).

The latent variable x ui corresponds to the variable λ ui = [WH T ] ui ∈ R perturbed by an additive noise ε ui , whose cumulative density function (c.d.f.) is denoted by F ε : R → [0, 1]. Thus, we obtain the following generative model, illustrated in Figure 1:

x ui = λ ui + ε ui , (2)
y ui = G b (x ui ).
(3)

10 0 10 2 10 4 10 6 x 0 2 4 6 8 10 v v = G b (x) Figure 2. Example of a quantization function x → G b (x).
The goal of MF models for ordinal data is therefore to jointly infer the latent variables W and H as well as the sequence of thresholds b.

Cumulative distribution function. The c.d.f. associated with the random variable y ui in Eqs. ( 2)-( 3) can be calculated as follows:

P[y ui ≤ v|λ ui ] = P[G b (x ui ) ≤ v|λ ui ] (4) = P[λ ui + ε ui < b v ] (5) = P[ε ui < b v -λ ui ] (6) = F ε (b v -λ ui ) . (7)
It follows that the function v → P[y ui ≤ v|λ ui ] is increasing since the sequence of thresholds is itself increasing. Moreover, the probability mass function (p.m.f.) associated to the ordinal data can be written as:

P[y ui = v|λ ui ] = P[y ui ≤ v|λ ui ] -P[y ui ≤ v -1|λ ui ] = F ε (b v -λ ui ) -F ε (b v-1 -λ ui ) . (8)
Some examples. If the c.d.f. is strictly increasing, we can rewrite Eq. ( 7) as:

F -1 ε (P[y ui ≤ v|λ ui ]) = b v -λ ui . (9)
Hence the name of CLM, since the factorization model is related to the c.d.f. of the ordinal data through a link function

F -1 ε : [0, 1] → R.
Various choices of noise (equivalently, of link function F -1 ε ) have been considered in the literature. We present some of these choices in what follows.

• Logit function. The use of the logit function was first proposed in [START_REF] Walker | Estimation of the probability of an event as a function of several independent variables[END_REF]. This model was popularized and renamed as "proportional odds model" by (Mc-Cullagh, 1980). The model can be rewritten as:

logit P[y ui ≤ v|λ ui ] = log P[y ui ≤ v|λ ui ] P[y ui > v|λ ui ] = b v -λ ui (10)
• Probit function. A common choice for the additive noise is ε ui ∼ N (0, σ 2 ) [START_REF] Chu | Gaussian processes for ordinal regression[END_REF][START_REF] Paquet | A hierarchical model for ordinal matrix factorization[END_REF][START_REF] Hernandez-Lobato | Probabilistic matrix factorization with non-random missing data[END_REF]. In that case the link function F -1 ε is the probit function. Inference can be carried out with an EM algorithm based on the latent variable x ui .

• Other choices like log-log or cauchit functions have also been considered [START_REF] Agresti | Categorical data analysis[END_REF]. The survey [START_REF] Ananth | Regression models for ordinal responses: A review of methods and applications[END_REF] recaps some of these choices.

Bernoulli-Poisson Factorization (BePoF)

In this section, we present BePoF [START_REF] Acharya | Nonparametric Bayesian factor analysis for dynamic count matrices[END_REF] which is a variant of PF for binary data (not directly related to the CLMs introduced above). It employs a model augmentation trick for inference that we will use in our own algorithm presented in Section 4.

The Poisson distribution can easily be "augmented" to fit binary data y ui ∈ {0, 1}. Indeed, it suffices to introduce a thresholding operation that binarizes the data. The corresponding generative hierarchical model is therefore given by:

n ui ∼ Poisson([WH T ] ui ), (11)
y ui = 1[n ui > 0], (12)
where n ui ∈ N is a latent variable and 1 is the indicator function. We denote by N ∈ N U ×I the matrix such that [N] ui = n ui . This variable can easily be marginalized by noting that P[y ui = 0] = Poisson(0|[WH T ] ui ) = e -[WH T ]ui . We obtain:

y ui ∼ Bern(1 -e -[WH T ]ui ) (13)
where Bern refers to the Bernoulli distribution. The conditional distribution of the latent variable n ui is given by:

n ui |y ui ∼ δ 0 , if y ui = 0, ZTP([WH T ] ui ), if y ui = 1. ( 14
)
where ZTP refers to the zero-truncated Poisson distribution and δ 0 to the Dirac distribution located in 0. The latent variable N can be useful to design Gibbs or variational inference (VI) algorithms for binary PF [START_REF] Acharya | Nonparametric Bayesian factor analysis for dynamic count matrices[END_REF] and we will employ a similar trick in Section 4.

Remark. The generative model presented in Eq. ( 13) is in the form

y ui ∼ Bern(G([WH T ] ui )) where G : R (or R + ) → [0, 1]. When [WH T ] ui ∈ R, the function
G can be the inverse of the probit [START_REF] Consonni | Mean-field variational approximate Bayesian inference for latent variable models[END_REF] or of the logit function for example. They are special cases of the model presented in Section 2.1 with V = 1. Mean-parametrized Bernoulli MF models have also been considered [START_REF] Lumbreras | Bayesian meanparameterized nonnegative binary matrix factorization[END_REF]. They correspond to G = Id and require additional constraints on the latent factors W and H in order to satisfy

[WH T ] ui ∈ [0, 1].

Ordinal NMF (OrdNMF)

In this section, we introduce OrdNMF which is a NMF model pecially designed for ordinal data. A difference with Section 2.1 is that we impose that both matrices W and H are non-negative. Thus, we now have

[WH T ] ui ∈ R + instead of [WH T ] ui ∈ R. We denote λ uik = w uk h ik so that λ ui = k λ uik = [WH T ] ui .
3.1. Quantization of the Non-negative Numbers

Our model works on the same principle as OrdMF (see Section 2.1) and seeks to quantize the non-negative real line R + . For this, we introduce the increasing sequence of thresholds b given by b

-1 = 0 < b 0 < • • • < b V -1 < b V = +∞ (
the thresholds are here non-negative). Moreover, we define the quantization function G b : R + → {0, . . . , V } like in Eq. ( 1) but with support R + .

As compared to Section 2.1, we now assume a non-negative multiplicative noise on x ui . This ensures the non-negativity of x ui and it seems well suited for modeling over-dispersion, a common feature of recommendation data. Let ε ui be a nonnegative random variable with c.d.f. F ε , we thus propose the following generative model:

x ui = λ ui • ε ui , (15)
y ui = G b (x ui ). (16)
Like before, our goal is to jointly infer the latent variables W and H as well as the sequence of the thresholds b. In our model, the c.d.f. associated to the ordinal random variable y ui becomes:

P[y ui ≤ v|λ ui ] = P[G b (x ui ) ≤ v|λ ui ] (17) = P[λ ui • ε ui < b v ] (18) = P ε ui < b v λ ui (19) = F ε b v λ ui . (20)
Therefore, we can deduce that the p.m.f. is given by:

P[y ui = v|λ ui ] = P[y ui ≤ v|λ ui ] -P[y ui ≤ v -1|λ ui ] (21) = F ε b v λ ui -F ε b v-1 λ ui . ( 22
)
Various functions F ε can be used which determine the exact nature of the multiplicative noise. Figure 3 displays the function λ → F ε (λ -1 ) for the examples considered next.

• Gamma noise: ε ui ∼ Gamma(α, 1). 1 The c.d.f. is given 1 The rate parameter β is fixed to 1 because of a scale invariance with λui. by

F ε (x) = γ(α,x) Γ(α) where γ(α, x) = x 0 t α-1 e -t dt is the lower incomplete gamma function. If α = 1, we recover an exponential noise ε ui ∼ Exp(1) whose c.d.f. is F ε (x) = 1 -e -x .
• Inverse-gamma (IG) noise:

ε ui ∼ IG(α, 1). 1 The c.d.f. is given by F ε (x) = Γ(α,x -1 ) Γ(α)
where Γ(α, x) = ∞ x t α-1 e -t dt is the upper incomplete gamma function. If α = 1, we obtain the c.d.f. F ε (x) = e -1/x .

• Any increasing function F ε : R + → [0, 1] defines a nonnegative random variable which can be used in OrdNMF.

OrdNMF with IG Noise (IG-OrdNMF)

In the rest of the paper, we focus on the special case where ε ui is a multiplicative IG noise with shape parameter α = 1, i.e., ε ui ∼ IG(1, 1)2 . We use the acronym IG-OrdNMF for this particular instance of OrdNMF.

For convenience we write θ v = b -1 v . The sequence θ corresponds to the inverse of the thresholds and is therefore decreasing, i.e., θ

-1 = +∞ > θ 0 > • • • > θ V -1 > θ V = 0.
Moreover, we denote by ∆ the positive sequence of decrements defined by

∆ v = θ v-1 -θ v for v ∈ {1, . . . , V }. We have θ v = V l=v+1 ∆ l and, in particular, θ V -1 = ∆ V .
Interpretation. In IG-OrdNMF model, the c.d.f. associated with an ordinal data y ui is given by:

P[y ui ≤ v|λ ui ] = e -λuiθv , (23)
or

P[y ui > v|λ ui ] = 1 -e -λuiθv , (24)
with v ∈ {0, . . . , V }. Therefore, BePoF (see Section 2.2) is a particular case of IG-OrdNMF with V = 1 and θ 0 = 1.

This formulation allows for a new interpretation of IG-OrdNMF. As a matter of fact, the event {y ui > v} is a binary random variable which follows a Bernoulli distribution: {y ui > v} ∼ Bern(1-e -λuiθv ). Then, we can see IG-OrdNMF as the aggregation of V dependent BePoF models for different thresholds of binarization v ∈ {0, . . . , V -1}.

Probability mass function. The p.m.f. of an observation is given by:

P[y ui = v|λ ui ] =      e -λuiθ0 , for v = 0, e -λuiθv -e -λuiθv-1 , for 1 ≤ v < V, 1 -e -λuiθ V -1 , for v = V. (25)
Then, the log-likelihood of λ ui can be written as:

log P[y ui = v|λ ui ] = -λ ui θ 0 , if v = 0, -λ ui θ v + log(1 -e -λui∆v ), else. ( 26
)
This expression brings up a linear term in λ ui and a nonlinear term of the form x → log(1 -e -x ), similar to the function used in Section 2.2.

Moreover, the expectation of the observations is welldefined and given by :

E(y ui |λ ui ) = V - V -1 v=0 e -λuiθv . (27)
Note that, in the context of ordinal data processing, the expectation is not a good statistic since it implicitly implies a notion of distance between classes. However, this quantity will be useful to build lists of recommendations. Indeed, the function λ ui → E(y ui |λ ui ) is increasing. Thus, the higher the λ ui = [WH T ] ui , the higher (in expectation) the level of interaction between the user and the item.

Bayesian Inference

We impose a gamma prior on the entries of both matrices W and H, i.e., w uk ∼ Gamma(α W , β W u ) and h ik ∼ Gamma(α H , β H i ). Gamma prior is known to induce sparsity which is a desirable property in NMF methods.

Augmented Model

As described in Section 3.2 , the log-likelihood for ordinal data such that v ∈ {1, . . . , V } brings up a non-linear term log(1 -e -x ) which is not conjugate with the gamma distribution, making the inference complicated. To solve this issue, we use the trick presented in Section 2.2 by augmenting our model with the latent variable:

n ui |y ui , λ ui ∼ δ 0 , if y ui = 0, ZTP(λ ui ∆ yui ), if y ui > 0. (28)
Moreover, as commonly done in the PF setting [START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF][START_REF] Gopalan | Scalable recommendation with hierarchical Poisson factorization[END_REF], we augment our model with the latent variable c ui |n ui , λ ui ∼ Mult(n ui , φ ui ), where Mult is the multinomial distribution and φ ui is a probability Var. Distribution

C q(c ui |n ui ) = Mult c ui ; n ui , φui N q(n ui ) = δ 0 , if y ui = 0 ZTP(n ui ; Λ ui ∆ yui ), if y ui > 0 W q(w uk ) = Gamma w uk ; αW uk , βW uk H q(h ik ) = Gamma h ik ; αH ik , βH ik
vector with entries λ uik λui . Therefore, for ordinal data y ui ∈ {1, . . . , V }, we obtain the following joint log-likelihood:

log p(y ui , n ui , c ui |λ ui ) = -λ ui θ yui-1 (29) + n ui log ∆ yui + k (c uik log λ uik -log c uik !) , s.t. n ui ∈ N * and n ui = k c uik .
Joint log-likelihood of IG-OrdNMF. We denote by Z = {N, C, W, H} the set of latent variables of the augmented model. Moreover, we define T v such that:

T v = θ 0 , if v = 0, θ v-1 , if v > 0. (30)
The joint log-likelihood of IG-OrdNMF is therefore given by:

log p(Y, N, C|W, H) = ui n ui log ∆ yui + k (c uik log λ uik -log c uik !) -λ ui T yui . (31)
It is important to note that n ui = 0 and c ui = 0 K when y ui = 0. Consequently, the variables N and C are partially observed and the inference take advantage of the sparsity of the observed matrix Y.

Variational Inference

The posterior distribution p(Z|Y) is intractable. We use VI to approximate this distribution by a simpler variational distribution q. Here, we assume that q belongs to the meanfield family and can be written in the following factorized form:

q(Z) = ui q(n ui , c ui ) uk q(w uk ) ik q(h ik ). ( 32
)
Note that the variables N and C remains coupled. We use a coordinate-ascent VI (CAVI) algorithm to optimize the parameters of q. The variational distributions are described in Table 1. The associated update rules are summarized in Algorithm 1.

Approximation and link with PF. Algorithm 1 can be simplified by assuming that q(n ui ) = δ 1 if y ui > 0. This amounts to replacing the non-linear term log(1 -e -x ) by log x in Eq. ( 26). However, this approximation will produce similar results only if x is very small, since log(1 -e -x ) = log x+o(x). In practice, this can only be verified a posteriori by observing that E q (n ui ) ≈ 1.

As mentioned above, BePoF is a special case of IG-OrdNMF for V = 1 and θ 0 = 1. Thus, we can notice that PF algorithm applied to binary data is an approximation of BePoF algorithm for q(n ui ) = δ 1 if y ui = 1.

Thresholds Estimation

A key element of threshold models is the learning of thresholds (corresponding here to θ parameters). For this, we use a VBEM algorithm. It aims to maximize the term E q (log p(Y, Z; θ)), w.r.t. the variables θ, which is given by:

E q (log p(Y, Z; θ)) = (33) ui E q (n ui ) log ∆ yui -E q (λ ui )T yui + cst, s.t. θ 0 > θ 1 > • • • > θ V -1 > θ V = 0.
Note that both terms T v (defined in Eq. ( 30)) and

∆ v = θ v-1 -θ v > 0 depend on the sequence θ.
Decrements optimization. We choose to work on the decrement sequence ∆ rather than on the threshold sequence θ. Indeed, by doing so, the decreasing constraint of θ becomes a non-negativity constraint of ∆. Moreover, we obtain only terms in x and log x in the function to be maximized. Thus, the problem can be solved analytically.

We can rewrite the term T v w.r.t. the sequence ∆ by noting that:

T v = V l=1 1[v ≤ l]∆ l , ∀v ∈ {0, . . . , V }.
Therefore, the optimization problem presented in Eq. ( 33) amounts to maximizing the following function:

E q (log p(Y, Z; ∆)) = ui V l=1 1[y ui = l]E q (n ui ) log ∆ l -1[y ui ≤ l]E q (λ ui )∆ l + cst, s.t. ∆ ≥ 0. ( 34
)
Algorithm 1 CAVI for IG-OrdNMF. Data: Matrix Y Result: Variational distribution q and thresholds θ 1 Initialization of variational parameters and thresholds θ; repeat 2 foreach couple (u, i) such as y ui > 0 do 3

Λ uik = exp (E q (log w uk ) + E q (log h ik )); Λ ui = k Λ uik ; E q (n ui ) = Λui∆y ui 1-e -Λ ui ∆y ui ; E q (c uik ) = E q (n ui ) Λ uik Λui ; 4 end 5 foreach user u ∈ {1, . . . , U } do 6 αW uk = α W + i E q (c uik ); βW uk = β W u + i T yui E q (h ik ); 7 end 8 foreach item i ∈ {1, . . . , I} do 9 αH ik = α H + u E q (c uik ); βH ik = β H i + u T yui E q (w uk ); 10 end 11
Update of thresholds: Eq. ( 35) and Eq. ( 36); Update of rate parameters β W u and β H i ; Calculate ELBO(q, θ); 12 until ELBO converge; Thus, we obtain the following update rules:

∆ l = ui 1[y ui = l]E q (n ui ) ui 1[y ui ≤ l]E q (λ ui )
, ∀l ∈ {1, . . . , V }, (35)

θ v = V l=v+1 ∆ l , ∀v ∈ {0, . . . , V -1}. (36)
Algorithm 1 scales with the number of non-zero values in the observation matrix Y. The complexity of OrdNMF is of the same order of magnitude as BePoF and PF. The only difference with these algorithms in terms of computational complexity is the update of the thresholds (Line 11 of Alg. 1).

Posterior Predictive Expectation.

The posterior predictive expectation E(Y * |Y) corresponds to the expectation of the distribution of new observations Y * given previously observed data Y. This quantity allows us to create the list of recommendations for each user. We can approximate it by using the variational distribution q:

E(Y * |Y) ≈ W,H E(Y * |W, H)q(W)q(H)dWdH. (37)
Unfortunately, this expression is not tractable. But for recommendation we are only interested in ordering items w.r.t. 37) the simpler score s ui = [E q (W)E q (H) T ] ui .

Experimental Results

Experimental Set Up

Datasets. We report experimental results for two datasets described below.

• MovieLens [START_REF] Harper | The movielens datasets: History and context[END_REF]. This dataset contains the ratings of users on movies on a scale from 1 to 10. These explicit feedbacks correspond to ordinal data. We consider that the class 0 corresponds to the absence of a rating for a couple user-movie. The histogram of the ordinal data is represented in blue on Figure 4. We pre-process a subset of the data as in [START_REF] Liang | Modeling user exposure in recommendation[END_REF], keeping only users and movies that have more than 20 interactions. We obtain U = 20k users and I = 12k movies.

• Taste Profile [START_REF] Bertin-Mahieux | The million song dataset[END_REF]. This dataset, provided by the Echo Nest, contains the play counts of users on a catalog of songs. As mentioned in the introduction, we choose to quantize these counts on a predefined scale in order to obtain ordinal data. We arbitrarily select the following quantization thresholds: [1,2,5,10,20,50,100,200,500]. For example, the class labeled 6 corresponds to a listening counts between 21 and 50. As for MovieLens, the class 0 corresponds to users who have not listen to a song. The histogram of the ordinal data are displayed in blue on Figure 4. We pre-process a subset of the data as before and obtain U = 16k users and I = 12k songs.

Be careful not to confuse the predefined quantization used to obtain ordinal data, with the quantization of the latent variable in OrdNMF model which is estimated during inference. Although we expect OrdNMF to recover a relevant scaling between the categories, there is no reason to get the same quantization function that was used for pre-processing.

Evaluation. Each dataset is split into a train set Y train and a test set Y test : the train set contains 80% of the non-zero values of the original dataset Y, the other values are set to the class 0; the test set contains the remaining 20%. All the compared methods are trained on the train set and then evaluated on the test set.

First, we evaluate the recommendations with a ranking metric. For each user, we propose a list of m = 100 items (movies or songs) ordered w.r.t. the prediction score presented in Section 4.4: s ui = [E q (W)E q (H T )] ui . The quality of these lists is then measured through the NDCG metric [START_REF] Järvelin | Cumulated gain-based evaluation of IR techniques[END_REF]. The NDCG rewards relevant items placed at the top of the list more strongly than those placed at the end. We use the relevance definition proposed in [START_REF] Gouvert | Recommendation from raw data with adaptive compound Poisson factorization[END_REF]:

rel(u, i) = 1[y test ui ≥ s]. (38)
In other words, an item is considered as relevant if it belongs at least to the class s in the test set. The NDCG metric is between 0 and 1, the higher the better.

Moreover, for the Taste Profile dataset, we calculate the log-likelihood of the non-zero entries on the test set, as it is done in [START_REF] Basbug | Hierarchical compound Poisson factorization[END_REF]:

L N Z = (u,i)∈Testset log p(y test ui |y test ui > 0, Ŵ, Ĥ), (39)
where y test ui is the quantized data, ŵuk = E[w uk ] and ĥik = E[h ik ] are the estimated latent factors.

Compared methods. We compare OrdNMF with three other models: PF, BePoF and discrete compound PF (dcPF) [START_REF] Basbug | Hierarchical compound Poisson factorization[END_REF] with a logarithmic element distribution as implemented in [START_REF] Gouvert | Recommendation from raw data with adaptive compound Poisson factorization[END_REF]. Each model is applied either to raw data (R), quantized data (Q) or binarized data (B). For the MovieLens dataset, two different binarizations are tested: one with a threshold at 1 (≥ 1) and one with a threshold at 8 (≥ 8). For the Taste Profile dataset, dcPF is applied to the count data (R) whereas OrdNMF is applied to the quantized data (Q).

For all models, we select the shape hyperparameters α W = α H = 0.3 among {0.1, 0.3, 1} [START_REF] Gopalan | Scalable recommendation with hierarchical Poisson factorization[END_REF]. The number of latent factors is chosen among K ∈ {25, 50, 100, 150, 200, 350} for the best NDCG score with threshold s = 8 for the MovieLens dataset, and s = 1 for the Taste Profile dataset. All the algorithms are run 5 times with random initializations and are stopped when the relative increment of the expected lower bound (ELBO) falls under τ = 10 -5 . The computer used for these experiments was a MacBook Pro with an Intel Core i5 processor (2,9 GHz) and 16 Go RAM. All the Python codes are available on https://github.com/Oligou/OrdNMF.

Prediction Results

Table 2 displays the results for the MovieLens dataset. First, we can compare BePoF with its approximation, i.e., PF applied to binarized data. BePoF is slightly better than PF for both binarizations, and requires less latent factors. Then, we observe that the choice of the binarization has a big impact on the NDCG scores. BePoF with data thresholded at 1 (≥ 1) perform well on small NDCG threshold s but has poor performance after. On the contray, with data thresholded at 8 (≥ 8), BePoF achieves best performances for NDCG s = 10 but poor performances for small s. OrdNMF does not exhibit such differences between NDCG scores and benefits from the additional information brought by the ordinal classes. Nevertheless, we can note a small decrease of the performance with s = 10 which is the hardest class to predict.

Table 3 displays the same kind of results for the Taste Profile dataset. Again, OrdNMF outperforms BePoF and PF which exploit less data information. OrdNMF is competitive with dcPF which gives the best results for the highest thresholds s. However, OrdNMF presents a higher log-likelihood score than dcPF. Thus, OrdNMF seems better suited to predict the feedback class of a user than dcPF. This observation is confirmed by the posterior predictive checks (PPC) presented below.

Posterior Predictive Check (PPC)

A PPC consists of generating new data based on the posterior predictive distribution p(Y * , W, H|Y) ≈ p(Y * |W, H)q(W)q(H), and then compare the structure of the original data Y with the artificial data Y * . Here, we focus on the distribution of the ordinal categories. Figure 4 presents the results of these PPCs for the Taste Profile dataset. The blue bars correspond to the empirical histogram of the data (Y train ), the red and orange bars correspond to the histograms of the simulated data obtained with OrdNMF and dcPF respectively. While dcPF fails to model the very large values present in the data (from the class 7, which corresponds to values greater than 50 plays), OrdNMF seems to precisely describe all the ordinal categories. This is also the case on the MovieLens dataset too. Even if the empirical histogram is here less regular, OrdNMF can adapt itself to all type of data through the inferred thresholds b.

Conclusion

We developed a new probabilistic NMF framework to process ordinal data. In particular, we presented IG-OrdNMF which is an extension of BePoF and conducted experiments on two different datasets. We show the ability of OrdNMF to process different kinds of ordinal data both explicit and implicit. This work opens up several exciting perspectives. As we described in Section 3.1, OrdNMF can be used for different choices of multiplicative noise. It would be interesting to develop OrdNMF for the exponential noise in a similar way than IG-OrdNMF. Finally, when applied to implicit data, it would be of particular interest to learn the pre-processing during the factorization, in order to automatically tune the level of pre-processing adapted to a given dataset. This is yet left for future investigations.

Figure 1 .

 1 Figure1. Graphical model of OrdMF. A latent variable X is introduced to make the link between the factorization term WH T and the ordinal data Y.

Figure 3 .

 3 Figure 3. Functions λ → Fε(λ -1 ) for gamma (left) and inversegamma (right) noises.

Figure 4 .

 4 Figure 4. PPC of the distribution of the classes in the Taste Profile dataset (left) and MovieLens dataset (right). The blue bars (Truth) represents the histogram of the classes in the train set. The colored bars represent the simulated histograms obtained from the different inferred OrdNMF or dcPF models. The percentages of non-zero values are written in parentheses.

Table 1 .

 1 Variational distributions for IG-OrdNMF.

Table 2 .

 2 Recommendation performance of OrdNMF using the MovieLens dataset. Bold: best NDCG score. R: raw data. B: binarized data.

					NDCG @100 with threshold s
	Model	Data	K	s = 1	s = 4	s = 6	s = 8	s = 10
	OrdNMF R	150 0.444 0.444 0.439 0.414 0.353
	BePoF	B (≥ 1) 50	0.433	0.430	0.421	0.383	0.310
	PF	B (≥ 1) 100 0.431	0.428	0.418	0.380	0.306
	BePoF	B (≥ 8) 50	0.389	0.393	0.399	0.408	0.369
	PF	B (≥ 8) 150 0.386	0.389	0.395	0.403	0.365

Table 3 .

 3 Recommendation performance of OrdNMF using the Taste Profile dataset. Bold: best NDCG or log-likelihood score. Q: quantized data. B: binarized data. R: raw data. The function λ ui → E(y * ui |λ ui ) being increasing, we can use instead of Eq. (

				NDCG @100 with threshold s
	Model	Data	K	s = 1 s = 3 s = 6 s = 11 s = 21 s = 51	log-lik
	OrdNMF Q dcPF R BePoF B (≥ 1) 250 0.210 0.170 0.149 0.131 250 0.213 0.174 0.153 0.135 150 0.209 0.173 0.154 0.137	0.123 0.128 0.120	0.117 0.121 0.115	-2.8 • 10 5 -3.0 • 10 5 N/A
	PF	B (≥ 1) 250 0.206 0.167 0.146 0.129	0.118	0.115	N/A
	this quantity.

IRIT, Université de Toulouse, CNRS, France

ISAE-SUPAERO, Université de Toulouse, France. Correspondence to: Olivier Gouvert <oliviergouvert@gmail.com>.

The expectation of a IG variable is not defined for α ≤ 1, however the model is still well-defined.

Acknowledgements

Supported by the European Research Council (ERC FACTORY-CoG-6681839) and the ANR-3IA (ANITI).