Marie Cottrell
email: marie.cottrell@univ-paris1.fr

Madalina Olteanu
email: madalina.olteanu@univ-paris1.fr

Fabrice Rossi
email: fabrice.rossi@univ-paris1.fr

Nathalie Vialaneix
email: nathalie.villa@toulouse.inra.fr

Nathalie Villa-Vialaneix

Theoretical and Applied Aspects of the Self-Organizing Maps

Keywords: SOM, batch SOM, Relational SOM, Stability of SOM

1 Brief history of the SOM Since its introduction by T. Kohonen in his seminal 1982 articles (Kohonen [1982a,b]), the self-organizing map (SOM) algorithm has encountered a very large success. This is due to its very simple definition, to the easiness of its practical development, to its clustering properties as well as its visualization ability. SOM appears to be a generalization of basic clustering algorithms and at the same time, provides nice visualization of multidimensional data.

The basic version of SOM is an on-line stochastic process which has been inspired by neuro-biological learning paradigms. Such paradigms had previously been used to model some sensory or cognitive processes where the learning is directed by the experience and the external inputs without supervision. For example, [START_REF] Ritter | On the stationary state of kohonen's self-organizing sensory mapping[END_REF] illustrate the somatosensory mapping property of SOM. However, quickly in the eighties, SOM was not restricted to neuro-biology modeling and has been used in a huge number of applications (see e.g. [START_REF] Kaski | Bibliography of self-organizing map (SOM) papers: 1981-1997[END_REF], [START_REF] Oja | Bibliography of self-organizing map (SOM) papers: 1998-2001 addendum[END_REF] for surveys), in very diverse fields as economy, sociology, text mining, process monitoring, etc.

Since then, several extensions of the algorithms have been proposed. For instance, for users who are not familiar with stochastic processes or for industrial applications, the variability of the equilibrium state was seen as a drawback because the learnt map is not always the same from one run to another. To address this issue, T. Kohonen introduces the batch SOM, [START_REF] Kohonen | Self-Organizing Maps[END_REF][START_REF] Kohonen | Comparison of som point densities based on different criteria[END_REF]) which is deterministic and thus leads 3 to reproducible results (for a given initialization). Also, the initial SOM (on-line or batch versions) was designed for real-valued multidimensional data, and it has been necessary to adapt its definition in order to deal with complex non vectorial data such as categorical data, abstract data, documents, similarity or dissimilarity indexes, as introduced in [START_REF] Kohonen | Median strings[END_REF], Kaski et al. [1998a], [START_REF] Kohohen | Self-organizing maps of symbol strings[END_REF]. One can find in [START_REF] Kohonen | Self-organization and associative memory[END_REF][START_REF] Kohonen | Self-Organizing Maps[END_REF][START_REF] Kohonen | Self-Organizing Maps, 3rd Edition[END_REF][START_REF] Kohonen | Essentials of self-organizing map[END_REF][START_REF] Kohonen | MATLAB Implementations and Applications of the Self-Organizing Map[END_REF] extensive lists of references related to SOM. At this moment more than 10 000 papers have been published on SOM or using SOM.

In this paper, we review a large selection of the numerous variants of the SOM. One of the main focus of this survey is the question of convergence of the SOM algorithms, viewed as stochastic processes. This departs significantly from the classical learning theory setting. In this setting, exemplified by the pioneering results of [START_REF] Pollard | Strong consistency of k-means clustering[END_REF], one generally assumes given an optimization problem whose solution is interesting: for instance, an optimal solution of the quantization problem associated to the k-means quality criterion. The optimization problem is studied with two points of view: the true problem which involves a mathematical expectation with respect to the (unknown) data distribution and its empirical counterpart where the expectation is approximated by an average on a finite sample. Then the question of convergence (or consistency) is whether the solution obtained on the finite sample converges to the true solution that would be obtained by solving the true problem.

We focus on a quite different problem. A specific stochastic algorithm such as the SOM one defines a series of intermediate configurations (or solutions). Does the series converge to something interesting? More precisely, as the algorithm maps the inputs (the data) to an output (the prototypes and their array), one can take this output as the result of the learning process and may ask the following questions, among others:

-How to be sure that the learning is over? -Do the prototypes extract a pertinent information from the data set? -Are the results stables? -Are the prototypes well organized?

In fact, many of these questions are without a complete answer, but in the following, we review parts of the questions for which theoretical results are known and summarize the main remaining difficulties. Section 2 is devoted to the definition of SOM for numerical data and to the presentation of the general methods used for studying the algorithm. Some theoretical results are described in the next sections: Section 3 explains the one dimensional case while in Section 4, the results available for the multidimensional case are presented. The batch SOM is studied in Section 5. In Section 6, we present the variants proposed by Heskes to get an energy function associated to the algorithm. Section 7 is dedicated to non numerical data. Finally, in Section 8, we focus on the use of the stochasticity of SOM to improve the interpretation of the maps. The article ends with a very short and provisional conclusion.

SOM for numerical data

Originally, (in Kohonen [1982a] and Kohonen [1982b]), the SOM algorithm was defined for vector numerical data which belong to a subset X of an Euclidean space (typically R p ). Many results in this paper additionally require that the subset is bounded and convex. There are two different settings from the theoretical point of view:

continuous setting: the input space X can be modeled by a probability distribution defined by a density function f ; -discrete setting: the data space X comprises N data points x 1 , . . . , x N in R p (In this paper, by discrete setting, we mean a finite subset of the input space).

The theoretical properties are not exactly the same in both cases, so we shall later have to separate these two settings.

Classical on-line SOM, continuous or discrete setting

In this section, let us consider that X ⊂ R p (continuous or discrete setting).

First we specify a regular lattice of K units (generally in a one-or twodimensional array). Then on the set K = {1, . . . , K}, a neighborhood structure is induced by a neighborhood function h defined on K × K. This function can be time dependent and, in this case, it will be denoted by h(t). Usually, h is symmetrical and depends only on the distance between units k and l on the lattice (denoted by dist(k, l) in the following)). It is common to set h kk = 1 and to have k kl decrease with increasing distance between k and l. A very common choice is the step function, with value 1 if the distance between k and l is less than a specific radius (this radius can decrease with time), and 0 otherwise. Another very classical choice is

h kl (t) = exp - dist 2 (k, l) 2σ 2 (t) ,
where σ 2 (t) can decrease over time to reduce the intensity and the scope of the neighborhood relations.

A prototype m k ∈ R p is attached to each unit k of the lattice. Prototypes are also called models, weight vectors, code-vectors, codebook vectors, centroids, etc. The goal of the SOM algorithm is to update these prototypes according to the presentation of the inputs in such a way that they represent the input space as accurately as possible (in a quantization point of view) while preserving the topology of the data by matching the regular lattice with the data structure. For each prototype m k , the set of inputs closer to m k than to any other one defines a cluster (also called a Voronoï cell) in the input space, denoted by C k , and the neighborhood structure on the lattice induces a neighborhood structure on the clusters. In other words, after running the SOM process, close inputs should belong to the same cluster (as in any clustering algorithm) or to neighbor clusters.

From any initial values of the prototypes, (m 1 (0), . . . , m K (0)), the SOM algorithm iterates the following steps:

1. At time t, if m(t) = (m 1 (t), . . . , m K (t))
denotes the current state of the prototypes, a data point x is drawn according to the density f in X (continuous setting) or at random in the finite set X (discrete setting). 2. Then c t (x) ∈ {1, . . . , K} is determined as the index of the best matching unit, that is

c t (x) = arg min k∈{1,...,K} x -m k (t) 2 , (1)
3. Finally, all prototypes are updated via

m k (t + 1) = m k (t) + ǫ(t)h kc t (x) (t)(x -m k (t)), (2)
where ǫ(t) is a learning rate (positive, less than 1, constant or decreasing).

Although this algorithm is very easy to define and to use, its main theoretical properties remain without complete proofs. Only some partial results are available, despite a large amount of works and empirical evidences. More precisely, (m k (t)) k=1,...,K are K stochastic processes in R p and when the number t of iterations of the algorithm grows, m k (t) could have different behaviors: oscillation, explosion to infinity, convergence in distribution to an equilibrium process, convergence in distribution or almost sure to a finite set of points in R p , etc. This is the type of convergence that we will discuss in the sequel. In particular, the following questions will be addressed:

-Is the algorithm convergent in distribution or almost surely, when t tends to +∞? -What happens when ǫ is constant? when it decreases? -If a limit state exists, is it stable? -How to characterize the organization? One can find in [START_REF] Cottrell | Theoretical aspects of the SOM algorithm[END_REF] and [START_REF] Fort | SOM's mathematics[END_REF] a summary of the main rigorous results with most references as well as the open problems without solutions until now.

Mathematical tools related to the convergence of stochastic processes

The main methods that have been used to analyze the SOM convergence are summarized below.

-The Markov Chain theory for constant learning rate and neighboring function, which is useful to study the convergence and the limit states. If the algorithm converges in distribution, this limit distribution has to be an invariant measure for the Markov Chain. If it is possible to prove some strong organization, it has to be associated to an absorbing class;

-The Ordinary Differential Equation method (ODE), which is a classical method to study the stochastic processes.

If we write down the equation (2) for each k ∈ K in a vector form, we get

m(t + 1) = m(t) -ǫ(t)Φ(x, m(t)), (3)
where Φ is a stochastic term. To study the behavior of such stochastic processes, it is often useful to study the solutions of the associated deterministic ordinary differential equation that describes the mean behavior of the process. This ODE is

dm dt = -φ(m), (4)
where φ(m) is the expectation of Φ(., m) with respect to the probability distribution of the inputs x (continuous setting) or the arithmetic mean (discrete setting).

Here the k th -component of φ is

φ k (m) = K j=1 h kj Cj (x -m k )f (x)dx, (5)
for the continuous setting or

φ k (m) = 1 N K j=1 h kj xi∈Cj (x i -m k ), (6)
that can be also written

φ k (m) = 1 N N i=1 h kc(xi) (x i -m k ), (7)
for the discrete setting.

The possible limit states of the stochastic process in Equation ( 2) would have to be solutions of the equation

φ(m) = 0.
Then if the zeros of this function were the minima of a function (most often called energy function), it would be useful to apply the gradient descent methods.

-The Robbins-Monro algorithm theory which is used when the learning rate decreases under conditions

t ε(t) = +∞ and t ε(t) 2 < +∞. (8)
Unfortunately some remarks explain why the original SOM algorithm is difficult to study. Firstly, for dimension p > 1, a problem arises: it is not possible to define any absorbing class which could be an organized state. Secondly, although the process m(t) can be written down as a classical stochastic process of Equation (3), one knows since the papers Erwin et al. [1992a,b], that it does not correspond to an energy function, that is it is not a gradient descent algorithm in the continuous setting. Finally, it must be emphasized that no demonstration takes into account the variation of the neighborhood function. All the existing results are valid for a fixed size and intensity of the function h.

The one-dimensional case

A very particular setting is the one-dimensional case: the inputs belong to R and the lattice is a one-dimensional array (a string). Even though this case is of a poor practical utility, it is interesting because the theoretical analysis can be fully conducted.

The simplest one-dimensional case

The simplest case was fully studied in the article [START_REF] Cottrell | Étude d'un processus d'auto-organisation[END_REF]. The inputs are supposed to be uniformly distributed in [0, 1], the lattice is a onedimensional array {1, 2, . . . , K}, the learning rate ǫ is a constant smaller than 1 2 , the neighborhood function is a constant step function h kl = 0 if |k -l| > 1 and 1 otherwise. In that case the process m(t) is a homogeneous Markov Chain with continuous state space. The organization we look for is simply the ordering (ascending or descending) and so is easy to characterize. Let us describe the main steps of the proof.

1. There exists a decreasing functional: the number of badly ordered triplets.

But this is not sufficient to prove the convergence, it has to be strictly decreasing with a strictly positive probability. 2. The set of ordered dispositions is an absorbing class, composed of two classes which do not communicate: the increasing sequences class and the decreasing sequences class. 3. One shows that ordering (topology preservation in this special case) takes place after a finite time with a probability which is greater than a positive bound, and that the hitting time of the absorbing class is almost surely finite.

0 1 j -1 j j + 1 j -1 j j + 1 j -1 j j + 1 j -1 j j + 1 Fig. 1.
Four examples of triplets of prototypes (mj-1, mj, mj+1). For each j, j -1 and j + 1 are its neighbors. The y-axis coordinates are the values of the prototypes that take values in [0, 1]. The first two triplets on the left are badly ordered. In the case under study, SOM will order them with a strictly positive probability. The last two triplets (on the right) are well ordered and SOM will never disorder them.

4. Then one shows that the Markov Chain has the Doeblin property: there exists an integer T , and a constant c > 0, such that, given that the process starts from any ordered state, and for all set E in [0, 1] n , with positive measure, the probability to enter in E with less than T steps is smaller than c vol(E). 5. This implies that the chain converges in distribution to a monotonous stationary distribution which depends on ǫ (which is a constant in that part). 6. If ǫ(t) tends towards 0 and satisfies the Robbins-Monro conditions (8), once the state is ordered, the Markov Chain almost surely converges towards a constant (monotonous) solution of an explicit linear system.

So in this very simple case, we could prove the convergence to a unique ordered solution such that

m 1 (+∞) < m 2 (+∞) < . . . < m K (+∞), or m 1 (+∞) > m 2 (+∞) > . . . > m K (+∞).
Unfortunately, it is not possible to find absorbing classes when the dimension is larger than 1. For example, in dimension 2, with 8 neighbors, if the x-and y-coordinates are ordered, it is possible (with positive probability) to disorder the prototypes as illustrated in Figure 2.

What we know about the general one-dimensional case

We summarize in this section the essential results that apply to the general one dimension case (with constant neighborhood function and in the continuous setting). References and precise statements can be found in [START_REF] Fort | SOM's mathematics[END_REF] and [START_REF] Cottrell | Theoretical aspects of the SOM algorithm[END_REF]. Compared to the previous section, hypothesis on the data distribution and/or the neighborhood function are relaxed. . This figure represents 2-dimensional prototypes (x-and y-axes are not shown but are the standard horizontal and vertical axes) which are linked as their corresponding units on the SOM grid. At this step of the algorithm, the x-and y-coordinates of the prototypes are well ordered. But contrarily to the one-dimensional case, this disposition can be disordered with a positive probability; in an 8-neighbors case, A is C's neighbor, but B is not a neighbor of C. If C is repeatedly the best matching unit, B is never updated, while A becomes closer and closer to C. Finally, the y coordinate of A becomes smaller than that of B and the disposition is disordered.

-The process m(t) is almost surely convergent to a unique stable equilibrium point in a very general case: ǫ(t) is supposed to satisfy the conditions ( 8), there are hypotheses on the density f and on the neighborhood function h. Even though these hypotheses are not very restrictive, some important distributions, such as the χ 2 or the power distribution, do not fulfill them. -For a constant ǫ, the ordering time is almost surely finite (and has a finite exponential moment). -With the same hypotheses as before to ensure the existence and uniqueness of a stable equilibrium x * , from any ordered state, for each constant ǫ, there exists an invariant probability measure π ǫ . When ǫ tends to 0, this measure concentrates on the Dirac measure on x * . -With the same hypotheses as before to ensure the existence and uniqueness of a stable equilibrium x * , from any ordered state, the algorithm converges to this equilibrium provided that ǫ(t) satisfies the conditions (8).

As the hypotheses are sufficiently general to be satisfied in most cases, one can say that the one-dimensional case is more or less well-known. However nothing is proved neither about the choice of a decreasing function for ǫ(t) to ensure simultaneously ordering and convergence, nor for the case of decreasing neighborhood function.

Multidimensional case

When the data are p-dimensional, one has to distinguish two cases, the continuous setting and the discrete one.

Continuous setting

In the p-dimensional case, we have only partial results proved by Sadeghi in [START_REF] Sadeghi | Convergence in distribution of the multi-dimensional kohonen algorithm[END_REF]). In this paper, the neighborhood function is supposed to have a finite range, the learning rate ǫ is a constant, the probability density function is positive on an interval (this excludes the discrete case). Then the algorithm weakly converges to a unique probability distribution which depends on ǫ.

Nothing is known about the possible topology preservation properties of this stationary distribution. This is a consequence of the difficulty of defining an absorbing organized state in a multidimensional setting. For example, two results of Flanagan and Fort-Pagès illustrate the complexity of the problem. These two apparently contradictory results hold. For p = 2, let us consider the set F ++ of simultaneously ordered coordinates (respectively x and y coordinates). We then have:

for a constant ǫ and very general hypotheses on the density f , the hitting time of F ++ is finite with a positive probability [START_REF] Flanagan | Self-organisation in kohonen's som[END_REF]), -but in the 8-neighbor setting, the exit time is also finite with positive probability [START_REF] Fort | About the kohonen algorithm: strong or weak self-organisation[END_REF]).

Discrete setting

In this setting, the stochastic process m(t) of Equations ( 2) and ( 4) derives from a potential function, which means that it is a gradient descent process associated to the energy. When the neighborhood function does not depend on time, [START_REF] Ritter | Neural computation and self-Organizing Maps, an Introduction[END_REF] have proven that the stochastic process m(t) of Equations ( 2) and

(3) derives from a potential, that is it can be written

m k (t + 1) = m k (t) + ǫ(t)h kc(x) (t)(x -m k (t)), = m k (t) -ǫ(t)Φ k (x, m(t)), = m k (t) -ǫ(t) ∂ ∂m k E(x, m(t)),
where E(x, m) is a sample function of E(m) with

E(m) = 1 2N K k=1 K j=1 h kj xi∈Cj m k -x i 2 , (9)
or in a shorter expression

E(m) = 1 2N N i=1 K k=1 h kc(xi) m k -x i 2 . ( 10
)
In other words the stochastic process m(t) is a stochastic gradient descent process associated to function E(m). Three interesting remarks can be made:

1. The energy function is a generalization of the distortion function (or intra classes variance function) associated to the Simple Competitive Learning process (SCL, also known as the Vector Quantization Process/Algorithm), which is the stochastic version of the deterministic Forgy algorithm. The SCL process is nothing else than the SOM process where the neighborhood function is degenerated, ie when h kl = 1 only for k = l and h kl = 0 elsewhere. In that case, E reduces to

E(m) = 1 2N N i=1 m c(xi) -x i 2 .
For that reason, E is called extended intra-classes variance. 2. The above result does not ensure the convergence of the process: in fact the gradient of the energy function is not continuous and the general hypotheses used to prove the convergence of the stochastic gradient descent processes are not valid. This comes from the fact that there are discontinuities when crossing the boundaries of the clusters associated to the prototypes, because the neighbors involved in the computation change from a side to another. However this energy gives an interesting insight on the process behavior. 3. In the 0-neighbor setting, the Vector Quantization algorithm converges, since there is no problem with the neighbors and the gradient is continuous. However there are a lot of local minima and the algorithm converges to one of these minima.

5 Deterministic batch SOM

As the possible limit states of the stochastic process (2) would have to be solutions of the ODE equation φ(m) = 0, it is natural to search how to directly get these solutions. The definition of the batch SOM algorithm can be found in [START_REF] Kohonen | Self-Organizing Maps[END_REF][START_REF] Kohonen | Comparison of som point densities based on different criteria[END_REF].

From Equation (5), in the continuous setting, the equilibrium m * must satisfy

∀k ∈ K, K j=1 h kj Cj (x -m * k )f (x)dx.
Hence, for the continuous setting, the solution complies with

m * k = K j=1 h kj Cj xf (x)dx K j=1 h kj Cj f (x)dx .
In the discrete setting, the analogous is

m * k = K j=1 h kj xi∈Cj x i K j=1 h kj |C j | = N i=1 h kc(xi) x i N i=1
h kc(xi) Thus, the limit prototypes m * k have to be the weighted means of all the inputs which belong to the cluster C k or to its neighboring clusters. The weights are given by the neighborhood function h.

Using this remark, it is possible to derive the definition of the batch algorithm.

m k (t + 1) = K j=1 h kj (t) C j(m k (t)) xf (x)dx K j=1 h kj (t) C j(m k (t)) f (x)dx . ( 11
)
for the continuous setting, and

m k (t + 1) = N i=1 h kc t (xi) (t)x i N i=1 h kc t (xi) (t) (12)
for the discrete case. This algorithm is deterministic, and one of its advantages is that the limit states of the prototypes depend only on the initial choices. When the neighborhood is reduced to the unit itself, this batch algorithm for the SOM is nothing else than the classical Forgy algorithm [START_REF] Forgy | Cluster analysis of multivariate data: efficiency versus interpretability of classifications[END_REF]) for clustering. Its theoretical basis is solid and a study of the convergence can be found in [START_REF] Cheng | Convergence and ordering of kohonen's batch map[END_REF]. One can prove [START_REF] Fort | Stochastic on-line algorithm versus batch algorithm for quantization and self organizing maps[END_REF][START_REF] Fort | Advantages and drawbacks of the batch kohonen algorithm[END_REF]) that it is exactly a quasi-Newtonian algorithm associated to the extended distortion (energy) E (see Equation ( 10)), when the probability to observe a x in the sample which is exactly positioned on the median hyperplanes (e.g., the boundaries of C k ) is equal to zero. This assumption is always true in the continuous setting but it is not relevant in the discrete setting since there is no guarantee that data points never belong to the boundaries which vary along the iterations.

The batch SOM algorithm is the extension of the Forgy algorithm with the introduction of the neighborhood between clusters, in the same way as the online SOM algorithm is for the Vector Quantization algorithm. It is not exactly a gradient descent algorithm, but it converges to a minimum of the energy E. Obviously there are many local minima. In conclusion, the relations between these clustering algorithms are summarized in Table 1.

Other algorithms related to SOM

As explained before, the on-line SOM is not a gradient algorithm in the continuous setting (Erwin et al. [1992a,b]). In the discrete setting, there exists an energy function, which is an extended intra-classes variance as in Equation ( 10), but this function is not continuously differentiable. To overcome these problems, [START_REF] Heskes | Energy functions for self-organizing maps[END_REF] proposes to slightly modify the on-line version of the SOM algorithm so it can be seen as a stochastic gradient descent on the same energy function. To do so, he introduces a new hard assignment of the winning unit and a soft version of this assignment.

Hard assignment in the Heskes's rule

In order to obtain an energy function for the on-line SOM algorithm, [START_REF] Heskes | Energy functions for self-organizing maps[END_REF] modifies the rule for computing the best matching unit (BMU). In his setting, Equation (1) becomes

c t (x) = arg min k∈{1,...,K} K j=1 h kj (t) x -m k (t) 2 (13)
The energy function considered here is

E(m) = 1 2 K j=1 K k=1 h kj (t) x∈Cj(m) x -m k (t) 2 f (x)dx, (14)
where C j (m) is the cluster (Voronoï cell) associated to the j-th prototype. The regularity properties of the energy function and of its gradient are summarized in Table 2, as discussed in [START_REF] Heskes | Energy functions for self-organizing maps[END_REF].

Discrete setting Continuous setting

Kohonen rule Energy : discontinuous (but finite on V ) Energy : continuous for computing BMU Gradient : discontinuous (infinite on V ) Gradient : discontinuous

Heskes rule

Energy : continuous Energy : continuous for computing BMU Gradient : discontinuous (finite on V ) Gradient : continuous Table 2. Smoothness of the energy function

Soft Topographic Mapping (STM)

The original SOM algorithm is based on a hard winner assignment. Generalizations based on soft assignments were derived in [START_REF] Heskes | Energy functions for self-organizing maps[END_REF] and [START_REF] Graepel | Self-organizing maps: generalizations and new optimization techniques[END_REF]. First, let us remark that the energy function in the discrete case can also be written as

E(m, c) = 1 2 K k=1 N i=1 c ik K j=1 h kj (t) m j (t) -x i 2
where c ik is equal to 1 if x i belongs to cluster k and zero otherwise. This crisp assignment may be smoothed by considering c ik ≥ 0 such that K k=1 c ik = 1. The soft assignments may be viewed as the probabilities of input x i to belong to class k.

Since the optimization of the energy function with gradient descent-like algorithms would get stuck into local minima, the problem is transformed into a deterministic annealing scheme. The energy function is smoothed by adding an entropy term and transforming it into a "free energy" cost function, parameterized by a parameter β :

F (m, c, β) = E(m, c) - 1 β S(c) ,
where S(c) is the entropy term associated to the full energy. For low values of β, only one global minimum remains and may be easily determined by gradient descent or EM schemes. For β → +∞, the free energy has exactly the same expression as the original energy function. When using deterministic annealing, one begins by computing the minimum of the free energy at low values of β and then attempts to compute the minimum for higher values of β (β may be chosen to grow exponentially), until the global minimum of the free energy for β → +∞ is equal to the global minimum of the original energy function.

For a fixed value of β, the minimization of the free energy leads to iterating over two steps given by Equations ( 15) and ( 16), in batch version, and very similar to the original SOM (the neighborhood function h is not varied during the optimization process) :

P(x i ∈ C k ) = exp(-βe ik ) K j=1 exp(-βe ij ) , (15)
where

e ik = 1 2 K j=1 h jk (t) x i -m j (t) 2 and m k (t) = N i=1 x i K j=1 h jk (t)P(x i ∈ C j ) N i=1 K j=1 h jk (t)P(x i ∈ C j ) (16)
The updated prototypes are written as weighted averages over the data vectors. For β → +∞, the classical batch SOM is retrieved.

Probabilistic views on the SOM

Several attempts have been made in order to recast the SOM algorithm (and its variants) into a probabilistic framework, namely the general idea of mixture models (see e.g. [START_REF] Mclachlan | Finite mixture models[END_REF]). The central idea of those approaches is to constrain a mixture of Gaussian distributions in a way that mimic the SOM grid. Due to the heuristic nature of the SOM, the resulting models depart quite significantly from the SOM algorithms and/or from standard mixture models. We describe below three important variants. Other variants are listed in e.g. [START_REF] Verbeek | Self-organizing mixture models[END_REF].

SOM and regularized EM One of the first attempt in this direction can be found in [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF]. Based on his work on energy functions for the SOM, Heskes shows in this paper that the batch SOM can be seen as a form of regularized Expectation Maximization (EM) algorithm1 .

As mentioned above, the starting point of the this analysis consists in introducing an isotropic Gaussian mixture with K components. The multivariate Gaussian distributions share a single precision parameter β, with the covariance matrix 1 β I, and are centered on the prototypes. However, up to some constant terms, the opposite of the log likelihood of such a mixture corresponds to the k-means quantization error. And therefore, maximizing the likelihood does not provide any topology preservation. Thus Heskes introduces a regularization term which penalizes prototypes that do not respect the prior structure (the term does not depend directly on the data points), see [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] for details. Then Heskes shows that applying the EM principle to the obtained regularized (log)likelihood leads to an algorithm that resembles the batch SOM one.

This interpretation has very interesting consequences, explored in the paper. It is easy for instance to leverage the probabilistic framework to handle missing values in a principled (non heuristic) way. It is also easy to use other mixtures e.g. for non numerical data (such as count data). However, the regularization itself is rather ad hoc (it cannot be easily interpreted as a prior distribution on the parameters, for instance). In addition, the final algorithm is significantly different from the batch SOM. Indeed, as in the case of the STM, crisp assignments are replaced by probabilistic ones (the crispness of the assignments is controlled by the precision parameter β). In addition, as in STM, the neighborhood function is fixed (as it is the core of the regularization term). To our knowledge, the practical consequences of those differences have not been studied in detail on real world data. While one can argue that β can be increased progressively and at the same time, one can modify the neighborhood function during the EM algorithm, this might also have consequences that remain untested.

SOM and variational EM Another take at this probabilistic interpretation can be found in [START_REF] Verbeek | Self-organizing mixture models[END_REF]. As in [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] the first step consists in assuming a standard mixture model (e.g. a K components Gaussian isotropic mixture for multivariate data). Then the paper leverages the variational principle (see e.g. [START_REF] Michael I Jordan | An introduction to variational methods for graphical models[END_REF]).

In summary, the variational principle is based on introducing an arbitrary distribution q on the latent (hidden) variables Z of the problem under study. In a standard mixture model, the hidden variables are the assignment ones, which map each data point to a component of the mixture (a cluster in the standard clustering language). One can show that the integrated log likelihood of a mixture model with Θ as parameters, log p(X|Θ), is equal to the sum of three components: the complete likelihood (knowing both the data points X and the hidden variables Z) integrated over the hidden variables with respect to q, E q log p(X, Z|Θ), the entropy of q, H(q), and the Kullback-Leibler divergence, KL(q|p(Z|X, Θ)), between q and the posterior distribution of the hidden variables knowing the data points p(Z|X, Θ). This equality allows one to derive the EM algorithm when the posterior distribution of the hidden variables knowing the data points can be calculated exactly. The variational approach consists in replacing this distribution by a simpler one when it cannot be calculated.

In standard mixture models (such as the multivariate Gaussian mixture), the variational approach is not useful as the posterior distribution of the hidden variables can be calculated. However [START_REF] Verbeek | Self-organizing mixture models[END_REF] propose nevertheless to use the variational approach as a way to enforce regularity in the mixture model. Rather than allowing p(Z|X, Θ) to take an arbitrary form, they constrain it to a subset of probability distributions on the hidden variables that fulfill topological constraint corresponding to the prior structure of the SOM. See [START_REF] Verbeek | Self-organizing mixture models[END_REF] for details.

This solution shares most of the advantages of the older proposal in Heskes [2001], with the added value of being based on a more general principle that can be applied to any mixture model (in practice, [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] makes sense only for the exponential family). In addition, [START_REF] Verbeek | Self-organizing mixture models[END_REF] study the effects of shrinking the neighborhood function during training and conclude that it improves the quality of the solutions. Notice that, in [START_REF] Verbeek | Self-organizing mixture models[END_REF], the shared precision of the Gaussian distributions (β) is not a meta-parameter as in [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] but a regular parameter that is learned from the data.

The Generative Topographic Mapping The Generative Topographic Mapping (GTM, [START_REF] Bishop | GTM: The generative topographic mapping[END_REF]) is frequently presented as a probabilistic version of the SOM. It is rather a mixture model inspired by the SOM rather than an adaptation. Indeed the aim of the GTM designers was not to recover a learning algorithm close to a SOM variant, but rather to introduce a mixture model that enforce topology preservation.

The GTM is based on uniform prior distribution on a fixed grid which is mapped via an explicit smooth nonlinear mapping to the data space (with some added isotropic Gaussian noise). It can be seen as a constrained Gaussian mixture, but with yet another point of view compared to [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] and [START_REF] Verbeek | Self-organizing mixture models[END_REF]. In [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF], the constraint is enforced by a regularization term on the data space distribution while in [START_REF] Verbeek | Self-organizing mixture models[END_REF] the constraint is induced at the latent variable level (via approximating p(Z|X, Θ) by a smooth distribution). In the GTM the constraint is induced on the data space distribution because it is computed via a smooth mapping. In other words, the centers of the Gaussian distributions are not freely chosen but rather obtained by mapping a fixed grid to the data space via the nonlinear mapping.

The nonlinear mapping is in principle arbitrary and can therefore implement various type of regularity (i.e. topology constraints). The use of Gaussian kernels lead to constraints that are quite similar to the SOM constraints. Notice that those Gaussian kernels are not to be confused with the isotropic Gaussian distributions used in the data space (the same confusion could arise in [START_REF] Verbeek | Self-organizing mixture models[END_REF] where Gaussian kernels can be used to specify the constraints on p(Z|X, Θ)).

Once the model has been specify (by choosing the nonlinear mapping), its parameters are estimated via an EM algorithm. The obtained algorithm is quite different from the SOM (see [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] for details), at least in its natural formulation. However the detailed analysis contained in [START_REF] Heskes | Self-organizing maps, vector quantization, and mixture modeling[END_REF] shows that the GTM can be reformulated in a way that is close to the batch SOM with probabilistic assignments (as in e.g. the STM). Once again, however, this is not exactly the same algorithm. In practice, the results on real world data can be quite different. Also, as all the probabilistic variants discussed in this section, the GTM benefits from the probabilistic setting that enables principled missing data analysis as well as easy extensions to the exponential family of distributions in order to deal with non numerical data.

Non numerical data

When the data are not numerical, the SOM algorithm has to be adapted. See for example [START_REF] Kohonen | Median strings[END_REF], [START_REF] Kohonen | Self-organizing maps of symbol strings[END_REF], Kaski et al. [1998a], [START_REF] Kohohen | Self-organizing maps of symbol strings[END_REF], [START_REF] Kohonen | Self-Organizing Maps, 3rd Edition[END_REF], [START_REF] Kohonen | How to make large self-organizing maps for nonvectorial data[END_REF], [START_REF] Kohonen | Essentials of self-organizing map[END_REF], [START_REF] Kohonen | MATLAB Implementations and Applications of the Self-Organizing Map[END_REF], [START_REF] Cottrell | Neural networks for complex data[END_REF], where some of these adaptations are presented. Here we deal with categorical data collected in surveys and with abstract data which are known only by a dissimilarity matrix or a kernel matrix.

Contingency table or complete disjunctive table

Surveys collect answers of the surveyed individuals who have to choose an answer to several questions among a finite set of possible answers. The data can consist in a simple contingency table, where there are only two questions, and where the entries are the numbers of individuals who choose a given pair of categories, -a Burt table, that is a full contingency table between all the pairs of categories of all the questions, -a complete disjunctive table that contains the answers of all the individuals, coded in 0/1 against dummy variables which represent all the categories of all the questions.

In all these settings, the data consist in a positive integer-valued matrix, which can be seen as a large "contingency table". In classical data analysis, one uses Multiple Correspondence Analysis (MCA) that are designed to deal with these tables. MCA is nothing else than two simultaneous weighted Principal Component Analysis (PCA) of the table and of its transposed, using the χ2 distance instead of the Euclidean distance. To use a SOM algorithm with such tables, it is therefore sufficient to apply a transformation to the data, in order to take into account the χ 2 distance and the weighting, in the same way that it is defined to use Multiple Correspondence Analysis. After transforming the data, two coupled SOM using the rows and the columns of the transformed table can thus be trained. In the final map, related categories belong to the same cluster or to neighboring clusters. The reader interested by a detailed explanation of the algorithm can refer to Bourgeois et al. [2015a] . More details and real-world examples can also be found in [START_REF] Cottrell | Som-based algorithms for qualitative variables[END_REF], [START_REF] Cottrell | How to use the Kohonen algorithm to simultaneously analyse individuals in a survey[END_REF]. Notice that the transformed tables are numerical data tables, and so there is no particular theoretical results to comment on. All the results that we presented for numerical data still hold.

Dissimilarity Data

In some cases, complex data such as graphs (social networks) or sequences (DNA sequences) are described through relational measures of resemblance or dissemblance, such as kernels or dissimilarity matrices. For these general situations, several extensions of the original algorithm, both in on-line and batch versions, were proposed during the last two decades. A detailed review of these algorithms is available in [START_REF] Rossi | How many dissimilarity/kernel self organizing map variants do we need?[END_REF].

More precisely, these extensions consider the case where the data are valued in an arbitrary space X , which is not necessarily Euclidean. The observations are described either by a pairwise dissimilarity D = (δ(x i , x j )) i,j=1,...,N , or by a kernel matrix K = (K(x i , x j )) i,j=1,...,N 2 . The kernel matrix K naturally induces an Euclidean distance matrix, but the dissimilarity matrix D may not necessarily be transformed into a kernel matrix.

The first class of algorithms designed for handling relational data is based on the median principle (median SOM ): prototypes are forced to be equal to an observation, or to a fixed number of observations. Hence, optimal prototypes are computed by searching through (x i ) i=1,...,N , instead of X , as in [START_REF] Kohohen | Self-organizing maps of symbol strings[END_REF], [START_REF] Kohonen | How to make large self-organizing maps for nonvectorial data[END_REF], El Golli et al. [2004a] and El Golli et al. [2004b]. The original steps of the algorithm are thus transformed in a discrete optimization scheme, which is performed in batch mode:

1. Equation ( 1) is replaced by the affectation of all data to their best matching units: c(x i ) = arg min k=1,...,K δ(x i , m k (t)); 2. Equation ( 2) is replaced by the update of all prototypes within the dataset (x i ) i=1,...,N : m k (t) = min xi : i=1,...,N N j=1 h c(xi)j (t)δ(x i , x j ). Since the algorithm explores a finite set when updating the prototypes, it is necessarily convergent to a local minimum of the energy function. However, this class of algorithms exhibits strong limitations, mainly due to the restriction of the prototypes to the dataset, in particular, a large computational cost (despite efficient implementations such as in [START_REF] Conan-Guez | Fast algorithm and implementation of dissimilarity self-organizing maps[END_REF]) and no interpolation effect which yields to a deterioration of the quality of the map organization.

The second class of algorithms, kernel SOM and relational/dissimilarity SOM, rely on expressing prototypes as convex combinations of the input data. Although these convex combinations do not usually have sense in X (consider, for instance, that input data are various texts), a convenient embedding in an Euclidean or a pseudo-Euclidean space gives a sound theoretical framework and gives sense to linear combinations of inputs.

For kernel SOM, it is enough to use the kernel trick as given by Aronszajn [1950] which prove that there exists a Hilbert space H, also called feature space, and a mapping ψ : X → H, called feature map, such that K(x, x ′ ) = ψ(x), ψ(x ′ ) H . In the case where data are described by a symmetric dissimilarity measure, they may be embedded in a pseudo-Euclidean space ψ : x ∈ X → ψ(x) = (ψ + (x), ψ -(x)) ∈ E, as suggested in [START_REF] Goldfarb | A unified approach to pattern recognition[END_REF]. E may be written as the direct decomposition of two Euclidean spaces, E + and E -, with a non-degenerate and indefinite inner product defined as

ψ(x), ψ(y) E = ψ + (x), ψ + (u) E + -ψ -(x), ψ -(u) E -
The distance naturally induced by the pseudo-Euclidean inner product is not necessarily positive.

For both kernel and relational/dissimilarity SOM, the input data are embedded in H or E and prototypes are expressed as convex combinations of the images of the data by the feature maps. For example, in the kernel case,

m k (t) = N i=1 γ t ki ψ(x i ) , with γ t ki ≥ 0 and i γ t ki = 1.
The above writing of the prototypes allows the computation of the distance from an input x i to a prototype m k (t) in terms of the coefficients γ t ki and the kernel/dissimilarity matrix only. For kernel SOM, one has

ψ(x i ) -m k (t) 2 = γ t k T Kγ t k -2K i γ t k + K ii ,
where K i is the ith row of K and (γ t k ) T = γ t k,1 , ..., γ t k,N . For relational/dissimilarity SOM, one obtains a similar expression

ψ(x i ) -m k (t) 2 = D i γ t k - 1 2 γ t k T Dγ t k .
The first step of the algorithm, finding the best matching unit of an observation, as introduced in Equation ( 1), can thus be directly generalized to kernels and dissimilarities, both for on-line and batch settings.

In the batch framework, the updates of the prototypes are identical to the original algorithm (see Equation ( 12)), by simply noting that only the coefficients of the x i 's (or of their images by the feature maps) are updated :

m k (t + 1) = N i=1 h kc t (xi) (t) N j=1 h kc t (xj ) (t) ψ(x i ) ⇔ γ t+1 ki = h kc t (xi) (t) N j=1 h kc t (xj ) (t) (17)
This step is the same, both for batch kernel SOM, [START_REF] Villa | A comparison between dissimilarity SOM and kernel SOM for clustering the vertices of a graph[END_REF], and for batch relational SOM, [START_REF] Hammer | Topographic processing of relational data[END_REF].

In the on-line framework, updating the prototypes is similar to the original algorithm, as in Equation ( 2). Here also, the update rule concerns the coefficients γ t ki only, and the linear combination of them remains convex :

γ t+1 k = γ t k + ε(t)h kc t (xi) (t) 1 i -γ t k , (18)
where x i is the current observation and 1 i is a vector in R N , with a single non-null coefficient, equal to 1, on the i-th position. As previously, this step is identical for on-line kernel SOM, Mac Donald and [START_REF] Donald | The kernel self organising map[END_REF], and for on-line relational SOM, [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF].

In the case where the dissimilarity is the squared distance induced by the kernel, kernel SOM and relational SOM are strictly equivalent. Moreover, in this case, they are also fully equivalent to the original SOM algorithm for numerical data in the feature (implicit) Euclidean space induced by the dissimilarity or the kernel, as long as the prototypes are initialized in the convex hull of the input data. The latter assertion induces that the theoretical limitations of the original algorithm also exist for the general kernel/relational versions. Furthermore, these may worsen for the relational version since the non-positivity of the dissimilarity measure adds numerical instability when using a gradient-descent like scheme for updating the prototypes.

The third class of algorithms uses the soft topographic maps setting introduced in section 6.2. Indeed, in the algorithm described in equations ( 15) and ( 16), the soft assignments depend on the distances between input data and prototypes only, while prototypes update consists in making an update if the coefficients of the input data. Using a mean-field approach and similarly to the previous framework for kernel and dissimilarity/relational SOM, [START_REF] Graepel | Self-organizing maps: generalizations and new optimization techniques[END_REF] obtain the extensions of soft topographic mapping (STM) algorithm for kernels and dissimilarities. The updates for the prototype coefficients are then expressed as

γ ki (t + 1) = K j=1 h jk (t)P(x i ∈ A j ) N l=1 K j=1 h jk (t)P(x l ∈ A j ) , (19)
where m k (t) = N i=1 γ t ki ψ(x i ) and ψ is the feature map.

8 Stochasticity of the Kohonen maps for the on-line algorithm

Starting from a given initialization and a given size of the map, different runs of the on-line stochastic SOM algorithm provide different resulting maps. On the contrary, the batch version of the algorithm is a deterministic algorithm with always provides the same results for a given initialization. For this reason, the batch SOM algorithm is often preferred over the stochastic one because its results are reproducible. However, this hides the fact that all the pairs of observations which are associated in a given cluster do not have the same significance. More precisely, interpreting a SOM result, we can use the fact that close input data belong to close clusters, i.e. their best matching units are identical or adjacent. But if two given observations are classified in the same or in neighboring units of the map, then they may not be close in the input space. This drawback comes from the fact that there is not perfect matching between a a multidimensional space and a one-or two-dimensional map. More precisely, given a pair of observations (data), {x i , x j }, three cases can be distinguished, depending on the way their respective mapping on the map can be described:

significant association: the pair is classified in the same cluster or in neighboring clusters because x i and x j are close in the input space. The observations are said to attract each other; -significant non-association: the pair is never classified in neighboring clusters and x i and x j are remote in the input space. The observations are said to repulse each other; -fickle pair : the pair is sometimes classified in the same cluster or in neighboring clusters but x i and x j are not close in the input space: their proximity on the map is due to randomness.

The stochasticity of the on-line SOM results can be used to precisely qualify every pairs of observations by performing several runs of the algorithm. The question is addressed in a bootstrap framework in de [START_REF] De Bodt | Statistical tools to assess the reliability of self-organizing maps[END_REF] and used for text mining applications in Bourgeois et al. [2015a,b]. The idea is simple: since the on-line SOM algorithm is stochastic, its repetitive use may allow to identify the pairs of data in each case. More precisely, if L is the number of different and independent runs of the on-line SOM algorithm and if Y i,j denotes the number of times x i and x j are neighbors on the resulting map in the L runs, a stability index can be defined: for the pair (x i , x j ), this index is equal to:

M i,j = Y i,j L .
Using an approximation of the binomial distribution that would hold if the data were neighbors by chance in a pure random way, and a test level of 5%, for a K-units map, the following quantities are introduced

A = 9 K and B = 1.96 9 KL (1 - 9 K ). ( 20
)
These values give the following decision rule to qualify every pair {x i , x j }:

if M i,j > A+B, the association between the two observations is significantly frequent; -if A -B ≤ M i,j ≤ A + B, the association between the two observations is due to randomness. {x i , x j } is called a fickle pair; -if M i,j < A -B, the non-association between the two observations is significantly frequent.

In de [START_REF] De Bodt | Statistical tools to assess the reliability of self-organizing maps[END_REF], the method is used in order to qualify the stability and the reliability of the global Kohonen map, while both other papers (Bourgeois et al. [2015a] and Bourgeois et al. [2015b]) study the fickle data pairs for themselves. In these last work, the authors also introduce the notion of fickle word which is defined as an observation which belongs to a huge number of fickle pairs by choosing a threshold.

These fickle pairs and fickle words can be useful in various way: first, fickle pairs can be used to obtain more robust maps, by distinguishing stable neighboring and non neighboring pairs from fickle pairs. Also, once identified, fickle words can be removed from further studies and representations: for instance, Factorial Analysis visualization is improved. In a text mining setting, Bourgeois et al. [2015b] have shown that a graph of co-occurrences between words can be simplified by removing fickle words and Bourgeois et al. [2015a] have used the fickle words for interpretation: they have shown that the fickle words form a lexicon shared between the studied texts.

Conclusion

We have reviewed some of the variants of the SOM, for numerical and non numerical data, in their stochastic (on-line) and batch versions. Even if a lot of theoretical properties are not rigorously proven, the SOM algorithms are very useful tools for data analysis in different contexts. Since the Heskes's variants of SOM have a more solid theoretical background, SOM can appear as an easyto-develop approximation of these well-founded algorithms. This remark should ease the concern that one might have about it.

On a practical point of view, SOM is used as a statistical tool which has to be combined with other techniques, for the purpose of visualization, of vector quantization acceleration, graph construction, etc. Moreover, in a big data context, SOM-derived algorithms seem to have a great future ahead since the computational complexity of SOM is low (proportional to the number of data). In addition, it is always possible to train the model with a sample randomly extracted from the database and then to continue the training in order to adapt the prototypes and the map to the whole database. As most of the stochastic algorithm, SOM is particularly well suited for stream data (see [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF] which proposed a "patch SOM" to handle this kind of data). Finally, it would also be interesting to have a look at the robust associations revealed by SOM, to improve the representation and the interpretation of too verbose and complex information.

  Fig.2. This figure represents 2-dimensional prototypes (x-and y-axes are not shown but are the standard horizontal and vertical axes) which are linked as their corresponding units on the SOM grid. At this step of the algorithm, the x-and y-coordinates of the prototypes are well ordered. But contrarily to the one-dimensional case, this disposition can be disordered with a positive probability; in an 8-neighbors case, A is C's neighbor, but B is not a neighbor of C. If C is repeatedly the best matching unit, B is never updated, while A becomes closer and closer to C. Finally, the y coordinate of A becomes smaller than that of B and the disposition is disordered.

Table 1 .

 1 Comparison summary.

		on-line stochastic batch deterministic
	No neighbor VQ, SCL, k-means Forgy, moving centers
	With neighbors	SOM	batch SOM

EM is the standard algorithm for mixture models.

A kernel is a particular case of symmetric similarity such that K is a symmetric matrix, semi-definite positive with K(xi, xi) = 0 and satisfies the following positive constraint∀ M > 0, ∀ (xi)i=1,...,M ∈ X , ∀ (αi)i=1,...,M , i,j αiαj K(xi, xj) ≥ 0.


Madalina Olteanu
email: madalina.olteanu@univ-paris1.fr

Nathalie Villa-Vialaneix
email: nathalie.villa@toulouse.inra.fr

Sparse online self-organizing maps for large relational data

Keywords: relational data, online relational SOM, sparse approximations

During the last decades, self-organizing maps were proven to be useful tools for exploring data. While the original algorithm was designed for numerical vectors, the data became more and more complex, being frequently too rich to be described by a fixed set of numerical attributes. Several extensions of the original SOM were proposed in the literature for handling kernel or dissimilarity data. Most of them use the entire kernel/dissimilarity matrix, which requires at least quadratic complexity and becomes rapidly unfeasible for 100 000 inputs, for instance. In the present manuscript, we propose a sparse version of the online relational SOM, which sequentially increases the composition of the prototypes.

Introduction

The self-organizing map (SOM) algorithm, [START_REF] Kohonen | Self-Organizing Maps, 3rd Edition[END_REF], was proven, over the years, to be a powerful and convenient tool for clustering and visualizing data. While the original algorithm was 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 ressemblance/dissemblance (professional trajectories). During the past twenty years, the SOM algorithm was extended for handling relational data, either described by kernels (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), or by dissimilarities (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). All these extensions are based on the same underlying principle : the dissimilarity or the kernel implicitly define an Euclidean (or pseudo-Euclidean) space in which the prototypes can be expressed as convex combinations of the embedded input data. However, when the goal is to explore large data sets, the relational approaches may become rapidly unfeasible. Indeed, complex relational data often have a large dimensionality. Moreover, kernel and relational SOM rely on the knowledge of the dissimilarity matrix for the entire data set, which generates at least quadratic complexity for the algorithms. As stressed in [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF], algorithms will be slow for data sets with 10,000 observations and impossible to run on a normal computer for 100,000 input data. In addition to the complexity issue, expressing prototypes as convex combinations of the entire data set has a second drawback, as emphasized in [START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF]: the interpretability of the prototypes and of the model is lost. In order to tackle these two issues, several approaches were introduced for relational data, all of them seeking for a sparse representation of the prototypes and a linear (in the number of observations) computational cost. [START_REF] Rossi | Accelerating relational clustering algorithms with sparse prototype representation[END_REF] use the natural sparsity of the prototypes in batch relational k-means in order to reduce the complexity. The natural sparsity is enhanced by selecting the K (K fixed) closest inputs to each prototype. In [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF], the complexity is reduced using iterative "patch clustering". First, the data are split into P patches of size n P (P fixed). A prototype-based clustering algorithm in batch version (neural gas or SOM) is then run on a patch P t and the resulting prototypes, which may be viewed as compressed representations of the data already seen, are added as new data points to the next patch, P t+1 . Moreover, the full vector of coefficients is replaced by the K closest input data (K fixed). With this method, linear time and constant space representation are obtained. Another technique consists in using the Nyström approximation ( [START_REF] Gisbrecht | The nystrom approximation for relational generative topographic mappings[END_REF]) for the dissimilarity matrix. This technique also leads to a linear computational cost in the number of input data, but is strongly dependent on the intrinsic dimensionality of the given dissimilarity matrix, which has to be of low rank and entirely known in advance. All these cited approaches are batch algorithms. In the online framework, [START_REF] Mariette | Bagged kernel SOM[END_REF] propose a bagging approach for kernel SOM. Data is split into B subsamples of size n B (B fixed), the online kernel SOM is trained on each subsample and, after training, the most representative K observations are chosen for each prototype (K fixed). Eventually, a final map is trained on the resulting most representative observations. The algorithm has the advantage of being parallelizable, although it does not consider all the advantages of an online implementation. In the present paper, we propose a sparse version of the online relational SOM algorithm, which takes further advantage of the online setting. Instead of expressing prototypes as convex combinations of the entire data set from the beginning, the size and the composition of the prototypes are sequentially increased with each new input fed to the algorithm. When the size of the prototypes becomes too large, prototypes are made sparse by deleting all the insignificant coefficients. Different approaches for selecting the most interesting observations are reported in [START_REF] Hofmann | Learning interpretable kernelized prototype-based models[END_REF]. In this manuscript, we use a slightly different technique, by interpreting the coefficients as a probability distribution and by selecting the most probable observation: a global probability mass ν is fixed and the largest coefficients summing to ν are kept. In this way, more flexibility is allowed to the prototypes which are no longer represented by a fixed number K of observations, but by the necessary number of observations allowing an "almost complete" knowledge of the composition of the prototypes (if ν is chosen close to 1). The rest of the paper is organized as follows: Section 2 recalls the online relational SOM, while Section 3 introduces the sparse version of the online relational SOM. The equivalent algorithm for kernels is briefly described in Section 4, while Section 5 contains some examples on real data-sets.

Online relational SOM

In this section we shall briefly recall the principles of the online relational SOM (RSOM) algorithm, as introduced in [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF]. Throughout the rest of the paper, let us suppose that the input data, x 1 , ..., x N , belong to some arbitrary space G and can be described through a dissimilarity measure δ, such that δ ij = δ (x i , x j ). The dissimilarity measure is supposed to verify some basic assumptions: symmetry (δ ij = δ ji ) and non-negativity (δ ij ≥ 0), for all i, j = 1, ..., N , and also δ ii = 0, for all i = 1, ..., N . The online RSOM 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 linked together by a neighborhood relationship H, expressed as a function of the distance between the units on the grid, d (u, u ′ ). The distance on the grid, d, may be chosen, for example, as the length of the shortest path between the units. The U prototypes are initialized either as random convex combinations of the input data or randomly among the input data. The extension of the original SOM algorithm is based on two key ideas:

-First, prototypes are written as (symbolic) convex combinations of the input data, p u = N i=1 β u,i x i , with β u,i ≥ 0 and N i=1 β u,i = 1, for all u = 1, ..., U . This definition is justified by the fact that, when a dissimilarity is given, it can be viewed as the dot product of the images by a mapping function φ into a pseudo-Euclidean space [START_REF] Goldfarb | A unified approach to pattern recognition[END_REF]: the prototypes are thus truly the convex combinations of (φ(x i )) i in this space (see [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF][START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF] for further explanations).

-Second, the distance between an input data x i and a prototype p u can be written only in terms of the dissimilarity matrix of the input data and the coefficients β u,i as follows:

x i -p u 2 = ∆ i β T u - 1 2 β u ∆β T u , (1)
where ∆ = (δ ij ) i,j=1,...,N , ∆ i represents the i-th row of the matrix ∆ and β u = (β u,1 , ..., β u,N ) is the vector of coefficients for the prototype p u .

Expressing the prototypes as convex combinations of the input data and computing the distances between observations and prototypes as in Equation (1) consists, in fact, in a generalization of the original SOM algorithm. Indeed, one can easily see that the two are equivalent if the dissimilarity δ is the squared Euclidean distance and if the prototypes of the original SOM are initialized within the convex hull of the input data. This general framework allowing an elegant writing of the algorithm for complex data described by dissimilarities was introduced initially for the online version of kernel SOM (KSOM) in [START_REF] Donald | The kernel self organising map[END_REF]. Afterwards, extensions and rediscoveries were described for batch relational SOM in [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF], batch kernel SOM in [START_REF] Boulet | Batch kernel SOM and related Laplacian methods for social network analysis[END_REF] and online relational SOM in [START_REF] Olteanu | On-line relational and multiple relational SOM[END_REF]. A detailed and complete comparison of these methods and their equivalences may be found in [START_REF] Rossi | How many dissimilarity/kernel self organizing map variants do we need?[END_REF].

The distance computation in Equation ( 1) may be theoretically justified in the very general setting of dissimilarities by extending the Hilbert embedding for kernels to a pseudo-Euclidean embedding, as shown, for example, in [START_REF] Hammer | Topographic mapping of large dissimilarity data sets[END_REF].

The online relational SOM algorithm is summarized in Algorithm 1. The neighborhood function H is supposed to verify the following assumptions: H : R → R, H(0) = 1 and lim x→+∞ H(x) = 0. In the setting of Algorithm 1, H t decreases piecewise linearly, while µ(t) vanishes at the rate 1 t .

Algorithm 1 Online relational SOM

1: For all u = 1, . . . , U and i = 1, . . . , N, initialize β 0 ui such that β 0 u,i ≥ 0 and

N i β 0 u,i = 1. 2: for t = 1, . . . , T do 3:
Randomly choose an input xi 4:

Assignment step: find the unit of the closest prototype

f t (xi) ← arg min u=1,...,U ∆i β t-1 u T - 1 2 β t-1 u ∆(β t-1 u ) T

5:

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

β t u ← β t-1 u + µ(t)H t (d(f t (xi), u)) 1i -β t-1 u
where 1i is a vector with a single non null coefficient at the ith position, equal to one. 6: end for 3 Sparse online relational SOM Similarly to relational SOM, prototypes are 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 significant of them. In order to guarantee the sparsity of the writing as well as similar properties with the original online relational SOM algorithm, several issues have to be verified.

1. Prototypes have to be initialized at random among the input data. Hence, the observations have to be randomly presented to the algorithm. The first U observations will be then used as initial values for the U prototypes. 2. The dissimilarity between a new input data and a prototype, written as a convex combination of the most significant past observations, has to be computed. This can be achieved using the following formula x k -

p u 2 = j∈I(t) β u,j δ (x k , x j ) -1 2 i∈I(t)
j∈I(t) β u,i β u,j δ (x i , x j ), where p u = j∈I(t) β u,j x j and I(t) contains the indices of the most significant inputs already fed to the algorithm before x k is chosen. 3. Prototypes are sparse combinations of the input data. Hence, prototypes are periodically updated and the most coefficients only are selected. The updates may be performed throughout the iteration using either a deterministic design (the number of updates is fixed and updates are uniformly distributed during the learning of the map), or a random design (the updates are distributed according to some geometric distribution. The parameter of the geometric distribution may depend on the total number of iterations and on the size of the neighborhood). 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 and interpretability of the prototypes, the most significant coefficients are selected according to their value, by fixing a threshold: let 0 < ν ≤ 1 be the selected threshold (if ν = 1, the algorithm is no longer sparse, but the original one). For u = 1, ..., U , the coefficients are ordered in descending order for each prototype β u,(1) , ..., β u,(♯I(t)) , where β u,(1) = max i∈I(t) β u,i and β u,(♯I(t)) = min i∈I(t) β u,i . Consider N u such that N u = arg min n=1,...,♯I(t) n i=1 β u,(i) ≥ ν . The most significant coefficients are updated as follows

β u,(i) =      β u,(i) Nu j=1 β u,(j) , if (i) ≤ N u 0 , if (i) > N u
The sparse online relational SOM algorithm is summarized in Algorithm 2.

The kernel version

In some cases, data may be described by a kernel, K, instead of a dissimilarity. We shall recall that a kernel is a symmetric similarity such that K (x i , x i ) = 0 and which satisfies the following positive constraint:

∀M > 0, ∀ (x i ) i=1,...,M ∈ G, ∀ (α i ) i=1,...,M ∈ R, M i,j=1 α i α j K (x i , x j ) ≥ 0.
According to [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], there exists a Hilbert space H, also called feature space, as well as a feature map ψ : G → H, such that K(x, x ′ ) = ψ(x), ψ(x ′ ) H . Similarly to the dissimilarity case, the prototypes are defined as convex combinations of (the images by ψ of) (x i ) i . The distance between an input data x k and some prototype p u is then computed as the squared distance induced by the kernel Algorithm 2 Sparse online RSOM 1: For all u = 1, . . . , U , initialize p 0 u among the first U input data: β 0 u = 1 U u , where 1 U u is a vector of length U with a single non-null coefficient on the u-th position, equal to 1. Initialize I(0) = {1, ..., U }. 2: for t = 1, . . . , T do 3:

Randomly choose an input x k , k ∈ {1, ..., N }.

4:

Assignment step: find the unit of the closest prototype

f t (x k ) ← arg min u=1,...,U   j∈I(t-1) β t-1 u,j δ (x k , xj) - 1 2 β t-1 u ∆ I(t-1) β t-1 u T   ,
where ∆ I(t-1) = (δ (xi, xj )) i,j∈I(t-1) . 5:

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

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

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

β t u ← 1 -µ(t)H t (d(f t (x k ), u)) β t-1 u , 0 + µ(t)H t (d(f t (x k ), u))(0, ..., 0 ♯I(t-1)
, 1)

11:

I(t) = I(t -1) ∪ {k}. 12:
end if 13:

Sparse representation: 14:

if t is an update instant (deterministic or random design) then 15:

Sparsely update the prototypes: ∀ u = 1, . . . , U ,

β t u,(1) ≥ ... ≥ β t u,♯I(t) , Nt,u = arg min n=1,...,♯I(t) n i=1 β t u,(i) ≥ ν β t u,(i) =        β t u,(i) N t,u j=1 β t u,(j) , if (i) ≤ Nt,u 0 , if (i) > Nt,u 16: end if 17: end for x k -p u 2 = K (x k , x k ) -2 i∈I(t) β u,i K (x k , x i ) + i,j∈I(t) β u,i β u,j K (x i , x j ).
The sparse online relational SOM can thus be immediately adapted for kernels. Algorithm 2 has to be modified only in the assignment step which becomes 1: Assignment step: find the unit of the closest prototype

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

Examples

The sparse version introduced in the present manuscript was compared to the online relational SOM on two real data sets. For the sparse version, several values were considered for the threshold ν. The sparse updates were performed either in a uniform deterministic design (fixed number of updates), or at random, according to a geometric distribution. The performances of the sparse RSOM and the online RSOM were then compared in terms of average computational time (in seconds), quantization and topographic errors and sparsity (number of non-zero coefficients). Scripts were all implemented under the free statistical software environement R. 1. Average results for Astraptes fulgerator (100 random initializations). The first column contains the number of updates (deterministic design). The third column is the computational time (provided in seconds). The last column is the average number of non zero coefficients in the prototypes. The bolded values correspond to the results at least as good as the online RSOM.

Astraptes fulgerator. The first data set was introduced in [START_REF] Hebert | Ten species in one: DNA barcoding reveals cryptic species in the neotropical skipper butterfly astraptes fulgerator[END_REF]. In contains information on 465 Amazonian butterflies, each of them described by a sample of their DNA. Each input data is a DNA sequence of length 350. The Kimura distance for genetical sequences, as introduced in [START_REF] Kimura | A simple method for estimating evolutionary rates of base substitutions through comparative studies of nucleotide sequences[END_REF], was computed and the resulting distance matrix was used as input for relational and sparse relational SOM. For both algorithms, 100 different initializations with 2 500 iterations each were performed on a square grid of size 5 × 5. The results are summarized in Tables 1 and2 for the deterministic and random designs respectively. Professional trajectories. The second example comes from [START_REF] Rousset | Classifying qualitative time series with som: The typology of career paths in france[END_REF]. It contains information about 2 000 people having graduated high-school in 1998 and monitored during 94 months afterwards. For each individual, a categorical sequence of length 94, giving his monthly professional status is available. In all, there are nine possible situations, from permanent contracts to unemployment. The dissimilarity used for these data is the optimal matching (OM) distance, as introduced in [START_REF] Abbott | Optimal matching methods for historical sequences[END_REF]. Here, 100 different initializations with 10 000 iterations each were performed on a square grid of size 10 × 10. The sparse version was compared to the standard online relational SOM (itself run from 100 different initializations and 10 000 iterations). The results for the deterministic design are summarized in Table 3 (due to the lack of space, we do not report here the results with a random design, which are quite similar). It is interesting to note that the sparsity has a strong influence on the computational time: increasing the number of updates tends to decrease the computational time since the prototypes are regularly cleared from unnecessary coeffi-cients. The computational time compared to the standard version is at least 10 times smaller in the sparse version for this large dataset. On the contrary, the performances, measured in terms of quantization and topographic errors, can be affected by a too large sparsity but the best ones remain close to those of the standard version. 3. Average results for "professional trajectories" (100 random initializations, updates were made with a deterministic design). Simulations were all performed on a server with OS Debian 8 Jessie, 8 processors AMD Opteron 8384 with 4 cores each and 256 Go RAM.

Conclusion and future work

A sparse version of the online relational SOM algorithm was proposed, by sequentially increasing the composition of the prototypes and sparsely updating them. The algorithm was compared with the online ROM on two real data sets and the sparse version appeared to achieve very similar performances as compared to the original algorithm, while improving computational time and prototype representation.

Table 2 .

 2 Average results for Astraptes fulgerator (100 random initializations, updates were made with a random design).

	nb. updates	ν comp. time (s) quantization err. topographic err. nb. coefs
	50	0.80	1.92	0.00093	0.0353	5.44
	50	0.85	2.09	0.00078	0.0176	7.35
	50	0.90	2.37	0.00069	0.0145	11.02
	50	0.95	2.92	0.00067	0.0102	21.75
	50	0.99	4.02	0.00068	0.0068	45.51
	25	0.80	2.50	0.00067	0.0210	9.92
	25	0.85	2.88	0.00066	0.0114	14.09
	25	0.90	2.94	0.00066	0.0107	20.41
	25	0.95	3.56	0.00064	0.0057	29.63
	25	0.99	4.66	0.00066	0.0053	51.93
	10	0.80	4.23	0.00062	0.0132	22.48
	10	0.85	4.69	0.00065	0.0072	28.41
	10	0.90	5.18	0.00065	0.0098	33.97
	10	0.95	5.14	0.00065	0.0051	43.34
	10	0.99	6.30	0.00067	0.0033	59.95
	Online RSOM		12.18	0.00067	0.0051