Chabane Chenane

Dalil Ichalal

Hichem Arioui

Said Mammar

Lateral Dynamics Reconstruction for Sharp'71 Motorcycle Model with P2I Observer

The main objective of this paper is the reconstruction of lateral dynamics and both roll angle and steering torque of single track vehicles (motorcycle, scooter, etc.). For that purpose, the well-known motorcycle model developed by Sharp in 1971 is used. This model characterizes the lateral dynamics of a motorcycle [16]. The roll angle is not observable in the obtained structure, for overcoming this problem, the model is transformed in order to take into account the roll angle as an unknown input as well as the steering torque. A Proportional two integrals (P2I) observer is then proposed for estimating simultaneously all the variable states, the lateral forces and both roll angle and steering torque. This study is a part of the ongoing work of the research team on the design of preventive safety systems for motorcycles users. Simulation results and discussions are given in order to illustrate the effectiveness of the proposed observer.

I. INTRODUCTION

Recently, the park of single track vehicles is constantly increasing, upsetting driving practices and road traffic. Unfortunately, this expansion resulted in a growth of traffic fatalities. The statistics endorse this statement and riders are considered as the most vulnerable road users. In 2010, the French Agency of Road Safety made a finding of around 1000 deaths (25% of traffic fatalities), while the traffic volumes of motorcycles does not exceed 1%, [START_REF]characteristics and causes of road accidents[END_REF]. Many research projects are initiated to fulfill this issue in order to propose an enhancement in term of security through preventive and / or active safety systems, [START_REF] Amans | Rider project . research on accidents involving powered two-wheelers[END_REF], [START_REF] Bekiaris | Saferider european project[END_REF].

The achievement of safety systems depends on the proper knowledge of: 1) the dynamics of single track vehicle, and 2) the evolution of its states (observation/estimation) , and to lesser extent, 3) the road geometry. Regarding the first point, several studies were carried out in order to understand the motorcycle dynamics [START_REF] Sharp | Advances in the modelling of motorcycle dynamics[END_REF], [START_REF] Cossalter | Motorcycle Dynamics[END_REF], stability analysis (eigenmodes) of PTW [START_REF] Åström | Bicycle dynamics and control[END_REF], optimal and safe trajectories [START_REF] Bobbo | Application of the optimal maneuver method for enhancing racing motorcycle performance[END_REF] and the proposition of risk functions [START_REF] Slimi | Motorcycle speed profile in cornering situation[END_REF], [START_REF] Evangelou | Influence of road camber on motorcycle stability[END_REF] to control loss or equilibrium margins. These works are little sustainable if they are not supported by an efficient sensor system helping in the estimation of some dynamic states.

The use of sensors is not always possible for two main reasons: 1) instrumentation can be very expensive and leading inevitably to expensive new bikes, and 2) according to used technologies, the measurement noise can seriously compromise the future safety systems. Accordingly, we suggest the use of robust observation techniques to overcome {chabane.chenane, dalil.ichalal, hichem.arioui, said.mammar}@ibisc.univ-evry.fr the aforementioned shortcomings. Nevertheless, including all methodologies (Luenberger observer, Takagi-Seguno based observer, Extended Kalman filter), very few studies have been conducted on the estimation of motorcycle dynamic states [START_REF] Ichalal | Observer design for two-wheeled vehicle: A takagi-sugeno approach with unmeasurable premise variables[END_REF], [START_REF] Van Daal | Design and automatic tuning of a motorcycle state estimator[END_REF]. The present paper proposes a robust proportional two integral (P2I) approach, [START_REF] Koenig | Unknown input proportional multiple-integral observer design for linear descriptor systems: application to state and fault estimation[END_REF], [START_REF] Ichalal | Simultaneous state and unknown inputs estimation with PI and PMI observers for Takagi-Sugeno model with unmeasurable premise variables[END_REF], [START_REF] Linder | Estimating and accommodating unknown actuator faults with pi observers[END_REF], [START_REF] Shafai | Simultaneous disturbance attenuation and fault detection using proportional integral observers[END_REF], helping in states observation of linear motorcycle model and the reconstruction of rider's action (steering torque). Disturbance and the rider action are assumed to be almost affine. This assumption is quite realistic according to common variation of roads profile. An H ∞ performance index is included during the design in order to attenuate the effect of the second derivative of the Unknown Input (UI) on the estimation error.

This paper is organized as follows: section III is dedicated to the motorcycle model description. Section IV presents the robust estimation of the motorcycle states and UI reconstruction. Finally, simulation results and conclusions complete the paper.

II. OVERVIEW ON SAFETY QUANTIFICATION

One of the long-term objective of our studies is concern with the problem of quantifying the risk of loss-of-control over a motorcycle when cornering. To make a Safe Cornering, riders should observe some wariness : 1) maintain an appropriate speed when entering the corner, 2)choose a good position on the road, 3) use Counter steering phenomenon to lean the bike and 4) in curve, maintain a constant speed and smoothly accelerate when exit.

These principles must be inflected with other factors that may make the situation critical. Indeed, a road with limited adhesion, 2) poor weather conditions are dangerous.

Regarding early warning systems, loss of control, for standard cars, is expressed by the maximum speed at which a vehicle can be kept on the road while moving at a constant speed on a circular section giving by:

v max = g • µ lat ρ (1)
where, g is the acceleration due to gravity and µ lat is the maximum available side friction, v m ax represents the authorized longitudinal velocity in cornering situation and ρ is the curvature of the road.

Other risk function are proposed by the National Highway Traffic Safety Administration (NHTSA) which recommend a maximum safe speed governed by the following equation:

v max = g ρ φ r + µ lat 1 -φ r µ lat (2)
where φ r is the road super-elevation angle. Please refer to [START_REF] Slimi | Motorcycle speed profile in cornering situation[END_REF] for more precise models adapted for motorcycles.

In general, computing the lateral friction µ lat involves all the dynamic states of the bike and a good interpretation of the tire-road contact. This makes the success of such warning system strongly dependent on the availability of dynamic states of the motorcycle (efficient sensors or robust observers).

III. DYNAMIC DESCRIPTION AND MODEL OF MOTORCYCLE

To study the dynamics of Two Wheels vehicle, often it is assumed that overall vehicle and driver consists of interconnected rigid bodies with the possibility of rotation of each body around predefined axes [START_REF] Cossalter | A motorcycle multi-body model for real time simulations based on the natural coordinates approach[END_REF]. According to the assumed number of rigid bodies, a number of corresponding degrees of freedom which characterizes the system motion is considered. Then, motorcycle's motion can be characterized by two main modes: in-plane mode representing movements in its plane of symmetry (longitudinal and pitch displacements) and the out-of-plane mode represented by the lateral dynamics when cornering [START_REF] Sharp | Advances in the modelling of motorcycle dynamics[END_REF].

In this work we exploited the reference model Sharp'71 [START_REF] Sharp | The stability and control of motorcycle[END_REF], which is a linear one when assuming both a constant forward speed and Pacejka's linear form of the lateral forces. The study concerns the normal driving situation with where the driver do not make high maneuvers (small roll and steering angle) and not sudden movement (small variations of the steering torque) which justifies all the linearizations leading to a linear model.

The considered model of the motorcycle contains two bodies: the rear body G r which includes the chassis, engine and the rear wheel, the front body G f which represents the steering assembly and the front wheel. This configuration allows 4 degrees of freedom: the lateral displacement y, the yaw rotation ψ, the roll inclination φ and the handlebar steering angle δ with respect to the rider torque input τ applied on the motorcycle's handlebar (Figure 1).

According to all these assumptions, The resulting set of four second-order differential equations of motion is compactly written as:

M         vy ψ Ḟy f Ḟyr δ φ         + H         v y ψ F y f F yr δ φ         + P δ φ + G τ = 0 (3)
A full description of all terms of matrices is given in the Appendix.

G r G f v x v y F y f F yr δ , τ φ ψ τ d -τ b ε ε Fig. 1. Geometrical representation of the Sharp's motorcycle model
Now we recast the equation of motion in state space formulation as follows:

ẋ(t) = A x(t) + Bτ (4)
The components of the state variable vector x(t) are motion coordinates lateral forces and velocities, such that

x(t) = [δ φ v y ψ F y f F yr δ φ ] T .
The Linear Time Invariant matrices A and B are given by:

A = 0 I -M -1 P -M -1 H , B = 0 -M -1 G
Measuring the roll angle rate φ (t), one can note that it is not possible to obtain directly φ (t) by integration. Indeed, computing φ (t) from its derivative requires the knowledge of the initial conditions which are not necessarily known. Then the second component φ (t) of the state vector is not observable. To overcome this limitation, the system's equations are rewritten by separating the observable and non observable state variables. In the obtained new model, the dimension of the state vector is reduced and it does not contain the roll angle φ (t). The roll angle together with the steering torque are considered as unknown inputs. One obtains the following model:

ẋ(t) = Ax(t) + B 1 φ + B 2 τ (5) y(t) = Cx(t) (6)
where the new state vector x(t) is the state vector given by

[δ v y ψ F y f F yr δ φ ] T . A = 0 I -M -1 P(:,1) -M -1 H
, B 1 = P(:,1)

B 2 = 0 -M -1 G C =     1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1    
As a conclusion, the motorcycle model, is expressed as a state space linear model with two unknown inputs: the roll angle and the steering torque. In the next section, a Proportional two Integrals observer (P2I) is designed for estimating simultaneously the state and unknown inputs.

IV. PROPORTIONAL TWO INTEGRALS OBSERVER DESIGN

Observers design for linear models is a problem which has been dealt with intensively since Luenberger's work [START_REF] Luenberger | An introduction to observers[END_REF]. And Proportional-Integral (PI) observer is one of the extensions based on adding an integral action of the estimation error which makes it effective in estimating constant and slowly time varying unknown inputs. The PI observer can be extended for the estimation of non constant UI having a polynomial form by the use of Proportional Multiple Integral loops (PMI). Here it is assumed that the two inputs (roll angle and steering torque) are unknown. The proposed observer considers two integral actions which leads to a proportional two-Integrals (P2I) observer which is able to estimate the two inputs and their first derivatives. Consequently, The form of the unknown inputs for which the asymptotic convergence can be obtained are in first order polynomial form. This correspond to the condition that the second derivatives of the unknown inputs are zero. In order to extend the class of unknown inputs which can be estimated by the P2I observer, it is possible to relax the last condition by assuming that the second derivatives of the UI are not zero but bounded. The observer is then designed in such a way to estimate the unknown inputs, their first derivatives and to minimize the effect of the bounded second derivatives on the state and unknown input estimation errors.

For that purpose the previous system ( 5) is augmented by considering additional state variables ξ1 = φ , ξ2 = ξ 1 = φ , ξ3 = τ, ξ4 = ξ 3 = τ which leads to the augmented system described by state vector, x = x T ξ 1 ξ 2 ξ 3 ξ 4

T . The new state space system is given by:

ẋ(t) = Ā x(t) + B ū(t) (7) y(t) = C x(t) (8)
with

Ā =       A 0 B 1 0 B 2 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0       , B =       0 0 1 0 0 0 0 1 0 0       ū(t) = χ(t) α(t) and C = C 0
where χ(t) = φ (t) and α(t) = τ(t).

The P2I observer has the following form:

               ẋ(t) = A x(t) + B 1 ξ2 + B 2 ξ4 + K p (y(t) -ŷ(t)) ξ1 (t) = K i1 (y(t) -ŷ(t)) ξ2 (t) = ξ1 + K i2 (y(t) -ŷ(t)) ξ3 (t) = K i3 (y(t) -ŷ(t)) ξ4 (t) = ξ3 + K i4 (y(t) -ŷ(t)) (9)
It is clear that ξ1 = φ is estimated in the initial state vector x(t), so we can use it to compute directly φ without need to minimize φ .

Then the new obtained P2I observer can be expressed as follows:

           ẋ(t) = A x(t) + B 1 ξ2 + B 2 ξ4 + K p (y(t) -ŷ(t)) ξ2 (t) = φ + K i1 (y(t) -ŷ(t)) ξ3 (t) = K i2 (y(t) -ŷ(t)) ξ4 (t) = ξ3 + K i3 (y(t) -ŷ(t)) (10)
where K p , K i1 , K i2 and K i3 are the observer gain matrices. By defining the new augmented state vector x = x T ξ 2 ξ 3 ξ 4 T , the state space observer becomes:

ẋ(t) = Āx (t) + L(y(t) -Cx (t)) (11)
where

Ā =     A B 1 0 B 2 Ψ 0 0 0 0 0 0 0 0 0 1 0     , B =     0 0 1 0    
and C = C 0 and Ψ = 0 0 0 0 0 0 1

The gain of the observer is given by the matrix K to be designed which includes all the gains K p , K i1 , K i2 and K i3 . Let us define the state estimation error e(t) = x -x. Its time derivative is given by:

ė(t) = ẋ(t) -ẋ(t) (12) = ( Ā -K C)e(t) + Bα(t) (13)
Consider the Lyapunov function:

V (t) = e T (t)Pe(t) (14)
where P is a symmetric and positive definite matrix. The time derivative of the Lyapunov function is given by:

V = e T (t)[P( Ā -K C) + ( Ā -K C) T P]e(t) + α T BT Pe(t) + e(t) T P Bα(t) (15)
The observer gain K is computed in order to stabilize the system generating the state estimation error and also to attenuate the effect of the second derivative α(t) of the Unknown Input on the estimation errors, namely T eα ∞ < γ, γ > 0, where T eα represents the transfer from α(t) the error

e(t).

Let us define the L 2 -gain of the system as the quantity:

sup e(t) 2 α(t) 2 , α(t) 2 = 0 ( 16
)
where the L 2 -norm of e(t) and α(t) are defined by:

e(t) 2 = ∞ 0 e(t) T e(t)dt 1/2 (17) α(t) 2 = ( ∞ 0 α(t) T α(t)dt) 1/2 (18)
Then, if the inequality

V (t) + e T e -γ 2 α T α < 0 (19)
holds, the state estimation error dynamics is stable and the transfer from α(t) to e(t) is bounded by γ. By replacing [START_REF] Shafai | Simultaneous disturbance attenuation and fault detection using proportional integral observers[END_REF] in [START_REF] Slimi | Motorcycle speed profile in cornering situation[END_REF] one obtains:

e T (t)[P( Ā -K C) + ( Ā -K C) T P + I]e(t) + α T (t) BT Pe(t) + e T (t)P Bα(t) -γ 2 α T (t)α(t) < 0 (20)
Since the final objective is to derive LMI conditions, let us consider the change of variable Z = PK. The inequality ( 20) is equivalent to:

e T (t)[ ĀT P + P Ā -CT Z T -Z C + I]e(t) + α T (t) BT Pe(t) + e T (t)P Bα(t) -γ 2 α T (t)α(t) < 0 (21)
In matrix formulation, ( 21) is equivalent to: by choosing γ as a variable and using the change of variables γ = γ 2 . After solving this optimization problem, the gains of the observer are obtained by

e(t) α(t) T ĀT P + P Ā -CT Z T -Z C + I P B BT P -γ 2 I e(
K = P -1 Z (25)
and the attenuation gain is given by γ = γ.

V. SIMULATION ANALYSIS The main objective of this work is to reconstruct the lateral dynamics, especially, the roll angle, the lateral forces and the steering torque. To limit the effect of the oscillatory phenomenon in the transit phase the poles of the matrix A -LC are assigned in a LMI region S defined by S = {z ∈ C | Re(z) < -a, |z| < R}, which is an intersection between the left plan defined by the Re(z) < -a and the disc with center (0, 0) and radius R. Thus, we solve simultaneously the proposed optimization problem and the LMI constraints, corresponding to the LMI region, given as follows (for more details see [START_REF] Chilali | H-infinity design with pole placement constraints : an LMI approach[END_REF]).

ĀT P + P Ā -ZC -C T ZT + 2aP < 0 (26)

-RP P Ā -ZC ĀT P -C T ZT -RP < 0 (27)
In all the simulations we fixed a = 10 and R = 50. Simulation results prove the efficiency of the motorcycle model used here as illustrated in figure 2 where we show that the equilibrium condition [START_REF] Slimi | Assistance system security for Driving Vehicles of Two-Wheelers[END_REF] is verified. The torque has been chosen in such a way to excite the lateral vehicle dynamics which corresponds to a cornering maneuver as shown in the same figure.

A. Results without measurement noise

The two figures (3 and 4) show the evolution of the state variables and the estimate of each state. We see that the observer converges quickly and the state estimation errors converge accurately to real states. Figure 5 represents the roll angle estimation of the vehicle and the reconstruction of the steering torque, satisfactory results are then obtained.

B. Results with measurement noise

Now, consider the same observer in the presence of measurement noise of order 10% of maximal value of each output. The obtained results are depicted in figures 6, 7 and 8. One can conclude that even if the the measurements are affected by noises, an acceptable state and unknown input estimations are obtained.

VI. CONCLUSION AND FUTURE WORK

In this paper, a synthesis of P2I observer was presented to estimate the dynamic states of a two-wheeled vehicle, the roll angle that is often inaccessible to measure and reconstruct the steering torque which is the control input (rider's action on the handlebar). The stability of the observer is studied with Lyapunov theory and LMI conditions are established. Finally, the efficiency of the observer are illustrated by some simulation results. For Future works, it will be interesting to extend the approach for a nonlinear model of the motorcycle by using Takagi-Sugeno fuzzy structure and Linear Parameter Varying (LPV) model in order to take into account some nonlinear behaviors and the variation of the longitudinal velocity. In addition, it remains to achieve these results through a validation of the prototype test available in the laboratory, this aims to test some risk functions developed for strong risks related at preventive security systems development.

M =        M f + M r M f k 0 M f k M f k 2 + I rz + I f x sin ε 2 + I f z cos ε 2 0 0 0 σ f /V 0 0 0 M f e M f ek + I f z cos ε 0 M f j + M r h M f jk -C rxz + (I f z -I f x) sin ε cos ε 0 0 M f e 0 M f ek + I f z cos ε 0 0 σ r /V 0 0 I f z + M f e 2 0 M f e j + I f z sin ε M f j + M r h (M f jk -C rxz + (I f z -I f x ) sin ε cos ε 0 0 M f e j + I f z sin ε (M f j 2 + M r h 2 + I rx + I f x (cos ε 2 + I f z sin ε 2 )        H =        0 (M f + M r )V -1 -1 0 0 M f kV -l b -(i f y /R f ) sin εV C f 1 /V C f 1 l/V 1 0 -C f 1 t/V C r1 /V -C r1 b/V 0 1 0 0 (M f e + (i f y /R f ) sin εV t 0 K 0 (M f j + M r h + (i f y /R f ) + (L/R r ))V 0 0 (i f y /R f ) cos εV
/R f ) + (L/R r ))V 0 0 -(i f y /R f ) cos εV 0        P =        0 0 0 0 -(C f 1 cos ε +C f 2 sin ε -C f 2 0 -C r2 (tZ f -M f eg) sin ε (tZ f -M f eg) tZ f -M f eg -(M f j + M r h)g        , G =        0 0 0 0 -1 0       
where the parameters are defined as follows:

M f , M

Fig. 2 .

 2 Fig. 2. Vehicle trajectory (top) and Comparison between guiding force and lateral balance force (bottom)

Fig. 3 .Fig. 4 .

 34 Fig.3. mesurments (blue line) and her estimations (dashed red line)

Fig. 5 .Fig. 6 .

 56 Fig. 5. Roll angle and its estimation (top)and Steering torque and its estimation (bottom)

Fig. 7 .Fig. 8 .

 78 Fig.7. no mesured states (blue line) and her estimations (dashed red line)

  Furthermore, in order to enhance the performances of the observer, it is possible to minimize the transfer gain γ subject to LMI constraints.

				t) α(t)	< 0
				(22)
	The quadratic form (22) is negative definite if and only if
	the linear matrix inequality:
	ĀT P + P Ā -CT Z T -Z C + I P BT P -γ 2 I B	< 0	(23)
			The following
	optimization problem is then stated
	min	γ
	s.t.
	ĀT P + P Ā -CT Z T -Z C + I P BT P -γI B	< 0	(24)

holds. Finally, given a scalar γ, if there exists a symmetric and positive definite matrix P and a matrix Z such that the LMI (23) is satisfied, then the system generating the state estimation error is stable and the transfer from α(t) to e(t) is bounded by the L 2 gain γ.

  r Mass of front/rear frame (M = M f + M r ) j Distance between the center of gravity of the front frame and ground k Distance between the centers of gravity of each frame L f , L r Distance between the center of gravity e Distance between the fork and the center of gravity h Height of the center of gravity i f y , i ry Polar moment of inertia of front/rear wheel I rz Camber inertia of rear wheel. R f , R r radius of the front/rear wheel g Acceleration due to gravity. σ f , σ r Front and rear tire relaxation lengths respectively C f 1 C f 2 Front and rear tire cornering stiffnesses respectively C r1 C r2 Front and rear tire camber stiffnesses respectively t, V Pneumatic trail and Forward speed Z f , ε Front wheel load and Steering head angle


Rémi Bardenet
email: bardenet@lri.fr

Olivier Cappé
email: cappe@telecom-paristech.fr

Balázs Kégl
email: balazs.kegl@gmail.com

Adaptive Metropolis with Online Relabeling

We propose a novel adaptive MCMC algorithm named AMOR (Adaptive Metropolis with Online Relabeling) for efficiently simulating from permutation-invariant targets occurring in, for example, Bayesian analysis of mixture models. An important feature of the algorithm is to tie the adaptation of the proposal distribution to the choice of a particular restriction of the target to a domain where label switching cannot occur. The algorithm relies on a stochastic approximation procedure for which we exhibit a Lyapunov function that formally defines the criterion used for selecting the relabeling rule. This criterion reveals an interesting connection with the problem of optimal quantifier design in vector quantization which was only implicit in previous works on the label switching problem. In benchmark examples, the algorithm turns out to be fastconverging and efficient at selecting meaningful non-trivial relabeling rules to allow accurate parameter inference.

Introduction

Adaptive Metropolis (AM) [1] is a powerful recent algorithmic tool in numerical Bayesian data analysis. AM builds on a well-known Markov Chain Monte Carlo (MCMC) algorithm but optimizes the rate of convergence

to the target distribution by automatically tuning the design parameters of the algorithm on the fly. In the case of AM, the adaptive parameter tuning rule relies on the availability of an online estimate of the covariance of the target distribution. AM is considered a pioneering contribution in adaptive MCMC methodology, and it represents a significant step towards developing self-tuning, turn-key sampling algorithms which are necessary for using MCMC approaches in high throughput data analysis. The optimality criterion considered in AM relies on theoretical results derived when the target distribution is multivariate Gaussian [START_REF] Roberts | Weak convergence of optimal scaling of random walk Metropolis algorithms[END_REF][START_REF] Roberts | Optimal scaling for various Metropolis-Hastings algorithms[END_REF]. When the target is multi-modal or concentrated around a non-linear manifold, the algorithm still applies but its optimality is no longer guaranteed, and its practical performance is often suboptimal.

An important case where AM usually fails is when the target is the posterior distribution in Bayesian inference on a mixture model. In this case the mixture likelihood is invariant to permuting some of the mixture components, and the chosen prior often does not favor any permutation either. Although we usually have "nice" genuinely unimodal posteriors of the parameters of well-identified components, the posterior is highly multimodal with exponentially many peaks, one for each permutation of the components. Running a well-mixing MCMC on such a model results in useless marginal estimates for the parameters of the individual components due to the confusion between all possible labelings. This phenomenon is called label switching. Several approaches have been proposed to deal with this problem, usually in a post-processing step after the posterior sample has been produced [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components[END_REF][START_REF] Celeux | Bayesian inference for mixtures: The label-switching problem[END_REF][START_REF] Stephens | Dealing with label switching in mixture models[END_REF][START_REF] Marin | Bayesian modelling and inference on mixtures of distributions[END_REF][START_REF] Jasra | Markov Chain Monte Carlo methods and the label switching problem in Bayesian mixture modelling[END_REF][START_REF] Jasra | Bayesian inference for mixture models via Monte Carlo[END_REF][START_REF] Sperrin | Probabilistic relabelling strategies for the label switching problem in Bayesian mixture models[END_REF][START_REF] Papastamoulis | An artificial allocations based solution to the label switching problem in Bayesian analysis of mixtures of distribution[END_REF]. They all aim to solve the identifiability problem in order to produce meaningful marginal posteriors for the mixture parameters. Running vanilla AM on such a model has two pitfalls. If the chain does not mix well, we get stuck in one of the modes. The adaptive proposal can then be efficient but it is unclear what the empirical distribution of the chain is and it is quite likely that the results of inference will be ultimately biased with respect to the posterior distribution. On the other hand, if the MCMC chain does switch between components, the online sample covariance estimate will be too broad, resulting in poor adaptive proposals and slow convergence. Note that the latter situation is usually prevalent when using trans-dimensional samplers based on the Reversible Jump approach [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components[END_REF]: in this case, the dimensionchanging moves force label switching and running an embedded AM algorithm is, in our experience, almost useless.

To solve these challenging issues, we develop an Adaptive Metropolis algorithm with Online Relabeling (AMOR). The main idea is to combine the online relabeling strategy of [START_REF] Celeux | Bayesian inference for mixtures: The label-switching problem[END_REF] with the AM algorithm of [START_REF] Haario | An adaptive Metropolis algorithm[END_REF]. This results in a doubly adaptive algorithm that uses the sample statistics both for (i) adapting its proposal and for (ii) redefining the region onto which it constrains the chain by attempting relabeling in each iteration. We prove two important facts about the AMOR algorithm. First, we show that when the adaptation is frozen, the target of the algorithm is indeed a restriction of the original target distribution to one of its symmetric "modes". Thus, except for the automated selection of a particular relabeling strategy, the proposed algorithm does not modify the target distribution, a feature that is missing in some of the procedures commonly used to remedy label-switching. We also establish the existence of a Lyapunov function for AMOR, that is, a quantitative criterion which defines the set of possible limiting relabeling rules. This allows us to combine relabeling with adaptive proposals, something that is impossible in the usual setup where relabeling is done on the MCMC outcome in a postprocessing step. The proofs of these results -to be found in the supplementary material -also unveil interesting connections with the problem of optimal quantifier design in vector quantization for which we prove an extension of known results regarding quantification using Mahalanobis divergence.

The rest of the paper is organized as follows. In Section 2 we formalize the problem and motivate AMOR on a realworld example. In Section 3, we derive the new online relabeling procedure based on AM and present our main convergence results. We show experimental results in Section 4 and conclude in Section 5.

Adaptive Metropolis and label switching

In this section we briefly provide some more background regarding AM and the label switching problem. Readers familiar with these notions may skip this section.

We consider using MCMC sampling to explore the posterior distribution

π(x) p(x|y) ∝ p(y|x)p(x)
of the parameters x ∈X given the observation y. The posterior π(x), whose normalization constant is usually unknown, can be explored by running a Markov chain (X t ) with stationary distribution π. In this context, π is also said to be the target distribution of the MCMC chain. The Symmetric Random Walk Metropolis algorithm (SRWM; [START_REF] Metropolis | Equations of state calculations by fast computing machines[END_REF]; corresponding to the blue steps in Figure 2 in Section 3) is one of the most popular techniques for simulating such a chain (X t ). In SRWM the user has to provide a symmetric proposal kernel that will be used to propose a new sample X given the previous sample X t-1 . When the posterior is a distribution over a continuous space X = R d , the most common proposal kernel is a multivariate Gaussian N •|X t-1 , Σ .

The goal of Adaptive Metropolis (AM) (corresponding to the blue and green steps in Figure 2) is to automatically calibrate the design parameter Σ of SRWM. When the target π(x) is multivariate Gaussian with covariance Σ π , the optimal choice of Σ is of the order of (2.38) 2 Σ π /d [START_REF] Roberts | Weak convergence of optimal scaling of random walk Metropolis algorithms[END_REF][START_REF] Roberts | Optimal scaling for various Metropolis-Hastings algorithms[END_REF]. In practice, Σ π is unknown thus motivating the use of an estimate of the covariance of the posterior based on samples (X 1 , . . . , X t-1 ) generated so far. From a theoretical point of view, the conditional distribution of X t given the past then depends on the whole past, rendering the analysis of AM more challenging. The convergence of AM has been recently addressed under quite general conditions (see, e.g., [START_REF] Saksman | On the ergodicity of the adaptive Metropolis algorithm on unbounded domains[END_REF][START_REF] Fort | Convergence of adaptive and interacting Markov chain Monte Carlo algorithms[END_REF] and references therein).

The optimal choice for Σ is appropriate only when the target distribution is strongly unimodal [START_REF] Roberts | Optimal scaling for various Metropolis-Hastings algorithms[END_REF]. If the data y is assumed to be drawn independently from a mixture model, its likelihood is of the form

p(y|x)= i M m=1 α m f (y i |φ (m) ),
where i α i =1and α i ≥ 0, φ (m) denotes the ndimensional parameter vector of the mth component, and the parameter space is a subset of (R + × R n ) M . The likelihood p(y|x) in this case is invariant under any permutation of the mixture components. If the prior p(x) is exchangeable, that is, it also does not favor a particular permutation, then the posterior π(x) inherits the permutation invariance.

To illustrate the challenges of inference in this model, we present an example motivated by a signal processing problem of the water Cherenkov signals of the Pierre Auger Observatory [START_REF] Auger | Pierre Auger project design report[END_REF]. Figure 1(a)-1(b) display the MCMC sample corresponding to a single run of AM on an exponential mixture where the rates are known: only the location parameters and the mixture weights are estimated. A flat prior is taken. The red variable gets stuck in one of the mixture components, whereas the blue, green, and brown variables visit all the three remaining components. Marginal estimates computed for the blue, green, and brown variables are then mostly identical as seen on Figure 1(b). In addition, the shaded ellipses, depicting the marginal posterior covariances of each component's parameters, indicate that the resulting empirical covariance estimate is very broad, resulting in poor efficiency of the adaptive algorithm.

Several approaches have been proposed to deal with the label switching problem. The first solution consists in modifying the prior in order to make it select a single permutation of the variables, introducing an identifiability constraint [START_REF] Richardson | On Bayesian analysis of mixtures with an unknown number of components[END_REF]. This solution is known to cause artificial biases with respect to the posterior by not respecting its topology [START_REF] Marin | Bayesian modelling and inference on mixtures of distributions[END_REF]. An effort has then been made to adapt to the estimated posterior surface through the design of relabeling algorithms [START_REF] Stephens | Dealing with label switching in mixture models[END_REF][START_REF] Marin | Bayesian modelling and inference on mixtures of distributions[END_REF] that process the MCMC sample after the completion of the simulation run. These techniques look for a permutation P t of each individual sample point X t so as to minimize a posterior-based criterion depending on the whole chain history. [START_REF] Celeux | Bayesian inference for mixtures: The label-switching problem[END_REF] proposed an online version of the relabeling procedure in which the simulation of each X t is followed by a permutation P t of its components. The permutation P t is chosen to minimize a user-defined criterion that depends only on the past history of the chain up to time t. The major advantage of this online approach is that it is compatible with our objective of solving label switching "on the fly" in order to optimize AM for permutationinvariant models.

In both batch and online relabeling algorithms, inference is carried out by using relabeled samples. Since the permutation steps in the MCMC procedure modify the distribution of the chain (X t ), the target distribution is no longer the posterior π(x). [START_REF] Jasra | Markov Chain Monte Carlo methods and the label switching problem in Bayesian mixture modelling[END_REF] showed that, empirically, relabeling induces the learning of an appropriate identifiability constraint, but the existence of a target distribution and its relation with the original target π(x) has been an open problem, meaning that these relabeling techniques have remained mostly heuristics. [START_REF] Papastamoulis | An artificial allocations based solution to the label switching problem in Bayesian analysis of mixtures of distribution[END_REF] recently proved a convergence result for a non-adaptive, batch, identifiability constraint-based relabeling procedure, where the constraint however depends on both the unrelabeled sample and the user. Our ultimate goal is to prove a similar but sample-and user-independent result, tightening adaptivity of the MCMC procedure with adaptivity of the identifiability constraint.

An adaptive online relabeling algorithm

From now on, we consider the general case of MCMC sampling from a target distribution with density π(x) on X⊆ R d with respect to (w.r.t.) the Lebesgue measure, with d = qM. Let P be a finite group of d×d block permutation matrices, indexed by some permutations ν of {1, ..., M }, acting on X as follows: for x =(x 1 , ..., x M ) ∈ R d ,P ν x = (x ν -1 (1) , ..., x ν -1 (M ) ). As an example, let us take d =6 dimensions divided in three blocks of size 2. This would correspond to having a mixture model with M =3components and q =2parameters per component. The 6 × 6 matrix associated to the permutation ν that sends 1 onto 2, 2 onto 3 and 3 onto 1 is

P ν =   00I 2 I 2 00 0 I 2 0   ,
and for x = (1, 2, ..., 6) T , P ν x = (5, 6, 1, 2, 3, 4) T .

Let us now assume that π is invariant under the action of P, i.e. π(Px)=π(x) for any P ∈P. Our goal is to isolate a single mode out of the many identically repeated symmetric modes of the posterior. Formally, we are interested in restricting the target π to a fundamental domain D of the action of P, that is, finding a subset D ⊂ X which is minimal for the inclusion and for which {Px : x ∈D,P ∈ P} = X , up to a set of Lebesgue measure zero. Following [START_REF] Celeux | Bayesian inference for mixtures: The label-switching problem[END_REF][START_REF] Stephens | Dealing with label switching in mixture models[END_REF], we will select D so that the sample looks as Gaussian as possible, since we want to select a single mode of the symmetric modes of the target π. For this purpose, the domain D will be defined through the minimization of a Mahalanobis-type criterion.

For a d × d invertible covariance matrix Σ, let

N (x|µ, Σ) = 1 (2π) d/2 1 √ detΣ exp - 1 2 L (µ,Σ) (x)
be the density of a Gaussian distribution on R d with mean µ and covariance matrix Σ, where

L (µ,Σ) (x) = (x -µ) T Σ -1 (x -µ).

Derivation of the algorithm

Let C + d be the set of real d × d symmetric positive definite matrices, and let θ ∈ R d ×C + d . θ will later on be taken to be the concatenation of the mean and covariance of the chain (X t ), but for now, it is fixed to an arbitrary value. AMOR combines two actions: (i) sample a chain with target proportional to π D where D is a fundamental domain of the action of P and (ii) learn the domain D on the fly (see Figure 2). First assume that the adaptation is frozen, that is, consider AMOR when steps 12 and 13 in the pseudocode in Figure 2 are removed. In this case, we prove in Proposition 1 that our algorithm is a MCMC sampler with target distribution

π θ (x)=Z -1 θ π(x) V θ (x), where Z θ = V θ π(x)dx,
(1) and V θ is defined by

V θ = {x : L θ (x) = min P ∈P L θ (Px)}.
In other words, π θ (x) is the initial target π(x) restricted to the nearest-neighbor (Voronoï) cell V θ defined by the distortion measure L θ . . Define a sequence {X t ,t ≥ 0} as prescribed by Figure 2: set X 0 ∈ V θ and for t ≥ 0, (i) sample X ∼N (•|X t ,cΣ); (ii) conditionally on X, draw P uniformly over the set argmin P ∈P L θ (P X); (iii) set X t+1 = P X with probability α(X t , P X) and X t+1 = X t otherwise, where the acceptance ratio α is given by

α(x, x) = 1 ∧ π(x) π(x) Q∈P N (Qx|x, cΣ) Q∈P N (Qx|x, cΣ) . ( 2
)
Then {X t ,t ≥ 0} is a Metropolis-Hastings Markov chain with invariant distribution π θ .

Proof We first prove by induction that X t ∈ V θ for any t ≥ 0. This holds true for t =0. Assume that X t ∈ V θ . By construction, P X ∈ V θ ; therefore, X t+1 ∈ V θ .

Recall that for a permutation matrix P , P -1 = P T . Observe that given the current state X t , P X is sampled under the proposal distribution qθ (x t ,x)=

P ∈P p θ (P |P T x) N (P T x|x t ,cΣ) ,
where the conditional distribution p θ (•|X) is the uniform distribution over argmin P ∈P L θ (PX). As P is a group, p θ (•|PX)=p θ (•|X) for any P ∈P. Furthermore, for any x ∈ V θ , the support of the distribution qθ (x, •) is in V θ . This implies that for any x,

x ′ ∈ V θ , qθ (x, x ′ ) ∝ P ∈P N (Px ′ |x, cΣ) . AMOR π(x), X 0 , T , µ 0 , Σ 0 , c 1 S ← ∅ 2 for t ← 1 to T 3 Σ ← cΣ t-1 ⊲ scaled adaptive covariance 4 X ∼N •|X t-1 , Σ ⊲ proposal 5 P ∼ arg min P ∈P
L (µt-1,Σt-1) P X ⊲ pick an optimal permutation 6 X ← P X⊲ permute

7 if π(X) P N PX t-1 |X, Σ π(X t-1 ) P N PX|X t-1 , Σ > U [0, 1] then 8 X t ← X⊲ accept 9 else 10 X t ← X t-1 ⊲ reject 11 S ← S ∪ {X t } ⊲ update posterior sample 12 µ t ← µ t-1 + 1 t X t -µ t-1
⊲ update running mean and covariance

13 Σ t ← Σ t-1 + 1 t ((X t -µ t-1 )(X t -µ t-1 ) ⊺ -Σ t-1 ) 14 return S Figure 2:
The pseudocode of the AMOR algorithm. The steps of the classical SRWM algorithm are in blue, the adaptive MH algorithm adds the green steps, and the new online relabeling steps are in red. Notice the adaptation of both the proposal (line 4) and the selection mechanism through the dependence of L (µ,Σ) on (µ, Σ). Note that for practical reasons, a small εI d is often added to the covariance matrix in line 3, but [START_REF] Vihola | Can the adaptive Metropolis algorithm collapse without the covariance lower bound?[END_REF] recently confirmed that core AM does not lead to degenerate covariances. Note also that line 5 is usually a simple assignment of the optimal permutation. In case of ties, we draw uniformly from the finite set of optimal permutations. Therefore, for any x, x ′ ∈ V θ ,

α(x, x ′ )=1 ∧ π θ (x) π θ (x) P ∈P N (Px|x ′ ,cΣ) P ∈P N (Px ′ |x, cΣ) =1 ∧ π θ (x ′ ) π θ (x) qθ (x ′ ,x) qθ (x, x ′ ) .
V θ is a fundamental domain of the action of P on X . It contains only one copy of each genuine mode of π and it is sufficient to consider the restriction π θ of π to V θ : ∀x ∈ X , ∃P ∈P and y ∈ V θ s.t. x = Py and π(x)=π(y) (by the permutation invariance of π). The sets (PV θ ) P ∈P , where PV θ = {Px : x ∈ V θ }, cover X . If θ is further taken such that no P exists with µ = Pµ or P ΣP T =Σ, then the sets (PV θ ) P ∈P are pairwise disjoint. Proofs of these claims can be found in the supplementary material.

To sum up, we have a family of cells (PV θ ) P ∈P that are permutations of each other, and where each cell PV θ is the support of a permuted copy of the landscape of π (Figure 3). In view of the definition of V θ and L θ , forcing the chain to stay in V θ means that we want to make the sample look "as unimodal as possible".

The second step now is to find a convenient θ in such a way that MH is optimized: based on [START_REF] Roberts | Optimal scaling for various Metropolis-Hastings algorithms[END_REF], we want to choose the covariance matrix Σ such that Σ is proportional to the covariance matrix of π θ . This implies that θ =( µ, Σ) solves the fixed point equations

µ = xπ θ (x)dx, Σ= (x -µ)(x -µ) T π θ (x)dx.
(3) To achieve this goal, the MH-online relabeling (steps 3 to 11 in Figure 2) is combined with the adaptive steps 12 and 13. AMOR can thus be described as follows: let {K θ ,θ ∈ R d ×C + d } be the family of MH kernels described by Proposition 1. The posterior sample S in AMOR is obtained by running adaptive MCMC: given an initial value θ 0 =( µ 0 , Σ 0 ), define by induction {(X t ,θ t ),t ≥ 0} as follows:

• Sample X t ∼ K θt-1 (X t-1 , •),
• Update the parameter:

θ t = θ t-1 + 1 t H(X t ,θ t-1 ), (4)
where

H(x, θ) = (x -µ, (x -µ)(x -µ) T -Σ). (5)
Recent results on convergence of adaptive MCMC show that when this algorithm converges, there exists θ ⋆ ∈ R d ×C + d such that the distribution of the posterior sample {X t ,t ≥ 0} converges to π θ⋆ and n -1 n t=1 f (X t )π θ⋆ (f ) converges a.s. to zero for a wide family of functions f [START_REF] Fort | Convergence of adaptive and interacting Markov chain Monte Carlo algorithms[END_REF]. Therefore, the empirical mean and covariance of {X t ,t ≥ 0} converges to the expectation and covariance matrix of π θ⋆ . Hence θ ⋆ solves the fixed point equations [START_REF] Roberts | Optimal scaling for various Metropolis-Hastings algorithms[END_REF]. In Section 3.2 we will discuss sufficient conditions for the convergence of AMOR.

Allowing only scalar matrices (Σ t = λ t I) for selection and using fixed-covariance proposals, our relabeling mechanism coincides with the one of [START_REF] Celeux | Bayesian inference for mixtures: The label-switching problem[END_REF]. Beside the obvious generalization to full covariance matrices and the adaptivity of the proposal, the breakthrough of our algorithm lies in its modification of the acceptance ratio that allows to identify the limiting distribution of the sample {X t ,t ≥ 0} and to prove limit theorems (e.g., ergodicity and law of large numbers).

Convergence analysis

As discussed in section 3.1, AMOR is an adaptive MCMC algorithm such that for any θ ∈ R d ×C + d , K θ has its own invariant distribution π θ . Sufficient conditions for the convergence of such algorithms have been recently derived [START_REF] Fort | Convergence of adaptive and interacting Markov chain Monte Carlo algorithms[END_REF]. Three main conditions are required. The first condition is that adaptation has to be diminishing; this means that some suitable distance between two consecutive kernels K θt-1 and K θt vanishes. When adaptation relies on stochastic approximation as in the case of AMOR, this condition is easy to check [START_REF] Fort | Convergence of adaptive and interacting Markov chain Monte Carlo algorithms[END_REF][START_REF] Saksman | On the ergodicity of the adaptive Metropolis algorithm on unbounded domains[END_REF]. The second condition is the containment condition: roughly speaking, a kind of uniformin-θ ergodic behavior of the kernels

{K θ ,θ ∈ R d ×C + d } is required.
As shown in [START_REF] Fort | Convergence of adaptive and interacting Markov chain Monte Carlo algorithms[END_REF][START_REF] Saksman | On the ergodicity of the adaptive Metropolis algorithm on unbounded domains[END_REF], this property is related to the stability of the random sequence {θ t ,t ≥ 0}. The third condition is to control the weak convergence of the random invariant distributions {π θt ,t ≥ 0} to the limiting target π θ⋆ . In the case of AMOR, this amounts to control both the stability and the almost sure convergence of the sequence {θ t ,t ≥ 0}.

Therefore, a key ingredient for the proof of the convergence of AMOR is to show the stability and the convergence of the stochastic approximation sequence {θ t ,t ≥ 0}. Existing results cannot apply since the draws are neither i.i.d. nor Markovian. However, since the draws X t are obtained by applying a new kernel at each iteration, our idea is to adapt the existing results for the Markovian case. While the full proof of convergence of AMOR is out of the scope of this paper, we address here one of its most important aspects which unveils a novel relation between relabeling and vector quantization.

The main tool for the stability and convergence of stochastic approximation algorithms is the existence of a Lyapunov function [START_REF] Borkar | Stochastic Approximation: A Dynamical Systems Viewpoint[END_REF]. The space R d ×C + d is endowed with the scalar product θ 1 ,θ 4) converges to F [START_REF] Andrieu | Stability of stochastic approximation under verifiable conditions[END_REF]. Proposition 2 shows that the function w given by w(θ)=-log N (x|θ) π θ (x)dx is a natural candidate for the Lyapunov function. It is also proven that F is the set of points (µ, Σ) satisfying the fixed point equation (3). Proposition 2. Let us consider the AMOR algorithm (see Figure 2) with quadratic loss

L θ (x) = (x -µ) T Σ -1 (x -µ).

Define the mean field

h on R d ×C + d by h(θ)=E [H(X, θ)] = µ π θ -µ, (µ π θ -µ)(µ π θ -µ) T +Σ π θ -Σ) ,
and let

Θ= {θ ∈ R d ×C + d : ∀P ∈P, Σ = P ΣP T or µ = Pµ}.
Then ∇w(θ),h(θ) ≤0 for any θ ∈ Θ and ∇w(θ),h(θ) =0iff θ =( µ, Σ) solves the fixed point equation [START_REF] Roberts | Optimal scaling for various Metropolis-Hastings algorithms[END_REF].

The proof of Proposition 2 is given in the supplementary material (Corollary 1). We also show that for any θ ∈ Θ, w(θ) is, up to an additive constant, the Kullback-Leibler divergence between N (•|θ) and π θ . Equivalently w(θ) is the Kullback-Leibler divergence between the original unrelabeled posterior π and a mixture of Gaussians, whose component parameters are images of each other through the symmetries that leave π invariant. Finally, Proposition 3 establishes that when θ ∈ Θ, w(θ) is a distortion measure in vector quantization [START_REF] Graf | Foundations of Quantization for Probability Distributions[END_REF], which is the key property for the proof of Proposition 2, and unveils a novel link between clustering and relabeling. Notice that w(θ) is exactly a distortion measure, as the first term of the right-hand side of Proposition 3 is constant over the set {P ΣP T ,P ∈ P}.

Experiments

We benchmarked the AMOR algorithm on two Bayesian inference tasks. The first one is the problem of estimating the nine parameters ψ =( α i ,µ i ,σ i ) i=1,2,3 of a mixture of three one-dimensional Gaussians 3 i=1 α i N (.|µ i ,σ i ), taking wide flat priors over each parameter. We compared four algorithms: (i) an SRWM with an ordering constraint on the three means µ 1 ≤ µ 2 ≤ µ 3 , (ii) the original online relabeling algorithm of [START_REF] Celeux | Bayesian inference for mixtures: The label-switching problem[END_REF], (iii) the same online algorithm with modified acceptance ratio according to (2) (henceforth denoted as MC for Modified Celeux), and (iv) the AMOR algorithm. To quantify the performance after T iterations, we first selected the permutation of the running posterior mean components (μ (T ) i ) i=1,2,3 which minimized the sum of the ℓ 2 errors on the three estimates of the means µ i ,i = 1, 2, 3, and we considered the latter sum taken at this best permutation of the posterior mean:

S T = arg min τ ∈S3 3 i=1 (μ (T ) τ (i) -µ i ) 2 .
We repeated this experiment 100 times on 100 different datasets coming from parameters generated as follows: draw (α i ) ∼D(1), µ i ∼U (0,1) i.i.d., and σ i ∼U (0,0.05) i.i.d. This choice of generative distribution ensured a reasonable number of datasets containing overlapping Gaussians, thus provoking switching. Figure 5(a) depicts the performance measure S T averaged over the 100 datasets of this 9D experiment, versus time T . We use this averaging as a way to estimate the expected performance measure on a class of problems given by the generative distribution. AMOR significantly outperforms other approaches as it converges faster on average and to a better solution. As expected, imposing an ordering constraint on the means (RWM+OC) reveals a poor strategy leading to artificial additional bias. Note finally that the modification of the acceptance ratio did not increase drastically the performance of the Online Relabeling algorithm of [START_REF] Celeux | Bayesian inference for mixtures: The label-switching problem[END_REF] ("Original RWM+OR" vs. "Modified RWM+OR"), which is not a surprise since the additional factor in the ratio (2) is often close to 1. Figure 4 provides insight into how the two best methods (AMOR and MC) behaved after T 1 =1 K and T 2 = 30K iterations, presenting scatter plots of performances S 1,000 and S 30,000 . Each point corresponds to one of the 100 datasets. Clearly, starting from a rather random distribution of the errors, AMOR took the advantage after 30k iterations, while a few cases were still better treated by MC.

We also performed a higher dimensional experiment to further investigate the comparison between AMOR and CC. This time, the goal was to estimate the three means of a 10dimensional Gaussian mixture 3 i=1 1/3 N (.|µ i , 0.1I 10 ). Again, 100 datasets of 100 points each were generated with µ i ∼U (0,1) i.i.d. Again, as seen Figure 5(b), AMOR stabilizes earlier and selects a better region of R 30 than MC, thus illustrating the interest of combining adaptive selection and proposal mechanisms.

Conclusion

We derived AMOR, a novel MCMC algorithm for simulating from a permutation invariant target. AMOR combines Adaptive Metropolis with Online Relabeling. We justified its use by showing that it relies on a stochastic approximation procedure for which we exhibited a Lyapunov function. Experiments show that AMOR outperforms existing approaches for Bayesian inference in Gaussian mixture models. Besides, the theoretical framework of its analysis is appealing because (i) it generalizes previous approaches [5], (ii) it paves the way towards future work on the asymptotic behavior of relabeling algorithms and convergence of the samples {X t ,t ≥ 0}, and (iii) the Lyapunov function we derived exhibits an elegant relation with vector quantization techniques.

However, these features come at a price: future work should try to remove the need to sweep over all permutations in P, which is prohibitive when using AMOR with large |P|. Future work could also study variants of AMOR, using, e.g., a soft clustering selection mechanism, replacing the indicator appearing in the definition of π θ by a logistic indicator, or equivalently relabel each new Gaussian sample in AMOR by sampling from a multinomial distribution over permutations conditionally on the history statistics, a technique similar in spirit to probabilistic relabeling algorithms developed in [START_REF] Jasra | Bayesian inference for mixture models via Monte Carlo[END_REF] and [START_REF] Sperrin | Probabilistic relabelling strategies for the label switching problem in Bayesian mixture models[END_REF].

Figure 1 :Proposition 1 .

 11 Figure 1: The results of AM (top row) and AMOR (bottom row) algorithms on a mixture of exponential components.The right panels show the parameters of the four components (black dots: the x-coordinates are the location parameters φ (m) and the y-coordinates are the mixture weights α (m) ), the (unnormalized) mixture distribution (blue curve), and the marginal posteriors of the four location parameters (colored histograms). Colored ellipses are exp(1/2)-level sets of Gaussian distributions: the means are the Bayesian estimates for the location and weight parameters of each component, and the covariance is the marginal posterior covariance of each location/weight couple. The left panels show the four chains of the location parameters φ (m) (light colors), the running means (dark colors), and the mean of the running means (black curve). The AM algorithm shows heavy label switching among the three rightmost components whereas the AMOR algorithm separates the components nicely.

Figure 3 :

 3 Figure 3: Examples of tessellations (PV θ ) P ∈P . (a) The plane is always cut into two parts in a "butterfly" shape. The central intersection point can move along the first diagonal, eventually until ∞. (b) Borders between two cells in 3D are quadratic hypersurfaces.

Proposition 3 .L

 3 For any θ ∈ Θ, (P µ,P ΣP T ) (x) π(x)dx.

Figure 4 :

 4 Figure 4: Experimental comparison of the performance measure S for AMOR (x-axis) versus Modified Celeux (y-axis) in the 9D experiment.

Figure 5 :

 5 Figure 5: Experimental comparison of several relabeling approaches. The plots show the performance measure S T vs. T averaged over 100 datasets drawn from a common generative model.

  2 = µ T 1 µ 2 + Trace(Σ 1 Σ 2 ).A continuously differentiable function w : R d ×C + d → R + is a Lyapunov function for the mean field h if (i) for any θ ∈ R d ×C + d , ∇w(θ),h(θ) ≤0;(ii) the level sets {w ≤ C}, C>0, are compact subsets of R d ×C + d . When w(F) has empty interior, where F = {θ ∈ R d ×C + d : ∇w(θ),h(θ) =0 }, then the sequence {θ t ,t ≥ 0} defined by (

Acknowledgements

This work was supported by the ANR-2010-COSI-002 grant of the French National Research Agency.


Thomas Piton
email: thomas.piton@univ-nantes.fr

Julien Blanchard
email: julien.blanchard@univ-nantes.fr

Fabrice Guillet
email: fabrice.guillet@univ-nantes.fr

Keywords: customer data mining, data mining for marketing, recommender system, actionable knowledge discovery, industrial case-based application

Recommender systems can apply knowledge discovery techniques to the problem of making product recommendations. This aims to establish a customer loyalty strategy and thus to optimize the customer life time value. In this paper we propose CAPRE, a data-mining based methodology for recommender systems based on the analysis of turnover for customers of specific products. Contrary to classical recommender systems, CAPRE does not aspire to predict a customer's behavior but to influence that behavior. By measuring the actionability and profitability of customers, we have the ability to focus on customers that can afford to spend larger sums of money in the target business. CAPRE aggregates rules to extract characteristic purchasing behaviors, and then analyzes the counter-examples to detect the most actionable and profitable customers. We measure the effectiveness of CAPRE by performing a cross-validation on the MovieLens benchmark. The methodology is applied to over 10,000 individual customers and 100,000 products for the customer relationship management of VM Matériaux company, thus assisting the salespersons' objective to increase the customer value.

I. INTRODUCTION

In a world where the number of product and service offers can be overwhelming, Recommender Systems (RS) help customers find and evaluate products or services. RS are software tools and techniques providing suggestions for products to be use by a customer [START_REF] Mahmood | Improving Recommender Systems with Adaptive Conversational Strategies[END_REF] [START_REF] Resnick | Recommender Systems[END_REF]. They connect customers with potential products to purchase by associating either the content of products or the opinions of customers, or both. In recent years, RS have found numerous applications in industry, especially in e-commerce. They help to develop and maintain profitable relationships with customers [START_REF] Ren | An Evolving Information System Based on Data Mining Knowledge to Support Customer Relationship Management[END_REF], for example through cross or up-selling.

With most RS, the recommendations depend on the products that are currently bought by the customer [START_REF] Adomavicius | Toward the Next Generation of Recommender Systems: A survey of the Stateof-the-Art and Possible Extensions[END_REF]. This principle is perfectly appropriate to e-commerce websites, where the customer is alone for driving the product exploration. However, it is hardly adaptable to salespersons' visits to customers. We propose to adopt a different strategy: when salespersons visit customers, which product(s) can they recommend to trigger purchasing? Contrary to classical RS, our goal is not to predict customer behaviors but to influence customer behaviors. Making such recommendations requires (i) to identify characteristic purchasing behaviors, and (ii) to apply them on customers while respecting their own behaviors. This last point is crucial because the cost of an inappropriate recommendation is higher in the case of salespersons' visits than on e-commerce websites [START_REF] Linden | Amazon.com Recommendations: Item-to-Item Collaborative Filtering[END_REF]. The salespeople may even refuse to use the RS if they consider the recommendations not sufficiently relevant [START_REF] Cosley | Is Seeing Believing?: How Recommender System Interfaces Affect Users' Opinions[END_REF].

In this paper, we propose CAPRE (Customer Actionability and Profitability Recommendation), a new methodology for recommender systems based on the analysis of turnover for customers of specific products. More precisely, CAPRE aggregates rules to extract characteristic purchasing behaviors, and then analyzes the counter-examples to detect the most actionable and profitable customers. Recommendations are made by targeting the actionable counter-examples with the most profitable rules. Actually, we do not use rules as naive implications of purchasing behaviors. We use the rules as a model to detect counter-examples, i.e. customers with deviant purchasing behaviors that have the potential to spend larger sums of money in the target business. Our methodology is applied over 10,000 customers and 100,000 products for the customer relationship management of VM Matériaux.

A. Contributions of the Paper

The main contributions of this paper can be resumed as follows:

• We develop a new methodology for recommending products from a turnover database. • We propose a measure of actionability of the recommendations, based on the similarity between the examples and counter-examples. • We offer an original measure of the profitability of the recommendations, based on criteria defined by experts. • We apply CAPRE on an actual VM Matériaux trading group database composed of more than 10,000 customers and 100,000 products. • We measure the effectiveness of our system on the MovieLens benchmark with a cross-validation.

B. Organization

The rest of the paper is organized as follows: section II presents a state of the art in the recommender systems and rule mining fields. Then, section III discusses our methodology CAPRE with the two measures of actionability and profitability. We apply CAPRE in order to make recommendations on an actual dataset of VM Matériaux trading group 1 in section IV. We analyze the effectiveness of CAPRE in section V by performing a cross-validation on the MovieLens benchmark. The final section concludes and describes our future work.

II. RELATED WORK

The growth of Internet and the emergence of e-commerce has led the development of RS [START_REF] Resnick | GroupLens: an Open Architecture for Collaborative Filtering of Netnews[END_REF]. These allow companies to filter information, and recommend products to customers according to their preferences. Recommend products and services can enhance the relationship between customer and seller, and thus increase profits [START_REF] Zhang | An Associative Classification-Based Recommendation System for Personalization in B2C e-commerce Applications[END_REF]. Actually recommendation is not a new concept, it is already used on the Web through movies, books and music. RS have become an important research area since the appearance of the first papers on collaborative filtering [START_REF] Resnick | GroupLens: an Open Architecture for Collaborative Filtering of Netnews[END_REF]. RS are doubly useful: on one hand, they help customers filter through enormous numbers of products, and focus on the few ones that match their preferences; on the other hand, recommender systems help companies increase their sales on the long tail.

Many techniques have been developed, and several classifications have been proposed [START_REF] Resnick | Recommender Systems[END_REF]. Montaner et al. [START_REF] Montaner | A Taxonomy of Recommender Agents on the Internet[END_REF] provide a taxonomy and classify existing RS applications to specific application domains: entertainment (movies, music and IPTV), content (personalized newspapers, Web pages, elearning applications and e-mail filters), e-commerce (products to purchase such as books, cameras, PCs, etc.) and services (travel services, experts for consultation, houses to rent, matchmaking services). Examples of such applications include recommending books at Amazon.com, movies by MovieLens or music by Yahoo! Music. Schafer et al. [START_REF] Schafer | E-Commerce Recommendation Applications[END_REF] present a detailed classification of recommender systems in e-commerce, and elucidate how they can be used to provide personalized service customer loyalty. RS are usually classified into the following categories:

• Content-based recommendations: the customer will be recommended products similar to the ones the customers preferred in the past; • Collaborative recommendations: the customer will be recommended products that customers with similar preferences and behaviors liked in the past; • Hybrid approaches: this approach combine collaborative and content-based methods. In the following subsections, we will focus on collaborative filtering and rule mining for RS.

A. Recommender Systems Based on Collaborative Filtering

Collaborative recommender systems [START_REF] Resnick | GroupLens: an Open Architecture for Collaborative Filtering of Netnews[END_REF] produce customer specific recommendations of products based on patterns of ratings or purchases without need for exogenous information about either products or customers. In other words, the main idea is to exploit information about the past behavior of an existing customer for predicting which products the current customer will most probably like or be interested in. These types of systems are in widespread industrial use today, in particular as a tool in online retail sites.

There have been many collaborative systems developed in the academia and the industry. It can be argued that the Grundy system [START_REF] Rich | Readings in Intelligent User Interfaces[END_REF] was the first recommender system, which proposed using stereotypes as a mechanism for building models. Later on, the Tapestry system relied on each customer to identify like-minded customers manually [START_REF] Goldberg | Using Collaborative Filtering to Weave an Information Tapestry[END_REF]. GroupLens, Video Recommender and Ringo use collaborative filtering algorithms to automate prediction. In the recommendation field, the difficulty to collect descriptions customers and content is a recurrent problem.

B. Rule Mining for Recommender Systems

In the context of RS, a transaction could be viewed as the set of all the previous ratings or purchases of a customer. Association rules can be used to develop Top-N recommendations for each customer depending on past purchases. Frequently, products are sorted by the confidence of the rules [START_REF] Sarwar | Analysis of Recommendation Algorithms for E-Commerce[END_REF], using multi-level rules [START_REF] Kim | A Recommendation Algorithm Using Multi-Level Association Rules[END_REF] or depending on the antecedent length of rules [START_REF] Nakagawa | Impact of Site Characteristics on Recommendation Models Based On Association Rules and Sequential Patterns[END_REF].

Recently, rule mining has become a current application of RS. Leung et al. [START_REF] -K. Leung | A Collaborative Filtering Framework Based on Fuzzy Association Rules and Multiple-Level Similarity[END_REF] present a model based on fuzzy rules by exploiting similarities in product taxonomies. Liu et al. [START_REF] Liu | A Hybrid of Sequential Rules and Collaborative Filtering for Product Recommendation[END_REF] propose to extract customer purchasing behaviors using sequential rules. Finally, Smyth et al. [START_REF] Smyth | Case Studies in Association Rule Mining for Recommender Systems[END_REF] present two different case studies using association rules for RS.

III. CAPRE METHODOLOGY

In this section, we propose CAPRE, a new methodology to apply rules to recommendation systems. A real-world application is presented in section IV. CAPRE is divided into five steps (see figure 1).

A. Purchasing Behavior Rule Mining

1) Data: Let C be the set of all customers {c 1 , c 2 , . . . c n } and P be the set of all products {p 1 , p 2 , . . . p m } that can be recommended. Let u be a utility function that measures the usefulness of product p for customer c:

u : P × C → R p, c → u(p, c) (1)
where R is a part of R, e.g. positive integers or numbers.

In our methodology, u concerns a target business and has a financial meaning. For instance, it can be the turnover of product p for customer c. The greater u(p, c), the more the purchase of p by c is profitable for the company. The function u depends on customer transaction history and can be materialized by a matrix customers × products. This context is directly applicable in the usual framework of rating datasets (see section V). Then, we apply two transformations (see figure 2) on the values of the utility function, to make the customers commensurable and to extract rule behaviors on discretized data.

u(p, c) u % (p, c) u * (p, c)

Utility fonction Normalization Discretization Figure 2. Normalization and discretization of the utility function

• The first step is to make the turnover of a customer c comparable with another, regardless of their financial potential. For each customer c, the sales are transformed as a percentage of sales, representing the usefulness ratio of a product p to a customer c.

u % : P × C → [0; 100 %] p, c → u % (p, c) (2)
with,

u % (p, c) = u(p,c) p∈P u(p,c)
• The second step is to normalize the proportions u % (p, c), that we denote u * (p, c), by taking into account their distribution for each product. This normalization results in a discretization of proportions in several intervals

µ 1 < µ 2 < . . . < µ t . u * : P × C → U = {µ 1 , µ 2 , . . . µ t } p, c → u * (p, c) (3)
The normalization and discretization depend on the application. Each customer c ∈ C is described by various characteristic variables (numerical or categorical), such as age, income, etc. For example, the utility function u is the turnover of customers for product families such as plaster, concrete, etc. u * (p, c) = µ 1 (respectively µ 2 and µ 3 ), means that the turnover of customer c for product p is low (respectively medium and large). Each building contractor customer c is described by characteristic variables composed of the number of employees, the professional category, etc.

2) Mining Cohorts of Rules (step 1 on figure 1):

We apply an association rule mining algorithm [START_REF] Han | Frequent Pattern Mining: Current Status and Future Directions[END_REF] to extract purchasing behavior rules of the form:

u * (p X , ) = µ X , u * (p Y , ) = µ Y , . . . → u * (p Z , ) = µ Z ( 4
) where ∀ i p i ∈ P and ∀ j µ j ∈ U . The examples of the rule are the customers who satisfy the antecedent and the consequent (corresponding to the dots " " in ( 4)). By contrast, counter-example customers satisfy the antecedent and not the consequent 2 . Such a rule means that the customers who satisfy the antecedent tend to satisfy the consequent. In our methodology, as our goal is to make recommendations, we use rules with only one product in the consequent. For simplicity, we will denote u * (p X , c) = µ X by p X = µ X (where p X ∈ P and µ X ∈ U ). So, we denote the rules in (4) as follows:

p X = µ X , p Y = µ Y . . . → p Z = µ Z (5)
In the classical association rule terminology, p X = µ X is called an item. The set of rules is filtered by support and confidence thresholds. The resulting ruleset is denoted R. Given an item p X = µ X , a cohort is the set of all the rules in R concluding on this item: The interest of the cohorts is to merge all the rules which provide the same recommendation. By exploiting rule redundancy, we are able to extract more robust purchasing behavior models (with more examples). Given a cohort Ct, a Ct-Recommendation is to propose the consequent of Ct to a counter-example.

Ct(p X = µ X ) = {r ∈ R | consequent(r) = (p X = µ X )} (6)
For example, let us consider the cohort Ct(Gravel = c) composed of two rules: The Ct-recommendation consists in proposing some gravel to customers who are counter-examples of at least one of the rules.

B. Actionability of Counter-Examples

Let us now consider a cohort. With their partially unexploited purchase potential, the counter-examples can theoretically help the company to develop its turnover. However, all the counter-examples cannot be equally receptive to our recommendations. This is why our methodology focuses on the most actionable counter-examples, i.e. the counterexamples that are similar to examples according to the customer characteristic variables (age, gender, income, etc.). The more a counter-example is similar to the examples, the more it is likely to behave like an example, i.e. to purchase more of the consequent of the cohort.

1) Pre-actionability on Purchasing Variables (step 2 on figure 1): To be close to the examples Ct + of a cohort, a counter-example e -should not display an extreme purchasing behavior on the antecedents of the rules of the cohort Ct. Given the purchases of example distribution on the antecedents of the cohort, we define the median M and (Q 3 -Q 1 ) the interquartile range of the distribution. A counter-example e -∈ Ct -is pre-actionable if and only if:

M -(Q 3 -Q 1 ) ≤ p ∈ antecedents of the cohort Ct u % (p, e -) ≤ M + (Q 3 -Q 1 ) (7)
For example, given the cohort Ct(Gravel = c), the purchases of examples on the antecedents of the cohort are characterized by M = 40 % and 3 -Q 1 = 20 %. To be pre-actionable, a counter-example e -∈ Ct -must have spent between 20 % and 60 % of its turnover on cement, sand and concrete mixer.

2) Actionability on Characteristic Variables (step 3 on figure 1): In the multidimensional space of the customer characteristics, we need a distance measure to assess similarities between customers. As customer characteristics are numerical and categorical variables, we consider the framework of Factor Analysis of Mixed Data (FAMD) introduced by Escofier in 1979 (see [START_REF] Bécue-Bertaut | Multiple Factor Analysis and Clustering of a Mixture of Quantitative, Categorical and Frequency Data[END_REF] and [START_REF] Saporta | Simultaneous Analysis of Qualitative and Quantitative Data[END_REF] for further details). This framework provides the distance measure for mixed data that we use in CAPRE. The customer characteristics consist of numerical and categorical variables:

• Let N be the set of centered and reduced numerical variables.

• By applying a complete disjunctive coding on categorical variables, we generate a set I of indicators 3 (a binary variable by modality). The FAMD distance between customers c 1 and c 2 is defined by: 3 The values of these indicators are divided by the square root of frequency of each indicator.

d 2 (c 1 , c 2 ) = v i=1 (N i,c1 -N i,c2 ) 2 + w j=1 (I j,c1 -I j,c2 ) 2 (8)
For example, roofers and tilers purchase many common products. They can be close by their purchases (see section III-B1) but remote by their characteristic variables. For the rule Cement Bag → Tile, most examples are tilers. However, roofers tend to regularly purchase 25kg bags of cement to sink a fork or install a ridge roof. Even if the counter-examples are tilers and roofers, you have to recommend purchasing tiles only to tilers. Computing actionability allows to avoid recommend purchasing tiles to roofers.

To sum-up, a counter-example e -is actionable if it is close enough to the examples. Given two thresholds minDist and δ, e -is considered as δ-actionable for a cohort Ct if and only if:

1) e -is pre-actionable (see section III-B1);

2) e -respects the following inequality:

{e + ∈ Ct + | d(e -, e + ) minDist} Ct + δ (9)
with minDist a threshold defining a neighborhood around the counter-example e -, and δ a threshold defining a minimum ratio of neighbor examples to be actionable. Given the distribution of the examples in the space of customer characteristics, the distance from a counterexample to the barycenter of the examples is not necessarily representative. That is why we do not use the barycenter in our measure.

C. Profitability of the Cohorts (step 4 on figure 1)

The profitability of a counter-example corresponds to the value in Dollars of the Ct-Recommendation. For example, for each actionable counter-example (see section III-B), we compute the total turnover in Dollars that should have been spent if the counter-example had behaved as an example of the cohort. We denote by u( , c) the sum p∈P u(p, c). Given a threshold θ ∈ R + , a counter-example e -is defined as θ-profitable for the product p Z , i.e. for the recommendation Ct(p Z ) if and only if:

P rof it(Ct(p Z ), e -) = (u( , e -) * α) -u(p Z , e -) θ with the mean purchase ratio α = c ∈ Ct + u(p Z , c) c ∈ Ct + u( , c) (10
) Therefore, the profitability of a cohort of rules corresponds to the sum of the profitabilities of the θ-profitable counterexamples of the cohort. From this value, we deduce fixed and variable costs incurred by a business approach. Moreover, we suggest using a scoring to refine the profitability expected for each customer. The learning phase is guided by the examples of the cohort and the application phase by the profitable counter-examples.

D. Triggering the Most Actionable and Profitable Cohorts (step 5 on figure 1) 1) Marketing Recommendation Campaigns: Several scenarios exist for decision makers to profitably action the recommender system:

• Sorting cohorts in order to identify the most profitable Ct-Recommendations. Therefore, commercial actions can be triggered to promote products. • For each customer, sorting the cohorts in which the customer is a θ-profitable counter-example. Salespersons can trigger to one or more recommendations, creating a personalized relationship with the target customer. 2) Return On Investment: The methodology presented above allows us to propose an equation to calculate the expected ROI before triggering the most actionable and profitable cohorts, assuming that for each θ-profitable counterexample e -, the profitability P rof it(Ct(p), e -) is computed for the most profitable cohorts.

ROI = c∈C p∈P P rof it(Ct(p), c) -ω with ω = α + β -γ (11)
• α: operation fixed costs (salespersons' visits, communication, gifts, meals, advertising, etc.). • β: data mining costs: time for preprocessing, modeling and reporting. • γ: cost reduction due to time spent gain in the marketing department to establish recommendations.

IV. APPLICATION ON REAL DATA

A. Context

The VM Matériaux trading group offers more than 100,000 product references for professional building contractors whose have different characteristic variables and often atypical purchasing behaviors. Therefore, the use of our methodology is appropriate to profitably recommend products for target customers, creating a personalized relationship and assisting the salespersons' objective to increase the customer value.

B. Purchasing Behavior Rule Mining 1) Data:

To address the problem of over-specialization, we use the product taxonomy of VM Matériaux. The advantage of this approach is to obtain relevant rules, in which the most common products do not hide the less common products by their frequency. For example, the rarest products, such as photovoltaic panels, will be codified in a more general level, while the most common products, such as plaster boards, will be codified at a finer level. We decided to collect data from our datawarehouse composed of approximately 91 tera-bytes and 180 tables. 9,575 customers having achieved a turnover threshold (defined by decision makers) are selected and described through variables divided into three groups:

• 20 internal variables (in the datawarehouse): address, loyalty, main salesperson, authorized outstanding, etc. • 9 external variables (obtained from a Coface4 file): number of employees, professional category, etc. • 414 aggregates (computed from levels of the taxonomy): sum of turnover for each product family (cement, insulation, etc.) within the year before a datetime. This data preprocessing produces a model with 9,575 observations and 443 variables. In agreement with business requirements, the purchases are made commensurable and are discretized (see section III-A1) into three intervals a, b and c of the same frequency 1 3 , to segment the small, medium and large customer purchases, respectively.

2

) Mining Cohorts of Rules:

The extraction of association rules is performed using the algorithm CHARM based on indepth research of closed itemsets [START_REF] Zaki | CHARM: An Efficient Algorithm for Closed Itemset Mining[END_REF] with support threshold of 0.005 (i.e. about 50 customers) and a confidence threshold of 75 %. A set of 7,788 rules is extracted (see table I) and partitioned into 57 cohorts of rules.

C. Actionability of Counter-Examples

To select the actionable counter-examples of Ct(Plaster = c), we compute the parameter minDist = 28.43 which is the average distance between all customers (see section III-B2) and we define with the commercial direction (trading director, purchasing director and sales force team) the threshold δ = 10 %. Therefore, a counter-example will be actionable if and only if:

• It is pre-actionable, i.e. its purchasing behavior for the antecedents of the rules of the cohort are in the interquartile range defined in section III-B1. • It has 10 % of the examples Ct + of the cohort in its neighborhood at a maximal distance of 28.43. To apply these two steps, let us visualize on figure 3 Secondly, to be actionable, the 46 counter-examples must have 10 % of Ct + in their neighborhood at a distance less than or equal to 28.43 (see section III-B2). This second filter allows to prune 4 counter-examples. In the end, we select 42 actionable counter-example customers.

D. Profitability of the Cohorts

With the help of decision makers of VM Matériaux, we considered that the solicitation of a customer by a salesman were approximately $2 000. That is why, to run at a healthy profit, decision makers decide to set the threshold θ = $2, 000 from which a counter-example is profitable for a cohort. When applying [START_REF] Schafer | E-Commerce Recommendation Applications[END_REF] For the 38 profitable counter-examples, the expected profit for the cohort is estimated at about $420,000 of additional turnover, i.e. over $11,000 per customer. To refine the profitability expected for each customer, we decide to use a scoring with the data mining software program KXEN. To predict the continuous target of turnover per customer for the consequent plaster = c, we choose to calculate a tendency score with a ridge regression. The learning phase is guided by the 150 examples of the cohort and the application phase by the 38 $2000-profitable counterexamples for the purchase of plaster. The learning dataset consists of a few customers. However, they represent a very homogeneous population, i.e. with a strong purchasing potential on the consequent of the cohort. Therefore, the accuracy of the model rises. For the 38 profitable counterexamples, the expected profit for the cohort is estimated at about $264,373 of additional turnover, i.e. over $7,000 per customer.

E. Triggering the Most Actionable and Profitable Cohorts

Our methodology has enabled the sales team to analyze the behavior of their customers in order to understand why some of them do not purchase certain products (price?, competition?, etc.). Therefore, other marketing measures will be implemented:

1) To promote products on targeted customers (free samples) by selecting the most profitable cohorts. 2) To target customers in a region or a store and recommend some products based on cohorts in which these customers are the most profitable. Finally, based on profitable recommendations, supplier partnerships may be agreed on to highlight products.

V. EXPERIMENTAL VALIDATION WITH MovieLens

In this section, we present the experimental results on the movie-rating benchmark. MovieLens benchmark was collected by GroupLens Research Project at the University of Minnesota and contains 100,000 movie-rating data from 943 users on 1,682 movies. Each user has rated at least 20 movies from 1 (strongly unsatisfactory) to 5 (strongly satisfactory). The dataset was converted into a binary usermovie matrix C × P that has 943 customers and 1,682 movies that were rated by at least one of the users. We apply our methodology on the five usual learning datasets named U [START_REF] Mahmood | Improving Recommender Systems with Adaptive Conversational Strategies[END_REF][START_REF] Resnick | Recommender Systems[END_REF][START_REF] Ren | An Evolving Information System Based on Data Mining Knowledge to Support Customer Relationship Management[END_REF][START_REF] Adomavicius | Toward the Next Generation of Recommender Systems: A survey of the Stateof-the-Art and Possible Extensions[END_REF][START_REF] Linden | Amazon.com Recommendations: Item-to-Item Collaborative Filtering[END_REF].base and the five usual application datasets named U [START_REF] Mahmood | Improving Recommender Systems with Adaptive Conversational Strategies[END_REF][START_REF] Resnick | Recommender Systems[END_REF][START_REF] Ren | An Evolving Information System Based on Data Mining Knowledge to Support Customer Relationship Management[END_REF][START_REF] Adomavicius | Toward the Next Generation of Recommender Systems: A survey of the Stateof-the-Art and Possible Extensions[END_REF][START_REF] Linden | Amazon.com Recommendations: Item-to-Item Collaborative Filtering[END_REF]

A. Cross-Validation and Comparison

With a support threshold of 0.05 and a confidence threshold of 50 %, we extract rules and generate the cohorts summarized in table II. For robustness considerations, only the cohorts with at least 20 rules are retained. To evaluate the Top-N recommendations, we use the two classical metrics widely using in the recommender system community namely, precision and recall. We privilege the precision measure to minimize the number of false positive recommendations. We present in table III the Top-10 and Bottom-10 recommendations sorted by average precision P r. We notice the good results reached with CAPRE on average precision P r on the Top-10 recommendations (even the Bottom-10 results are acceptable). The errors measured by precision and recall are similar to those generally measured on MovieLens with RMSE (Root Mean Square Error). Aggregating rules into cohorts allows to build a more robust model than single rules, this reduces the errors of recommendation. The pre-actionability, actionability and profitability steps allow to prune the recommendations and target users: • Profitability: we could easily assign a monetary value on each movie, considering that the recommendation of a movie like Independence Day yields $5 per streaming viewing. Therefore, the expected profitability of the cohort could be $465. We can note that CAPRE methodology allows (i) to prune sharply the number of recommendations and (ii) to filter efficiently the number of potential customers. Moreover, we observe customer rating behaviors throughout other cohorts: However, in our methodology, we can note that for the cohort recommending Jaws, composed of 7 rules, we prune the number of counter-examples from 144 to 75 actionable customers. By this way, we target more finely the recommendation to the most receptive customers.

VI. CONCLUSIONS AND FUTURE WORKS

In this paper we propose CAPRE, a new methodology for recommender systems based on the analysis of turnover of customers for specific products. By measuring the actionability of rules based on the similarity between the examples and counter-examples, we present an original measure of the profitability of rules. Thus, the methodology allows not only to make recommendations, but also to identify and sort the most receptive customers for these recommendations, and estimate the related profit. CAPRE is applied to over 10,000 individual customers and 100,000 products in the framework of the customer relationship management of VM Matériaux company, thus assisting the salespersons' objective to increase the customer value. The most actionable and profitable recommendations are triggered and an expected ROI is computed, with 7,788 rules and 57 cohorts, representing about $264,373 of additional turnover for plaster recommendation. Thus, recommender systems applied to the analysis of the turnover of customers would improve customer relationships and generate qualitative and quantitative profits for companies. The cross-validation performed on the MovieLens benchmark shows that CAPRE makes good recommendations with little errors (even on the Bottom-10 results).

  1 http://www.vm-materiaux.com

Figure 1 .

 1 Figure 1. The Five Key Steps of CAPRE

  R is partitioned into as many cohorts as consequents of rules. For the cohorts to be robust, we retain only the cohorts with a sufficient number of rules. The sets of examples and counter-examples of the cohort, noted Ct + (p X = µ X ) and Ct -(p X = µ X ) are the union of the sets of examples and counter-examples of the rules respectively. Detecting counter-examples provides us with a precise group of customers who can theoretically allow the company to run at a healthy profit.

•

  Cement = c AND Sand = c → Gravel = c has 190 examples and 70 counter-examples. • Sand = c AND Concrete Mixer = c → Gravel = c has 135 examples and 40 counter-examples.

  c → Plaster = c To illustrate CAPRE, we focus on the cohort to recommend the purchase of plaster Ct(Plaster = c). The cohort is composed of 12 rules, 150 distinct examples and 50 distinct counter-examples.

  the distributions of examples Ct + and counter-examples Ct - on purchases from the antecedents of the cohort. Antecedent purchases of Ct + are characterized by a median M = 37.32 % and an interquartile range Q3-Q1 = 28.05 %. Therefore, the counter-examples outside the interval [9.27 %; 65.37 %] are considered non pre-actionable. This filtering allows to prune 4 counter-examples.

Figure 3 .

 3 Figure 3. Distributions of Examples and Counter-Examples on Purchases from the Antecedents of a rule of the Cohort Ct(Plaster = c)

Figure 4 .

 4 Figure 4. Summary of Counter-Example Filtering for Ct(Plaster = c)

  .test. Training and validation datasets contain 80 % and 20 % of all ratings respectively. Customer characteristic variables, i.e. N ∪ I, are the age, the gender and the occupation. We gathered the ratings into three classes a = [1, 2], b = [3] and c = [4, 5].

Table I SOME

 I RULES EXTRACTED BY THE CHARM ALGORITHM

  of our methodology, 38 actionable counter-examples are $2000-profitable for the purchase of plaster.

	50	pre-actionability	46	actionability	42	profitability	38

  Independence Day = c, 102 customers who watched the action movie Men in Black and voted positively (4 or 5) have also watched the action movie Independence Day by voting positively. However, 43 counter-examples have watched and voted positively Men in Black but did not appreciate Independence Day (vote 1 or 2 or 3). CAPRE allows to obtain 99 cohorts of rules. Let us focus on the cohort recommending the movie Independence Day (see tableIV). The cohort is composed of four rules.

	Table III
	PRECISION (PR) AND RECALL (RC) OF RECOMMENDATIONS
	Top-10	Average Pr Rc
	Rear Window	91,50 78,77
	The Shawshank Redemption	91,07 87,67
	Casablanca	90,02 88,16
	Schindlers List	89,97 90,97
	To Kill a Mockingbird	89,03 78,13
	One Flew Over Cuckoos Nest 88,28 90,39
	The Silence of the Lambs	87,98 89,90
	The Usual Suspects	86,76 87,32
	Star Wars	86,58 89,70
	The Godfather	85,01 86,32
	Bottom-10
	The Terminator	72,87 90,57
	2001: A Space Odyssey	72,32 85,15
	Toy Story	71,48 87,61
	Die Hard	71,36 78,05
	Contact	70,85 66,68
	Indiana Jones	69,81 89,40
	Forrest Gump	69,73 89,14
	Twelve Monkeys	68,04 84,35
	Back to the Future	67,83 92,48
	Jaws	67,29 84,40

B. Example and Discussion

Let us consider the rule Men in Black = c →

Table IV A

 IV COHORT OF FOUR RULES RECOMMENDING Independence Day

	Sup.(%) Conf.(%)	Rules
	0.13	54.12	The Rock = c → Ind. Day = c
	0.12	52.28	Star Trek = c → Ind. Day = c
	0.11	54.12	Men in Black = c → Ind. Day = c
	0.10	56.54	Star Wars = c & The Rock → Ind. Day = c

•

  Some cohorts have scattered examples because their consequents are a very popular movie, such as Star Wars. The behavior of the examples is difficult to characterize through ratings or characteristic variables. That is why, we little prune the counter-examples. • However, for horror movies, examples show much more typical behaviors. Counter-example customers who view and appreciate at least one horror movie of the cohort are generally actionable and profitable because of a more restricted customer rating behavior. That is why, we widely prune the counter-examples. Finally, in [23], Onuma et al. showed that for Horror Movie Fan, their system recommends A Nightmare on Elm Street, The Shining and Jaws.

In CAPRE, the eligible counter-examples are only the customers who purchase p Z at a level lower than µ Z .

http://www.coface.com

We are currently using CAPRE in a real business application. The goal is to use extensively CAPRE for any databases of the type customers × products. This article offers the prospect for further work. First, the measure of profitability can be improved: the actionability of cohorts could be modulated by the number of actionable counter-examples and the net margin of recommended products. Secondly, an application of recommendations can be integrated in the corporate information system for decision makers. Finally, an OLAP-based approach can be adapted to address decision makers by trying to provide aggregated recommendations (e.g. for brands or product families) to certain segments of customers.

ACKNOWLEDGMENT

The authors would like to thank Pierrick Richard and Gaëtan Blain (VM Matériaux) for supporting this work.