Francesco Cambi
email: fcambi@bridgeconsulting.it

Pierluigi Crescenzi
email: pierluigi.crescenzi@unifi.it

Linda Pagli
email: linda.pagli@unipi.it

Analyzing and Comparing On-Line News Sources via (Two-Layer) Incremental Clustering

Keywords: 1998 ACM Subject Classification H.2.8 Database Applications: Data mining, H.3.3 Information Search and Retrieval: Clustering, H.3.5 Online Information Services: Web-based services Keywords and phrases text mining, incremental clustering, on-line news Digital Object Identifier 10.4230/LIPIcs.FUN.2016

In this paper, we analyse the contents of the web site of two Italian news agencies and of four of the most popular Italian newspapers, in order to answer questions such as what are the most relevant news, what is the average life of news, and how much different are different sites. To this aim, we have developed a web-based application which hourly collects the articles in the main column of the six web sites, implements an incremental clustering algorithm for grouping the articles into news, and finally allows the user to see the answer to the above questions. We have also designed and implemented a two-layer modification of the incremental clustering algorithm and executed some preliminary experimental evaluation of this modification: it turns out that the two-layer clustering is extremely efficient in terms of time performances, and it has quite good performances in terms of precision and recall.

Introduction

The web is a huge source of data, which are produced by companies, institutions and individuals, and, most of the times, are available for free. The use we can make of all this information is limited just by our imagination. In the last few years, for example, there has been a quite significant amount of research devoted to the analysis of the so-called on-line social networks. One of the most recent examples of such analysis is the adaptation of the well-known "six degrees of separations" phenomenon to the Facebook network: it has, indeed, has been observed that, in the case of this on-line social network, the degrees of separation are less than four [START_REF] Edunov | Three and a half degrees of separation[END_REF]. In this paper, instead, we show how the information available on the web sites of news agencies and newspapers can be used in order to answer several questions about the news on-line system, such as the following ones.

What are the most relevant news? Clearly, answering to this question in an automatic way implies defining the notion of relevance, that is, a news score function. In this paper, we propose one possible definition (based on the position occupied by the news in the web sites) and apply it to all the news collected from mid October until mid December of the last year. It turns out that the news with the maximum score is the shooting down of an American drone in Turkey. Among the news with the highest score, there is one concerning a case of homicide. The score pattern of this latter news is quite interesting since it has three main peaks: one peak corresponds to the disappearance of the victims, the second peak corresponds to the rising of the hypothesis of homicide, and the third peak corresponds to the identification of possible suspects. We believe that it is interesting to analyse this kind of patterns because they help us to better understand the evolution of a news. What is the average lifespan of a news? Similar questions have already been posed and answered in different contexts. For instance, a network of hundreds of media sites have been recently looked into, and it has been found that most posts had a shelf life of 2 days, with the median stretching out to about 2.5 days [START_REF] Parse | What is the Lifespan of an Article?[END_REF]. Closer to our analysis are the results presented in [START_REF] Bhattacharya | Sharing News Articles Using 140 Characters: A Diffusion Analysis on Twitter[END_REF], where the authors study the diffusion patterns of news articles from several popular news sources, determine the lifespan of a news article on Twitter by the time difference between the last and first tweet posted containing the URL to that article, and show that for most of the news media companies, about 45% of their articles survive beyond 18 hours. In this paper, we analyse almost 10000 articles, grouped into approximately 7000 news, and we observe that the average lifespan of a news is a little bit less than one day. The news with the longest lifespan concerns vaccination campaigns in Italy (with a lifespan greater than 19 days), while one of the many news with the shortest lifespan concerns the rescue of two bear cubs in Laos (with a lifespan of at most one hour).

Are different news sources really different? Once again, answering this question implies defining a notion of similarity among the on-line news sources. In this paper, we propose such a definition, and we apply it to the six on-line sources we have analysed. Interestingly enough, it turns out that the two news agencies are quite similar, one newspaper seems to echo the two news agencies, while the other thre newspapers are equally and significantly distant from the two news agencies, but in opposite directions.

In order to answer the above questions, we have developed a web-based tool which allows us to perform the following operations.

Web scraping. The tool can download the articles (formed by the title, the abstract, and the text) from the web site of two Italian news agencies, and of four popular Italian on-line newspapers. Almost 10000 articles have been downloaded from mid October 2015 to mid December 2015. These articles have been appropriately processed in order to be subsequently analysed in terms of their similarity.

Clustering. As we said, we have defined a similarity measure between two articles, which is based on quite standard text analysis methods. By using these similarity values, a simple incremental clustering algorithm has been implemented: this algorithm basically insert a new article into the cluster with the highest average similarity, if this average similarity is greater than a specific threshold (otherwise a new cluster is created containing only the new article). This incremental clustering algorithm (which is different from previous incremental clustering approaches applied to news detection such as the one analysed in [START_REF] Azzopardi | Incremental Clustering of News Reports[END_REF][START_REF] Borglund | Event-Centric Clustering of News Articles[END_REF]), performs very well in terms of precision and recall results. However, the time performance of the algorithm degrades as the number of articles (and, hence, of clusters) increases. To improve the time performances, we have designed and implemented a two-layer variation of the incremental clustering algorithm (which is also different from similar two-layer approaches such as the one used while analysing Twitter messages in [START_REF] Petkos | Two-level Message Clustering for Topic Detection in Twitter[END_REF]). In this variation, a centroid is associated with each cluster, and the centroids are F. Cambi, P. Crescenzi, and L. Pagli

9:3

clustered by making use of the same approach (that is, by using the average similarity values). Once the best cluster of centroids has been determined, the original incremental clustering is applied to the corresponding clusters. As far as we know, this two-layer approach is different from all previous techniques, which have been used to cluster news. Even if its performances in terms of precision degrade of about 10% the time performance of the new approach are drastically better: indeed, clustering approximately 3000 articles requires less than two seconds, while the original approach required more than 2 minutes.

Statistical analysis.

In order to answer the first two questions listed above, we have defined a score function of news, which is based on the position of the articles, included in the same news, in the front page of the web site of the on-line news sources. The position in the front page (along with variation of the Kendall tau measure) has been used also to define similarity measure between different newspapers, and, hence, to answer the third question. Apart from the fact that our clustering algorithms are different from previous approaches used in news detection (even if, due to the huge literature in the field of incremental clustering, they are similar to several other clustering approaches), as far as we know, this is the first time that such techniques are used to analyse and explicitly compare on-line news sources, on the ground of the order in which news appear in the web sites of the news sources themselves.

The paper is structured as follows. In Section 2, we describe the web scraping performed by our tool, and we present some preliminary descriptive statistics. In Section 3 we introduce our incremental clustering algorithm in order to detect news from articles. In Section 4 we present the results about the score and the lifespan of news, and about the comparison of the six analysed news on-line sources. In Section 5 we describe the two-layer variation of the clustering algorithm. Finally, in Section 6 we conclude and propose some possible directions for future research.

The Web Scraping Process

Web scraping is a technique to extract data from web sites, by making use of software libraries that simulates surfing on the web. Fortunately, most of the on-line news sources use well-defined XML schemata and tag systems, which makes the content independent of the graphical rendering of the web site. This feature (which is usually part of what is called the "web semantic") allowed us to quite easily download the contents of the articles which were present on the main web page of the analysed news sources. For instance, in the left part of Figure 1 it is shown part of the main page of one of the most popular Italian on-line newspaper. The XML schema used by this page easily allowed us to detect the main article column (rounded in green), and, within this column, the main article sub-column (rounded in red). In the figure, we have three articles, whose content can then be downloaded by using the hypertext link associated with the title of the article. In the current implementation of our tool, we have then decided to ignore the articles appearing in other sub-columns, since in the vast majority of the cases these articles seem to be of relatively smaller importance. However, it would not be too difficult to change the web scraping module of our tool in order to let it download also these articles.

In the vast majority of the cases, the content of an article consists of three components: the title, the abstract, and the text. In the example shown in the right part of Figure 1, the title is rounded in orange, the abstract is rounded in red, and the text is rounded in green. Once again, the XML schema used by the web site of the on-line news source allows us to easily detect these three components (by also eliminating all types of content we are not interested in, such as video and pictures). After having downloaded the title, the abstract, and the text of an article, this information is stored in a MySQL database, in order to be used for the news detection process that we are going to describe in the next section. This web scraping process has been started at the beginning of October 2015, and it is still in action: currently, we have downloaded more than 23000 articles.

Some descriptive statistics

As we just said, we have downloaded 23335 articles (during approximately five months) from the web sites of two news agencies NA1 and NA2, and of four newspapers NP1, NP2, NP3, and NP4. These articles are quite uniformly distributed among the six news sources, as it is shown in Figure 2. Note that the download of these articles has been done hourly (almost every day), but, clearly, we have avoided to download again articles that had been previously downloaded. In other words, the number of daily downloaded articles is not just 24 times the number of articles included in the home page of the web site of a news source, since the same articles can remain on the home page itself for several hours. Indeed, in Figure 3 we show the minimum, maximum, and average numbers of distinct daily downloaded articles from each of the six news sources. Once again, these numbers seem to uniformly distributed: this suggest that, on the average, approximately 40 new articles are produced every day in the main area of each home page. The activity of a news source will be better analysed in Section 4.1, in which we will try to measure the degree of variability of the home page of a news source in terms of news (that is, clusters of articles), instead that in terms of articles.

Another simple statistics that can be immediately computed, once a sufficiently large dataset of articles has been downloaded, is the word length of the articles themselves. This number does not really differs among the six news sources (apart form the second news agency that seems to produce, on average, shorter articles). This, indeed, can be deduced from Figure 4, where the average lengths of the titles, the abstracts, and the bodies of the six news sources are shown. More precisely, it seems that, on average, a title (respectively, an abstract, and a body) contains 10 (respectively, 20, and 400) words. The news agency NA2, however, uses shorter body lengths and no abstract at all. This uniformity among the article lengths can be, maybe, justified because of the standards normally used within the printed news system (mainly due to the limited number of pages of a newspaper): however, it is quite surprising that similar standards are also used within the on-line news system, where, in theory, no limit is a-priori existing on the number of pages an article should utilise. It is also worth noting that, if we consider the maximum word length (instead of the average one), the first newspaper NP1 is clearly producing much longer articles: indeed, its maximum article word length is 7216 which is almost twice the second longest maximum length (that is, 3730 in the case of the fourth newspaper NP4).

We have also computed the frequency of the words used in the title or in the abstract of an article. In particular, in Figure 5 of Nuovo Centro-Destra. More surprisingly, Grillo (129 occurrences), a popular comedian and blogger who is the leader of Movimento 5 Stelle, is only fifth in this ranking, preceded also by Salvini (242 occurrences), leader of Lega Nord. Interestingly enough, the second newspaper NP2 is systematically the one that most frequently includes the names of the first five leaders into the title or the abstract of an article, apart from the case of Grillo (as it is shown in Figure 6). On the other hand, the second news agency NA2 is the news source which almost always includes the name of the leaders less than the other news sources: in this case, however, we should again observe that this agency does not use abstracts at all. Finally, we also tried to verify somehow whether the following statement taken from Wikipedia corresponds to reality: "the major news agencies generally prepare hard news stories and feature articles that can be used by other news organisations with little or no modification" [13]. To this aim, we made use of the Levenshtein distance between two strings (that is, the minimum number of insertions, deletions and substitutions required to change one string into the other) [START_REF] Levenshtein | Binary codes capable of correcting deletions, insertions, and reversals[END_REF] in order to define a similarity measure between two string as follows: the difference between the length of the longest string and the Levenshtein distance, divided by the length of the longest string. We then computed for a (small) sample set of articles the similarity of their texts with the texts of the articles of the first news agency NA1. As a result of this comparison, we found a quite high similarity between the articles. The maximum found similarity has been, indeed, equal to 0.96. The two articles reaching this value reported about the first American democratic party debate, and their length was 4527 and 4555, respectively. This implies that the two articles were almost identical. The Clustering Algorithm As we said in the introduction, we focused our attention on news, instead of articles. Indeed, several different articles can refer to the same news, either because they are present on different web sites, or because the news evolved and new articles concerning the news itself have been produced. For this reason, we had to design a clustering of the articles into clusters corresponding to news. To this aim, we first have to define a similarity measure between articles, and subsequently apply a clustering algorithm on the ground of the similarity values.

The similarity measure

In order to define a similarity measure between articles, we made use of quite standard text mining techniques (see, for example, the first chapter of [START_REF] Leskovec | Mining of Massive Datasets[END_REF]). First of all, we cleaned the contents of the articles by eliminating the so-called "stop words", that is, the several hundred most common words in Italian that do not carry any significance by themselves (such as, for example, the articles). Successively, we applied a stemming algorithm for reducing derived words to their word root so that related words map to the same root, even if this root is not in itself a valid word (note that algorithms for stemming have been studied in computer science since the late sixties [START_REF] Lovins | Development of a Stemming Algorithm[END_REF]). Finally, we identified the keywords of an article a, that is, the words to be used in order to measure the similarity of a with other articles, as the words appearing in its title title(a) and in its abstract abstract(a), and the ones appearing in its text text(a) starting with a capital letter. In the following, we will denote by key(a) the set of keywords of article a. For each set of articles A and for each word w, the inverse document frequency of w with respect to A is defined as

idf(w, A) = log |A| |{a ∈ A : w ∈ text(a)}| .
Moreover, the term frequency of w with respect to an article a ∈ A is defined as tf(w, a) = occ(w, a) |a| where occ(w, a) denotes the number of occurrences of w in text(a), and |a| denotes the number of words in text(a). For each pair of articles a 1 and a 2 in the set A, we can then define their TF-IDF vectors as follows. Let key(a 1 , a 2 ) = key(a 1 ) ∪ key(a 2 ) = {k 1 , . . . , k m } be the ordered set of keywords of either a 1 or a 2 . Then, for any i with 1 ≤ i ≤ m, we define

tfidf a1,a2 [i] = tf(k i , a 1 ) • idf(k i , A)
and

tfidf a2,a1 [i] = tf(k i , a 2 ) • idf(k i , A).
Finally, the cosine similarity cosim(a 1 , a 2 ) between a 1 and a 2 is the cosine value of the angle formed by the two vectors tfidf a1,a2 and tfidf a2,a1 . Formally, it is defined as

cosim(a 1 , a 2 ) = m i=1 tfidf a1,a2 [i]tfidf a2,a1 [i] m i=1 tfidf a1,a2 [i] 2 m i=1 tfidf a2,a1 [i] 2 .
Note that this value is as close to 1 as the two vectors are similar, that is, as the two articles have similar TF-IDF vectors.

F U N 2 0 1 6

9:8

Analyzing and Comparing On-Line News Sources

The incremental clustering

By using the above defined similarity measure, we have designed and implemented the following quite simple incremental clustering algorithm. Assume that a set A of articles has already been clustered into a set N of n news or clusters. Indeed, N can be seen as a function N : A → [n], where [n] denotes the set of integers between 1 and n. Assume also that a new article x has been downloaded and, hence, has to be classified in either one of the already existing n clusters or in a new one. For each cluster c with c ∈ [n], let A c = {a ∈ A : N (a) = c}. We then define the similarity between x and c as the average of the similarities between x and all articles in c. More formally,

sim(x, c) = 1 |A c | a∈Ac cosim(x, a).
Let c * be the cluster c for which sim(x, c) is maximum: if sim(x, c * ) is greater than or equal to a given threshold τ , then N (x) = c * (and x is assigned to cluster c * ), otherwise N (x) = n + 1 (and a new cluster is created).

Threshold estimation

In order to apply the above clustering algorithm, we have to determine the threshold τ . To this aim, we have manually clustered the set T containing the first 3000 downloaded articles, and we then have determined the threshold value that produced the best results. In particular, given two clustering N 1 and N 2 with n 1 and n 2 clusters, respectively, we can define the matching function M N1,N2 : [n 1 ] → [n 2 ] as follows:

M N1,N2 (c 1 ) = max c2∈[n2] jac(c 1 , c 2 )
where jac(c 1 , c 2 ) denotes the Jaccard index of c 1 and c 2 , defined as |N1(c1)∩N2(c2)| |N1(c1)∪N2(c2)| . In order to evaluate the clustering N τ obtained by using threshold τ and applied to T , we have first computed the function M Nτ ,N * , where N * is the manual cluster, and subsequently computed the recall and precision values of each cluster in N τ . The precision of a cluster c is defined as

pre(c) = |{a ∈ T : a ∈ M Nτ ,N * (c)} ∩ T c | |T c | ,
while the recall of a cluster c is defined as

rec(c) = |{a ∈ T : a ∈ M Nτ ,N * (c)} ∩ T c | |{a ∈ T : a ∈ M Nτ ,N * (c)}| .
The recall and precision values can be combined by obtaining the F -measure, which is defined as follows:

F(c) = 2 p(c) • r(c) p(c) + r(c) .
The average weighted F -measure of a clustering N τ is then defined as the average value of the F -measure weighted with respect to the sizes of the clusters. Formally,

F(N τ ) = nτ i=1 F(c i )|c i | nτ i=1 |c i |
where n τ denotes the number of clusters in N τ and |c i | denotes the size of the i-th cluster (that is, |c i | = |{a ∈ T : N τ (a) = i}|). The threshold τ has then been determined by selecting the one who produced the higher average weighted F -measure. It turned out that the best value of τ is equal to 0.35: with this value of τ the average weighted F -measure is equal to 0.89.

Analyzing and Comparing the On-Line News Sources

Once the articles have been clustered, and, hence, the news have been detected, we can now answer to several interesting questions concerning the on-line news system.

Activity of a News Source

A front page f = (t, L) of a news source is a time stamped ordered list L of news, which are the news appearing in the home page of the web site of the news source at time t. We would like to measure what is the degree of variability of the front pages of the same news source during the day. To this aim, we have to compare two ordered lists of news, which is equivalent to comparing two top k lists for similarity/dissimilarity. We then decided to make use of the averaging Kendall distance analysed in [START_REF] Fagin | Comparing Top K Lists[END_REF], which is a modification of the Kendall's tau metric between permutations (see the textbook [START_REF] Kendall | Rank Correlation Methods[END_REF]) for the case when we only have the top k members of the ordering. Given two ordered lists L 1 and L 2 of distinct numbers taken from a domain D, for each pair (x, y) ∈ D × D, the contribution K x,y (L 1 , L 2 ) of x and y to the averaging Kendall distance of L 1 and L 2 is defined as follows (in the following we denote by a < Li b the fact that a precedes b in L i ).

1. If x and y belong to both L 1 and L 2 , then K x,y (L 1 , L 2 ) = 0 if they are in the same order in both lists, otherwise K x,y (L 1 , L 2 ) = 1. 2. Else if x and y belong to L i , for i ∈ {1, 2}, and either x or y belongs to L 3-i , then

K x,y (L 1 , L 2 ) = 0 if x < Li y and x belongs to L 3-i or if y < Li x and y belongs to L 3-i , otherwise K x,y (L 1 , L 2 ) = 1. 3. Else if x belongs to L i , for i ∈ {1, 2}, and y belongs to L 3-i , then K x,y (L 1 , L 2 ) = 1.

4.

Else if x and y belong to L i , for i ∈ {1, 2}, and neither x nor y belongs to L 3-i , then

K x,y (L 1 , L 2 ) = 1 2 . 5. Else K x,y (L 1 , L 2 ) = 0.
The averaging Kendall distance K(L 1 , L 2 ) of L 1 and L 2 is then equal to the sum of the contributions of all possible pairs (x, y) ∈ D × D. Given two front pages f 1 = (t 1 , L 1 ) and f 2 = (t 2 , L 2 ) we then define their distance as the averaging Kendall distance of the corresponding lists L 1 and L 2 of news. Since we have downloaded articles every hour, we have then computed the distance between the front page taken at a given time and the front page taken one hour later (that is, t 2 is equal to t 1 plus 3600 seconds). In Figure 7 we show the plot of these distances computed for the front pages of one news agency (upper part of the figure) and of one newspaper (lower part of the figure). In both cases, it is quite evident that there is a periodic behaviour in updating the web site: however, it is also clear that while the web site of the news agency is updated quite uniformly, the updating process of the web site of the newspaper is significantly concentrated in a given moment of the day (which is, quite obviously, the morning).

Distances between Different News Sources

By using the averaging Kendall distance, we can also define a distance between two different news sources (note that, in the previous case, we have used the averaging Kendall distance to define a distance between two different front pages of the same news source). Indeed, let t 1 , . . . , t n be the time instants in which we have downloaded the articles and let f S i = (t i , L S i ) be the front page of the news source S taken at time t i . The distance between two news F U N 2 0 1 6 Figure 7 Analysing and comparing the activity of two different news sources.

sources S 1 and S 2 at time t i is then defined as K(L S1 i , L S2 i ), and the distance d(S 1 , S 2 ) the two news sources is the average of all such distances, that is,

d(S 1 , S 2 ) = 1 n n i=1 K(L S1 i , L S2 i ).
Once we have computed these distances, we have then applied the well-known multidimensional technique [START_REF] Cox | Multidimensional Scaling[END_REF] in order to plot the news sources on a two-dimensional plane. The result is shown in Figure 8. As it can be seen from the figure, the two news agencies (in black and grey) are in the centre of the plot and quite close each to the other. One of the four newspapers (in green) is also in the centre of the plot and quite close to both the two news agencies: this suggest that this newspaper acts more as an echo of the news agencies. The other three newspapers, instead, are quite far from the two news agencies and almost equidistant from them (even if at opposite sides): this might imply that these newspapers do indeed elaborate the new produced by the news agencies and present them in the front pages in different ways.

News Score

The front pages can also be used to determine the score of a news. Indeed, it is well known that a web reader usually reads a web page from top to bottom and from left to right. Hence, we can assume that an article, which appears first in the list of a front page, should receive a higher score than an article which appears after. We have decided to assign to each article a penalty proportional to its position in the front page. More formally, if an article a appears in position i of a front page f formed by n articles, then its penalty is equal to i-1 n . In other words, the score of this article is

rel(a, f ) = n -i + 1 n .
We can agglomerate the score of the single articles and obtain the score of a news c with respect to a front page, as follows:

rel(c, f ) = a:N (a)=c
rel(a, f ) . Figure 9 The score of the news concerning the shooting dwon of an American drone.

The total Score of a news is just the sum of its score with respect to all front pages. In Figure 9 we show the score values, during the last four months, of the most relevant news among the almost 7000 news obtained by clustering approximately 10000 articles. This news concerns the shooting down of an American drone in Turkey. As we can see from the figure, the score of this news has two main peaks: one peak corresponds to the news of a shooting down, while the second peak corresponds to the discover that the object shot down was an American drone. It is interesting to analyse this kind of figures because they help us to better understand the evolution of a news. For example, in Figure 10 the evolution of a news concerning a case of homicide in Italy is shown: in this case we have several peaks corresponding, more or less, to the disappearance of the victims, to the rising of the hypothesis of homicide, and to the identification of possible suspects.

Lifespan of a News

Finally, the news detection, performed by means of the clustering algorithm, allows us to estimate the lifespan of a news. In this case, we have first to define what the lifespan is. Indeed, it does not seem to be correct to consider the lifespan as the temporal interval between the first time the news has appeared and the last time in which it is still present in some news source, since during this interval there might have been long sub-intervals in which the news was not present at all. A typical example of this phenomenon is given by a news concerning a homicide. Typically, at the beginning the news is present for few hours, and then disappears. When some evolution in the investigation takes place, the news appears again and very fast disappears. Finally, when the killer is found, the news appears again and definitively disappears. All this can happen during a very long interval time, which cannot really be interpreted as the lifespan of the news, since during the vast majority of this interval the news was not present at all. We then decided to consider as the lifespan of a news as the effective time in which at least one article included in the news was present in at least one source of news. By using this definition, we computed the lifespan of the 7000 news produced by the algorithm with input the 10000 downloaded articles. It turned out that on average the lifespan of a news is approximately twenty hours: the maximum lifespan is a little bit more than 19 days, and it is reached by a news concerning vaccination campaigns inItaly, while the minimum lifespan is one hour, and it is reached, for example, by a news concerning the rescue of two bear cubs in Laos.

Improving the Clustering Time Performance

Despite of the very good performances in terms of precision and recall values, the incremental clustering algorithm described above has the disadvantage of degrading its performances while increasing the number of articles to be clustered. Indeed, any new article has to be compared with any other article already clustered in order to compute its average similarity with all the existing clusters. In order to improve the time performances of the clustering algorithm, we have then designed and implemented a two-layer variation of the previously described incremental algorithm. In this variation, each cluster of the first level is "represented" by its centroid, which is a non-existing article whose keyword set is the union of all the keyword sets of the articles in the cluster, and whose occurrence function is the average of all the occurrence functions. The centroids are themselves clustered by still using the same incremental approach (with maybe a different threshold value). More formally, the centroid of a set of articles A is defined as a "dummy" article m such that ) is smaller than a given threshold τ 2 , then N 2 (x) = n 2 + 1 and N 1 (x) = n 1 + 1 (that is, a new cluster is created both in the first and in the second level). Otherwise, the incremental clustering algorithm described in the previous section is applied to the sub-clustering of N 1 determined by the centroids included in c * 2 . If x is inserted in a cluster c * 1 , then a new centroid of c * 1 is computed and substituted inside the cluster in C 2 containing the previous centroid. Otherwise, a new cluster is created both in the first and in the second level (that is, N 2 (x) = n 2 + 1 and N 1 (x) = n 1 + 1).

key(m) = a∈A k(a)

Threshold estimation

In order to apply the above two-layer incremental clustering algorithm, we have to determine the two thresholds τ 1 and τ 2 . To this aim, we have used again the manual clustering T of the first 3000 downloaded articles, and we have chosen the two values of τ 1 and τ 2 which produced the higher average F -measure. It turned out that the best value of τ 1 and τ 2 is equal to 0.19 and 0.48, respectively: with these values of τ 1 and τ 2 the average F -measure is equal to 0.8. As expected, the obtained average F -measure is lower than the one obtained in the one-layer clustering algorithm. However, the time performance improvement is quite impressive: the clustering of the 3000 articles was executed in less than 2 seconds, while the original approach required more than 2 minutes. Since the number of articles which are downloaded is increasing quite rapidly (approximately 200 articles every day), we believe that a little loss in precision is a reasonable price to be paid in order to make the clustering really efficient.

Conclusions

In this paper we have proposed a new tool for analysing and comparing on-line news sources.

The main ingredients of this tool are a web scraping module, a news detection module based on a (two-layer) incremental clustering algorithm, and a statistical analysis module, which allowed us to answer several questions concerning the Italian on-line news system. Although the study is at the beginning and it is limited to the Italian news system, the results are quite encouraging: some of them are surprising (for instance, the average lifespan of a news), others can be seen just as curiosities, as for instance the news with minimal lifespan (that is, the rescue of two bear cubs in Laos) or which are the most relevant news according to their position in the site, compared with the truly most relevant ones. Another curiosity is to know that some of the online newspapers add very little information to their news, that result almost identical to that of the news agencies. Apart from integrating our tool with other statistical analysis, we think that the most interesting possible future research directions concerns the possibility of improving our incremental clustering algorithm. Indeed, this can be done either by considering a k-level clustering algorithm with k > 2, or by using the estimated average lifespan of the news F U N 2 0 1 6

9:14

Analyzing and Comparing On-Line News Sources in order to eliminate from the clustering execution all clusters (that is, news) which are sufficiently old.

FFigure 1 Figure 2

 12 Figure 1 The structure of the main page of a news source (left) and of an article (right).

Figure 3 Figure 4

 34 Figure 3The minimum, maximum, and average number of daily downloaded articles from mid October 2015 until mid January 2016.

FFigure 5 Figure 6

 56 Figure[START_REF] Edunov | Three and a half degrees of separation[END_REF] The number of occurrences of the name of the leaders of the Italian parties in the title or in the abstract of an article.

FFigure 8

 8 Figure 8 Plotting the news sources on the plane (the two news agencies are in black and in grey).

FFigure 10

 10 Figure 10The score of the news concerning a case of homicide in Italy.

  A two-layer clustering is defined as a pair (N 1 , N 2 ) of two functions N 1 : A → [n 1 ], where n 1 is the number of clusters in the first level, andN 2 : C 1 → [n 2 ], where C 1 is the set of centroids of the n 1 clusters at the first level and n 2 is the cardinality of the set C 2 of clusters in the second level. Given a two-layer clustering (N 1 , N 2 ) and given a new article x ∈ A, let c * 2 be the cluster c 2 ∈ C 2 for which sim(x, c 2 ) is maximum: if sim(x, c * 2

	and, for each w ∈ key(m),
	tf(w, m) =	1 |A|	a∈A	tf(w, a).

Acknowledgements. We would like to thank Alice Marchetti and Andrea Marino to suggest us to include some descriptive statistics of the collection of downloaded articles. The second and third authors received additional support from the Italian Ministry of Education, University and Research, under PRIN 2012C4E3KT National research project AMANDA.


Nicolas Perrin
email: perrin@isir.upmc.fr

Philipp Schlehuber-Caissier

Fast diffeomorphic matching to learn stable nonlinear dynamical systems

Keywords: Nonlinear dynamical systems, Diffeomorphic mapping, Imitation learning, Lyapunov stability, Dynamical movement primitives

We propose a new diffeomorphic matching algorithm and use it to learn nonlinear dynamical systems with the guarantee that the learned systems have global asymptotic stability. For a given set of demonstration trajectories, we compute a diffeomorphism that maps forward orbits of a reference stable time-invariant system onto the demonstrations, thereby deforming the whole reference system into one that reproduces the demonstrations, and is still stable.

Introduction

We consider the problem of learning dynamical systems (DS) from demonstrations. More precisely, given a list of trajectories (x i (t)) observed as timed sequences of points in R d , the objective is to build a smooth autonomous system ẋ = f (x) (i.e. a vector field) that reproduces the demonstrations as closely as possible.

The ability to construct such DS is an important skill in imitation learning (see for example [START_REF] Schaal | Computational approaches to motor learning by imitation[END_REF]). The learned systems can be used as dynamical movement primitives generating goaldirected behaviors [START_REF] Ijspeert | Dynamical movement primitives: learning attractor models for motor behaviors[END_REF].

Modeling motion primitives with DS is convenient for closed loop implementations, and their generalization to unseen parts of the state space provides robustness to spatial perturbations. Moreover, the choice of autonomous (i.e. time-invariant) systems, while not always suitable or preferable, is interesting in many situations as they are inherently robust to temporal perturbations.

The most common motion primitives consist of point-topoint motions, i.e. movements in space that stop at a given target. They correspond to globally asymptotically stable DS. But classical learning algorithms cannot provide the guarantee that their output is always stable. They might produce DS with instabilities or spurious attractors. This issue has recently been studied by Khansari-Zadeh and Billard [START_REF] Khansari-Zadeh | Learning stable nonlinear dynamical systems with gaussian mixture models[END_REF][START_REF] Khansari-Zadeh | Learning control lyapunov function to ensure stability of dynamical system-based robot reaching motions[END_REF] who proposed several approaches to learn stable nonlinear DS. One of the main ideas they investigated consists in learning a Lyapunov function candidate 1 L that is highly compatible with the demonstrations in the following sense: at almost every point x i (t j ), the estimated or measured velocity v i (t j ) is such that its scalar product with the gradient of L is negative: v i (t j )•∇L(x i (t j )) < 0.

Once L is found, a learning algorithm optimizes a weighted sum of DS that admit L as a common Lyapunov function, therefore ensuring the global stability of the resulting DS. Alternatively, L can be used to modify motion primitives by correcting trajectories whenever they would violate the compatibility condition.

The main limitation of this method comes from the difficulty to find good Lyapunov candidates. In SEDS (Stable Estimator of Dynamical Systems), one of the first algorithms proposed by Khansari-Zadeh and Billard, the Lyapunov function is set to be the l 2 -norm squared ( • 2 ), which means that all trajectories produced by the learned DS can only monotonically decrease in norm. In their more recent algorithm CLF-DM (Control Lyapunov Function-based Dynamic Movements), the search for a Lyapunov function candidate is done among a set called Weighted Sums of Asymmetric Quadratic Functions (WSAQF). It highly increases the set of DS that can be learned, but the restrictions remain significant (for instance, no WSAQF can take the same value on two distinct points x and λx with λ > 0). Neumann et al. [START_REF] Neumann | Neural learning of stable dynamical systems based on data-driven lyapunov candidates[END_REF] compute a "Neurally Imprinted Lyapunov Candidate" (NILC) via quadratic programming, obtaining similar restrictions to that of WSAQF. In short, a way to efficiently find Lyapunov candidates is to fix a simple stable DS with which they must be compatible (ẋ = -x in the case of WSAQF and NILC). An issue with this approach is that it necessarily restricts the search to a small convex subset of the whole set of Lyapunov candidates.

To go further, Neumann and Steil [START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF] suggested to initially compute a Lyapunov candidate with one of the above methods, and then apply a simple diffeomorphism (of the form x → η(x)x, with η(x) ∈ R ≥0 ) that deforms the space and transforms the Lyapunov candidate into the function x → x 2 , thus simplifying the trajectories of the demonstrations. In the deformed space, an algorithm like SEDS is then more likely to learn a stable DS that reproduces faithfully the demonstrations.

In this paper, we propose a more direct diffeomorphismbased approach. Our contribution is twofold.

• First, we introduce a new algorithm for diffeomorphic matching (Sections 2 and 3) and show from experimental comparisons that it tends to be one or two orders of magnitude faster than one of the state-of-the-art algorithms.

• Then, we explain how it can be used to directly map simple trajectories of a DS like ẋ = -x onto the trajectories of the training data (Section 4). This gives a new way to generate Lyapunov candidates as well as globally asymptotically stable smooth autonomous systems reproducing the demonstrations.

The most direct applications of this work are in motor control and robotics, but we believe that learning stable nonlinear DS and computing Lyapunov candidates can be useful for various types of systems and control design problems.

Diffeomorphic locally weighted translations

Given a smooth (symmetric positive definite) kernel function k ρ (x, y), depending on some parameter ρ, such that ∀x, k ρ (x, x) = 1 and k ρ (x, y) → 0 when yx → ∞, given a "direction" v ∈ R d and a "center" c ∈ R d , we consider the following locally weighted translation:

ψ ρ,c,v (x) = x + k ρ (x, c)v. Theorem 1. If ∀(x, y) ∈ R d × R d , ∂k ρ ∂x (x, y) • v > -1, then ψ ρ,c,v is a diffeomorphism.
Proof. For x ∈ R d , let us define:

h x : r ∈ R → r + k ρ (x + rv, c) ∈ R. If ∂k ρ ∂x (x, y) • v > -1
, we get: ∀r ∈ R, dh x dr (r) > 0. Since h x (r) tends to -∞ when r tends to -∞, and to +∞ when r tends to +∞, we deduce that there exists a unique s(x) ∈ R such that h x (s(x)) = 0. Moreover, ψ ρ,c,v (x

+ rv) = x + h x (r)v, thus ψ ρ,c,v (x + rv) = x if and only if r = s(x) or v = 0. Besides, any y ∈ R d verifying ψ ρ,c,v (y) = x is necessarily of the form x + rv, so ψ ρ,c,v (y) = x implies y = x + s(x)v.
We conclude that ψ ρ,c,v is invertible, and:

ψ -1 ρ,c,v (x) = x + s(x)v.
Finally, the inverse function theorem can be applied to prove that ψ ρ,c,v is a diffeomorphism.

With Gaussian Radial Basis Function (RBF) kernel:

We now consider the following kernel (with ρ ∈ R >0 ):

k ρ (x, y) = exp -ρ 2 x -y 2 .
We have:

∂k ρ ∂x (x, y) • v = -2ρ 2 exp -ρ 2 x -y 2 (x -y) • v,
with the lower bound:

∂k ρ ∂x (x, y) • v ≥ -2ρ 2 exp -ρ 2 x -y 2 x -y . v .
The expression on the right takes its minimum for x-y = 1 √ 2ρ , which yields:

∂k ρ ∂x (x, y) • v ≥ - √ 2 v ρ exp - 1 2 . Applying Theorem 1, v = 0 or ρ < ρ max (v) = 1 √ 2 v exp 1 2 implies that ψ ρ,c,v is a diffeomorphism.
In that case, s(x), and as a result ψ -1 ρ,c,v (x), can be very efficiently computed with Newton's method.

A diffeomorphic matching algorithm

In this section we are interested in the following problem: given two sequences of distinct points X = (x i ) i∈{0,...,N} and Y = (y i ) i∈{0,...,N} , compute a diffeomorphism Φ that maps each x i onto y i , either exactly or approximately. More formally, defin-

ing dist(A, B) = 1 N+1 i a i -b i
2 for two sequences A and B

of N + 1 points, and denoting by Φ(X) the sequence of points (Φ(x i )) i∈{0,...,N} , we want to find a diffeomorphism Φ that minimizes dist(Φ(X), Y).

Since the sequences X and Y can be very different in shape, to the best of our knowledge the state-of-the-art existing techniques to solve this problem are based on the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework introduced in the seminal article by Joshi and Miller [START_REF] Joshi | Landmark matching via large deformation diffeomorphisms[END_REF]. Its core idea is to work with a time dependent vector field v(x, t) ∈ R d (t ∈ [0, 1]), and define a flow φ(x, t) via the transport equation:

dφ(x, t) dt = v(φ(x, t), t),
with φ(x, 0) = x. With a few regularity conditions on v (see [START_REF] Dupuis | Variational problems on flows of diffeomorphisms for image matching[END_REF] for specific requirements), x → φ(x, t) is a diffeomorphism. The resulting diffeomorphism Φ(x) = φ(x, 1) is given by:

Φ(x) = x + 1 0 v(φ(x, t), t)dt.
Using an appropriate Hilbert space for the vector fields x → v(x, t), they can be associated with an infinitesimal cost whose integration is interpreted as a deformation energy. Various gradient descent algorithms have been proposed to optimize v with respect to a cost that depends both on the deformation energy and on the accuracy of the mapping, whether the objective is to map curves [START_REF] Glaunès | Large deformation diffeomorphic metric curve mapping[END_REF], surfaces [START_REF] Vaillant | Surface matching via currents[END_REF], or, as in our case, points [START_REF] Guo | Diffeomorphic point matching[END_REF].

The LDDMM framework has several advantages. For example, it tries to minimize the deformation, and allows the computation of similarity measures between diffeomorphic geometrical objects. However, Φ is not in closed-form, so once obtained, evaluating it requires an integration that can be slightly timecostly. In our context, Φ can be used inside a control law, so its evaluation (and that of Φ -1 ) needs to be very fast.

We propose a completely different approach to diffeomorphic matching, based on the diffeomorphic locally weighted translations presented in the previous section, which are functions that can be evaluated extremely quickly. The proposed algorithm is based on simple principles, yet it does not seem to have been studied in the past.

The algorithm

We fix a number of iterations K, and two parameters 0 < µ < 1 and 0 < β ≤ 1. Typically, on the examples presented in this paper, we use K = 150, µ ≈ 0.9 and β ≈ 0.5.

Initially, we define Z := X. Every iteration j updates Z, and can be briefly described by the three following steps:

1. we select the point p j in Z that is the furthest from its corresponding target q in Y; 2. we consider the locally weighted translation ψ ρ j ,p j ,v j of direction v j = β(qp j ), center p j , and Gaussian RBF kernel k ρ j , optimizing ρ j ∈ [0, µρ max (v j )] to minimize the error between ψ ρ j ,p j ,v j (Z) and Y; 3. we perform the update: Z := ψ ρ j ,p j ,v j (Z).

The resulting diffeomorphism is:

Φ = ψ ρ K ,p K ,v K • • • • • ψ ρ 2 ,p 2 ,v 2 • ψ ρ 1 ,p 1 ,v 1
Here is a description of the algorithm in pseudo-code:

1: Input: X = (x i ) i∈{0,...,N} and Y = (y i ) i∈{0,...,N} 2: Parameters: Z := ψ ρ j ,p j ,v j (Z) 13: end for 14: return (ρ j ) j∈{1,...,K} , (p j ) j∈{1,...,K} , (v j ) j∈{1,...,K} Remarks:

K ∈ N >0 , 0 < µ < 1, 0 < β ≤ 1 3: 4: Z = (z i ) i∈{0,...,N} 5: Z := X 6: for j = 1 to K do
• Here we have presented a version of the algorithm in which the parameter β is constant, but we can also make it vary iteration after iteration, for example by slowly increasing towards 1.

• The line 11 of the algorithm performs a nonlinear optimization, but since it depends only on one bounded real variable, a minimum can be found very quickly and precisely.

• We can add a fixed upper bound ρ M > 0 for all ρ j , and a regularization term in the cost of the optimization problem of line 11, to prevent the diffeomorphism from overly deforming the space to get a perfect matching. Simply using sequences with a large number of points has a similar effect (and it almost does not slow the algorithm down).

• Again in line 11, dist can be replaced by any other distance. In practice, writing X and Y as (N+1)-by-2 matrices and taking the largest singular value norm of (X-Y) seems to give slightly better results than with the Euclidean distance. • Nothing prevents the algorithm from getting stuck in a local minimum, so a general proof of convergence cannot be found. However, as we show in the next sections, experimental results give empirical evidence that the algorithm is efficient and converges quickly in practice, even on difficult matching problems. In future work, we will try to further improve the algorithm and get convergence proofs under realistic assumptions.

Experimental evaluation

We compare our algorithm to an implementation of diffeomorphic matching in the LDDMM framework developed by J. Glaunès (the "Matchine" software [START_REF]Matchine" software by J. A. Glaunès[END_REF]). Given a sequence of points Y = (y i ) i∈{0,...,N} representing a trajectory, we set X = (x i ) i∈{0,...,N} = y 0 + i N (y Ny 0 ) i∈{0,...,N} and apply our algorithm or the LDDMM algorithm to construct a diffeomorphism Φ such that Φ(X) and Y match. Figure 1 displays the result of our algorithm on four 2D trajectories, and Table 1 shows a comparison of the results obtained on these trajectories with our algorithm and the LDDMM algorithm. For each trajectory, we try with representations as sequences of 21, 51 and 101 points (i.e. N = 20, N = 50, N = 100). For both algorithms, the same parameters are kept across all the trials. In all cases, our algorithm provides a substantial speedup. For example, with 51 points, Φ is learned in average 58 times faster and evaluated 240 times faster, while the error dist(Φ(X), Y) is 2.67 times smaller. The tests were made on an Intel(R) Core(TM) i7-4700MQ @ 2.4 GHz with 4GB of RAM.

N

Learning stable nonlinear dynamical systems

Overview of the method

Two autonomous systems ẋ = f (x) and ẋ = g(x) are said diffeomorphic, or smoothly equivalent, if there exists a diffeomorphism Φ : R d → R d such that:

g(Φ(x)) = J Φ (x) f (x),
where J Φ (x) is the Jacobian matrix: J Φ (x) = ∂Φ ∂x (x). When this equation is verified, Φ maps the orbits of the DS ẋ = f (x) onto the orbits of the DS ẋ = g(x). In particular, if one of the two diffeomorphic DS is globally asymptotically stable, then both are. And if L is a Lyapunov function for the DS ẋ = f (x) , then L • Φ -1 is a Lyapunov function for the DS ẋ = g(x) .

The objective of our approach is to learn a diffeomorphism Φ that maps forward orbits of the DS ẋ = -x (i.e. line segments) onto the observed trajectories. As a consequence, x → Φ -1 (x) will be a Lyapunov candidate. Additionally, it will be possible to directly get a stable DS of the form ẋ = -γ(Φ -1 (x))J Φ (Φ -1 (x))Φ -1 (x), with γ : R d → R >0 , that reproduces approximately the demonstrations. Remark: this The diffeomorphism Φ, that maps the trajectory X onto Y, transforms a stable DS with X as a forward orbit into a stable DS with Y as a forward orbit (top row). Using this property, we can modulate the initial vector field while keeping X unchanged to obtain systems with different behaviors that all reproduce the demonstration Y.

system can be written x = Φ(z), with ż = -γ(z)z (the DS ż = -γ(z)z and ż = -z have the same orbits).

Trajectories being represented as sequences of points, this problem of orbits mapping can be cast as diffeomorphic matching. In the case of a unique demonstration Y = (y(t i )) i∈{0,...,N} , with t i = i∆t, we want to find a diffeomorphism that maps X = y(0) + i N 0y(0) i∈{0,...,N} onto Y (the trajectory is assumed to arrive at the target: y(t N ) = 0). To do so, we simply use the algorithm presented in the previous section. The diffeomorphism Φ K obtained after K iterations can be such that Φ K (0) 0, so we add an extra iteration that picks p K+1 = Φ K (0) and v K+1 = 0 -Φ K (0). This ensures that the final diffeomorphism Φ verifies Φ K (0) = 0. Remark: the structure of Φ makes it easy to efficiently compute J Φ (x) at any given point.

Results

The top row of Figure 2 shows the result of mapping the straight trajectory X (on the left) onto the trajectory Y (on the right). The diffeomorphism Φ that realizes this matching also transforms the entire DS ẋ = -x into a nonlinear stable DS that reproduces the trajectory Y (as the forward orbit of y(0)).

Modifying the initial DS without changing the forward orbit of x(0) leads, by application of Φ, to another DS that still reproduces Y. On the bottom row of Figure 2, we use a linear system that keeps x(0) as an eigenvector associated with eigenvalue -1 (ensuring that its orbit is not modified), but has an eigenvalue smaller than -1 in the orthogonal direction (unlike the DS ẋ = -x). This results in a transformed DS that "tracks" more agressively the trajectory Y. We evaluated our approach on the LASA Handwriting Dataset, similarly to [START_REF] Khansari-Zadeh | Learning stable nonlinear dynamical systems with gaussian mixture models[END_REF][START_REF] Khansari-Zadeh | Learning control lyapunov function to ensure stability of dynamical system-based robot reaching motions[END_REF][START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF]. On all cases shown in Figure 4, a set of 7 trajectories ending at the same point demonstrate a single pattern of handwriting motion. For each of these patterns, we create an average timed sequence of points Y = (y(i∆t)) i∈{0,...,N} based on the 7 trajectories, and apply our matching algorithm to construct a diffeomorphism Φ that maps X = N-i N y(0) i∈{0,...,N} onto Y. This gives a Lyapunov candidate x → Φ -1 (x) . We compare it to the optimized WSAQF Lyapunov candidates obtained with the method of Khansari-Zadeh and Billard [START_REF] Khansari-Zadeh | Learning control lyapunov function to ensure stability of dynamical system-based robot reaching motions[END_REF] also used in [START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF]. On the 1st column are displayed level sets of the WSAQF Lyapunov candidates, and on the 2nd column level sets of our Lyapunov candidates. We can observe that the level sets of the Lyapunov candidates produced by our method have a richer geometry and exhibit variations that are more suitably adapted to the training data.

Of course, following the method of Khansari-Zadeh and Billard [START_REF] Khansari-Zadeh | Learning control lyapunov function to ensure stability of dynamical system-based robot reaching motions[END_REF], these Lyapunov candidates can be used to correct any learned DS to ensure global asymptotic stability. But as mentioned above, the diffeomorphism also provides a way to directly get a stable DS that reproduces the motion pattern. We define a function γ : R d → R >0 such that, starting at x(0) = y(0) with t = 0, the DS ẋ = -γ(x)x produces a trajectory that passes by the points N-i N y(0) at times i∆t, for i ∈ {1, . . . , N -1}, and converges asymptotically towards 0 for t > (N -1)∆t. A simple choice for γ is γ(x) = y(0) N∆t x for x ≥ y(0) N and γ(x) = y(0)

N otherwise (but it is easy to design a smoother function with the same desired properties). Φ transforms the DS ẋ = -γ(x)x into one that reproduces the demonstrations and their velocity profiles, as shown in Figure 3. The eigenvalue in the direction orthogonal to y(0) can be adjusted according to the variability of the 7 demonstrations. The vector fields we obtained are shown on the 4th column of Figure 4, and the vector fields obtained with the τ-SEDS (WSAQF) method of Neumann and Steil [START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF] are shown on the 3rd column.

Conclusion

In this paper, we presented a new algorithm for diffeomorphic matching and a way to use it to learn stable nonlinear au-tonomous dynamical systems from demonstrations.

While the demonstrations were 2D single motion patterns in the results we presented, our algorithm scales well to higher dimensions (because all its parameters are dimensionindependent) and can handle multiple motion patterns, although in some cases topological issues may prevent the convergence of the matching. It should also be noted that we can only produce vector fields that are diffeomorphic to the DS ẋ = -x, which is not true for all globally asymptotically stable smooth autonomous systems. A related limitation concerns 2D spiral trajectories, which cannot be reproduced. In future work, we will try to combine our approach with existing methods to extend its possibilities. For instance, Kronander et al. [START_REF] Kronander | Incremental motion learning with locally modulated dynamical systems[END_REF] suggest to iteratively reshape DS by locally applying full-rank modulations such as scalings and rotations. It does not guarantee the global asymptotic stability of the resulting DS, but we see in this work an approach that could be complementary to ours.

The main advantages of our method are:

1. speed: unlike [START_REF] Khansari-Zadeh | Learning control lyapunov function to ensure stability of dynamical system-based robot reaching motions[END_REF] and [START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF], we do not need a second learning phase once the Lyapunov candidate has been found, and we do not rely on numerical optimization of parameters whose number rapidly increases with the dimensionality; instead, the simple iterative algorithm we use has a constant number of parameters, and we have empirically verified its quick convergence for many difficult matching problems. 2. simplicity: of implementation because the algorithm is very short, but also of use thanks to the small number of parameters to adjust.

For these reasons, we believe it can be applied with ease to efficiently learn a large variety of stable autonomous systems, with applications in dynamic movement primitives construction or more generally in control design. If instead of ẋ = -x, we start with an initial DS that has a stable limit cycle, our approach can be adapted to learn limit cycle systems.

Finally, being significantly faster than a state-of-the-art algorithm, our diffeomorphic matching algorithm itself might be of interest for completely different applications, such as for example image registration. • 3rd column: streamlines of the stable DS produced by the τ-SEDS (WSAQF) approach [START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF]. • 4th column: streamlines of the stable DS produced by our approach.

• All columns: the demonstrations are displayed in black.

• 3rd & 4th columns: the reproduced trajectories are in red.

  j := β(qp j ) 11:ρ j := arg min ρ∈[0,µρ max (v j )]dist(ψ ρ,p j ,v j (Z), Y)12:

Figure 1 :

 1 Figure 1: On the left, the dashed curve is a trajectory represented by a sequence of points Y = (y i ) i∈{0,...,N} . The solid line is X = y 0 + i N (y Ny 0 ) i∈{0,...,N} . The right side shows the result of the application of the diffeomorphism Φ constructed by our algorithm to map X onto Y.

Figure 2 :

 2 Figure2: The diffeomorphism Φ, that maps the trajectory X onto Y, transforms a stable DS with X as a forward orbit into a stable DS with Y as a forward orbit (top row). Using this property, we can modulate the initial vector field while keeping X unchanged to obtain systems with different behaviors that all reproduce the demonstration Y.

Figure 3 :

 3 Figure3: On the left: a smooth autonomous systems, learned with our method, that reproduces a motion pattern (demonstrations are in black, reproduced trajectories in red). The trajectories on the right show that the velocity profiles are quite accurately reproduced as well (again, demonstrations in black and reproductions in red).

Figure 4 :

 4 Figure 4: • 1st column: level sets of the WSAQF Lyapunov candidates [4, 6]. • 2nd column: level sets of the Lyapunov candidates obtained with our approach.• 3rd column: streamlines of the stable DS produced by the τ-SEDS (WSAQF) approach[START_REF] Neumann | Learning robot motions with stable dynamical systems under diffeomorphic transformations[END_REF]. • 4th column: streamlines of the stable DS produced by our approach.• All columns: the demonstrations are displayed in black.• 3rd & 4th columns: the reproduced trajectories are in red.

Table 1 :

 1 Comparison of experimental results

			our algorithm LDDMM
	Learning: average	20	0.25 s	2.78 s
	duration of the	50	0.25 s	14.5 s
	construction of Φ	100	0.26 s	53.3 s
	Forward evaluation:	20	3.05 ms	157 ms
	average duration of the 50	3.35 ms	804 ms
	computation of Φ(X)	100	3.72 ms	3130 ms
	Backward evaluation:	20	29.8 ms	145 ms
	average duration of the 50	35 ms	798 ms
	computation of Φ -1 (Y) 100	38.5 ms	3110 ms
	Accuracy: average value of dist(Φ(X), Y)	20 50 100	3.49×10 -3 8.32×10 -3 9.51×10 -3	18.2×10 -3 22.2×10 -3 22×10 -3

In this paper, a Lyapunov function candidate (or simply Lyapunov candidate) is a continuously differentiable function from R d to R ≥0 taking the value 0 at the target point 0 and with no other local minimum.

Preprint submitted to Systems & Control Letters December 18, 2015