Marco Corneli

Charles Bouveyron

Pierre Latouche

Co-Clustering of ordinal data via latent continuous random variables and a classification EM algorithm

Keywords: categorical data, model based clustering, ICL. 1

come

Co-clustering of ordinal data via latent continuous random variables and not missing at random entries

Introduction

Data clustering plays a central role in all scientific and industrial fields where there is a need of data analysis. The goal of data clustering is to group similar data together, to provide a synthetic view of a dataset. Thus, clusters are homogeneous groups of data that can be interpreted in a common way. The nature of the data requires different types of clustering techniques to be applied. This paper focuses on peculiar categorical data in which categories are ordered: the ordinal data [START_REF] Agresti | Analysis of ordinal categorical data[END_REF]. Such data is very common in marketing researches (see for instance [START_REF] Dillon | Marketing research in a marketing environment[END_REF], where people are asked to evaluate products and/or services on an ordinal scale. As an example, the last section of this paper focuses on a related dataset, where Amazom customers review fine foods. In more details, the dataset consists of an ordinal data matrix whose entry (i, j) is the note attributed to the j-th products by the i-th customer.

In the last decades, several clustering algorithms for ordinal data have been introduced in the literature [START_REF] D'elia | A mixture model for preferences data analysis[END_REF][START_REF] Podani | Braun-blanquet's legacy and data analysis in vegetation science[END_REF][START_REF] Gouget | Utilisation des modèles de mélange pour la classification automatique de données ordinales[END_REF][START_REF] Jollois | Classification de données ordinales: modèles et algorithmes[END_REF][START_REF] Giordan | A clustering method for categorical ordinal data[END_REF][START_REF] Fernández | Mixture-based clustering for the ordered stereotype model[END_REF][START_REF] Ranalli | Mixture models for ordinal data: a pairwise likelihood approach[END_REF][START_REF] Biernacki | Model-based clustering of multivariate ordinal data relying on a stochastic binary search algorithm[END_REF]. Moreover, a recent work of McParland and Gormley (2016) adopted a model-based clustering approach for mixed data (including ordinal data) based on latent Gaussian random variables. In the example of the Amazon dataset, the clustering algorithms listed above can group the rows of the ordinal data matrix such that two customers share the same cluster if they tend to rate products similarly. However, when the number of columns of the matrix is high, uncovering relevant row clusters is a particularly hard problem. Moreover, the number of variables makes the interpretation of the row clusters challenging. A solution to this issue is provided by coclustering, which aims at simultaneously clustering the rows and the columns of a data matrix, thus providing a partition for rows and another one for columns. Several model-based co-clustering methods rely on the latent block model (LBM, [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF] who deals with matrices of binary data. Other extensions of LBM can tackle counting [START_REF] Govaert | Latent block model for contingency table[END_REF], real [START_REF] Lomet | Sélection de modèle pour la classification croisée de données continues[END_REF], categorical [START_REF] Keribin | Estimation and selection for the latent block model on categorical data[END_REF], functional [START_REF] Bouveyron | The functional latent block model for the co-clustering of electricity consumption curves[END_REF] and ordinal data (Jacques and Biernacki, 2018). Up to our knowledge, the The present section describes a statistical model to generate arrays of ordinal data. In order to properly take into account data sparsity (see Section 2.2), the model extends the binary LBM of the original paper of [START_REF] Govaert | Block clustering with bernoulli mixture models: Comparison of different approaches[END_REF]. Thus, the first part of this section is devoted to a description of the binary LBM, while the second part focuses on the latent framework adopted to model ordinal data.

Latent block model

Consider an M ×P incidence matrix A such that A ij is equal to 1 if one interaction between i and j is observed (e.g. the i-th user assigns a score to the j-th product), 0 otherwise.

Rows are assumed to be clustered into Q row clusters and columns into L column clusters.

An hidden vector R, of length M , is such that R i = q if the i-th row of A is in the q-th row cluster. Moreover, the i-th row of A is assumed to be assigned to its row cluster according to a multinomial distribution M(1, ρ := {ρ 1 , . . . , ρ Q }), where ρ q > 0, for all q, and Q q=1 ρ q = 1. Thence, being the cluster of the i-th row recorded into R i , it holds that P(R i = q) = ρ q , ∀i.

Similarly, an hidden vector C, of length P , is such that C j = l iff the j-th column of A is in the l-th column cluster. The j-th column of A is assigned to its column cluster according to a multinomial distribution

M(1, δ := {δ 1 , . . . , δ L }),
where δ l > 0, for all l and L l=1 δ l = 1, so that

P(C j = l) = δ l , ∀j.
The two vectors R and C are further assumed to be independent. In the following, when no confusion arises, the equivalent 0-1 notation will be employed. In that case, R will denote a binary M × Q matrix and if the i-th row of A is in the q-th row cluster, R iq = 1 and R iv = 0, for all v = q. Similarly, C will denote a binary P × L matrix.

Conditionally on R and C, the entries of A are all independent and such that

A ij |{R, C} ∼ B(π R i C j ), ∀i ≤ M, j ≤ P
where B(p) denotes a Bernoulli distribution of parameter p, π ql ∈ [0, 1], for all q, l and π := {π ql } q,l . According to this model, the complete data likelihood is

p(A, R, C|π, ρ, δ) = M i=1 P j=1 π A ij R i C j (1 -π R i C j ) 1-A ij M i=1 ρ R i P j=1 δ C j .
(1)

Modeling ordinal data

Let us now consider an M × P matrix Y , whose ordinal entry Y ij , conditionally on A ij = 1, is a random variable taking values in {1, . . . , K}, for some K ∈ N * , not depending on the pair (i, j). For the pairs (i, j) such that A ij = 0, we assume that Y ij = 0. Note that the matrix Y contains both observed and missing data. The observed data are the values Y ij corresponding to A ij = 1 and in real applications these values could be scores that users assign to some products. However, one user could (and generally will) rate only a subset of products. Thus, the unrated products are seen as missing values and coded as 0 in Y . Now, the sparsity of A is modelled by LBM. For instance, users densely ranking a single class of products are more likely to be clustered together (via π) when fitting the model to the data. Hence, the link between Y and A has an important consequence: the missing data in Y are not missing at random (see both [START_REF] Little | Statistical analysis with missing data[END_REF]Jacques and Biernacki, 2018). Before going further, two assumptions should be made.

Assumption 1. Henceforth, we assume that an ordinal scale is consistently defined. For instance, in the example of customers evaluating products, 1 always means "very poor" and K always means "excellent". The assumption is necessary, otherwise the results obtained when fitting the model to the data would be completely misleading. The analyst should therefore take this into account when designing the data collection.

1. If it was not the case, a scale conversion pre-processing algorithm (see for instance [START_REF] Gilula | A study into mechanisms of attitudinal scale conversion: A stochastic ordering approach[END_REF] should be employed to normalize the number of levels.

The model that we assume to be generating Y relies on hidden Gaussian random variables Z ij such that

Z ij |{A ij = 1, R, C} ∼ N (µ R i C j , σ 2 R i C j ). (2)
Henceforth, µ := {µ ql } q,l and σ 2 := {σ 2 ql } q,l will denote the sets of the Gaussian parameters. Similarly to Y ,

Z ij |{A ij = 0, R, C} = 0 a.s. ( 3
)
and all the random variables Z ij are collected into an hidden M × P matrix denoted by Z.

Assume that K -1 unknown real numbers (thresholds) γ := (γ 1 , . . . , γ K-1 ) are such that

-∞ =: γ 0 < γ 1 < • • • < γ K-1 < γ K := ∞.
Then, conditionally on the event

{A ij = 1} Y ij := K k=1 k1 ]γ k-1 ,γ k [ (Z ij ), (4)
where 1 Ω (•) is the indicator function over the set Ω ⊂ R. Finally, conditionally on {A ij = 1}

as well as R and C, we assume that the pairs (Y ij , Z ij ) are mutually independent. Hence, the joint density of (Y, Z) can be written as

p(Y, Z|A, R, C, µ, σ 2 ) = M i=1 P j=1 φ(Z ij ; µ R i C j , σ 2 R i C j )1 ]γ Y ij -1 ,γ Y ij [ (Z ij ) A ij , (5)
where φ(•; µ ql , σ 2 ql ) is the probability density function of a Gaussian distribution N (µ ql , σ 2 ql ) and we used that, conditionally on {A ij = 0}, the pairs (Y ij , Z ij ) are equal to (0, 0) a.s.

Eqs. 1-5 can be combined to obtain the complete data likelihood

p(Y, Z, A, R, C|θ) = p(Y, Z|A, R, C, µ, σ 2 )p(A, R, C|π, ρ, δ), (6)
where θ := {µ, σ 2 , π, ρ, δ} denotes the set of the model parameters.

Inference

In the first part of this section, the numbers Q of row clusters and L of column clusters are assumed to be given. A model selection criterion will be detailed later. Now, we aim at estimating the model parameters θ as well as the most likely posterior values of R and C. Let us start with a remark.

Remark 1 (Thresholds). It is immediate to see that either γ or (µ, σ 2 ) need to be fixed in order for the model parameters to be identifiable and, from a generative point of view, it seems reasonable to fix γ, as it can be seen in Eq. (4). However, notice that once γ is fixed to some value (for instance by randomly selecting K -1 Gaussian quantiles and sorting them) and the model fitted to the data, the estimated (parameters of the) random variables in Z lie in a space which is in general not related with the range of the ordinal entries in Y . Thus, in order to have easily interpretable results, γ is fixed as γ = (1.5, 2.5, . . . , (K -0.5)).

In order to illustrate the estimation strategy in detail, we need the following proposition.

Proposition 1. Conditionally on the event {A ij = 1, R iq = 1, C jl = 1}, the random variable Y ij has probability mass function

P(Y ij = k|A ij = 1, R iq = 1, C jl = 1) = η (q,l) k 1 {1,...,K} (k), (7)
where

η (q,l) k := Φ γ k -µ ql σ ql -Φ γ k-1 -µ ql σ ql . ( 8
)
Proof. By the definition of marginal probability density function, it follows that

p(Y ij |A ij = 1, R, C, θ) = R p(Y ij , z|A ij = 1, R, C, θ)dz = R φ(z; µ R i C j , σ 2 R i C j )1 ]γ Y ij -1 ,γ Y ij [ (z)dz = Φ γ Y ij -µ R i C j σ R i C j -Φ γ Y ij -1 -µ R i C j σ R i C j ,
where Φ(•) denotes the cumulative density function of the Gaussian distribution N (0, 1).

If we denote

η (q,l) k := Φ γ k -µ ql σ ql -Φ γ k-1 -µ ql σ ql , it is immediate to verify that K k=1 η (q,l) k = 1.
Propostion 1 has two important consequences. First, the posterior density function of

Z ij can be obtained by p(Z ij |Y ij , A ij = 1, R, C, θ) = φ(Z ij ; µ R i C j , σ 2 R i C j ) η (R i ,C j ) Y ij 1 ]γ Y ij -1 ,γ Y ij [ (Z ij ). ( 9
)
The above probability density function defines a truncated Gaussian distribution and it is fully tractable. Second, due to independence arguments, the marginal likelihood

p(Y |A, R, C, θ) can be computed as p(Y |A, R, C, θ) = M i=1 P j=1 η (R i ,C j ) Y ij A ij , (10)
where we used again that Y ij = 0 when A ij = 0, a.s.

In the light of these results, it is possible to design an estimation strategy, called C-EM, consisting of the following two steps:

1. C step. The model parameters being fixed to a local optimum, log p(Y, A, R, C|θ)

is maximized with respect to R and C in a greedy fashion. This classification step replaces the Expectation step in the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. See also [START_REF] Celeux | A classification em algorithm for clustering and two stochastic versions[END_REF] for a description of Classification EM algorithms.

2. M step. R and C being fixed, the likelihood

log p(Y, A, R, C|θ) = log p(Y |A, R, C, µ, σ 2 ) + log p(A, R, C|π, ρ, δ) (11)
is maximised with respect to the model parameters θ. As we will see in the next sections, the maximization with respect to (π, ρ, δ) is straightforward. On the contrary, the first term on the right hand side of the above equation (detailed in Eq. 10) cannot be directly maximized with respect to µ and σ 2 and no close formulas can be derived for these Gaussian parameters. Therefore, we will rely on Q × L independent EM algorithms to maximize this term with respect to (µ, σ 2 ).

The above two steps are alternatively repeated until convergence of log p(Y, A, R, C|θ).

Each step is detailed in the following sections.

C step

Let us focus on the log-likelihood in Eq. ( 11) and assume that the model parameters are fixed to the (local) optima θ, obtained in the M step (see Section 3.2). The goal of the C step is to maximize the left hand side of Eq. ( 11) with respect to R and C . No closed formula exists for such combinatorial maximization problem and testing all possible combinations (Q M L P ) would be computational prohibitive. Thus, we rely on a greedy search strategy not looking for all the possible solutions. Greedy strategies are quite popular in network and bipartite network analysis (see for instance [START_REF] Côme | Model selection and clustering in stochastic block models based on the exact integrated complete data likelihood[END_REF][START_REF] Wyse | Inferring structure in bipartite networks using the latent blockmodel and exact icl[END_REF].

The basic idea is to swap each row (column) of A to the row (column) cluster leading to the highest increase of the log-likelihood log p(Y, A, R, C| θ). Of course, if no swap increases the log-likelihood, the row (column) is not moved.

Assume that the i-th row of A formerly belonging to the cluster q is now moved to the cluster q . Furthermore, let us denote by R * the row label vector R after that the swap occurred. Thus

∆ i:q →q := log p(Y, A, R * , C| θ) -log p(Y, A, R, C| θ) = log p(Y |A, R * , C, η) p(Y |A, R, C, η) + log p(A|R * , C, π) p(A|R, C, π) + log p(R * |ρ) p(R|ρ)
is the increase (possibly null or negative) of the log-likelihood in Eq. 11 when moving row i from cluster q to cluster q . Moreover, we used η := {η

(q,l) k } q,l,k , where η (q,l) k
is defined in (8). By using Eqs. 1 and 10, we can further obtain

∆ i:q →q = P j=1 A ij log η(q ,C j ) Y ij η(q ,C j ) Y ij + P j=1 A ij log πq C j πq C j + (1 -A ij ) log 1 -πq C j 1 -πq C j + log ρq ρq . (12)
This quantity can be computed in O(P ) and can be used to rank the possible swaps of the i-th row of A to all row clusters. An equivalent formula can be obtained to assess the contribution of a column swap into a column cluster.

It is important to notice that, since Q and L are fixed, one row (column) alone in its current row (column) cluster is not allowed to move. In case one row (column) remains alone in its group, another criterion (that will be introduced in Section 3.4) will decide whether that group is suppressed or not.

M step

The label vectors R and C are fixed throughout this section. Notice that the maximization of the right hand side of Eq. ( 11) with respect to (π, ρ, δ) only involves the second term.

Moreover, this maximization is straightforward. Taking the logarithm in Eq. 1, differentiating with respect to π, ρ and δ and setting the derivatives equal to zero leads to the following stationary points πql : =

M i=1 P j=1 R iq C jl A ij M i=1 P j=1 R iq C jl , (13) ρq :
= M i=1 R iq M , (14) δl :
= P j=1 C jl P , (15)
for all q, l.

As anticipated in Section 3.1, the maximization of the first term on the right hand side of Eq. ( 11) is more challenging. Since no close formula for such optimization does exist, let us consider the following inequality

log p(Y |A, R, C, µ, σ 2 ) ≥ E Z log p(Y, Z|A, R, C, µ, σ 2 ) p(Z|Y, A, R, C, µ 0 , σ 2 0 ) , ( 16
)
where the expectation is taken with respect to Z following the posterior probability density function p(•|Y, A, R, C, µ 0 , σ 2 0 ). The inequality comes from a standard variational decomposition (see for instance Ch.10, [START_REF] Bishop | Pattern Recognition and Machine Learning (Information Science and Statistics)[END_REF], it holds for all (µ, σ 2 ) and it turns into an equality when (µ, σ 2 ) is equal to (µ 0 , σ 2 0 ). Since the posterior distribution of Z is known and tractable, the EM algorithm can be used to provide numerical estimates of μ and σ2 . M-Expectation. We now focus on the right hand side of the inequality in Eq. ( 16). By taking the logarithm of Eq. 5, it holds that

E Z log p(Y, Z|A, R, C, µ, σ 2 ) = - 1 2 M i=1 P j=1 A ij log σ 2 R i C j + E Z ij (Z ij -µ R i C j ) 2 σ 2 R i C j + c = - 1 2 Q q=1 L l=1 M i=1 P j=1 A ij R iq C jl m (2) ij + µ 2 ql -2µ ql m (1) ij σ 2 ql + log σ 2 ql + c, (17)
where c regroups the constant terms not depending on (µ, σ 2 ). We assumed that

Z ij ∈ ]γ Y ij -1 , γ Y ij [ a.s. for all i, j such that A ij = 1 and m (1) ij : = E Z ij [Z ij ] = µ R i C j -σ R i C j φ(β ij ) -φ(α ij ) Φ (β ij ) -Φ (α ij ) , (18) m
(2)

ij : = E Z ij [(Z ij ) 2 ] = (µ R i C j ) 2 -2σ R i C j µ R i C j φ(β ij ) -φ(α ij ) Φ(β ij ) -Φ(α ij ) -σ 2 R i C j β ij φ(β ij ) -α ij φ(α ij ) Φ(β ij ) -Φ(α ij ) -1 , (19)
with

α ij = γ Y ij -1 -µ R i C j σ R i C j , β ij = γ Y ij -µ R i C j σ R i C j ,
and E Z ij is the expectation taken with respect to Z ij following the probability density function in Eq. 9. Notice that both m

(1) ij and m

(2) ij only depend on the pair (i, j) via R i , C j and Y ij . Thus, m

(1) ij and m

(2) ij are the same for all pairs (i, j) in clusters (q, l), respectively, associated with the score k.

M-Maximization.

Once the expectation in Eq. 17 is computed, the right hand side of the equality can be maximized with respect to (µ, σ 2 ). The maximization of µ ql can be performed independently of σ 2 ql , but the opposite is not true. Thus, we first differentiate the right hand side of Eq. 19 with respect to µ ql and set the derivative equal to zero to obtain the following stationary point μql :=

M i=1 P j=1 A ij R iq C jl m (1) ij M i=1 P j=1 A ij R iq C jl , (20)
for all q, l. In order to compute the optimal σ 2 ql , μ is plugged into Eq. 17 in place of µ and then differentiating with respect to σ 2 ql and setting the partial derivative to zero leads to

σ2 ql := M i=1 P j=1 A ij R iq C jl m (2) ij M i=1 P j=1 A ij R iq C jl -μ2 ql , (21)
for all q, l.

The two steps of the EM algorithm described so far are part of the M step. Therefore they are called M-Expectation (Eqs. 18-19) and M-Maximization (Eqs. 20-21). They are alternatively applied up to convergence. We stress that all the equations involving M-Expectation and M-Maximization factorized over q and l. Thus, Q × L independent EM algorithms are used, one for each pair (q, l) of clusters and this task can be done in parallel.

Initialization

Assuming that Q and L are momentarily fixed, the C-EM algorithm described in the previous sections needs some initial values of R and C to be provided. Then, a first M step can be implemented, followed by a greedy search C step and so on. In this paper two different initialization strategies are considered:

1. Multiple random initializations. Both the initial R and C are independently sampled from multinomial distributions with uniform parameters. Since the C-EM algorithm is not guaranteed to converge toward a global optimum, the algorithm is provided in the applications with several independent initializations. The estimates R and Ĉ leading to the highest log-likelihood are finally retained. This initialization strategy is assessed in Section 4.4.

2. K-means initialization. Two k-means algorithms are independently run on the rows and the columns of the matrix Y . The C-EM algorithm is then initialised with the estimates R and Ĉ provided by the two k-means.

Notice that missing values could be present in Y . This is not a problem when adopting random initializations, but it can be one when using k-means. Indeed, when the proportion of missing data (i.e. zeros) in Y is very large, the k-means algorithm will provide very poor initial estimates of R and C. In a similar scenario, it is preferable to opt for multiple random initializations.

The pseudocode in Algorithm 1 summarizes the estimation routine detailed so far.

Algorithm 1 Pseudocode

1: function Estim(Y ,A,Q,L, type) 2:
(R, C) ← Init(Y , type) type is "multiple random" or "k-means" 3:

while log p(Y, A, R, C|θ) increases do

4: θ ← M step Including Q × L M-EM algorithms 5: (R, C) ← C step 6:
end while 7:

return ( R, Ĉ, θ)

8: end function

Model selection

So far, the numbers Q and L of row and column clusters were assumed to be known. Of course, in real applications, this assumption is too restrictive. Thus, we now detail a model selection criterion we propose to select the numbers of row and column clusters.

In clustering contexts, the integrated classification likelihood (ICL, [START_REF] Biernacki | Choosing starting values for the em algorithm for getting the highest likelihood in multivariate gaussian mixture models[END_REF] criterion is often used to approximate a complete data integrated log-likelihood and to select the number of components. In our case

ICL(Q, L) ≈ log p(Y, A, R, C) = D θ log p(Y, A, R, C|θ)ν(θ)dθ,
where, in a Bayesian framework, the model parameters θ are seen as random variables and the integral is taken over the support D θ of any prior probability density function ν(•). The following proposition details the functional form of the ICL for our model.

Proposition 2. An ICL(Q, L) criterion for the generative model described in Section 2 is

ICL(Q, L) = log p(Y, A, R, Ĉ| θ) -QL log D - QL 2 log(M P ) - Q -1 2 log M - L -1 2 log P, (22)
where θ, R, Ĉ are the stationary points obtained after convergence of the algorithm described in Sections 3.2-3.1.

Proof. By definition of conditional probability density function, it follows that

ICL(Q, L) ≈ log p(Y |A, R, C) + log p(A, R, C).
The last term on the right hand side is the complete data integrated log-likelihood of a binary LBM. This log-likelihood can be approximated as follows

log p(A, R, C) ≈ max π,ρ,δ log p(A, R, C|π, ρ, δ) - QL 2 log(M P ) - Q -1 2 log M - L -1 2 log P, (23)
see for instance [START_REF] Keribin | Model selection for the binary latent block model[END_REF]. In order to approximate log p(Y |A, R, C) we propose a BIC-like approximation [START_REF] Schwarz | Estimating the dimension of a model[END_REF] log

p(Y |A, R, C) ≈ max µ,σ 2 log p(Y |A, R, C, µ, σ 2 ) -QL log D, ( 24
)
where D is the number of ordinal (i.e. not zero) entries in Y and 2QL accounts for the number of parameters in µ and σ 2 . Combining Eqs. 23-24, the proposition is proven.

When fitting to data the model presented in Section 2, the ICL criterion in Eq. ( 22) Notice that, as long as either Q or L lie on the boundary of the grid, Q max and/or L max can be increased to obtain a solution which is interior. In this sense, the grid search described so far is exhaustive.

is
However, when dealing with massive datasets, this strategy could be computationally prohibitive since the number of elements in the grid is Q×L. Alternatively, a greedy search algorithm can be employed to select Q and L. The ICL criterion is still used, but it is not computed for all values of Q and L in the grid. We propose here a greedy search algorithm inspired by the one introduced in [START_REF] Keribin | The latent block model: a useful model for high dimensional data[END_REF]. Initially, both Q and L are set to 1. In a second step, Q is increased by one and the value of ICL(2, 1) is computed via the C-EM algorithm and recorded. Then, ICL(2, 1) is compared with ICL(1, 2), obtained by setting Q = 1 and L = 2. Thus, if

ICL(2, 1) > ICL(1, 2)
Q is definitely set to 2 and L to 1. The opposite otherwise. This routine is recursively applied in such a way that, if Q * and L * are the current values of Q and L, then ICL(Q * + 1, L * ) is compared with ICL(Q * , L * + 1) and Q * and L * are updated accordingly. The algorithm stops when no further increase in the ICL criterion is possible. A pseudo-code illustrating the greedy search detailed so far is reported in Appendix A.1.

Experiments on synthetic data

The estimation procedure detailed in Section 3 is now tested on simulated data. Two main scenarios are considered: the former adopting the generative model described in Section 2 and the latter adopting the generative model described in Jacques and Biernacki (2018).

Henceforth, these two generative models will be referred to as OLBM (Ordinal Latent Block Model, our proposal) Co-Clustering and BOS (Binary Ordinal Search) Co-Clustering. For each simulated scenario, both OLBM-CC and BOS-CC are fitted to the data to provide estimates of R and C, and the results are compared. The parameters of BOS-CC are estimated via the R package ordinalClust 1 . A third competitor approach is considered in this section, namely LBM for continuous data (cLBM) as described in [START_REF] Bhatia | blockcluster: An r package for model-based co-clustering[END_REF].

This model is fitted to the data via the R package blockcluster 2 . Unfortunately, this package does not support missing data in Y . Therefore, comparisons will not always be possible.

Data simulated according to OLBM-CC

For this first experiment, the data are simulated according to the generative model described in Section 2. We consider incidence matrices with M = 150 rows, grouped into Q = 3 row U 1 , . . . , U K-1 independent random variables, uniformly distributed in [0, 1]. Without loss of of generality, let us assumed that they are sorted. A last variable U K = 1 is introduced.

Then:

γ i := K + 1 2 + Φ -1 (U i ) ,
for all i ∈ {1, . . . , 5} and γ 0 = -∞. Then, the Gaussian parameters µ and σ 2 are set to

µ =      0 0 ζ -ζ -ζ ζ      , σ 2 =      1.2 1.2 1.4 1.4 1.0 1.0     
, where ζ ≥ 0 is a real parameter controlling how the block distributions differ: as long as

ζ is far enough from 0, we expect that the estimation algorithms would correctly estimate the row and the column clusters. When ζ approaches to 0, the column clusters become indistinguishable and the row groups are only separated via the matrix σ 2 (see Figure 1).

The last parameter to set is π, defining the probability of observing an ordinal entry in Y . Three different setups are considered.

No missing data. In this framework, π ql = 1, for all q, l and only the Gaussian parameters (µ, σ 2 ) induce a block structure. For each value of ζ, fifty datasets Y are independently simulated according to the setup described so far and the three approaches (OLBM-CC, cLBM, BOS-CC) are fitted to each dataset. At first, the selection of Q and L is not consid- Missing at random data. In this framework, the probability of a missing data in Y is independent of the cluster assignments. In other words, it is the same for each entry Y ij . We set π ql = 0.7, for all q, l (30% of missing data in Y , on average) and repeat the experiment.

As previously mentioned, the blockcluster package does not support missing data. Thus, in the reminder of this section, OLBM-CC will only be compared with BOS-CC. The results can be observed in Figure 3, where boxplots of the ARIs can be observed in two sub-figures. The one on the left hand side reports the ARIs for the row labels estimates. The ARIs for the column labels estimates are reported in Figure 3b. These results are coherent with the ones in Figure 2. Since missing at random data does not bring any information, this experiment reduces to the previous one with less ordinal data.

This explains the slightly worse performance of OLBM-CC and BOS-CC in Figure 3 with respect to Figure 2.

Not missing at random data. In the reminder of this section, the assumption of missing at random data is relaxed. An alternative setup can be obtained by adopting the following connectivity matrix

π =      0.5 0.7 0.7 0.5 0.7 0.5     
, meaning that the probability of a missing data in Y is no longer independent on the cluster assignments. It is lower (30%) in cluster pairs (1, 2), (2, 1) and (3, 1) and higher (50%) for the remaining block pairs. Other settings being unchanged, the experiment is repeated (50 simulated data matrices Y for each value of ζ and Q, L known). Results can be seen in via OLBM-CC (blue bars) are more accurate (higher ARIs) for both rows and columns. In particular, the average column ARI is around 1 even when ζ = 0. Indeed, the π matrix itself induces an additional block structure that discriminates the two column clusters although ζ = 0. Such a framework, in which the block structure is accentuated by the connectivity patterns, is very common in real data (see also Section 5). In contrast with OLBM-CC, BOS-CC has very similar performances with or without missing at random data. Since it cannot deal with block dependent missing data, the additional information carried by π cannot be exploited by the model.

Data simulated according to BOS-CC

In order to present fair results regarding methods other than ours, the datasets are now generated according to the BOS-CC model. The reader is referred to Jacques and Biernacki

µ =   2 3 1 2 1 3   , =   (0.4 -ζ) (0.4 -ζ) (0.4 -ζ) (0.4 -ζ) (0.4 -ζ) (0.4 -ζ)   ,
where, as long as ql = 0, µ ql can be seen as the mode of the ordinal entries associated with the pair (q, l). The parameter ql measures the dispersion around the mode, which is minimal when ql = 1 and maximal when ql = 0 (in this case the ordinal entries of the pair (q, l) are uniformly distributed in {1, . . . , K}). As in the previous section, ζ is a real parameter which controls how different the block distributions are. Here, however, As it can be seen, 50% of the ordinal entries corresponding to the block pairs (2, 1), (2, 2) and (3, 2) are randomly replaced by 0 according to π. Instead, only 30% of the ordinal entries corresponding to the remaining block pairs are replaced by 0. Other settings being unchanged, the experiment is repeated and results can be seen in Figure 6. As in the previous section, OLBM-CC exploits the information carried by π to produce better estimates of R and C, whereas BOS-CC performs slightly worse than in the missing at random framework. Notice that, when ζ approaches to zero, the only parameter allowing to discriminate clusters is π. However, while the row clusters are perfectly discriminated, the second and the third column clusters are indistinguishable when looking at π. This explains why the median columns ARI for OLBM-CC is around 0.6 on the right hand side of Figure 6.

Scalability

When dealing with massive datasets, the question of the scalability of the algorithm is of course of great interest. A deep understanding of the computational complexity of the C-EM algorithm (see Sections 3.2 and 3.1) is outside the scope of this paper. Nonetheless, this section aims at providing some insights about the scalability of the OLBM-CC estimation algorithm. In particular, we aim at assessing how the algorithm behaves when either the number of rows/columns of Y or the number of ordinal levels K increase. At first, M is set equal to P and both vary between 50 and 500. Data are missing at random, in such a way that the (mean) number of ordinal entries in Y is equal to 10 × M . Then, the data are simulated according to the OLBM-CC generative model and all the remaining parameters (Q, L, K, µ, σ 2 , γ) are as in Section 4.1. For each simulated matrix Y , the three estimation algorithms considered so far (OLBM-CC, cLBM, BOS-CC) are provided with Figure 8 highlights another feature of the C-EM estimation algorithm: its scalability with respect to number of ordinal levels K. The previous experiment was repeated with M = P = 100, but now K ranges is {3, . . . , 13}. For each value of K, 10 OLBM-CC estimation algorithms are independently run. In Figure 8, the average running times are plotted versus the number of ordinal levels K. As it can be seen, the computing time of OLBM-CC does not seem to be dependent on the number of ordinal levels.

Initialization

In the previous experiments, for each simulated matrix Y , the OLBM-CC estimation algorithm was provided with a single k-means initialization. This section aims at comparing k-means and random initializations (see Section 3.3). Here, the data is simulated according to the BOS-CC generative model described in Section 4.2, with ζ = 0.125 and missing at random values. As in can be seen in Figure 5, our approach (blue bars) works quite well when ζ = 0.125. However, some outliers can be observed both in row and column ARIs.

Thus, 50 matrices Y are independently sampled according to the setup described in Section 4.2 and the OLMB-CC estimation algorithm is run on each matrix, provided with two different initializations: a k-means initialization and a purely random one. Not surprisingly, as in can be seen in Figure 9a, both the row and the column label estimates are more accurate when a k-means initialization is provided. The experiment is now repeated but the multiple random initializations detailed in Section 3.3 is adopted in place of a single k-means. More in details, 10 independent random initializations are used for each dataset.

The results can be seen in Figure 9b. Two remarks can be made: first, 10 random initializations are enough to sensibly reduce the gap between the multiple random initializations and k-means initializations; second, the outliers ARIs in Figure 5 are no more present, due to the increased number of k-means initializations. We stress that Figures 9a and9b refer to the very same simulated dataset.

Model selection

So far, the numbers Q of row clusters and L of column clusters were assumed to be known.

However, in real applications, the pair (Q, L) needs to be estimated from the data. This can be done via the ICL criterion in Eq. ( 22). In order to assess the criterion, the BOS-CC 2), the number of times ICL correctly estimates Q and L is lower and it is not surprising. As it can be seen in Eq. 22, the estimated complete data log-likelihood plays a central role in the computation of ICL.

Thus, a less accurate estimate of R and/or C leads to a lower value of the log-likelihood and hence of the ICL.

Amazon fine foods

This section focuses on a real dataset consisting of reviews of fine foods from Amazon. The dataset can be freely downloaded at https://snap.stanford.edu/data/web-FineFoods.

html. A time horizon of 10 years is considered, up to October 2012. The number of reviews reported is 568,464 and in the original dataset, each row corresponds to one review. Some additional information is reported for each review: the user/product numerical identifiers, a summary of the review and a rating attributed to the product by the user. The rating is expressed via an integer number spanning from 1 (very bad) to 5 (very good). are reported (without hat, to keep the plot uncluttered). We present hereafter some (not exhaustive) remarks about the results.

1. As it can be seen by looking at both Figures 10 and 11, the users in row cluster q = 4 have a peculiar behaviour, both in terms of missing values and score assignments.

They assign most of the scores to the products in column cluster l = 5 and they do not review goods in clusters l = 2 and l = 3. Notice also that the co-cluster (q = 4, l = 5) is the only one not containing missing values (π 45 = 1). The products in column cluster l = 5 are all herbal teas of the same brand "alvita" and the most common rating is 5 (over 60% of scores). When taking a look to the texts associated with the scores on the Amazon website, we noticed that most texts are similar to each other. One might think that users in cluster q = 4 are paid to review.

2. Still on column cluster l = 5. As it can be seen in Figure 11, products in this cluster are only rated by users in row clusters q = 2, q = 4 and q = 6. The score distributions are very different from one row cluster to another. Users in row cluster q = 6 only rate 2 and the underlying Gaussian distribution is peaked in 2. Users in cluster q = 2 all note 1, but in that case the mode of the underlying Gaussian distribution is slightly shifted toward 0 with a higher variance. It reflects the fact that 1 is the worst note that one user can assign (some users would even rate worse if they could).

3. Products in column cluster l = 1 are a mix of food and beverages, including (e.g.) cat food, teas, coffees and chips. Together with products in cluster l = 4 they are scored different shapes from one row cluster to another: very positive ratings for q = 1, negative ratings for q = 2, neutral ratings for q = 3 and positive but (on average) not excellent ratings for q = 5. The column cluster q = 2 is certainly the one exhibiting the highest variety of preferences. Notice also that, the amount of missing data is higher in cluster pairs (q = 3, l = 2) and (q = 5, l = 2) than in cluster pairs (q = 1, l = 2) and (q = 2, l = 2).

5. In Figure 10, row clusters q = 3 and q = 5 look very similar (except for column cluster l = 1, where users q = 3 note more frequently than those q = 5). However, when looking at Figure 11, we see that the way they note is quite different. Users q = 5 are more demanding than users q = 3 on goods l = 1 and less demanding on goods l = 3. More important, on column cluster l = 3, the one where they both are more active (π 33 = 14.48% and π53 = 20.71%), users q = 3 are more enthusiastic than users q = 5.

6. Users in row cluster q = 2 are definitely the more demanding ones: the means of their underlying Gaussian distributions are constantly the leftmost ones. For instance, products in column cluster l = 6 are globally well rated except for users in q = 2.

Similar analyses can be done for the remaining blocks. This experiment demonstrated that OLBM-CC can be fitted to large and very sparse datasets to provide a synthetic and comprehensive view.

Conclusion and perspectives

A new method for the co-clustering of ordinal data has been introduced in this paper. This method relies on the binary LBM to manage data sparsity and adopts latent Gaussian random variables to generate ordinal entries in a data matrix. In our view, the reduced computational burden of the estimation procedure, the modeling of missing data and the easy interpretation of the latent distributions are the main advantages of the outlined approach.

A Appendix

A.1 Greedy model selection Calculate T 2 = ICL(Q * , L * + 1)

6:
if T 1 < T 2 then 7:

Q * = Q * + 1 8:
else if T 2 ≤ T 1 then

  computed for several values of Q and L. The pair ( Q, L) leading to the highest value of ICL(Q, L) is finally retained. An exhaustive strategy would consist into fixing some sufficiently high values Q max and L max and computing ICL(Q, L) for all (Q, L) in the grid {1, . . . , Q max } × {1, . . . , L max }.

  clusters and P = 100 columns, grouped into L = 2 column clusters. Rows and columns are randomly assigned to their clusters in uniform proportions. The observed ordinal entries in Y take values in {1, . . . , K} and K = 5. Recalling that Φ(•) denotes the standard Gaussian cumulative distribution function, the thresholds γ 1 , . . . , γ K are set as follows. We sample 1 https://cran.r-project.org/web/packages/ordinalClust/index.html 2 https://cran.r-project.org/web/packages/blockcluster/index.html

Figure 1 :

 1 Figure 1: Histograms of the ordinal entries of Y organized by block pairs. On the left hand side figure, a high value of ζ induces asymmetries in both the histograms and the underlying Gaussian distributions. See in particular the block pairs (2, 1),(2, 2) and (3, 1),(3, 2). On the right hand side figure, ζ = 0 and the two columns are indistinguishable.

Figure 2 :

 2 Figure 2: Results on the dataset simulated according to OLBM-CC -Not missing data.

Figure 3 :

 3 Figure 3: Results on the dataset simulated according to OLBM-CC -Missing at random data.

Figure 4 .

 4 Figure 4. As expected, with respect to the missing at random setup, the estimates obtained

Figure 4 :

 4 Figure 4: Results on the dataset simulated according to OLBM-CC -Not missing at random data.

(

  2018) for a full description of the model. The simulated incidence matrices have now M = 100 rows and P = 150 columns. The rows are clustered in Q = 2 groups and the columns in L = 3 groups. The model parameters are

ζ

  ∈ [0, 0.4] and the contrast is maximum when ζ = 0 whereas the row and column clusters are indistinguishable when ζ = 0.4. Missing values are injected into Y in two different ways. Missing at random data. For some values of ζ in [0, 4] fifty data matrices Y are independently sampled according to the setup detailed so far. Then, 30% of the entries of each matrix Y is randomly selected and replaced by missing values. The true values of Q and L are assumed to be known and both OLBM-CC and BOS-CC are fitted to each dataset. Results can be seen in Figure 5. Not surprisingly, BOS-CC globally provides

Figure 5 :

 5 Figure 5: Results on the dataset simulated according to BOS-CC -Missing at random data.

Figure 6 :

 6 Figure 6: Results on the dataset simulated according to BOS-CC -Not missing at random data.

  times (without BOS-CC).

Figure 7 :

 7 Figure 7: Figure 7a reports the running times of the three competitor algorithms versus the number of rows of Y , with (M = P ). Figure 7b zooms on the running times of OLBM-CC and cLBM.

Figure 8 :

 8 Figure 8: The number K of ordinal levels (on the horizontal axis) varies between 3 and 13. The running times of the C-EM estimation algorithm for OLBM-CC are computed and box-plotted (10 runs for each K).

Figure 9 :

 9 Figure 9: Fifty data matrices Y are simulated according to the BOS-CC generative model, with ζ = 0.125. The OLBM-CC estimation algorithm is run on each dataset provided with one k-means initialization and one random initialization. Boxplots of the resulting ARIs are plotted in Figure 9a. Figure 9b reports the results of the same experiment but the number of provided initializations of each type is 10.

  Easier setup, 50 simulated datasets. In bold, the number of times the actual values of Q and L are recovered by ICL. Harder setup, 50 simulated datasets. In bold, the number of times the actual values of Q and L are recovered by ICL. correctly estimated by the ICL criterion. The criterion succeeds 49 times over 50 and only fails once, by selecting L = 4. In the harder scenario (Table

Figure 10

 10 Figure 10 reports the reorganised incidence matrix A. Darker regions correspond to lower portions of missing data in the corresponding co-clusters. As in can be seen, the ICL criterion selected Q = L row clusters and L = 6 column clusters. The way the scores are assigned on each co-cluster can be assessed by looking at Figure 11. The score frequencies of each co-cluster are plotted as histograms and, in grey, one can see the estimated underlying Gaussian distributions. On the top of each histogram, the estimated parameters ρ, δ and π

Figure 10 :

 10 Figure 10: The incidence matrix A of the Amazon fine food notes reorganised according to the estimates R and Ĉ, provided by the C-EM algorithm. The ICL criterion selected Q = 6 row clusters and L = 6 column clusters.

Table 3 :

 3 To focus on the most meaningful part of the data, we only considered the users reviewing more than 26 By doing that, an ordinal data matrix Y with M = 1, 644 rows and P = 1, 733 columns was obtained by neglecting all the information but the ratings. The entry Y ij was either an ordinal entry (a score) or a missing value. The number of observed ordinal entries in Y The score frequencies in Y for the Amazon fine food data.

	(the scores) is 32, 836, corresponding to 98.85% of missing data. The score frequencies are
	reported are reported in Table 3.
	Scores	1	2	3	4	5
	Frequencies 1849 2126 4174 7912 16775

Due to the dimensions of Y , adopting a grid search to select Q and L via the ICL criterion (see Section 3.4) would be very long. Thus, we opted for the greedy search scheme described in Section 3.4. Moreover, as pointed out in Section 3.3, a k-means initialization is useless when the rows/columns of Y contain a majority of missing values and here it is the case. Therefore, for each value of (Q, L), the C-EM algorithm was initialized with multiple (25 times) random initializations. The highest ICL criterion (for each pair) was finally retained. The co-clustering of Y provided by our method can be observed in Figures 10 and 11

.

by users of all row clusters. However, users seem to be more satisfied by the products l = 4 than by those l = 1. It can clearly be seen in Figure 11 for row clusters q = 2, q = 5 and q = 6. Indeed, the products l = 4 are the best rated in the dataset: the mode is 5 in all row clusters.

4. In Figure 11 we see that the note distributions in column cluster l = 2 have very First, it could be useful to assess the advantages/disadvantages of using not Gaussian latent random variables to model ordinal data. Indeed, other distributions could be employed to capture some features in the data, e.g. asymmetric frequencies. Second, Section 3.4 describes a greedy search algorithm to select the number of row/column clusters when an exhaustive grid search is computationally prohibitive. Alternative greedy schemes could certainly be implemented, for instance based on the genetic algorithms described in Scrucca (2016).


Farah Haidar

Arnaud Kaiser

Brigitte Lonc

Pascal Urien

C-ITS PKI protocol: Performance Evaluation in a Real Environment

Keywords: C-ITS security, PKI, performance evaluation

In the near future, vehicles and roadside units (RSU) will communicate and cooperate by broadcasting V2X messages over the vehicular network (IEEE 802.11p). These messages are used by safety applications to improve road safety and traffic efficiency. However, those messages could also be used in a malicious way to track vehicles.

Therefore, to guarantee drivers privacy, vehicles use pseudonym identities (or certificates) provided by a Public Key Infrastructure (PKI). During a trip, vehicles frequently change of certificates to make tracking much more difficult. They thus need to reload their certificates pool by requesting new ones to the PKI.

In this paper, we evaluate the performance of the PKI protocol regarding the reloading of certificates. We ran several tests while driving in order to quantify the number of certificates that can be reloaded from the PKI at different speeds. The obtained results show that 1) the end-to-end latency between a requesting vehicle and the PKI is non-negligible and 2) as speed increases, the number of successfully reloaded certificates decreases.

I. Introduction

Tomorrow's vehicles will communicate and cooperate by exchanging V2X messages in order to improve road safety and traffic efficiency.

The exchanged messages contain critical data like geographic position, speed, heading, etc. These data are directly linked to the driver's privacy because they enable vehicles tracking by eavesdropping the exchanged V2X messages. The protection of driver's privacy is a crucial element in vehicular networks as it is one of the important conditions for user's acceptance of connected vehicles. Therefore, in order to protect driver's privacy, vehicles use pseudonym certificates and frequently change them (i.e they change their digital identity) in such a way that it becomes much harder to track them. To enable frequent change of pseudonyms, vehicles use a pool of available pseudonyms (60 per week for C2C and maximum 100 per week for the European commission [START_REF]Intelligent Transport Systems (ITS); Security; Prestandardization study on pseudonym change management[END_REF]) that has to be reloaded over time. Consequently, vehicles need to communicate with the PKI to reload their pseudonym pool.

The question of when vehicles should request new pseudonyms can be raised. For instance, is it possible for a vehicle to reload pseudonyms while driving? If yes, how many pseudonyms can be reloaded and at which speed? What is the end-to-end latency to reload one pseudonym? These are the questions we address in this paper through real experimentations.

We ran several experimentations at different driving speeds and compare the performance of two communication profiles (with and without V2X security). The obtained results show that the number of pseudonyms successfully reloaded decreases when speed increases for both profiles. The end-to-end latency is quite high and remains constant versus speed.

The remaining sections of this paper is organized as follow: section II presents the related works. Section III briefly describes the PKI while section IV details our considered use case. Section V analyses the obtained results and section VI concludes this paper.

II. Related works

To the best of our knowledge only a few papers in the literature evaluate the performance of certificates reloading in real environment.

The latency of pseudonym certificates reloading can be divided in two parts: 1) the latency of the embedded operations (message creation, encryption and signature) and 2) the network latency. Authors of [START_REF] Moalla | Experimentation with the PRESERVE VSS and the Score@F system[END_REF] focus on the first part. They evaluate the performance of their security system using a set of performance indicators (signature/verification delays, packets signature/verification per second, pseudonym change latency). In this paper we go a step further by including the network latency and comparing two communication profiles.

In our previous work [START_REF] Haidar | On the performance evaluation of vehicular PKI protocol for V2X communications security[END_REF], we evaluated the performance of pseudonym certificates reloading by deploying an inlab testbed and comparing two communication profiles in terms of network overhead and end-to-end latency. In this paper we extend our evaluation by conducting real experimentations while driving to get additional performance indicators such as the number of pseudonyms that can be reloaded versus speed.

III. C-ITS PKI Architecture

The PKI is a set of entities that create, manage and distribute digital certificates. We implemented a PKI • ITS-S: end-entity of the system that requests certificates to the PKI and communicates with other endentities. Basically speaking, ITS-S requests pseudonym certificates to the AA. In order to prevent that an AA links differents pseudonyms to a same requesting ITS-S (thus breaking the ITS-S privacy), ITS-S are allowed to request only one pseudonym per request. That is, in order to get 10 new pseudonyms, an ITS-S must send 10 independant requests to the AA.

IV. Use case and experimentation setup

Our considered use case is depicted in figure 2. A vehicle is within range of a RSU that provides Internet access. The vehicle has very few pseudonyms and needs to refill its pool. It thus uses the RSU as gateway forwarder to reach the PKI and requests new pseudonyms.

We evaluate the amount of pseudonyms reloaded at the following speeds (km/h): 30, 50, 70 and 90. The "in-RSUrange" detection mechanism works as follow: the vehicle sends ICMPv6 Echo Requests (or "ping") to the PKI. When the vehicle enters the range of the RSU, it receives back an ICMPv6 Echo Reply from the PKI. This response indicates that the vehicle can communicate with the PKI through the RSU. The vehicle then starts requesting new pseudonyms. When the vehicle moves out of range of the RSU (i.e. when no more new pseudonyms are received), the requesting of new pseudonyms is stopped. We configured the vehicle in such a way that it continuously requests new pseudonyms as long as it is under the RSU radio coverage.

A. Communication profiles

According to ETSI standards [START_REF]Intelligent Transport Systems (ITS); Security; ITS communications security architecture and security management[END_REF], a vehicle can send a pseudonym request to the PKI using two communication profiles presented in figure 3: 1) TIG profile: TCP over IPv6 over G5 2) TI3G profile: TCP over IPv6 over GN6ASL over GN over G5 The difference between the two communication profiles is that TI3G provides security in the ITS G5 network between the vehicle and the RSU, whereas TIG does not provide any security. Figure 3 shows the message flow inside the vehicle to create and send a request for TIG (1→ 2→ 3→ 4→ 9) and TI3G ( 1→ 2→ 3→ 4→ 5→ 6→ 7→ 8→ 9) profiles.

On the other hand, a message sent using TI3G contains overhead and the RSU has to verify the request from the vehicle and sign the response from the PKI, which adds processing latency. In this paper, we evaluate these two communication profiles in term of end-to-end latency and number of pseudonyms reloaded at different speeds.

B. Protocol

Figure 4 shows the sequence diagram of the messages that are exchanged between the vehicle and the PKI through the RSU. The dotted arrows depict the additionnal operations generated when using TI3G communication profile. Basically speaking, the request is sent by the vehicle to the RSU that forwards it to the PKI and viceversa for the response. On the PKI side, a verification is Fig. 3. Data path within ETSI ITS-S architecture when requesting a pseudonym Fig. 4. Pseudonym (or AT) request/response using TIG/TI3G profile made to ensure that the requesting vehicle is registered and authorized to request new pseudonyms, as specified by the ETSI PKI protocol [START_REF]Intelligent Transport Systems (ITS); Security; Trust and privacy management[END_REF].

C. Test environment

The experimentations were conducted on the test track "Val d or" located at Versailles-Satory (France) and using our vehicle prototype. The vehicle is provided by RENAULT and the model is a MEGANE COUPE as depicted in figure 5. Figure 6 depicts a satellite view of the test track. We equipped the test track with one RSU depicted by the red antenna. The green line shows the radio coverage of the RSU (about 940 meters) whereas the red line shows the road segment which is out of the coverage of the RSU. We equipped our prototype vehicle with an on-board unit (OBU), a Samsung tablet, a Raspberry Pi that provides an in-vehicle Wi-Fi access point (to connect the tablet with the OBU), and two IEEE 802.11p antennas that we placed on the vehicle's roof. The list of hardwares/softwares and the configuration of the antennas are presented in table I.

V. Results and analysis

We conducted the experimentations as follow: we did two track laps for each speed (30, 50, 70 and 90 km/h) and for each communication profile (TIG and TI3G), resulting in a total of 16 laps. The main objectives of these experimentations are 1) to demonstrate the correct behaviour of the system in a realistic environment and 2) to evaluate the system performance in terms of number of pseudonyms reloaded and the end-to-end latency.

A. Number of pseudonyms reloaded

Figure 7 shows the number of pseudonyms successfully reloaded from the PKI versus the speed of the vehicle. The presented values are the sum of the two laps for each speed. We first observe that the number of pseudonyms successfully decreases as speed following an exponential shape. This result was expected as the faster the vehicle moves, the less time it remains under the RSU radio coverage, thus the less time it has to reload new pseudonyms. Second, we also observe that both curves have the same shape and the gap between them remains constant when speed increases. The gap is explained by the additional security processing required by the TI3G profile. Apart from the gap, the similar shape of both curves shows that the speed has the same impact on the performance, no matter what communication profile is used.

B. End-to-end latency

Figure 8 depicts the median end-to-end latency of the pseudonym request/response for both profiles versus speed. We observe that the end-to-end latency is always shorter when using TIG profile. This was also observed during our in-lab experimentations and is explained by the additional security processing required by the TI3G profile. Also, the median end-to-end values remains roughly constant versus speed, i.e. the speed has no impact on the time required to successfully reload a pseudonym. Finally, the median end-to-end latency is roughly half a second to reload one pseudonym, which is quite high, especially for a highly mobile network.

VI. Conclusion and future works

In this paper we evaluate the performance of pseudonyms reloading in a realistic environment. We com- The obtained results in terms of end-to-end latency show that requesting pseudonyms is always shorter when using TIG profile but remains still high (about half a second). As we expected, the speed has an impact on the number of reloaded pseudonyms as the faster a vehicle moves, the less time it has to reload new pseudonyms.

Future works include a performance evaluation and comparison with the integration of a Hardware Security Module (HSM) that accelerates cryptographic operations (and thus decrease the end-to-end latency), and evaluating the protocol under a congested ITS G5 network (high density of vehicles).

Fig. 1 .

 1 Fig. 1. PKI architecture

  Figure 1 depicts our PKI architecture that consists of the following entities: • Root Certificate Authority (RCA): RCA is the root of trust for all certificates within the PKI hierarchy. It delivers certificates to the Enrolment Authority (EA) and the Authorization Authority (AA) to authorize them to issue certificates to Intelligent Transport System-Station (ITS-S). • Enrolment Authority (EA): EA issues one Enrolment certificate (EC) per ITS-S. EC are considered as a proof of identity, thus used to identify and authenticate ITS-S within the PKI. • Authorization Authority (AA): AA issues pseudonyms to ITS-S that are used in V2X communications. • Distribution Center (DC): provides to ITS-S upto-date trust information necessary to validate that received information come from a legitimate and authorized ITS-S/PKI authority. • Operator: registers ITS-S and updates necessary information in the EA.

Fig. 2 .

 2 Fig. 2. Use case: pseudonyms (or AT) reloading through RSU

Fig. 5 .

 5 Fig. 5. In-vehicle equipments

Fig. 7 .

 7 Fig. 7. Number of pseudonyms (or AT) reloaded versus speed for both communication profiles

Fig. 8 .

 8 Fig. 8. Pseudonyms reload median end-to-end latency versus speed

Acknowledgments

This research work has been carried out in the framework of the Technological Research Institute SystemX, and therefore granted with public funds within the scope of the French Program Investissements d'avenir.