Jérôme J Mariette

Madalina Olteanu

Nathalie Vialaneix

Nathalie Villa-Vialaneix

Efficient interpretable variants of online SOM for large dissimilarity data

SOM, Sparse methods, Kernel, dissimilarity, K-PCA, Nyström Over the years, the self-organizing map (SOM) algorithm [START_REF] Kohonen | Self-Organizing Maps[END_REF] was proved to be a powerful and convenient tool for clustering and visualizing data [START_REF] Penn | Using self-organizing maps to visualize high-dimensional data[END_REF][START_REF] Pölzlbauer | Advanced visualization of self-organizing maps with vector fields[END_REF][START_REF] Sarlin | Cluster coloring of the self-organizing map: an information visualization perspective[END_REF][START_REF] Neme | Stylistics analysis and authorship attribution algorithms based on self-organizing maps[END_REF][START_REF] Yu | Visual query processing for efficient image retrieval using a SOM-based filter-refinement scheme[END_REF]. While the original algorithm had been designed for numerical vectors, the available data in the applications became more and more complex, being frequently too rich to be described by a fixed set of numerical attributes only. This is the case, for example, when data are described by relations between objects (individuals involved in a social network) or by measures of resemblance/dissemblance which are context specific (see [START_REF] Abbott | Optimal matching methods for historical sequences[END_REF][START_REF] Elzinga | Sequence similarity: a nonaligning technique[END_REF] for similarities between categorical sequences, [START_REF] Lozupone | Quantitative and qualitative β eiversity measures lead to different insights into factors that structure microbial communities[END_REF] for similarities between microbial diversity distributions, [START_REF] Yu | Representative distance: a new similarity measure for class discovery from gene expression data[END_REF] for similarities in gene expression data).

During the past twenty years, the SOM algorithm was extended to handle non numerical data. For example, SOM was adapted to categorical data analysis, by using a method similar to Multiple Correspondence Analysis in [START_REF] Cottrell | How to use the Kohonen algorithm to simultaneously analyse individuals in a survey[END_REF]. Another solution, called median SOM [START_REF] Kohohen | Self-organizing maps of symbol strings[END_REF], addressed the issue of data described by pairwise relations (similarities or dissimilarities): in this solution, the standard computation of the prototypes is replaced by an approximation within the original dataset. However, as prototypes are chosen among the data, their representation is very restrictive. In order to increase the flexibility of the prototypes, [START_REF] Conan-Guez | Fast algorithm and implementation of dissimilarity self-organizing maps[END_REF] proposed to represent a class by several prototypes, all chosen among the original dataset. But, this method seriously increases the computational time, while prototypes remain restricted to the original dataset and may generate possible sampling or sparsity issues.

A very different approach to handle relational data consists in relying on a (pseudo-)Euclidean framework, following the results of [START_REF] Aronszajn | Theory of reproducing kernels[END_REF] (for data described by a kernel) or of [START_REF] Goldfarb | A unified approach to pattern recognition[END_REF] (for dissimilarity data). This approach was developed in the framework of kernel SOM (see [START_REF] Donald | The kernel self organising map[END_REF] for the online version and [START_REF] Boulet | Batch kernel SOM and related Laplacian methods for social network analysis[END_REF] for the batch version), and in the framework of relational SOM (see [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF] for the online version and [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF] for the batch version). Kernel SOM and relational SOM are equivalent if the dissimilarity in relational SOM is the squared distance induced by the kernel. The key idea of this approach is to express prototypes as convex combinations of the images of the original data (x i ) i=1,...,n in a (pseudo-)Euclidean space in which the data are (implicitly) embedded by the kernel (or the dissimilarity): as stated in [START_REF] Rossi | How many dissimilarity/kernel self organizing map variants do we need?[END_REF][START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF], this solution yields several drawbacks due to the large dimensionality of the embedding space (which is equal to the number of observations, n). Firstly, the complexity (in n) is strongly increased and becomes at least quadratic. As stressed in [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF], algorithms will be slow for datasets with 10,000 observations and impossible to run on a normal computer for 100,000 input data. Secondly, the results are highly unstable: especially in the online (also called stochastic) version of the algorithm, two different runs of the method can provide very different results. Thirdly, one of the most important features of the SOM algorithm is lost: in standard numerical SOM, clusters are represented by a single prototype valued in the data space. These prototypes help to interpret the obtained clusters and thus the overall map organization. In kernel/relational SOM, prototypes are given as n coefficients that correspond to a resemblance with each of the n observations: they do not correspond themselves to an observation in the original data space and as such, prototypes are not much more informative than the clustering itself.

In conclusion, kernel and relational extensions of the standard SOM algorithm are hardly practicable when the dataset is large. This is due partly to the number of observations, but also to the dimensionality of the (embedded) data which is directly related to this number. To address this issue, strategies usually used to handle large datasets or datasets with a high dimensionality are useful and they can even be combined.

Review of methods for large datasets and high dimensional datasets

Different strategies were developed and are available in the literature to handle large datasets (when the number of observations is large) and high-dimensional datasets (when the dimension of the dataset is large). For large datasets, standard approaches include i) divide and conquer approaches [START_REF] Chu | Map-Reduce for machine learning on multicore[END_REF][START_REF] Chen | A split-and-conquer approach for analysis of extraordinarily large data[END_REF][START_REF] Del Rio | On the use of MapReduce for imbalanced big data using random forest[END_REF] in which data are split into several bits of data which are processed separately. The results are aggregated afterwards to obtain a final solution which is supposed to well approximate the solution that would have been obtained if the entire dataset had been processed at once; ii) subsampling methods [START_REF] Bǎdoiu | Approximate clustering via coresets[END_REF][START_REF] Yan | Fast approximate spectral clustering[END_REF][START_REF] Kleiner | A scalable bootstrap for massive data[END_REF][START_REF] Laptev | Early accurate results for advanced analytics on mapreduce[END_REF][START_REF] Meng | Scalable simple random sampling and stratified sampling[END_REF], which consist in using a restricted subset (usually carefully designed) of the original data, in order to approximate the solution that could have been obtained with the entire dataset and iii) online updates [START_REF] Saffari | On-line random forests[END_REF][START_REF] Denil | Consistency of online random forests[END_REF], in which the results are updated with sequential steps, each having a low computational cost.

A particular case of the subsampling strategy is the Nyström approximation [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF], which consists in sampling a small number of rows/columns in square matrices and in obtaining an approximation of its eigendecomposition at a very reduced computational cost. The eigendecomposition is even exact when the matrix is of low rank (when the size of the subsample is larger than the rank of the matrix). This method is frequently used for kernel and dissimilarity-based algorithms.

For high-dimensional data, the strategies are a bit different and include i) sparse approaches [START_REF] Hochking | The analysis and selection of variables in linear regression[END_REF][START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], in which a subset of the variables is selected to build the final predictive model. This subset can be obtained from sequential exploration (stepwise strategies), from approximation heuristics or by using a sparse penalty term within the model (LASSO); ii) dimension reduction (DR) techniques, that can be linear (PCA for instance or random projections as in [START_REF] Yu | Incremental semi-supervised clustering ensemble for high dimensional data clustering[END_REF]) or nonlinear [START_REF] Lee | Nonlinear Dimensionality Reduction, Information Science and Statistics[END_REF]. DR methods embed the data in a small dimensional space and are usually mainly used for visualization and exploratory analysis. However, if the embedding can be obtained at a low cost, it can be used as an approximation of the high-dimensional dataset on which more costly algorithms may be applied. SOM itself is a dimension reduction method but, as stressed before, the computational complexity of its kernel and relational versions is high. Finally, a particular case of DR techniques is model-based clustering methods, which use mixture distributions and embed the data in a low-dimensional subspace that is the best suited for clustering (see [START_REF] Bouveyron | Model-based clustering of highdimensional data: a review[END_REF] for a review).

Kernel/relational extensions for large datasets

Several extensions for kernel and relational data of the standard SOM algorithm, or of related kernel/relational algorithms (such as, e.g., k-means, LVQ, topographic maps...) have already been proposed in the literature. They use ideas coming from the strategies handling large and/or highdimensional datasets cited above. Most of them seek a simplified/sparse representation of the prototypes and/or a reduced computational time.

In the relational k-means framework, [START_REF] Rossi | Accelerating relational clustering algorithms with sparse prototype representation[END_REF] proposed a sparse extension of the batch algorithm: every prototype is represented by at most K (K fixed) observations by cluster, that are selected at each step of the algorithm. In the supervised framework, [START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF] used a similar strategy for batch LVQ, by selecting the most representative observations (with different methods to obtain them, including approximation heuristics and L 1 penalty) in every cluster and at each step of the algorithm. A similar method was used in [START_REF] Hofmann | Efficient approximations of robust soft learning vector quantization for non-vectorial data[END_REF], combined with the Nyström approximation of the LVQ algorithm, in order to obtain sparse prototypes at very low computational cost. The Nyström approximation was also used for obtaining faster versions of topographic mapping methods [START_REF] Gisbrecht | The Nyström approximation for relational generative topographic mappings[END_REF][START_REF] Zhu | Approximation techniques for clustering dissimilarity data[END_REF] and for reducing the computational cost of the clustering. Another subsampling strategy was used in a nonlinear (kernel) DR framework to allow processing large datasets, in [START_REF] Gisbrecht | Parametric nonlinear dimensionnality reduction using kernel t-SNE[END_REF].

However, these approaches do not lead to a simplified (and thus interpretable) representation of the prototypes. Furthermore, all of them are restricted to the batch framework and most of them are performed after each iteration of a batch algorithm, i.e., after all observations have been processed at least once. An alternative to these methods consists in splitting the data into several subsets on which independent algorithms are trained: in [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF], the complexity is reduced using iterative "patch clustering" that mixes "divide and conquer" and "online updates" methods. First, the data are split into B patches of size n B ( n, B fixed). A prototype-based clustering algorithm in batch version (neural gas or SOM) is then run on a patch P t . The resulting prototypes, which may be viewed as compressed representations of the data already seen, are then added as new data points to the next patch, P t+1 . Moreover, the full vector of coefficients is replaced by the Q closest input data (Q fixed). With this method, linear time and constant space representation are obtained but the sequential training may influence the final result since all observations are not processed evenly.

In the same line of thoughts, [START_REF] Mariette | Bagged kernel SOM[END_REF] propose a bagging approach for kernel SOM. Data are split into B subsamples of size n B ( n, B fixed), the online kernel SOM is trained on each subsample and, after training, the most representative Q observations are chosen for each prototype (Q fixed). Eventually, a final map is trained on the resulting most representative observations. In this method, parallel computing techniques can be used for reducing the computational time. However, the results of the B trained SOMs are not used as such but only to select the most representative observations in the dataset.

Contributions of the article

In the present article, we propose and compare two methods to obtain sparse prototypes and a reduced computational cost in the online SOM algorithm. The first one uses a reduction dimension in the embedding space, which can be efficiently performed with Nyström technique. This method combines ideas coming from the high dimension and the large data problems. However, it is not specific to the online setting. We thus compared it with another approach which takes advantage of the online framework to provide sparse prototypes at all iteration steps of the algorithm: the coefficients are interpreted similarly to an amount of mass and the most important observations are selected by using a fixed global probability mass, ν, such that only the largest coefficients summing to at most ν are kept. This second proposal also takes ideas from large data problems (online updates) and from high dimension (variable selection with sparsity).

For the sake of simplicity, most of the paper is written in the framework of kernel SOM but its extension to relational SOM is straightforward and briefly explained in Section 5. The rest of the paper is organized as follows: Section 2 recalls the online kernel SOM (online K-SOM). Section 3 describes the dimension reduction approach for online K-SOM while Section 4 presents the direct sparse version of online K-SOM. Comparisons and numerical experiments are reported in Section 6.

Kernel SOM (K-SOM)

This section describes the theoretical background of the online version of the SOM algorithm and its extension to data described by a kernel. In the following, we consider a set of observations (x i ) i=1,...,n which take values in an arbitrary space X . X is equipped with a kernel K : X × X → R that provides pairwise similarity between the observations,

K ij := K(x i , x j ). K is assumed symmetric (K ij = K ji ) and positive (∀ N ∈ N, ∀ (α i ) i=1,...,N ⊂ R, ∀ (x i ) i=1,...,N ⊂ X , i,i α i α i K ii ≥ 0).
According to [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], K is the dot product in a uniquely defined Hilbert space (H, ., . ) of the images of x i and x j by a uniquely defined feature map φ : X → H:

K ij = φ(x i ), φ(x j ) .
The SOM algorithm aims at mapping the input data onto a low dimensional grid (usually a two-dimensional rectangle), composed of U units, each of them described by a prototype p u , u = 1, . . . , U . The units are related together by a neighborhood relationship H, expressed as a function of the distance between the units on the grid, d(u, u ), H : R → R. Classically, H is chosen such that H(0) = 1, lim x→+∞ H(x) = 0 and is decreased during the training. Also, the distance on the grid, d, may be chosen as the length of the shortest path between the units or as the Euclidean distance between the coordinates of the units positioned at (1, 1), (1, 2), . . . , (m, m ) (where m × m = U ).

In standard (numerical) SOM, the data take values in X = R d and the kernel is the standard dot product in this space K ij = x i x j . In this case, the prototypes are also valued in R d . In kernel SOM, X is arbitrary and prototypes are not easily defined in this space, which may not be Euclidean or equipped with standard operations such as the sum. To allow for a more flexible representation of the prototypes (i.e., a representation which is not restricted to the data already observed, (x i ) i ), the implicit feature space (H, ., . ) is used and the algorithm is re-defined in this space. The prototypes are expressed as convex combinations of the images by φ of the original data:

p u = n i=1 β ui φ(x i ), β ui ≥ 0 and i β ui = 1.
As said before, the feature map is uniquely defined from the kernel K. However, it is usually not explicitly given and thus, all calculations are based on the values of the coefficients

β u = (β u1 , . . . , β un ) ∈ R n only.
The kernel version of online SOM [START_REF] Donald | The kernel self organising map[END_REF][START_REF] Rossi | Accelerating relational clustering algorithms with sparse prototype representation[END_REF] is thus directly derived from the computations of the standard numerical SOM performed in the feature space H: the norm and the dot product in H can be obtained using the kernel, without the need of explicitly knowing H or φ (this is the so-called "kernel trick"). Online K-SOM is provided in Algorithm 1. In this algorithm, µ(t) usually vanishes at the rate 1 t and the final clustering is defined as (C u ) u=1,...,U , where

C u = {x i : f (x i ) = u} with f := f T .
The issue with kernel SOM is that, when n is large, the total number of coefficients β ui to learn is equal to n×U , which yields a complexity of O(n 2 U ) (for one iteration) instead of O(dU ) + O(ndU ) for the standard numerical SOM in R d . Hence, this algorithm cannot be used to analyze datasets with more than a few thousands observations. Also, the representation of the prototypes is so flexible that the results are highly unstable with different final clusterings and different prototypes for each run of the algorithm. Finally, as stated in [START_REF] Rossi | How many dissimilarity/kernel self organizing map variants do we need?[END_REF], one of the main challenges of this kind of algorithm is that the easiness in interpreting the resulting map is lost: in standard SOM, prototypes are easily interpreted because they are described by values of well known variables which are the variables also describing the observations in the input dataset. In kernel SOM, the prototypes are characterized by their proximity to all individuals in the dataset. Hence, interpreting the clustering requires to be able to understand the relationships of the prototypes with all individuals, which is probably as difficult as analyzing the entire dataset directly.

Algorithm 1 Online kernel SOM

1: For all u = 1, . . . , U and i = 1, . . . , n, initialize (β 0 ui ) u,i such that β 0 ui ≥ 0 and n i β 0 ui = 1. 2: for t = 1, . . . , T do

3:

Randomly choose an input x i 4:

Assignment step: find the unit of the prototype closest to φ(x i ) ∈ H

f t (x i ) ← arg min u=1,...,U p t-1 u -φ(x i ) 2
which is equivalent to minimize, over u, j,j β t-1 uj β t-1 uj K jj -2 j β t-1 uj K ij ;

5:

Representation step: ∀ u = 1, . . . , U ,

β t u ← β t-1 u + µ(t)H t (d(f t (x i ), u)) 1 i -β t-1 u
where 1 i is a vector with a single non null coefficient at the ith position, equal to one. 6: end for

A dimension reduction technique for K-SOM

This section describes a first approach to reduce the dimensionality of the prototypes, thus improving the computational complexity as well as the interpretability of the prototypes. Our proposal is to rely on a preprocessing of the data based on PCA in the feature space (Kernel PCA, denoted by K-PCA and described in Section 3.1) and then to express the SOM algorithm in the subspace of the feature space H which is spanned by the first eigenvectors of the K-PCA (Section 3.2). Finally, the computational complexity of the approach can be further reduced by performing the K-PCA thanks to the Nyström approximation (Section 3.3).

Kernel PCA (K-PCA)

K-PCA, introduced in [START_REF] Schölkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF], is a PCA analysis performed in the feature space (H, ., . ) induced by the kernel. A centered data matrix is first com-puted using:

K ij := φ(x i ) - 1 n n l=1 φ(x l ), φ(x j ) - 1 n n l=1 φ(x l ) = K ij - 1 n n l=1 (K il + K jl ) + 1 n 2 n l,l =1 K ll , (1)
which yields to the modified centered kernel

K = K -1 n I T n K -1 n KI n + 1 n 2 I T n KI n ,
in which I n is the vector with n entries equal to 1. The eigenvectors (α k ) k=1,...,n ∈ R n and the corresponding eigenvalues (λ k ) k=1,...,n are obtained from this centered matrix by solving the following eigenvalue problem:

Kα = λα. ( 2
)
This problem is equivalent to finding eigenvectors, in the feature space H, of the covariance matrix of the (centered) images of the original data by φ (the feature map associated to the centered kernel K). These vectors, (a k ) k=1,...,n ∈ H lie in the span of {φ(x i )} i=1,...,n and can be expressed as:

a k = n i=1 α ki φ(x i )
where α k = (α ki ) i=1,...,n . a k = (a ki ) i=1,...,n are orthonormal in H:

∀ k, k , a k , a k = α k Kα k = δ kk with δ kk = 0 if k = k 1 otherwise .
The principal components are the coordinates of the projections of the images of the original data, (φ(x i )) i onto the eigenvectors (a k ) k≥1 which can be expressed as:

a k , φ(x i ) = n j=1 α kj K ji = K i. α k ,
where K i. is the i-th row of the kernel K. According to Equation (2), we thus obtain that a k , φ(x i ) = λ k α ki . The original data are projected on these axes with:

P a k (φ(x i )) = a k , φ(x i ) a k = (λ k α ki )a k .
The result of K-PCA can be used to approximate the data in a reduced space by selecting p axes (a k ) k=1,...,p in the feature space H (with p n), on which the data can be projected.

K-PCA SOM

We suppose in this section that the kernel K is centered. Otherwise, without loss of generality, the centering procedure described in Equation ( 1) can be applied.

The main idea developed in this section is to define the U prototypes of the SOM in A = Span{a 1 , . . . , a p } instead of the entire feature space. A prototype can thus be written as:

p u = p k=1
β uk a k with β uk ≥ 0 and k β uk = 1 and the two steps of the SOM algorithm can thus be rewritten as:

• assignment step: when the observation x i is picked at random, it is affected to:

f (x i ) := arg min u p u -φ(x i ) 2
in which:

p u -φ(x i ) 2 = p k,k =1 β uk β uk a k , a k -2 p k=1 β uk a k , φ(x i ) + φ(x i ) 2 = p k,k =1 β uk β uk a k , a k -2 p k=1 β uk λ k α ki + φ(x i ) 2 = β u β u -2β u Λα .i + φ(x i ) 2 , (3)
where β u = (β uk ) k=1,...,p , Λ = Diag (λ 1 , . . . , λ p ) and α .i is the ith column of α = [α 1 , . . . , α p ] . This step is thus equivalent to minimize β u β u -2β u Λα .i over u ∈ {1, . . . , U } and is also equivalent to minimize p u -P A (φ(x i )) 2 over u. This result was expected since, by the definition of P A (φ(x i )), we have that

p u -φ(x i ) 2 = p u -P A (φ(x i )) 2 + P A (φ(x i )) -φ(x i ) 2
, in which the second term does not depend on u;

• representation step: the gradient descent like step is given by

p u = p u + µH(d(f (x i ), u)) (P A (φ(x i )) -p u ) = k β uk a k + µH(d(f (x i ), u)) k (λ k α ki )a k - k β uk a k ,
which is equivalent to update the coefficients as:

β u = β u + µH(d(f (x i ), u)) (Λα .i -β u ) .
This approach is simply the standard (numerical) SOM with entries the n×p matrix α Λ. The computational complexity of this approach is thus reduced from O(n 2 U T ) (online K-SOM, T iterations) to O(pU T ) + O(npU ) (online numeric SOM in R p , cost of the T iterations and cost of the final clustering computation), once the K-PCA is given. Hence, since T is generally of the order of n, this gives a linear complexity for the algorithm. However, K-PCA itself can have a large complexity. This issue is addressed in the next section.

The interpretation of the prototypes is done similarly as for a standard PCA: the axes a 1 , ..., a p of the prototypes are interpreted with respect to the observations (x i ) i which contribute most to their definition. Then, prototypes are interpreted by means of their coefficients on these axis.

K-PCA from Nyström approximation

Kernel-based methods, and especially K-PCA, do not scale well when the number of observations, n, is large. For instance, the computational complexity of the eigendecomposition of K is O(n 3 ). An effective approach for improving the scalability of kernel methods is the Nyström approximation [START_REF] Williams | Using the Nyström method to speed up kernel machines[END_REF] and a similar technique is used in [START_REF] Zhu | Approximation techniques for clustering dissimilarity data[END_REF] for improving the computational cost of topographic maps for dissimilarity data. In the latter reference, the authors use the Nyström approximation to reduce the computational cost of

p u -φ(x i ) 2 from O(n 2 ) to O(m 2 n),
in which m is a number of observations taken at random within the original dataset. m is supposed to be small compared to n and close to the kernel rank.

In the approach introduced in Section 3.2 and unlike article [START_REF] Zhu | Approximation techniques for clustering dissimilarity data[END_REF], it is the pre-processing step which is addressed by the Nyström approximation, while the computational cost of calculating the clustering of the map is handled by the use of the numerical version of the algorithm based on the K-PCA pre-processing. More precisely, the eigendecomposition of a kernel matrix K (which is supposed to be centered) is approximated by selecting m observations, T m among (x i ) i=1,...,n , and by using the eigendecomposition of the reduced matrix K (m) = (K(x i , x j )) i,j∈Tm . In practice, the selected observations T m are chosen at random, although more efficient sampling techniques such as the ones described and evaluated in [START_REF] Kumar | Sampling techniques for the Nyström method[END_REF] could also be used.

If the eigenvalues and the (orthonormal) eigenvectors of K (m) are denoted by (λ

(m) k ) k and (v (m)
k ) k respectively, then the eigenvalues and (orthonormal) eigenvectors of K are given by

λ k n m λ (m) k and v ki m n 1 λ (m) k K (n,m) i. v (m) k , with K (n,m) i.
the i-th row of the matrix K (n,m) = (K(x j , x j )) j=1,...,n, j ∈Tm . If the rank of K (m) is equal to the rank of the original matrix K, the approximation even becomes an equality. Then, assuming that the kernel K (which is supposed to be centered) is known or at least that any of the pairs K(x i , x j ) (i, j = 1, . . . , n) can be computed at low cost, the K-PCA requires to obtain the entries (λ k α ki ) k=1,...,p ∈ R p for all i = 1, . . . , n, where (α k ) k are the eigenvectors of K which are orthogonal with respect to the norm induced by the kernel. We can easily show that

α k = v k √ λ k (4)
because, by definition, as v k are eigenvectors of K, so are α k and, in addition, using the equality of Equation ( 4), we have that:

α k Kα k = 1 √ λ k √ λ k v T k Kv k = 1 √ λ k √ λ k v T k λ k v k = δ kk .
Therefore, K-PCA can be computed with entries the rows of the n×p matrix α Λ with ∀ i = 1, . . . , n and ∀ k = 1, . . . , p,

λ k α ki = λ k v ki = 1 λ (m) k K (n,m) i. v (m) k = K (n,m) i. α (m) k (5)
with α

(m) k = v (m) k λ (m) k
. This simplified K-PCA is a good approximation of the K-PCA with kernel K. The complexity of the preprocessing is reduced from O(n 3 ) (standard K-PCA) to O(m 3 )+O(nm 2 ) (respectively, K-PCA based on the reduced kernel and approximation). The complete algorithm is provided in Algorithm 2.

Algorithm 2 Online K-PCA SOM 1: Nyström approximation of K-PCA 2: Select at random m observations T m = {x i(1) , . . . , x i(m) } from the original dataset 3: Compute the first p eigenvalues and orthonormal eigenvectors of K (m) , (λ

(m) k ) k=1,...,p , (v (m)
k ) k=1,...,p and obtain, for k = 1, . . . , p, α Randomly choose an input i ∈ {1, . . . , n}

(m) k = v k √ λ k 4: Compute K (n,m) α (m) k k=1,...,p , which is an n × p matrix B = (b ik ) i=1,...,n, k=1,...,p 5: K-PCA SOM 6: Initialize prototypes randomly: ∀ u = 1, . . . , U , β 0 u ∈ [0, 1]

9:

Assignment step: find the unit of the prototype closest to i-th row of B, b i :

f t (x i ) ← arg min u=1,...,U β t-1 u -b i 2 R p 10:
Representation step: ∀ u = 1, . . . , U ,

β t u ← β t-1 u + µ(t)H t (d(f t (x i ), u)) b i -β t-1 u 11: end for 4. Direct sparse K-SOM
In this section, we present an alternative approach to obtain sparse prototypes while taking advantage of the online updates of the stochastic version of the K-SOM algorithm. Here, unlike in K-PCA SOM, the prototypes are directly written as convex combinations of the observations, but, in this case, they are restricted to the input data already fed to the algorithm and, more particularly, to the most important of them. This algorithm has already been described in the framework of dissimilarity data and tested on restricted datasets in [START_REF] Olteanu | Sparse online self-organizing maps for large relational data[END_REF].

Description of the approach

To ensure sparsity through the algorithm, the prototypes are initialized at random among the input data. Thus, the method first selects at random U observations that are used as initial values for the U prototypes. Then, at a given step t of the algorithm, a prototype is written as a convex combination of the most important past observations: if I u (t -1) denotes the set of the most important observations selected for prototype p u before step t, p u writes p u = j∈Iu(t) β uj φ(x i ). Finally, the distance in the feature space H between a new input, x i selected at random, and p u is given by:

φ(x i )-p u 2 = j,j ∈Iu(t-1)
β uj β uj K(x j , x j )-2 j∈Iu(t-1)

β uj K(x j , x i )+K(x i , x i ).
(6) In order to maintain prototypes as sparse combinations of the input data, they are periodically updated and the most important coefficients only are kept. The update instants may be performed throughout the iterations using various strategies: for instance, they can be uniformly distributed during the learning process or distributed according to some geometric distribution. The parameter of the geometric distribution may be fixed, during the whole training, or varying (in ascending or descending fashion) with the iterations. We ran several simulations on various data sets, using these four scenarios (results not shown). The results were globally similar, with a slight advantage for the "ascending random" strategy. In this case, the parameter of the geometric distribution was γ t = 1-µ(t) κ , where κ is a constant to be set. As suggested by [START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF] for a batch sparse LVQ method, sparsity could be achieved by selecting the first Q most important coefficients, where Q is a fixed integer. However, in order to allow for more flexibility in the expression of the prototypes, the most important coefficients are selected according to their value, by fixing a threshold: let 0 < ν ≤ 1 be the selected threshold. At time step t at which an update occurs and for every u = 1, ..., U , the coefficients of the prototype p u are first ordered in descending order, β u,(1) ≥ ... ≥ β u,( Iu(t)) . Then, the integer N u such that

N u = arg min k=1,..., Iu(t) k i=1 β u,(i) ≥ ν is introduced.
The most important coefficients are finally updated as follows

β u,(i) = β u,(i) Nu j=1 β u,(j) if (i) ≤ N u 0 if (i) > N u
, and I u (t) is updated accordingly afterwards by keeping the observations that correspond to non zero coefficients only.

The sparse relational SOM algorithm is entirely described in Algorithm 3.

Contrary to K-PCA SOM, in this version, the prototypes might directly be interpreted through the observations that are used in their representation. Also, the sparsity is updated during the training and the induced dimension reduction is thus not constrained to the efficiency of a given dimension reduction technique such as K-PCA. However, due to the sparse representation step, the algorithm can be computationally more expensive than K-PCA SOM and the amount of information preserved in the sparse representation is not as well controlled as in K-PCA projection. Finally, the complexity of the method, for T iterations, is not easily obtained: if a total mass equal to ν results in no more than Q(ν) observations for every prototype at each step, then the global complexity of the method has an upper bound of order

O ((Q(ν)) 2 U T ) + O (n(Q(ν)
2 ) (respectively for the iterations and the final clustering computation). However, the relation between ν and Q(ν) is hard to know in advance and can depend on the dataset distribution and on the training itself.

Variants of the sparse updates

In the line of [START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF], variants of the update step (step 15 in Algorithm 3) can be introduced. More precisely, some of the approaches introduced in [START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF] can be used almost directly in the online framework. The main difference is Algorithm 3 Sparse online K-SOM 1: For all u = 1, . . . , U , initialize p 0 u among U randomly selected observations in (x i ) i . Initialize I u (0) = {i(u)}, with i(u) ∈ {1, . . . , n} for all u and β 0 u = 1. 2: for t = 1, . . . , T do

3:

Randomly choose an input x i , i ∈ {1, ..., n}.

4:

Assignment step: find the unit with the prototype closest to φ(x i ):

f t (x i ) ← arg min u=1,...,U   (β t-1 u ) K Iu(t-1) β t-1 u -2 j∈Iu(t-1) β t-1 uj K(x j , x i )  
where K Iu(t-1) = (K(x j , x j ))) j,j ∈Iu(t-1) .

5:

Representation step: ∀ u = 1, . . . , U

6:

if i ∈ I u (t -1), then 7:

β t u ← β t-1 u + µ(t)H t (d(f t (x i ), u)) 1 i -β t-1 u 8: I u (t) = I u (t -1) 9:
else if i / ∈ I u (t -1), then 10:

β t u ← 1 -µ(t)H t (d(f t (x i ), u)) β t-1 u , 0 +µ(t)H t (d(f t (x i ), u))(0, ..., 0
Iu(t-1) , 1) 11:
I u (t) = I u (t -1) ∪ {i}.

12:

end if

13:

Sparse representation 14:

if t is an update instant then 15:

Sparsely update the prototypes: ∀ u = 1, . . . , U and with β t u,(1) ≥ ... ≥ β t u,( Iu(t)) , set

N t,u = arg min k=1,..., Iu(t) k i=1 β t u,(i) ≥ ν and ∀ (i), st i ∈ I u (t), β t u,(i) ←    β t u,(i) N t,u j=1 β t u,(j) if (i) ≤ N t,u 0 if (i) > N t,u
and I u (t) ← {i : (i) ≤ N t,u }.

16:

end if 17: end for that the number of iterations in the online version is much larger than in the batch version. In practice, this implies that the sparse approximation must have a low computational complexity because it is performed several times during the training. We thus restrict ourselves to the following two proposals which are the less computationally costly and can be used instead of step 15 in Algorithm 3:

• simple numerical heuristic approximation (N -num): instead of selecting the first coefficients which sum to a given amount of total mass, one can simply select the first N coefficients for a given N . This changes step 15 into

β t u,(i) = β t u,(i) N j=1 β t u,(j) if (i) ≤ N 0 if (i) > N ,
• simple geometric heuristic approximation (N -geom): as an alternative, one can search through the observations in I u the first N observations closest to prototype p u . Additionally to what is described in [START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF], we set the coefficients β u,i afterwards, in accordance to the new observations used to represent the prototype p u . Again, we used the previous values of the coefficients as a priori for the new ones.

Step 15 re-writes:

1. Ĩu (t) is the set of the first N observations, x i , in

I u (t) which min- imize p u -φ(x i ) 2 2. β t u,i ← β t u,i j∈ Ĩu(t) β t u,j if (i) ∈ Ĩu (t) 0 if i / ∈ Ĩu (t) , 3. I u (t) ← Ĩu (t)
Both approaches described above lead to a fixed number of observations for representing all the prototypes, contrary to what is proposed in Section 4.1. The advantage of our proposal is that prototypes may be represented by a number of observations which varies with the density of the data on the map: when a given cluster is isolated from its neighbors, this should yield to a more accurate way to represent the cluster. Also, for all three methods, the sparse approximation is maintained during all the training, as long as the sparse updates are performed frequently enough. This yields to a reduced computational cost of the assignment step (Step 4 in Algorithm 3) and, to a lesser extent, of the representation step (Step 10 in Algorithm 3).

As already stated in the previous section, the computational gain is not easily obtained and is a trade-off between the sparsity constraint imposed on the data (which itself depends on the total mass ν or on the number of observations N and on the number of update instants) and the complexity of the sparse representation step itself (Step 15 in Algorithm 3). Further discussions about computational time and efficiency of the different alternatives are provided in the experiments (Section 6).

The case of dissimilarity data

For the sake of clarity, the entire paper has been written in the framework of kernel data. However, its extension to dissimilarity data is almost straightforward, even in the case where the dissimilarity is not Euclidean. The present section briefly describes the modifications that must be made in this setting.

In the sequel, ∆ will denote a dissimilarity matrix with entries δ(x i , x j ), the dissimilarity between the observations x i and x j . ∆ is supposed to have some very basic properties, such as the positiveness of all entries (δ(x i , x j ) ≥ 0), symmetry (δ(x i , x j ) = δ(x j , x i )) and a null diagonal δ(x i , x i ) = 0. However, it may not necessarily be Euclidean (see, for instance [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF] for a discussion on the additional requirements which make this dissimilarity Euclidean). The extension of SOM to this case is called "relational SOM" and was described in [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF][START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF]. It is very similar to the kernel case, except that the assignment step of Algorithm 1 is replaced by

f t (x i ) ← arg min u=1,...,U β t-1 u ∆ i - 1 2 β t-1 u ∆β t-1 u , (7)
where ∆ i is the i-th row in matrix ∆. This equation is justified by the pseudo-Euclidean framework described in [START_REF] Goldfarb | A unified approach to pattern recognition[END_REF]: ∆ can be expressed as the difference between dot products of images of the original data by two mapping functions, ψ 1 and ψ 2 into two Euclidean spaces, E 1 and E 2 :

δ(x i , x j ) = ψ 1 (x i ) -ψ 1 (x j ) E 1 -ψ 2 (x i ) -ψ 2 (x j ) E 2 .
In this framework, Equation ( 7) is the exact distance calculation in the space E 1 ⊗ E 2 equipped with the pseudo dot product ., .

E 1 -., . E 2 . Moreover, if
∆ is Euclidean, then the similarity defined by

K(x i , x j ) = - 1 2 δ(x i , x j ) - 1 n n k=1 (δ(x i , x k ) + δ(x k , x j )) + 1 n 2 n k, k =1 δ(x k , x k ) , (8)
as suggested in [START_REF] Lee | Nonlinear Dimensionality Reduction, Information Science and Statistics[END_REF], is a kernel and the distances between observations in the feature space induced by this kernel are given by ∆. Thus, Equation ( 7) is exactly equivalent to the assignment step in standard SOM for the kernel K defined by Equation [START_REF] Elzinga | Sequence similarity: a nonaligning technique[END_REF].

Hence, sparse K-SOM can straightforwardly be extended to dissimilarity data using the assignment step restricted to selected observations in the sparse representation as given in Equation [START_REF] Abbott | Optimal matching methods for historical sequences[END_REF]. For SOM based on dimension reduction, the extension is also easy to obtain. For a non-Euclidean dissimilarity, ∆, the eigendecomposition of the similarity K defined by Equation ( 8) leads to positive and negative eigenvalues (because, in this case, K is not a kernel anymore). If n + ≤ n denotes the number of positive eigenvalues, the projection of the original data by K-PCA should be restricted to the first p components associated with the p largest eigenvalues such that p ≤ n + . This restriction is equivalent to performing K-PCA on a kernel pre-processed with the standard "clip" approach as suggested in [START_REF] Chen | Similarity-based classification: concepts and algorithm[END_REF].

Experimental results

Methodology

The present section is devoted to the evaluation of the proposed methods on several datasets and from several points of view. First, in Section 6.2, the methods presented in this article are assessed on both small and large datasets. This allowed us to compare the performances and computational efficiency of these methods with those of the standard K-SOM algorithm, in an extensive way. In the same section, the sparse approach is also compared, on the two small datasets, to the simple numerical and geometric heuristic approximation methods described in Section 4.2. Finally, in Section 6.3, the two new methods are used to produce a map on a very large dataset. In Section 6.3, computational times are reported and compared, the influence of the size of the prototypes is assessed and the Nyström approximation is tested. This example also serves as a use case to show how the results can be interpreted.

The datasets were chosen in the framework of this article, i.e., most of them are not standard numeric datasets. They include graphs, genomic sequences and categorical time series to illustrate the approach on a wide range of application and type of (dis)similarities. For every method and every set of parameters that we tested, the experiments were performed with 100 maps, so as to evaluate the variability of our conclusions.

Finally, performances are reported in terms of:

• standard quality measures used to evaluate SOM (see [START_REF] Pölzlbauer | Survey and comparison of quality measures for selforganizing maps[END_REF] for a description of quality measures in SOM): i) the quantization error (QE),

U u=1 i: x i ∈Cu φ(x i ) -p u 2 with φ(x i ) -p u 2
given by the kernel as in Equation (3) (K-PCA SOM) or as in Equation ( 6) (sparse K-SOM). This measure is a clustering quality measure, disregarding the map topology; ii) the topographic error (TE) which is the simplest of the topographic preservation measures: it computes the ratio over (x i ) i of the second best matching units that are in the direct neighborhood of the best matching units on the map. Since for the K-PCA SOM version, the distances φ(x i ) -p u 2 depend on the number of selected axes (the smaller the number of axes, the smaller the distances between observations in the corresponding projected subspace), we have normalized all QE by the average of the squared distances between all pairs of observations in the feature space, φ(x i ) -φ(x j ) 2 (for i = j) to leverage the impact of the feature space metric itself. However, the performances remain difficult to compare directly even with such a normalization, as shown in the results of the simulations. Thus, the average intra-cluster inertia (ICI) is also provided. This measure is equal to the average over the units u of i:

x i ∈Cu φ(x i ) -G Cu 2 in which G Cu = 1 Cu i: x i ∈Cu φ(x i
) is the center of gravity of C u . We can easily show that it is equal to

1 Cu i: x i ∈Cu K(x i , x i )-1 ( Cu) 2 i,i : x i ,x i ∈Cu K(x i , x i )
, for kernel data and to 1 2( Cu) 2 i,i : x i ,x i ∈Cu δ(x i , x j ) for dissimilarity data (see Appendix Appendix A for a proof);

• quality measures which take advantage of a prior external information: i) if an a priori classification of the observations is given, we use the normalized mutual information (NMI, [START_REF] Danon | Comparing community structure identification[END_REF]) between this a priori classification and the clustering induced by the SOM map. A perfect match between the two would correspond to an NMI equal to 1 whereas inde-pendent classifications would provide an NMI equal to 0; ii) when the studied dataset is a graph and no a priori classification is given, the quality of the clustering of the SOM map can nevertheless be evaluated by computing its modularity [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. A high modularity indicates a good clustering with respect to the graph structure;

• quality measures which aim at providing an indication of how much the results are stable between two runs of the algorithm. To do so, we provide the average NMI between any pair of clustering results obtained from the 100 runs, for a given method and with a given set of parameters. In the sequel, this measure is called stability.

Evaluation of the different sparse approaches on various datasets 6.2.1. Description of the datasets

This section aims at evaluating and comparing the two approaches described in Sections 3 and 4 on various datasets: two small datasets, for an intensive analysis, and three larger datasets, for a better evaluation of the gain in computational time.

More precisely, the two small datasets used for the experiments presented in this section are:

• a graph that gives the frequencies of co-appearance in a same chapter between characters in the novel "Les Misérables" from the French author Victor Hugo. The 77 vertices of the graph are the characters and the 254 edges are weighted by the number of co-appearances in a same chapter. For this graph, a kernel obtained from the shortestpath lengths is computed1 . The similarity described in Section 5 is used to perform K-PCA SOM: in this case, the similarity K defined in Equation ( 8) is not positive (only the first 67 components are positive);

• a dataset that contains 465 input data issued from ten unbalanced sampled species of Amazonian butterflies. This dataset was previously used by [START_REF] Hebert | Ten species in one: DNA barcoding reveals cryptic species in the neotropical skipper butterfly astraptes fulgerator[END_REF] to demonstrate the synergy between DNA barcoding and morphological-diversity studies. The notion of DNA barcoding comprises a wide family of molecular and bioinformatics methods aimed at identifying biological specimens and assigning them to a species. DNA barcoding data are composed of sequences of nucleotides, i.e., sequences of "A", "C", "G", "T" letters in high dimension (hundreds or thousands of sites). Specific distances and dissimilarities such as the Kimura-2P [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF] are usually computed and used in the experiments. A similarity was obtained from the Kimura-2P dissimilarity as described in Equation ( 8) to perform K-PCA. Again, this similarity is not positive. Only 246 eigenvalues are positive in the eigendecomposition of K.

These two datasets will be denoted, respectively, by "lesmis" and "astraptes" in the sequel. Additionally, three larger datasets are also used for comparison:

• a graph whose edges model hyper-links between blogs on US politics, recorded in 2005 by [START_REF] Adamic | The political blogosphere and the 2004 us election: divided they blog[END_REF] 2 . The 1,222 vertices of the graph represent the political blogs and are labeled by the political leaning (liberal or conservative). The graph contains 19,089 edges. The similarity used for this graph is again obtained from Equation ( 8), using the shortest-path lengths as the original dissimilarity measure. The undirected version of the shortest-path length was used to provide a symmetric dissimilarity, even though the original graph is directed. Only the first 779 eigenvalues of this similarity are positive;

• a DNA dataset similar to the "astraptes" data, except with a larger size. This dataset contains 614 sites of the CoI locus for 2,036 samples issued from the cowries family. The data were introduced in the DNA barcoding context in [START_REF] Meyer | DNA barcoding: error rates based on comprehensive sampling[END_REF]. In order to assess the ability for clustering of the proposed algorithms, the species with very few observations were removed. The final dataset used for simulations contained 1,414 samples and 47 species, the number of observations per species varying from 11 to 140 observations. The species were used as a priori classes for computing the NMI. Here also, the Kimura-2P dissimilarity [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF] was used for deriving the similarity matrix, for which only the first 476 eigenvalues were positive.

• a (standard) numerical dataset related to red variants of the Portuguese "Vinho Verde" wine [START_REF] Cortez | Modeling wine preferences by data mining from physicochemical properties[END_REF] 3 . The dataset contains 4 898 observations of 12 numeric variables based on physicochemical tests, such as the pH, the sulphates or the residual sugar. Additionally, the quality (a score between 0 and 10) of the wines is provided, which is used as an a priori class to compute the quality measure (NMI) described in Section 6.1 (and thus it is not used to compute the similarity). To stuck to the framework of the article, the Gaussian kernel was used to obtain a measure of similarity between wines: K ij = e -σ x i -x j 2 with σ equal to the median of

1 x i -x j 2 i<j .
In the sequel, these datasets will be refered as "polblogs", "cowrie" and "wines", respectively.

Evaluation of K-PCA SOM and direct sparse K-SOM

More precisely, the following methods are compared on these five datasets:

• a standard kernel SOM, K-SOM (or relational SOM if the dissimilarity is not Euclidean);

• K-PCA SOM, as described in Section 3. The datasets used in this section are not large enough to allow a relevant use of the Nyström approximation of the K-PCA;

• direct sparse K-SOM (or relational SOM) as described in Section 4, and denoted by sparse K-SOM in the sequel.

All methods were trained for a 5×5 grid ("lesmis"), a 8×8 grid ("astraptes") and a 10 × 10 grid ("polblogs", "cowrie" and "wines") with respectively 500 ("lesmis"), 2500 ("astraptes"), 6000 ("polblogs"), 7000 ("cowrie") and 8000 ("wines") iterations. These choices approximately correspond to default parameters in SOMbrero. All grids were equipped with a piecewise linear neighborhood with the distance between units calculated as the Euclidean distance between the unit coordinates in N * 2 . Furthermore, different sets of parameters, corresponding to different levels of sparsity, were tested:

• for K-PCA SOM, the only parameter to set was the dimension of the projection, p. This parameter was chosen with respect to a minimal ratio of entropy preserved in the projection and this ratio itself was varied in {20%, 40%, 60%, 80%};

• for the direct sparse K-SOM, two parameters had to be set: the first one is the mass parameter ν which was varied in {90%, 99%} for the two smallest datasets ("lesmis" and "astraptes") and in {95%, 99%} for the three largest datasets ("polblogs", "cowrie" and "wines") to ensure sparse results. The second parameter calibrated the update instants: random ascending updates were performed with κ ∈ {1, 50} for both datasets.

Tables 1-5 provide the results obtained over 100 maps (mean and standard deviation) for different quality criteria. For the smallest datasets, we focused on providing many different measures of the quality of the resulting map, as compared to the one obtained by the standard K-SOM: normalized QE, ICI, TE and modularity ("lesmis") or NMI ("astraptes"), as described in Section 6.1, are given and the last column of the results provides the final dimension of the prototypes (either the exact dimension for K-PCA SOM or the average over all prototypes for the sparse K-SOM). The computational time is not reported because, for these small datasets, it is not relevant.

For the three largest datasets, the computational time is reported but the results are less exhaustive on quality criteria (only the most important ones are reported: ICI, TE, modularity/NMI and stability). Also, for K-PCA SOM, only the clustering time is given. The computational time for the K-PCA itself is ∼2.59 seconds for "polblogs", ∼3.73 seconds for "cowrie" and ∼5.65 seconds for "wines", which is less than the computational time needed to train the map in all cases.

Discussion on the comparison

Results demonstrate a high efficiency of both the DR method and the sparse approach to decrease data dimensionality while producing accurate results in a reasonable computational time. For almost all quality criteria, both approaches introduced in the manuscript are in the range of or outperform the results obtained with the direct K-SOM at a very reduced computational time and a limited dimensionality of the resulting prototypes.

More specifically, QE is always better for K-PCA SOM with the smallest dimensionality. However, this is merely an effect of the dimensionality of Table 1: Performance results of K-PCA SOM and sparse K-SOM (average over 100 maps and standard deviation between parenthesis) for the "lesmis" dataset. Parameters for the methods are given between parenthesis after the method name (% of entropy preserved in the projection for K-PCA SOM and maximum mass, ν, and update parameter, κ, for random ascending updates in sparse K-SOM).

Methods the input data and is not a reliable criterion to compare results. ICI is a better way to measure the cluster homogeneity and shows that in almost all cases, K-PCA SOM and sparse K-SOM have comparable or even better ICI than the original approach (standard K-SOM) with a slight advantage for sparse K-SOM. Only, sparse K-SOM with the "wines" dataset exhibits a poor performance for this criterion. The quality of the organization of the map is measured with TE, which is often very comparable in both approaches to the direct K-SOM. However, it is again poor for the sparse K-SOM applied to the "wines" dataset and frequently tends to increase when ICI decreases. A good compromise between ICI and TE is reached for K-PCA SOM with an preserved entropy rate of 60% and for sparse K-SOM with ν equal to 90%/95% and κ equal to 50 for all datasets.

When comparing the different results with the a priori external information (through modularity or NMI), again the proposed approaches are comparable to or even better than the direct sparse K-SOM. However, in some cases ("lesmis", "astraptes", "cowrie"), it tends to highly deteriorate when the dimensionality is decreasing. This shows that a good tradeoff has to be found between interpretability (small dimensionality) and preservation of most of the information from the original dataset.

Finally, the stability of the results is often increased by K-PCA SOM (except for the "cowrie" dataset). This is an expected behavior: since the dimensionality of the input dataset is reduced, the prototypes lies in a reduced dimensionality space with less degrees of freedom. On the contrary, sparse K-SOM is not as appealing for this criterion because, even if the prototype representation is constrained, it is expressed in the original data space, which has a high dimension.

These good performance of our approaches come along with a high computational efficiency: K-PCA SOM and sparse K-SOM allows to highly reduce the computational times. Even if the K-PCA itself is a problem with a high complexity, it does not affect much the efficiency of K-PCA SOM on large datasets with a few thousands observations. The reason is that relational SOM is itself a time consuming algorithm, as explained in the introduction of this article. For a more challenging dataset (with more than ten thousands observations), a comparison and discussion is provided in Section 6.3.

In conclusion, both approaches introduced in this article are valid alternative to the direct K-SOM for large relational datasets. They both provide simplified prototype representation at a very reduced computational time. The dimensionality of the final prototypes is in all our examples 10 times or even 100 times smaller than in the direct sparse SOM version with no loss in accuracy of the resulting map. However, it should be noted that the easiness to interpret the prototypes is not exactly equivalent in K-PCA SOM and in sparse K-SOM: in the first approach, prototypes are expressed on the axes of the K-PCA which have to be interpreted in a previous step of the analysis, whereas, in sparse K-SOM, they are directly related to (a few number of) original observations. A detailed study is then presented in Sections 6.3.2 and 6.3.3.

Comparison with other approaches

In this section, the datasets "lesmis" and "astraptes" presented in Section 6.2.1 are used to compare the sparse approach to the simple numerical heuristic approximation (denoted N -num) and the simple geometric heuristic approximation (denoted N -geom) described in Section 4.2. Both strategies lead to a fixed number of observations for representing all the prototypes. Thus, for a fair comparison, this number was chosen considering the average final dimensionality of the prototypes obtained with the sparse method (last column in Tables 1 and2). The random ascending updates are performed using the same setting than in the original sparse version (i.e., random ascending updates κ ∈ {1, 50}).

Results obtained with N -num and N -geom over 100 maps are provided in Table 6 for "lesmis" and in Table 7 for "astraptes". In average, for both datasets, N -num gives better results than N -geom. As observed with the sparse method, both N -num and N -geom obtain best results with a final dimension equal to 8 for "lesmis". For all quality criteria except for the ICI, Table 6: Comparison between the different variants for the sparse updates (average over 100 maps and standard deviation between parenthesis) for the "lesmis" dataset. Parameters for the methods are given between parenthesis after the method name (maximum mass, ν, and update parameter, κ, for random ascending updates).

Methods QE (×100) ICI (×100) TE (%) Modularity (×100) Stability (%) Table 7: Comparison between the different variants for the sparse updates (average over 100 maps and standard deviation between parenthesis) for the "astraptes" dataset. Parameters for the methods are given between parenthesis after the method name (maximum mass, ν, and update parameter, κ, for random ascending updates). Methods QE (×100) ICI (×10 our variant of the sparse updates slightly improves the heuristic approximation approaches. This conclusion is supported by the results on "astraptes", except that N -num gives a better classification than our variant of the sparse updates.

N -num (4,
For the same level of sparsity, the results obtained with both N -num and N -geom slightly deteriorate the map quality compared to what can be observed with our version of the sparse updates.

Using K-PCA SOM and sparse K-SOM for mining job trajectories

This section presents the experiments performed on a more realistic dataset, with a larger sample size, obtained from the survey "Generation 98" [START_REF] Cottrell | How to use the Kohonen algorithm to simultaneously analyse individuals in a survey[END_REF][START_REF] Côme | Analysis of professional trajectories using disconnected self-organizing maps[END_REF]. The dataset contains information on career paths of 16,040 young people monitored during 94 months after having graduated in 1998. Nine categories are used to describe labor market statuses: permanent-labor contract, fixed-term contract, apprenticeship contract, public temporary-labor contract, on-call contract, unemployed, inactive, military service, education. Dissimilarities between career trajectories were computed using the optimal matching [START_REF] Needleman | A general method applicable to the search for similarities in the amino acid sequence of two proteins[END_REF][START_REF] Abbott | Optimal matching methods for historical sequences[END_REF] on the 12,560 unique career paths. This resulted in a non positive dissimilarity (6,651 eigenvalues out of 12 500 were found positive).

To assess the accuracy and the computational cost of both K-PCA SOM and sparse SOM, 100 maps were trained, using an R implementation of the methods, on the same 40-nodes computer without concurrent access. All maps were trained for a 10 × 10 grid, equipped with a piecewise linear neighborhood, with 60,000 iterations. The entropy ratio preserved by K-PCA SOM was varied in {20%, 40%, 60%, 80%}. For sparse K-SOM, the mass parameter ν was varied in {95%, 99%} and the update parameter κ was varied in {1, 50}. Only 10 maps were trained using the standard K-SOM due to its very high computational cost (more than ten days for each map).

Table 8 presents the results obtained in terms of QE, ICI, TE and CPU time (only the clustering time is reported, the K-PCA computational time is ∼8, 000 seconds. A detailed study to address this issue is made in Section 6.3.4). The last column provides the final dimension or number of coefficients of the prototypes.

As observed in Section 6.2.3, results demonstrate a high efficiency, in term of computational cost, of both K-PCA SOM and sparse K-SOM, while producing accurate results. Direct K-SOM takes more than ten days to train one map, whereas the slowest alternative strategy only requires sixteen minutes (i.e., the sparse K-SOM with ν = 99% and κ = 1). The results also confirm what was found with the toy datasets: K-PCA SOM provides a good trade-off between a good map quality and low dimensional prototypes and outperforms sparse K-SOM. The best results for K-PCA SOM are obtained with 20% entropy-rate preserved. However, this strategy selects only two dimensions, which increases the redundancy in the data and tends to produce clusters with few observations. Thus, the K-PCA SOM preserving 40% Table 8: Performance results of K-PCA SOM and sparse K-SOM (average over 100 maps and standard deviation between parenthesis) for the "trajectories" dataset. Parameters for the methods are given between parenthesis after the method name (% of entropy preserved in the projection for K-PCA SOM and maximum mass, ν, and update parameter, κ, for random ascending updates in sparse K-SOM).

Methods entropy should be preferred. The best results for the sparse K-SOM are obtained with a mass equal to 95%. Both K-PCA SOM and sparse K-SOM provide accurate results in a reasonable computational time. For sparse K-SOM, prototypes can be interpreted by inspecting the properties of the few observations used to represent them. For K-PCA SOM, the projection of the data on a subspace requires to interpret the axes of the K-PCA as an extra step in order to understand the meaning of the prototypes. A detailed study showing how the results of both approaches can be interpreted is performed in Sections 6.3.2 and 6.3.3.

Analysis of the influence of the prototype size

The job trajectory dataset was further used to assess the influence of the size of the prototypes on different characteristics. More precisely, we conducted an experiment to analyze the relation between the dimension of the projection (in K-PCA SOM) or the average number of coefficients per prototype (in sparse K-SOM) and some numerical characteristics of the algorithm: ICI, TE, CPU time and stability. A larger set of parameters, corresponding to varying dimensions of the prototypes, were tested by training 20 maps with each value of the parameters:

• for K-PCA SOM, the number of dimensions used for the projection • for sparse K-SOM, following the results given in Table 8, we set the random ascending update parameter κ to 1 and varied the maximum mass, ν in {0.9, 0.95, 0.975, 0.99, 0.995}. Figure 2 provides the evolution of different numerical characteristics of the maps versus the average number of coefficients per prototype.

For K-PCA SOM, up to approximately 50 (over a maximum of 12, 500), the dimension seems to have only a very limited effect on the quality of the map. All quality criteria stabilize after this value, with a slight tendency to improve and then to deteriorate again for very high dimensions. A strong computational time benefit can be observed when decreasing the dimensionality below 1, 000: this is a direct consequence of the quadratic complexity of K-SOM.

For sparse K-SOM, all characteristics, except for CPU time, tend to improve when the number of coefficients increases. However, TE seems to sta- .

bilize for approximately 50 coefficients per prototypes (ν = 0.99). However, since the computational time is not reduced from the same amount than in K-PCA SOM, testing the value of ν is not a good strategy. Moreover, ν = 0.995 gives prototypes with approximately 80 coefficients, which must not be increased to preserve interpretability.

Interpretability of K-PCA SOM

In this section, the map with the lowest ICI, among the 100 generated by the K-PCA SOM with 40% preserved entropy rate, is used to show how the results of K-PCA SOM can be interpreted, despite the K-PCA preprocessing. Its performances are equal to 5.74 (QE), 25.42 (ICI) and 7.66 (TE). The projection of the data on a subspace requires to interpret the K-PCA axes first. Figure 3 (left) presents the entropy supported by the first 15 axes and shows that the first two axes are enough to provide relevant information on the data. Figure 3 (right) displays the projections of the observations on the first two principal axes. The first axis represents 16.90% of the total entropy and opposes permanent-labor and fixed term contracts. This is supported by Figure 4 which shows the distribution of the 25 career paths with the smallest and the highest coordinates on the first two axes of the K-PCA. Stable job trajectories have the smallest coordinates on the first axis when fixed-term contract or unemployed have the highest. Figure 3 (right) also demonstrates that the second axis separates two kinds of precarious situations. Fixed-term contracts are opposed to highly precarious contracts such as unemployment, inactivity, on-call contract and public temporary-labor contract. The same observations can be made regarding the first 25 career paths with the smallest and highest coordinates on the second K-PCA axis, as shown in Figure 4.

The distribution of the job trajectories within each neuron of the map is represented in Figure 5. First note that the presented map is comparable in term of topology to the one described in [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF]. Different typologies can be highlighted: a fast access to permanent contracts (clear blue) on the bottom-left corner of the map, a transition through fixed-term contracts before obtaining stable ones (dark and then clear blue) on the map top-left corner, temporary jobs (dark blue) on the top-middle neurons, a long period of inactivity (yellow) or unemployment (red) on the map bottom-right corner.

The map organization is in accordance with the axis interpretation. Figure 6 (top) displays the average coordinates on the first and second axes in every cluster of the map. Results show a gradient of the observation coordinates on the first K-PCA axis between the bottom-left and the right side of the map. This gradient can also be seen on the map shown in Figure 5 and on the heatmaps presented on Figure 6 (bottom) which represent the cluster average of the career path modes 4 . This confirms that the first principal component (and corresponding diagonal on the map) separates permanent contracts from instable career paths. In Figure 6 (top), a gradient can also be observed for the second K-PCA axis between the top-left, where trajectories correspond to a fast access to permanent-labor contracts and the bottomleft corner of the map, where trajectories pertaining to precarious jobs are gathered.

Interpretability of sparse K-SOM

Similarly to the previous section, the present section provides a short discussion about one of the final results obtained from sparse K-SOM. The The resulting distribution of the job trajectories within the clusters of the map is provided in Figure 7. This distribution is fairly similar to the one obtained in Section 6.3.2: the left hand side of the map corresponds to a fast access to permanent contracts whereas the right hand side corresponds to different types of precarious situations. Two main differences can be highlighted: first, the class are more homogeneous in sparse K-SOM, especially at the border of the map. This is a direct effect of the dimension reduction in K-PCA SOM: since the dimension reduction increases redundancy in the dataset, some clusters (mostly located at the borders of the map) contain more observations and are thus less homogeneous. Second, the precarious situations (on the right hand side of the map) are organized a bit differently (with on-call contracts in the middle or the bottom of the map). However, both representations are realistic, with most of the clusters in the map being homogeneous.

A similar representation is provided in Figure 8 (left) but restricted to the observations which are involved in the prototype definition. The right part of this figure displays the number of such observations. Two conclusions can be derived from these graphics: the first one is that the observations involved in the prototype definition have a distribution very similar to the distribution of the entire set of observations included in the corresponding cluster. They are thus a selected subset of observations representative of their cluster. Moreover, as their number is very restricted compared to the total number of observations included in a cluster (approximately 14/15 observations as shown in the left part of Figure 8), they are a convenient way for the user to make sense of the prototypes and thus, of the corresponding cluster, since an exhaustive inspection of these observations becomes possible.

Nyström approximation

To evaluate the relevance of using a Nyström approximation of the K-PCA, the K-PCA SOM with 40% preserved entropy rate is used as a reference. Table 9 presents K-PCA SOM results using a Nyström approximation with different rates of observations sampled to perform the approximation. This rate was varied in {100%, 25%, 10%, 5%, 1%}, in which the 100%-results are reported from Table 8. The coefficients given to the K-PCA SOM are restricted to the first eight K-PCA axes everywhere, to avoid any bias related to data dimensionality. The computational time is reported in Table 9 and gives the time needed to perform the K-PCA only, excluding the training and clustering times. Results demonstrate a high efficiency in terms of computational time of the Nyström approximation while producing accurate results. In fact, none of the tested values lead to deteriorate the map quality in term of QE, ICI and TE, while the K-PCA is ∼ 1000 times faster when using 10% of observations. The best ICI is even obtained using only 1% of the observations. The clustering stability decreases with the number of observations used by the Nyström approximation, even if the stability is still high when using at least 10% of the observations.

The maps with the smallest ICI among the 100 maps generated from a Nyström approximation using 1% and 5% of the observations are displayed in Figure 9. Results show the ability of the Nyström approximation to preserve a realistic representation of the dataset while reducing the computational time.

The map obtained with 5% of the observations shows an organization similar to the one presented in Section 6.3.2, except for one fact: precarious situations, located on the right side of the map, are organized differently. Unemployment and public temporary-labor contracts are inverted between the top and the bottom of the map. With a Nyström approximation using 1% of the observations, trajectories mostly containing fixed-term contracts are located on the right hand side of the map, what is not observed on the map obtained with 5% of the observations and on maps presented in Section 6.3.2 and in Section 6.3.3. As expected, clusters obtained using 1% of the observations are less homogeneous than those obtained with 5%. The differences between these two maps might have different causes. Firstly, the instability of the SOM algorithm can explain the differences in terms of map organization: different runs of the algorithm give different results. This is particularly critical when the dataset to be analyzed is high dimensional as can be the "Generation 98" survey (even with a subsampling rate of 1%, the dimensionality of the problem is still larger than 100). This issue could be addressed by aggregating strategies, as described in [START_REF] Mariette | Aggregating self-organizing maps with topology preservation[END_REF]. Secondly, the differences between the two maps in Figure 9 might be explained by the high redundancy of trajectories with fixed-term contracts in the dataset: a very small subsampling might enforce the over-representation of these trajectories and affect the result. Such a problem could be addressed by using more efficient sampling techniques such as the ones described in [START_REF] Kumar | Sampling techniques for the Nyström method[END_REF].

The choice of the ratio m/n of observations to select in order to obtain accurate results highly depends on the quality of the kernel approximation provided by the Nyström technique. This quality is strongly influenced by the rank of the kernel, which can not easily be obtained when n is very large. Adaptative sampling technique for the Nyström algorithm, such as the one described in [START_REF] Drineas | Relative-error CUR matrix decompositions[END_REF][START_REF] Gittens | Revisiting the nystrom method for improved large-scale machine learning[END_REF] are based on an unequal probability sampling, which is performed iteratively and depends on the reconstruction error. [START_REF] Kumar | Sampling techniques for the Nyström method[END_REF] proposes an improved version in which the full kernel is not even needed to estimate the reconstruction error. Such methods could be relevant to assess Table 9: Performance results of the K-PCA SOM with K-PCA performed through a Nyström approximation (average over 100 maps and standard deviations between parenthesis) for the "trajectories" dataset. After the method name and between parenthesis, the percentage of observations used to perform the approximation is given. the evolution of the quality reconstruction in a growing sample and to stop the Nyström sampling when this quality is considered good enough.

Methods

Conclusion

The contributions of the present manuscript to the analysis of (dis)similarity data with topographic maps are twofolds: firstly, we have proposed a new version of the kernel and relational SOM algorithms, called sparse K-SOM, which ensures a sparse representation of the prototypes. Secondly, this approach has been compared to a preprocessing of the data by a dimension-reduction technique (K-PCA). We have also investigating the use of a Nyström approximation technique to ensure a better scalability of the method.

The experiments performed on several real datasets showed that both presented methods allow to strongly decrease the overall computational time on large datasets. The interpretability of the results is also improved since the prototypes have a much lower dimensionality. Moreover, the accuracy of the final map, in terms of cluster homogeneity, quality of the organization or adequacy to external a priori information is not deteriorated.

In average, K-PCA SOM gives better results in term of map and clustering quality than the sparse approach. However, prototypes returned by sparse K-SOM can be directly interpreted by inspecting the properties of the few observations used to represent them. In K-PCA SOM, the axes of the K-PCA have to be interpreted as an extra step to understand the prototypes meaning: this step is fairly standard in K-PCA. In conclusion, when selecting one or the other method, the user should also consider his/her need to easily interpret the result.

With the introduction of these methods, a further step is taken in allowing relational SOM to deal with massive data sets. Future works should investigate the multiple relational SOM, presented in [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF], with a view to integrate several sources of data of different types, while preserving the interpretability of the results.

1 .

 1 State-of-the art on SOM for (dis)similarity data

p and p k=1 β 0 uk = 1 7 :

 7 for t = 1, . . . , T do 8:

Figure 1 : 20 .

 120 Figure 1: K-PCA SOM. Average performances over 20 maps (ICI, TE, computational time in seconds and stability as measured by NMI) versus the dimensionality. Error bars correspond to the standard error standard deviation√ 20

Figure 2 :

 2 Figure 2: Sparse K-SOM. Average performances over 20 maps (ICI, TE, computational time in seconds and stability as measured by NMI) versus the average sparsity (average number of coefficients per prototypes) for different values of ν. Error bars correspond to the standard error standard deviation√ 20

Figure 3 :

 3 Figure 3: Entropy preserved by the 15 first axes on the left and projection of the observations on the first two principal components on the right. Colors represent the contract that appears the most often (mode) in the trajectory.

Figure 4 :

 4 Figure 4: Distribution of the 25 career paths with smallest and highest coordinates on the first two axes of the K-PCA.

Figure 5 :

 5 Figure 5: K-PCA SOM. For each neuron of the map, job trajectories distribution is represented using the observations classified in the corresponding unit. Colors represent the type of contract.

Figure 6 :

 6 Figure 6: K-PCA SOM. Representation of the SOM map with neurons filled using colors according to the average coordinate of the observations for the first (on the top-left) and the second (on the top-right) principal component. On the bottom, the map neurons are filled using colors according to the average of the career path modes.

Figure 7 :

 7 Figure 7: Sparse K-SOM. For each neuron of the map, job trajectories distribution is represented using the observations classified in the corresponding unit. Colors represent the type of contract.

Figure 8 :

 8 Figure 8: Sparse K-SOM. For each neuron of the map, job trajectories distribution for the observations with a positive coefficients for the corresponding neuron (left) and number of such observations (right).

Figure 9 :

 9 Figure 9: K-PCA SOM performed through a Nyström approximation using 1% (left) and 5% (right) of the observations: For each neuron of the map, job trajectories distribution is represented using the observations classified in the corresponding unit. Colors represent the type of contract.

Table 2 :

 2

		QE (×100)	ICI (×100)	TE (%) Modularity (×100) Stability (%) Dimension
	K-SOM	23.40 (0.44) 47.57 (2.53)	3.01 (2.52)	31.76 (2.80)	85.04 (2.13)	77
	K-PCA (80%)	18.30 (0.46) 47.31 (2.51)	3.30 (3.01)	32.02 (3.27)	86.20 (2.21)	27
	K-PCA (60%)	10.32 (0.55) 48.09 (2.41)	3.36 (3.65)	32.17 (2.63)	89.56 (2.17)	11
	K-PCA (40%)	2.55 (0.25)	52.62 (2.42) 2.64 (3.30)	25.96 (1.77)	91.85 (2.08)	4
	K-PCA (20%)	0.92 (0.18) 52.35 (2.21)	2.83 (3.32)	24.53 (1.92)	92.68 (2.02)	2
	sparse (90%, 1)	29.34 (1.20) 48.80 (2.88)	9.64 (7.00)	31.44 (4.09)	77.72 (2.94)	4
	sparse (90%, 50) 23.25 (5.61) 46.63 (2.70) 4.45 (3.62)	30.96 (3.03)	83.64 (2.34)	8
	sparse (99%, 1)	23.36 (4.06) 46.63 (2.71) 3.25 (2.85)	31.40 (3.03)	84.64 (2.18)	13
	sparse (99%, 50) 23.40 (4.49) 47.07 (2.57)	3.31 (2.72)	31.74 (3.41)	84.99 (2.35)	15
	Methods	QE (×100)	ICI (×10 4 )	TE (%)	NMI (%) Stability (%) Dimension
	K-SOM	1.02 (0.06)	3.37 (0.43)	1.39 (1.99) 82.19 (0.59) 94.57 (0.88)	459
	K-PCA (80%)	0.21 (0.02)	4.34 (0.70)	1.79 (2.20) 90.73 (2.86) 95.20 (1.05)	5
	K-PCA (60%)	0.10 (0.01)	6.40 (1.42)	2.71 (2.64) 80.99 (0.67) 94.73 (0.96)	3
	K-PCA (40%)	0.06 (0.01)	9.33 (1.52)	2.64 (2.27) 79.14 (0.59) 94.49 (0.95)	2
	K-PCA (20%)	0.01 (0.00) 23.15 (2.06) 24.55 (8.84) 73.02 (0.61) 94.70 (0.85)	1
	sparse (90%, 1)	2.08 (1.46)	3.59 (8.97)	2.94 (3.09) 86.34 (1.24) 91.87 (1.40)	5
	sparse (90%, 50) 0.99 (0.10)	3.08 (5.70)	3.72 (4.16) 82.60 (0.66) 93.62 (0.86)	8
	sparse (99%, 1)	0.98 (0.05)	3.30 (4.90)	1.51 (2.54) 82.22 (0.58) 94.66 (0.84)	25
	sparse (99%, 50) 20.99 (0.06) 3.38 (4.55)	1.39 (1.63) 82.16 (0.59) 94.73 (0.77)	25

Performance results of K-PCA SOM and sparse K-SOM (average over 100 maps and standard deviation between parenthesis) for the "astraptes" dataset. Parameters for the methods are given between parenthesis after the method name (% of entropy preserved in the projection for K-PCA SOM and maximum mass, ν, and update parameter, κ, for random ascending updates in sparse K-SOM).

Table 3 :

 3 Performance results of K-PCA SOM and sparse K-SOM (average over 100 maps and standard deviation between parenthesis) for the "polblogs" dataset. Parameters for the methods are given between parenthesis after the method name (% of entropy preserved in the projection for K-PCA SOM and maximum mass, ν, and update parameter, κ, for random ascending updates in sparse K-SOM).

	Methods	ICI (×100)	TE (%)	NMI (%) Stability (%) CPU time Dimension
	K-SOM	84.01 (0.85)	21.93 (1.51)	20.56 (0.30) 64.81 (0.73) 18269 (2378)	1599
	K-PCA (80%)	86.39 (0.84)	14.01 (1.42)	20.93 (0.27) 69.68 (0.79) 27.60 (4.45)	257
	K-PCA (60%)	86.61 (0.77)	14.87 (1.35)	20.96 (0.26) 70.62 (0.82) 18.88 (2.12)	121
	K-PCA (40%)	89.35 (2.28)	13.56 (2.12)	21.14 (0.32) 60.27 (1.73) 16.55 (2.45)	50
	K-PCA (20%)	91.46 (0.87) 12.06 (1.51) 21.30 (0.25) 77.19 (0.77) 13.35 (2.35)	13
	sparse (95%, 1)	85.20 (1.40)	34.03 (2.80)	20.64 (0.41) 55.14 (0.41) 60.13 (2.50)	8
	sparse (95%, 50) 84.00 (0.89) 32.61 (2.31)	20.66 (0.32) 60.39 (0.32) 67.89 (7.86)	12
	sparse (99%, 1)	84.18 (0.75)	23.95 (1.66)	20.46 (0.30) 63.98 (0.30) 88.19 (5.30)	34
	sparse (99%, 50) 84.06 (0.76)	24.15 (1.62)	20.51 (0.30) 64.30 (0.30) 96.73 (6.75)	34

Table 4 :

 4

	Methods	ICI (×100)	TE (%)	NMI (%) Stability (%) CPU time Dimension
	K-SOM	2.07 (0.24)	1.57 (1.29) 90.45 (0.86) 94.79 (1.04) 3771 (1534)	1414
	K-PCA (80%)	2.05 (0.19)	1.92 (1.39)	90.26 (0.73) 94.55 (0.89) 16.68 (1.45)	23
	K-PCA (60%)	1.94 (0.21)	2.33 (1.33)	89.11 (0.81) 94.04 (0.80) 15.19 (1.31)	10
	K-PCA (40%)	3.15 (0.30)	3.85 (1.96)	85.72 (0.80) 90.49 (0.93) 15.35 (2.41)	4
	K-PCA (20%)	6.26 (0.28)	5.11 (2.18)	74.47 (0.48) 87.99 (1.17) 15.04 (2.27)	2
	sparse (95%, 1)	2.09 (0.26)	2.04 (1.31)	89.64 (0.97) 91.78 (1.18) 70.03 (1.23)	8
	sparse (95%, 50) 1.84 (0.20) 2.38 (1.82)	90.37 (0.70) 93.78 (0.92) 62.66 (3.27)	13
	sparse (99%, 1)	20.36 (0.22) 1.55 (1.36)	90.34 (0.79) 94.45 (0.91) 90.73 (4.16)	37
	sparse (99%, 50) 20.00 (0.22) 1.54 (1.24) 90.37 (0.79) 94.50 (0.93) 81.14 (2.89)	37

Performance results of K-PCA SOM and sparse K-SOM (average over 100 maps and standard deviation between parenthesis) for the "cowrie" dataset. Parameters for the methods are given between parenthesis after the method name (% of entropy preserved in the projection for K-PCA SOM and maximum mass, ν, and update parameter, κ, for random ascending updates in sparse K-SOM).

Table 5 :

 5 Performance results of K-PCA SOM and sparse K-SOM (average over 100 maps and standard deviation between parenthesis) for the "wines" dataset. Parameters for the methods are given between parenthesis after the method name (% of entropy preserved in the projection for K-PCA SOM and maximum mass, ν, and update parameter, κ, for random ascending updates in sparse K-SOM).

	Methods	ICI (×100)	TE (%)	NMI (%) Stability (%) CPU time Dimension
	K-SOM	22.10 (0.50) 10.37 (1.34) 11.86 (0.31) 74.23 (0.81) 13480 (10575)	1222
	K-PCA (80%)	21.94 (0.51) 10.15 (1.31) 11.78 (0.30) 74.96 (0.77) 20.00 (3.29)	28
	K-PCA (60%)	22.43 (0.61)	8.80 (1.30)	11.72 (0.34) 75.37 (0.88) 17.95 (1.74)	9
	K-PCA (40%)	25.44 (0.74)	7.08 (1.33) 11.86 (0.26) 77.65 (0.88) 17.32 (2.33)	4
	K-PCA (20%)	34.99 (0.36) 0.13 (0.72) 11.01 (0.23) 85.75 (1.68) 15.99 (1.88)	2
	sparse (95%, 1)	47.62 (1.46) 14.29 (1.35) 11.75 (0.36) 66.21 (0.79) 81.48 (5.02)	8
	sparse (95%, 50) 45.90 (1.08) 11.95 (1.39) 11.76 (0.38) 70.72 (0.77) 97.60 (14.28)	13
	sparse (99%, 1)	44.14 (0.95) 11.26 (1.36) 11.80 (0.30) 74.00 (0.73) 141.17 (12.28)	40
	sparse (99%, 50) 44.17 (0.93) 11.58 (1.24) 11.83 (0.32) 74.13 (0.77) 146.93 (14.98)	39

The graph as well as the shortest path-lengths are included in the R package SOMbrero[START_REF] Boelaert | SOMbrero: an r package for numeric and non-numeric self-organizing maps[END_REF].

The graph is available at http://www-personal.umich.edu/ ~mejn/netdata/ polblogs.zip.

The dataset is available at https://archive.ics.uci.edu/ml/datasets/Wine+ Quality.

To compute mode averages, job market contracts have been converted to numerical labels from 1, for the permanent-labor contract, to 9 for education.

Acknowledgements

We thank the two anonymous reviewers and the editor for valuable comments and suggestions which helped to improve the quality of the paper. Part of this work was carried out using the resources from the INRA GENOTOUL bioinformatic platform http://bioinfo.genotoul.fr. We are grateful to the platform staff and especially to Marie-Stéphane Trotard and Didier Laborie.

Appendix A. Formula for the average intra-cluster intertia

The intra cluster inertia is the average over u ∈ {1, . . . , U } of the quantities

Expansion for kernel data

For kernel data, this quantity equals

Expansion for dissimilarity data

Using the result of Equation ( 2) in [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF], we get that

in which ν u = 1 Cu 1 Cu where the entries of 1 Cu are equal to 1 for indexes i such that x i ∈ C u and to 0 otherwise. Thus

( C u ) 2 j,j : x j ,x j ∈Cu δ(x j , x j ) = 1 2 (C u ) 2 i,i : x i ,x i ∈Cu δ(x i , x i ).