András Dobó
email: dobo@inf.u-szeged.hu

János Csirik
email: csirik@inf.u-szeged.hu

Comparison of the best parameter settings in the creation and comparison of feature vectors in distributional semantic models across multiple languages

Keywords:

Measuring the semantic similarity and relatedness of words is important for many natural language processing tasks. Although distributional semantic models designed for this task have many different parameters, such as vector similarity measures, weighting schemes and dimensionality reduction techniques, there is no truly comprehensive study simultaneously evaluating these parameters while also analysing the differences in the findings for multiple languages. We would like to address this gap with our systematic study by searching for the best combination of parameter settings in the creation and comparison of feature vectors in distributional semantic models for English, Spanish and Hungarian separately, and then comparing our findings across these languages. During our extensive analysis we test a large number of possible settings for all parameters, with more than a thousand novel variants in case of some of them. As a result of this we were able to find such combinations of parameter settings that significantly outperform conventional settings combinations and achieve state-ofthe-art results.

Introduction

In many NLP problems, including information retrieval [START_REF] Hliaoutakis | Information retrieval by semantic similarity[END_REF], spelling correction [START_REF] Budanitsky | Semantic distance in WordNet: An experimental, applicationoriented evaluation of five measures[END_REF] and noun compound interpretation [START_REF] Dobó | Interpreting noun compounds using paraphrases[END_REF] among many others, knowing the semantic similarity or semantic relatedness of words can be very useful. Although distributional semantic models (DSMs) [START_REF] Baroni | Don't count, predict! A systematic comparison of context-counting vs. context-predicting semantic vectors[END_REF] calculating these measures have many possible parameters, most research focuses only on one or two aspects of these models, while using some conventional settings for the rest of the parameters (e.g. cosine as vector similarity and (positive) pointwise mutual information as weighting). Therefore a truly comprehensive study evaluating the numerous parameters of DSMs for any language is still missing, and would be needed, as also suggested by [START_REF] Levy | Improving Distributional Similarity with Lessons Learned from Word Embeddings[END_REF]. Moreover, despite the fact that the best parameter settings for the parameters can differ for different languages, the vast majority of papers consider DSMs for only one language (mostly English), or consider multiple languages but without a real comparison of findings across languages. In this article we would like to address these gaps.

There are two distinct phases of DSMs in general: the extraction of statistical information from raw text, and the creation and comparison of feature vectors for words based on the extracted information. Within this study we focus on the parameters of the second phase, as the two phases are relatively distinct from each other, and the number of possible combinations of parameter settings (CPS) for the second phase is already well in the trillions. So we are searching for the best CPS of those 10 parameters during the creation and comparison of feature vectors in DSMs that we considered important, for English, Spanish and Hungarian separately, and then we compare our findings for the different languages.

A detailed description of our analysis can be found in [START_REF] Dobó | A comprehensive study of the parameters in the creation and comparison of feature vectors in distributional semantic models[END_REF], where tests were done only for English, without any comparison of findings across languages, and with far less settings tested for several parameters.

Background

Although there are a vast number of studies dealing with DSMs, most of them only consider one or two parameters of these models, and take the others granted with some standard setting. Most commonly, these models use cosine as vector similarity [START_REF] Yih | Measuring Word Relatedness Using Heterogeneous Vector Space Models University of Michigan[END_REF][START_REF] Bruni | Multimodal distributional semantics[END_REF][START_REF] Baroni | Don't count, predict! A systematic comparison of context-counting vs. context-predicting semantic vectors[END_REF][START_REF] Pennington | GloVe : Global Vectors for Word Representation[END_REF][START_REF] Levy | Improving Distributional Similarity with Lessons Learned from Word Embeddings[END_REF][START_REF] Iosif | Cognitively Motivated Distributional Representations of Meaning[END_REF][START_REF] Salle | Matrix Factorization using Window Sampling and Negative Sampling for Improved Word Representations[END_REF][START_REF] De Deyne | Predicting human similarity judgments with distributional models: The value of word associations[END_REF][START_REF] Speer | ConceptNet 5.5: An Open Multilingual Graph of General Knowledge[END_REF][START_REF] Salle | LexVec[END_REF][START_REF] Christopoulou | Mixture of Topic-Based Distributional Semantic and Affective Models[END_REF][START_REF] Vakulenko | Calculation of Semantic Distances Between Words: From Synonymy to Antonymy[END_REF] and (positive) pointwise mutual information as weighting [START_REF] Islam | Semantic text similarity using corpus-based word similarity and string similarity[END_REF][START_REF] Yih | Measuring Word Relatedness Using Heterogeneous Vector Space Models University of Michigan[END_REF][START_REF] Bruni | Multimodal distributional semantics[END_REF][START_REF] Levy | Improving Distributional Similarity with Lessons Learned from Word Embeddings[END_REF][START_REF] Iosif | Cognitively Motivated Distributional Representations of Meaning[END_REF][START_REF] Salle | Matrix Factorization using Window Sampling and Negative Sampling for Improved Word Representations[END_REF][START_REF] Salle | LexVec[END_REF]. Further, they also usually do not care for the interaction of these parameters, and experiment with the considered parameters one by one, and not simultaneously. Of course there are some studies that experiment with several parameters with multiple possible settings [START_REF] Curran | From Distributional to Semantic Similarity p[END_REF][START_REF] Lapesa | A Large Scale Evaluation of Distributional Semantic Models: Parameters , Interactions and Model Selection[END_REF][START_REF] Kiela | A systematic study of semantic vector space model parameters[END_REF], but even these are far from being truly comprehensive.

Moreover, most models were only tested for English and neglect any other languages despite the fact that DSMs might work differently across multiple languages. Of course, there are several studies in which results were presented for languages other than English, including Spanish [START_REF] Hassan | Cross-lingual semantic relatedness using encyclopedic knowledge[END_REF][START_REF] Moldovan | Semantic similarity: normative ratings for 185 Spanish noun triplets[END_REF][START_REF] Camacho-Collados | Semeval-2017 task 2: Multilingual and cross-lingual semantic word similarity[END_REF] and Hungarian [START_REF] Dobó | Magyar és angol szavak szemantikai hasonlóságának automatikus kiszámítása[END_REF][START_REF] Novák | Magyar szóbeágyazási modellek kézi kiértékelése[END_REF]. However, even those that include multiple languages usually only present some test results for the different languages separately, without any real analysis of the differences in the findings between the languages.

Data and evaluation methods

For input we used information extracted from the British National Corpus (BNC), the Spanish Wikicorpus (EsWiki) [START_REF] Reese | Wikicorpus: A Word-Sense Disambiguated Multilingual Wikipedia Corpus[END_REF], and the 23.01.2012 dump of the Hungarian Wikipedia (HuWiki) for English, Spanish and Hungarian, respectively, with the help of the information extraction methods of [START_REF] Dobó | Magyar és angol szavak szemantikai hasonlóságának automatikus kiszámítása[END_REF], [START_REF] Dobó | Computing semantic similarity using large static corpora[END_REF] and [START_REF] Salle | Matrix Factorization using Window Sampling and Negative Sampling for Improved Word Representations[END_REF].

Tests were done on parts of the MEN [START_REF] Bruni | Multimodal distributional semantics[END_REF] dataset for English, the Spanish WordSimilarity-353 [START_REF] Hassan | Cross-lingual semantic relatedness using encyclopedic knowledge[END_REF], the Moldovan [START_REF] Moldovan | Semantic similarity: normative ratings for 185 Spanish noun triplets[END_REF] and the Spanish Rubenstein-Goodenough [START_REF] Camacho-Collados | A Framework for the Construction of Monolingual and Cross-lingual Word Similarity Datasets[END_REF] datasets for Spanish, and parts of the Hungarian version of the TOEFL [START_REF] Dobó | Magyar és angol szavak szemantikai hasonlóságának automatikus kiszámítása[END_REF] and Rubenstein-Goodenough datasets for Hungarian. The last was constructed the same way as the Hungarian TOEFL and Miller-Charles datasets in [START_REF] Dobó | Magyar és angol szavak szemantikai hasonlóságának automatikus kiszámítása[END_REF].

Out of these datasets only the MEN dataset is truly reliable, as the others are rather small and except for the Moldovan dataset just translated from English datasets, during which they can be distorted. The Hungarian datasets are especially small, and the type of the TOEFL dataset also makes the results on it even less reliable compared to the other datasets. However, due to the lack of truly suitable resources, we had to settle for these.

In case of the TOEFL dataset, the accuracy (A) of the models on the questions were calculated, while in case of all the other datasets, the Pearson's (P) and Spearman's (S) correlations with the gold standard scores and their modified harmonic mean (H) were computed, as follows:

H(P, S) = 2 × P × S |P | + |S| (1)
For more information, please refer to [START_REF] Dobó | A comprehensive study of the parameters in the creation and comparison of feature vectors in distributional semantic models[END_REF].

Our heuristic analysis

As the number of possible CPSs are in the magnitude of trillions, we had to use a heuristic approach to find the best one in case of each language. First, each parameter was tested separately on a development dataset, where a candidate list of settings were selected for each parameter. Then all combinations of the selected settings of all parameters were tested on a different dataset, to find the best CPS. These were done for the three languages separately. The 10 parameters tested, together with the number of settings tried for them, are listed in Table 1. For several parameters, a large number of novel settings were tested. These were either brand new settings, modified versions of conventionally used settings, or combinations of multiple settings. A detailed description of the tested parameters and their settings can be found in [START_REF] Dobó | A comprehensive study of the parameters in the creation and comparison of feature vectors in distributional semantic models[END_REF].

We have to note that when using singular value decomposition (SVD) for dimensionality reduction or the various smoothing options, we usually did a smaller number of runs than in other cases due to our limited resources. Further, in case of Spanish we had to set MWFFreq to 3 instead of NoLimit when using SVD due to the too many features otherwise, which would have made running SVD unmanageable.

First we have done this two-step heuristic analysis for English and evaluated the results extensively in [START_REF] Dobó | A comprehensive study of the parameters in the creation and comparison of feature vectors in distributional semantic models[END_REF]. Then we have repeated the same analysis, with a greatly increased number of settings for several parameters, for English, Spanish and Hungarian, and compared the findings across the different languages in this article.

Results

Results of the first phase

During the first phase of our analysis multiple runs were done for each setting of every parameter, and the most promising ones in case of each parameter were selected to be included in the second phase. In case of English, we used half of the development part of the MEN dataset for evaluation, while for Spanish the Spanish WordSimilarity-353 dataset and for Hungarian half of the Hungarian TOEFL dataset was employed. The top 5 performing settings for each parameter are listed in Table 2 in case of each language.

Although presenting the definition of all settings for every parameter would be impossible within this article due to their large number, below we briefly define a couple of them to help interpreting our most important results.

In case of vector similarity measures, we have defined many new variants based on one or more conventional measures. For example, the best one for English is a combination of the Pearson, MarylandBridge [START_REF] Deza | Encyclopedia of distances[END_REF] and AdjCos [START_REF] Shalaby | Measuring Semantic Relatedness using Mined Semantic Analysis[END_REF] measures, with some additional transformations:

P earsM bAdjCosM od-3.Lb(u, v) = 1, d ≥ 0.1 d 0.1 , d < 0.1 d = 0.5 × n i=1 sgn(ui -ū) × lb(|ui -ū| + 1) × sgn(vi -v) × lb(|vi -v| + 1) lbinv n i=1 (lb(|ui -ū| + 1)) 2 + n i=1 sgn(ui -ū) × lb(|ui -ū| + 1) × sgn(vi -v) × lb(|vi -v| + 1) lbinv n i=1 (lb(|vi -v| + 1)) 2 lbinv(x) = min(max(sgn(x) × (2 |x| -1), -2 100 ), 2 100 ) (2)
On the other hand, the best vector similarity measure for Spanish is rather different. It is a modified and transformed version of the Hindle r measure [START_REF] Lin | Automatic retrieval and clustering of similar words[END_REF]:

LinHindleRM od-7.1.2.Cu(u, v) = 3 n i=1 lhr-1.Cu(ui, vi) n i=1 u 2 i × n i=1 v 2 i lhr-1.Cu(x, y) = min(x 3 , y 3 ), x = 0 ∧ y = 0 0, otherwise (3)
Weighting schemes were constructed similarly as vector similarity measures, and here too there were also numerous new variants. For example, the best one for English is a combination of PMIAplha [START_REF] Levy | Improving Distributional Similarity with Lessons Learned from Word Embeddings[END_REF], PMI with Laplace smoothing [START_REF] Turney | From Frequency to Meaning: Vector Space Models of Semantics[END_REF], Unisubtuples [START_REF] Pecina | Lexical association measures and collocation extraction[END_REF] and PMI with discounting factor [START_REF] Pantel | Discovering word senses from text[END_REF]: In case of feature transformation, we have experimented with transforming either raw frequencies or weights, both before and after normalization, and 7 different transformation functions were tried in all cases. For smoothing, we tried different versions of the Kneser-Ney smoothing [START_REF] Dobó | Multi-D Kneser-Ney Smoothing Preserving the Original Marginal Distributions[END_REF]. In case of dimensionality reduction, we tried a couple of different techniques, including SVD and the method of [START_REF] Islam | Semantic text similarity using corpus-based word similarity and string similarity[END_REF]. And for minimum limits on word-feature weights we have tried the following two novel variants with multiple limit values: limit(w, minV alue) = w if w ≥ minV alue minV alue otherwise [START_REF] Camacho-Collados | A Framework for the Construction of Monolingual and Cross-lingual Word Similarity Datasets[END_REF] zero(w, minV alue) = w if w ≥ minV alue

P miAl-T c3T w0S2P 4(x, y) = f xy f xy + 1 × min(f x , f y ) min(f x , f y ) + 1 × lb n α × f xy f x × f 0.75 y -3.29 × 1 a + 1 b + 1 c + 1 d a = f xy , b = f x -f xy , c = f y -f xy , d = n -f x -f y + f xy f x = fx + 1, f y = fy + 1, f xy = fxy + 1, n = n + 1, n α =   |V | i=1 f 0.75 i   +
0 otherwise (6)
The other parameters are much less complex and more commonly used in NLP, thus one should be able to have enough understanding of them from Table 1.

A more detailed description of the different parameters and settings can be found in [START_REF] Dobó | A comprehensive study of the parameters in the creation and comparison of feature vectors in distributional semantic models[END_REF] (although with far less settings for several parameters).

Results of the second phase

In the second phase all possible combinations of the settings of each parameter were tested in case of all three languages, in order to find the best CPS for all languages. The second half of the development part of the MEN dataset was used for testing in case of English, while the Moldovan dataset and the second part of the Hungarian TOEFL dataset were used for Spanish and Hungarian, respectively. The top 5 performing CPSs for each language are presented in Table 3.

Results on the test datasets

The best CPS for English was tested on the test part of the MEN dataset (MT), and the best CPSs for Spanish and Hungarian were tested on the respective version of the Rubenstein-Goodenough dataset (RG) to give us the final results. The best CPS of each language was also evaluated on the datasets of the other languages, to provide us a way of comparison. The results of these test can be found in Table 4.

Evaluation and discussion

In this section we evaluate our results presented in the previous sections. Please note that the scores are not fully comparable across languages, even when considering the same datasets on different languages, as except for the Moldovan dataset all of the used Spanish and Hungarian datasets were constructed by translating the English versions, and thus the results on them can be distorted and less reliable than on their English counterparts. Furthermore, the Spanish and Hungarian datasets, especially the latter ones, are rather small, which also makes them less reliable than the English ones.

As there are many differences in the syntax and morphology of the different languages, we anticipated from the beginning that there will be at least some small differences in our findings for the different languages. However, our intuition was that our findings for the different languages will be subtle, and we will be able to find good and rather language-independent CPSs. As English and Spanish belong to the family of Indo-European languages, while Hungarian does not, we expected that the results for English and Spanish will be similar due to this. Further, as both Spanish and Hungarian have very rich morphology, we expected that there will also be a higher similarity between our results for Spanish and Hungarian because of this. We anticipated that the least similarities will be between English and Hungarian, as these languages are the least similar to each other.

In the first phase of our analysis we could observe that some of the parameters worked exactly the same way or very similarly across languages. These parameters were the weighting scheme, feature transformation, vector normalization and minimum limits on word-feature frequencies. These finding are in line with our initial intuitions. Dimensionality reduction seemed to be similar for English and Spanish, while a bit different for Hungarian. Smoothing seemed to perform similarly for Spanish and Hungarian, while differently for English. Minimum limits on word-feature weights seem to behave a bit differently for all three languages. However, it was interesting to see that the results for vector similarity measures, stop-word filtering and minimum limits on feature frequencies were rather similar for English and Hungarian, but different for Spanish, which is contrary to what we anticipated.

In the second phase, although there were similarities in the found best CPSs across the different languages, one could also observe many differences. Here too, the weighting schemes, feature transformation and minimum limits on word-feature frequencies were mostly similar. Compared to the first phase vector similarity, smoothing and min-imum limits on feature frequencies were also alike for all languages. The other parameters showed a different behaviour for at least one language compared to the others.

We have to note here that there were actually two distinct CPSs with the same best score for English, and they were only different in their DimRed parameter setting. We have chosen the one with the "IslamInkpen 0.05" setting as BestEn, as that setting achieved better performance in the first phase than the "NoDimRed" setting in the other CPS. Furthermore, for Hungarian there were even more CPSs with the same best score. We have used a similar approach in selecting the BestHu version, as we have done in case of the BestEn version. However, as these different CPSs with the same best results have different settings in case of some parameters, one has to be careful drawing conclusions from the best CPSs of the different languages, and thus any conclusions drawn from them should be taken with some reservations.

The final conclusions for the parameters are the following:

-VecSim: for all languages measures based on cosine similarity achieve the best results -Weight: measures based on PMI dominate the top of the table by far in case of all languages -FeatTransf: no transformation and transforming the word-feature weights after normalization preforms best for all languages -DimRed: dimensionality reduction seems to help in most situations: while in case of English the IslamInkpen version performed the best alongside no dimensionality reduction, for Spanish and Hungarian SVD is superior to these options -Smooth: the no smoothing option clearly outperforms all others for all languages -VNorm: for English the L 1 option clearly seems to be the best, while for Spanish and Hungarian the best CPSs use either L 1 or L 2 normalization, and most CPSs achieve the same or very similar results with either -StopW: stop-word filtering seems to improve the results to some extent in case of Spanish, while it does not in case of English and Hungarian -MWFFreq: no limit is by far superior to the other options for all languages -MWFWeight: no limit seems to be the best option in case of Spanish and Hungarian, while the Zero option with different parameters seems to excel in case of English -MFFreq: a low limit or no limit seems to be best in case of all languages (as noted before, in case of SVD for Spanish we had to use a limit of 3 instead of no limit for computational reasons)

As we anticipated, there were parameters where the results for Spanish and Hungarian were similar, but different for English. However, it was interesting that we did not find any parameters that were alike for English and Spanish, but different for Hungarian. Further, to our surprise we found such a parameter, where the results were similar for English and Hungarian, but different for Spanish. These latter findings were in contrast to our initial intuition.

Although all Spanish scores in the second phase are much lower than the English and Hungarian ones, these are almost completely due to the dataset used, and do not mean that the found Spanish CPSs are worse than their English and Hungarian counterparts, as it was noted in the beginning of this section and can be seen from our results on Table 5. Comparison of our best English CPS with a conventional and a state-of-the-art CPS, using the information extracted from the BNC with the method of [START_REF] Salle | Matrix Factorization using Window Sampling and Negative Sampling for Improved Word Representations[END_REF] It was interesting to see that in the cross-language experiments on the test datasets the order of the Best CPSs of the different languages with respect to their performance is different in case of the datasets of the three languages. The best English CPS was always superior to its Spanish counterpart, but it has no absolute superiority over the best Hungarian CPS. Further, there is also no clear ranking between the best Spanish and Hungarian CPSs. It was also interesting to see that in case of the Spanish dataset, although the best Spanish CPS achieved rather good results, actually it achieved the lowest score out of the three best CPSs tested. It was the same for the best Hungarian CPS on the Hungarian dataset too.

All in all, there seems to be no clear ranking between the best CPSs of the different languages, and all of them achieved good results on the datasets of all languages. So, although we got different best CPSs for the different languages, all of them seem to be rather language-independent. These findings give us a strong intuition that our heuristic approach was good, and that our found best CPSs for all languages and their results are robust and reliable.

To further prove that our heuristic approach was successful, that our results are robust and reliable, and that our found best CPSs perform much better than conventional CPSs, we compared our best English CPS with the conventional cosine with positive pointwise mutual information setting (CosPPmi) and the state-of-the-art original settings combination (OSC) of [START_REF] Salle | Matrix Factorization using Window Sampling and Negative Sampling for Improved Word Representations[END_REF], using the information extracted from the BNC with the method of [START_REF] Salle | Matrix Factorization using Window Sampling and Negative Sampling for Improved Word Representations[END_REF] and the MT dataset for all tests. The results of these tests, presented in Table 5, clearly show that our found best CPS is robust, and is not just superior to conventional settings, but to a current state-of-the-art CPS too.

The fact that the best CPSs found in the second phase are not simply made up of the best parameter settings in the first phase proves that our intuition was correct, and the parameters of DSMs need to be tested simultaneously, rather than separately.

Conclusions

Within this article we have presented a systematic analysis of the parameters in the creation and comparison of feature vectors in distributional semantic models for English, Spanish and Hungarian, including novel parameters and novel parameter settings. To our best knowledge, we are the first to do such a detailed analysis for these parameters, and also to do such an extensive comparison of them across multiple languages.

With our heuristic approach we were searching for the best combination of parameter settings for all three languages. In accordance with our intuition, there were several parameters that worked very similarly in case of all three languages. We also found such parameters that were alike for Spanish and Hungarian, and different for English, which we also anticipated. However, it was interesting to see that there was such a parameter that worked similarly for English and Hungarian, but not for Spanish, and we did not find any parameters that worked similarly for the two Indo-European languages, but differently for Hungarian.

Although we have found that the very best results are produced by different settings combinations for the different languages, our cross-language tests showed that all of them work rather well for all languages. Based on this we think that our heuristic approach was successful, and we could find such combinations of parameter settings that are rather language-independent, and give robust and reliable results. Further, our best English CPS, incorporating multiple novel parameter settings, significantly outperformed both conventional and state-of-the-art parameter combinations.

Although our results seem rather robust and reliable for Spanish and Hungarian too, it would be interesting to redo our analysis on larger and more reliable Spanish and Hungarian datasets to check whether we could find even better CPSs for these languages, when such datasets will become available in the future.

Table 1 .

 1 The tested parameters, with the number of settings tested for each

	Parameter	Abbreviation Count
	Vector similarity	VecSim	1221
	Weighting scheme	Weight	2907
	Feature transformation	FeatTranf	22
	Dimensionality reduction DimRed	21
	Smoothing	Smooth	5
	Vector normalization	VNorm	3
	Stop-word filtering	StopW	2
	Minimum limits on word-feature frequencies	MWFFreq	6
	Minimum limits on word-feature weights	MWFWeight	26
	Minimum limits on feature frequencies	MFFreq	14

Table 2 .

 2 The top 5 performing setting for each parameter in case of all 3 languages, in descending order of maximum H scores

	Parameter	English		Spanish		Hungarian
		Setting	H	Setting	H	Setting	H
		PearsMbAdjCosMod-3.Lb 0.71	LinHindleRMod-7.1.2.Cu	0.37	PearsMbMod-1.Lb	0.80
		PearsMbAdjCosMod-4.Lb 0.71	LinHindleRMod-6.1.2.Cu	0.36	PearsMbMod-4.Lb	0.80
	VecSim	PearsMbAdjCosMod-2.Lb 0.71	LinHindleRMod-1.1.2.Cu	0.36	MbMod-6.Lb	0.80
		PearsMbAdjCosMod-6.Lb 0.71	LinHindleRMod-7.1.2.Sq	0.36	PearsMbMod-5.Lb	0.80
		PearsMbAdjCosMod-6.Sigm 0.71	LinHindleRMod-3.1.2.Sq	0.36	PearsMbMod-2.Lb	0.80
		PmiAl-Tc3Tw0S2P4	0.70	PmiAlUnis-Tc4Tw3S2P2	0.38	PmiAl-Tc4Tw2S1P1	0.85
		PmiAl-Tc3Tw0S2P0	0.70	PmiAlUnis-Tc4Tw3S2P1	0.38 PmiAlUnisAm-Tc0Tw3S2P1 0.85
	Weight	PmiAlUnis-Tc3Tw0S0P4	0.70	PmiAlUnis-Tc4Tw2S1P1	0.38 PmiAlUnisAm-Tc0Tw2S2P2 0.85
		PmiAlUnis-Tc3Tw0S0P0	0.70 PmiAlUnisAm-Tc4Tw3S2P5 0.38	PmiAl-Tc4Tw3S0P2	0.85
		PmiAl-Tc4Tw0S2P5	0.70	PmiAlUnis-Tc4Tw2S1P2	0.38 PmiAlUnisAm-Tc0Tw2S2P1 0.85
		Weight AftNorm Lb	0.67	Weight AftNorm Lb	0.34	Weight AftNorm Sqrt	0.85
		Freq Sq	0.67	Weight AftNorm Sigm	0.34	Weight BefNorm Sqrt	0.85
	FeatTransf	Weight BefNorm Sigm	0.67	NoTransf	0.34	Weight AftNorm Sigm	0.80
		Weight AftNorm Sigm	0.67	Weight BefNorm Lb	0.34	NoTransf	0.80
		NoTransf	0.67	Weight BefNorm Sigm	0.34	Weight AftNorm Lb	0.80
		SVD 200	0.70	SVD 100	0.37	IslamInkpen 0.025	0.80
		SVD 100	0.70	SVD 200	0.36	IslamInkpen 0.25	0.80
	DimRed	SVD 300	0.69	SVD 500	0.35	SVD 200	0.80
		SVD 500	0.68	SVD 300	0.34	IslamInkpen 0.005	0.78
		IslamInkpen 0.05	0.67	IslamInkpen 0.01	0.34	IslamInkpen 0.01	0.78
		NoSmooth	0.67	Freq KNS	0.34	Freq KNS	0.83
		Weight KNS	0.65	Freq MDKNSPOMD	0.34	Freq MDKNSPOMD	0.80
	Smooth	Freq KNS	0.62	NoSmooth	0.34	NoSmooth	0.78
		Freq MDKNSPOMD	0.62	Freq MKNS	0.33	Weight KNS	0.75
		Freq MKNS	0.57	Weight KNS	0.31	Freq MKNS	0.73
		L 2	0.67	L 2	0.34	L 2	0.80
	VNorm	L 1	0.67	L 1	0.34	NN	0.78
		NN	0.67	NN	0.34	L 1	0.75
	StopW	false true	0.67 0.67	true false	0.34 0.34	false true	0.80 0.75
		NoLimit	0.67	NoLimit	0.34	NoLimit	0.80
		2	0.60	2	0.30	3	0.68
	MWFFreq	3	0.57	3	0.29	2	0.63
		5	0.54	7	0.28	5	0.58
		7	0.49	5	0.27	7	0.55
		Zero 0.05	0.68	Zero -0.2	0.34	Zero 0	0.80
		Zero 0.1	0.68	Limit -0.1	0.34	Limit -0.02	0.80
	MWFWeight	Zero -0.05	0.68	Limit -0.2	0.34	Zero -0.05	0.80
		Limit -0.01	0.67	Limit -0.5	0.34	Limit -0.01	0.80
		Zero 0.02	0.67	NoLimit	0.34	Zero -0.01	0.80
		NoLimit	0.67	100	0.36	2	0.80
		2	0.67	50	0.36	NoLimit	0.80
	MFFreq	3	0.67	30	0.35	3	0.80
		5	0.67	20	0.35	20	0.80
		7	0.67	15	0.35	15	0.80

Table 3 .

 3 The top 5 performing CPSs for each language with their achieved scores, in descending order of maximum H scores

	Lang #					Parameter settings	P	S	H
			VecSim				Weight	FeatTransf
		1	PearsMbAdjCosMod-3.Lb DimRed Smooth VNorm StopW MWFFreq MWFWeight MFFreq PmiAl-Tc3Tw0S2P0 NoTransf	0.72 0.71 0.71
			IslamInkpen 0.05 NoSmooth	L1	false	NoLimit	Zero 0	NoLimit
			VecSim				Weight	FeatTransf
		2	PearsMbAdjCosMod-3.Lb DimRed Smooth VNorm StopW MWFFreq MWFWeight MFFreq PmiAl-Tc3Tw0S2P0 NoTransf	0.72 0.71 0.71
			NoDimRed	NoSmooth	L1	false	NoLimit	Zero 0	NoLimit
			VecSim				Weight	FeatTransf
	En	3	PearsMbAdjCosMod-3.Lb DimRed Smooth VNorm StopW MWFFreq MWFWeight MFFreq PmiAl-Tc3Tw0S2P0 NoTransf	0.72 0.71 0.71
			NoDimRed	NoSmooth	L1	false	NoLimit	Zero -0.05	NoLimit
			VecSim				Weight	FeatTransf
		4	PearsMbAdjCosMod-3.Lb DimRed Smooth VNorm StopW MWFFreq MWFWeight MFFreq PmiAl-Tc3Tw0S2P0 NoTransf	0.72 0.71 0.71
			IslamInkpen 0.05 NoSmooth	L1	false	NoLimit	Zero -0.05	NoLimit
			VecSim				Weight	FeatTransf
		5	PearsMbAdjCosMod-4.Lb DimRed Smooth VNorm StopW MWFFreq MWFWeight MFFreq PmiAl-Tc3Tw0S2P0 NoTransf	0.72 0.71 0.71
			NoDimRed	NoSmooth	L1	false	NoLimit	Zero -0.05	NoLimit
			VecSim				Weight	FeatTransf
		1	DimRed	Cos		Pmi-Tc1Tw3S2P0 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Lb	0.43 0.44 0.44
			SVD 100	NoSmooth	L2	true	NoLimit	NoLimit	3
			VecSim				Weight	FeatTransf
		2	DimRed	Cos		PmiAl-Tc3Tw3S2P0 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Lb	0.43 0.43 0.43
			SVD 100	NoSmooth	L2	true	NoLimit	NoLimit	3
			VecSim				Weight	FeatTransf
	Es	3	DimRed	Cos		Pmi-Tc1Tw3S2P0 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Lb	0.43 0.43 0.43
			SVD 100	NoSmooth	L2	true	NoLimit	NoLimit	100
			VecSim				Weight	FeatTransf
		4	DimRed	Cos		PmiAl-Tc3Tw3S2P0 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Lb	0.43 0.43 0.43
			SVD 100	NoSmooth	L2	false	NoLimit	NoLimit	3
			VecSim				Weight	FeatTransf
		5	DimRed	Cos		Pmi-Tc1Tw3S2P0 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Lb	0.43 0.43 0.43
			SVD 100	NoSmooth	L1	true	NoLimit	NoLimit	3
	Lang #					Parameter settings	A
			VecSim				Weight	FeatTransf
		1	MbCosAm DimRed	NPmiAlpha-Tc4Tw4S0P4 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Sqrt	0.65
			SVD 200	NoSmooth	L2	false	NoLimit	NoLimit	2
			VecSim				Weight	FeatTransf
		2	MbCosAm DimRed	NPmiAlpha-Tc4Tw4S0P4 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Sqrt	0.65
			SVD 200	NoSmooth	L1	false	NoLimit	NoLimit	2
			VecSim				Weight	FeatTransf
	Hu	3	DimRed	Cos		NPmiAlpha-Tc4Tw4S0P4 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Sqrt	0.65
			SVD 200	NoSmooth	L2	false	NoLimit	NoLimit	2
			VecSim				Weight	FeatTransf
		4	DimRed	Cos		NPmiAlpha-Tc4Tw4S0P4 Smooth VNorm StopW MWFFreq MWFWeight MFFreq Weight AftNorm Sqrt	0.65
			SVD 200	NoSmooth	L1	false	NoLimit	NoLimit	2
			VecSim				Weight	FeatTransf
		5	PearsMbMod-1.Lb DimRed Smooth VNorm StopW MWFFreq MWFWeight MFFreq NPmiAlpha-Tc4Tw4S0P4 Weight AftNorm Sqrt	0.63
			SVD 200	NoSmooth	L2	false	NoLimit	NoLimit	2

Table 4 .

 4 Results on the test datasets, in descending order of maximum H scores

	Lang Test set CPS	P	S	H
			BestEn 0.71 0.71 0.71
	En	MT	BestHu 0.67 0.68 0.67
			BestEs 0.63 0.63 0.63
			BestHu 0.83 0.83 0.83
	Es	RG	BestEn 0.82 0.80 0.81
			BestEs 0.80 0.79 0.80
			BestEn 0.73 0.72 0.72
	Hu	RG	BestEs 0.65 0.61 0.63
			BestHu 0.58 0.68 0.62

Table 4 )

 4 and the MT dataset for all tests too. It simply suggests that the dataset used for Spanish in this phase is considerably tougher than the ones used for English and Hungarian.

	CPS	P	S	H
	BestEn	0.67 0.67 0.67
	OSC of [29] 0.59 0.58 0.58
	CosPPmi 0.56 0.58 0.57
	the test datasets (see


Konstantinos Demertzis
email: kdemertz@fmenr.duth.gr1

Lazaros Iliadis
email: liliadis@civil.duth.gr2

Panayiotis Kikiras
email: kikirasp@uth.gr3

Nikos Tziritas

Cyber-Typhon: An Online Multi-task Anomaly Detection Framework

Keywords: Deep Content Inspection, Anomaly Detection, Multi-Task Learning, Online Learning, Restricted Boltzmann Machine, Critical Infrastructure Protection

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Information generated by complex environments such as the Internet of Things (IOT) ecosystem, has increased exponentially. The result is the inefficient management and storage of the total volume of generated information. This requires the adoption of complex data mining and analysis architectures [START_REF] Dedić | Towards Differentiating Business Intelligence, Big Data, Data Analytics and Knowledge Discovery[END_REF]. These architectures should incorporate specialized processing algorithms, that dynamically adapt to new standards or data, or to scaled data production as a equation of time [START_REF] Kiran | Lambda architecture for cost-effective batch and speed big data processing[END_REF]. Although mining Data Streams (DAST) is an emerging area, it poses enormous challenges to the data mining community. High transmission speed, change of data distribution, and high volume, raise the following issues that should be addressed [START_REF] Lin | The Lambda and the Kappa[END_REF]: ➢ High Velocity: Online DAST, arrive at a very high speed. Thus, it has become almost practically infeasible to scan all of them. This is also the case for the offline ones. ➢ Concept Drift: Frequent Patterns (FREP) keep changing, as data streams are time varying in nature. During the mining process, as new incoming data is added to the existing ones, some FREP may change their status to become infrequent and vice versa. This issue is known as the Concept Drift (CDR) problem. ➢ Unbounded Size (UNS): DAST are unbounded in size. Their size is unknown to the user in advance unlike the static data.

➢ Enormous Space Requirement (ESR): Huge amount of data are generated both in online and offline applications. There might not be enough space to store the data stream before processing. ➢ Unsteady Analysis Results (UAR): High speed as well as varying data distribution, may affect the analyzed results, due to which the mining outcome may be declined. In order to cope with this, DAST mining must be an incremental process. Anomaly Detection [START_REF] Demertzis | A Hybrid Network Anomaly and Intrusion Detection Approach Based on Evolving Spiking Neural Network Classification. e-Democracy[END_REF] over multiple data streams, is initially determined by the observation of a single (multi-variable) time series frequency, which constitutes the systems' quantitative performance parameters. A data stream 𝑠𝑖 may be coming from a system of sensors, and it consists of numerical values, where 𝑠𝑖(𝑡) stands for the data stream flow value at the time 𝑡, and 𝑡 𝜖 [0, +∞]. If 𝑛 sensors' flows are synchronized periodically to report their values, the whole of the multi-variable information at each time 𝑡 is represented with the following frame vector shown in equation 1, [START_REF] Krawczyk | Online ensemble learning with abstaining classifiers for drifting and noisy data streams[END_REF].

ΔΠ 𝑡 = (𝑠 1 (𝑡), 𝑠 2 (𝑡), … , 𝑠 𝑛 (𝑡) ∈ 𝑅 𝑛 ) (1) In practice, each flow forms a one-dimensional time series, while the frame vector flow (FVF) represents a multi-variable time series. Event detection on DAST is intended to determine the values 𝑠𝑖 (𝑡), which represent abrupt changes within an FVF. Particularly, each FVF of length 𝑛 is converted to a binary vector of the same length, where each value represents a possible change in the corresponding sensor flux. Such deviations from the normal behavior are called events, and binary vectors are called event vectors. An event may be an observation that does not conform to an expected standard in the data set (anomaly). Incidents may have been caused by a variety of reasons, like sensor failure or malequation, or deviations and substantial changes that may affect the system's behavior, such as Cyber Attacks (CYA) [START_REF] Zhang | Parallel multi-task learning[END_REF]. Therefore, a vector of events in time 𝑡 is represented by equation 2 below [START_REF] Lin | The Lambda and the Kappa[END_REF], [START_REF] Krawczyk | Online ensemble learning with abstaining classifiers for drifting and noisy data streams[END_REF].

ΔΣ 𝑡 = (𝑒 1 𝑡 , 𝑒 2 𝑡 , … , 𝑒 𝑛 𝑡 ) ∈[0,1] (2)
where 𝑒 𝑖 𝑡 = 𝑒 𝑖 (𝑡) is the binary value which represents the occurrence of an abnormal flow behavior, and its value is 𝑠𝑖(𝑡)=1 in time 𝑡. The values of 𝑠𝑖(𝑡) must be in an interval {0,1}. The error is calculated at each iteration, as data characteristics can change drastically and in an unpredictable way changing the typical, normal behavior. An object that may be considered abnormal, can then be included in the set of normal observations due to rapid developments in the data stream. Due to the fact that the data volume is unlimited, data mining is performed on a subset of the flow, called a sliding window, which obviously contains a small but recent percentage of the observations. The goal of the data flow processing algorithms is to minimize the cumulative error for all iterations, that can be calculated by the following equation 3 [START_REF] Lin | The Lambda and the Kappa[END_REF], [START_REF] Krawczyk | Online ensemble learning with abstaining classifiers for drifting and noisy data streams[END_REF].

𝐼 𝑛 [𝑤] = ∑ 𝑉(〈𝑤, 𝑥 𝑗 〉, 𝑦 𝑗 ) 𝑛 𝑗=1 = ∑(𝑥 𝑗 𝑇 𝑤 -𝑦 𝑗 ) 2 𝑛 𝑗=1 (3)
where 𝑥 𝑗 ∈ 𝑅 𝑑 , 𝑤 ∈ 𝑅 𝑑 and 𝑦 𝑗 ∈ 𝑅. We consider that Χ𝑖 × 𝑑 is a data matrix and Υ𝑖 × 1 is a matrix with target values, after the arrival of the first i data points. Assuming that the covariance matrix 𝛴 𝑖 = 𝑋 𝑇 𝑋 is reversible, the optimal solution 𝑓 * (𝑥) = 〈𝑤 * , 𝑥〉 is given by the following equation 4 [START_REF] Lin | The Lambda and the Kappa[END_REF], [START_REF] Krawczyk | Online ensemble learning with abstaining classifiers for drifting and noisy data streams[END_REF].

𝑤 * = (𝑋 𝑇 𝑋) -1 𝑋 𝑇 𝛶 = 𝛴 𝑖 -1 ∑ 𝑥 𝑗 𝑦 𝑗 𝑖 𝑗=1 (4)
First the covariance table is calculated by using the following equation 5

𝛴 𝑖 = ∑ 𝑥 𝑗 𝑖 𝑗=1
𝑥 𝑗 𝑇 [START_REF] Krawczyk | Online ensemble learning with abstaining classifiers for drifting and noisy data streams[END_REF].

The initial time complexity (CM) is calculated to be of 𝛰(𝑖𝑑 2 ) order, but after we inverse the 𝑋 𝑇 𝑋 (𝑑 × 𝑑) table, it increases to 𝛰(𝑑 3 ), while the rest of the required multiplications have an 𝛰(𝑑 2 ) 𝐶𝑀. This produces an overall CM of 𝛰(𝑖𝑑 2 + 𝑑 3 ) [3] [START_REF] Krawczyk | Online ensemble learning with abstaining classifiers for drifting and noisy data streams[END_REF] (d is the size of the window).

It is conceivable that robust systems ensuring reliability and high accuracy rates without requiring high availability of resources are required to safely approach problems stemming from knowledge mining processes. The above argument is further supported as follows: Let's suppose that the size of data points is equal to n and after the arrival of each new data point i = 1,2, ..., n, the recalculation of the solution is required. In this case the total time complexity would be equal to 𝛰(𝑛 2 𝑑 2 + 𝑛𝑑 3 ) [START_REF] Krawczyk | Online ensemble learning with abstaining classifiers for drifting and noisy data streams[END_REF]. The Multi-Τask Learning (MTL) is a robust method in order to face some of the most challenge of the Big Data Streams processing with Online Learning algorithms. MTL is a subfield of Machine Learning (ML) in which multiple learning tasks are solved at the same time, utilizing common elements or differences arising from the multiple tasks included in the case study [6][7]. More general MΤL is an inductive transfer method which generates many generalization features. The common features or differences that arise between distributed tasks during the training process are transferred or shared as guaranteed and unambiguous knowledge in subsequent relevant or unrelated tasks, maximizing the accuracy of the model. MTL is efficient because regularization induced by requiring an algorithm to perform well on a related task can be superior to regularization that prevents overfitting by penalizing all complexity uniformly. The following approaches are characteristic cases of MTL [START_REF] Baytas | Asynchronous multi-task learning ICDM[END_REF][7]: ➢ Τask grouping and overlapping, where tasks are grouped or provided in an overlapping way so that there is relevance, capable to lead in the use of cognitive or learning relationships. ➢ Exploiting unrelated tasks, where the common learning of non-relevant tasks using the same input data, can be beneficial for learning main tasks of an application field. ➢ Transfer of knowledge, where transfer of relevant knowledge is carried out to achieve learning from correspondingly trained models. ➢ Group online adaptive learning, where transfer of previous experience or knowledge into continually changing environments takes place.

Literature Review

Qian and Sherif [START_REF] Chen | A model-based approach to self-protection in computing system[END_REF] have applied autonomous computing technology, to monitor SCADA system performance. Their approach proactively estimates upcoming attacks for a given system model of a physical infrastructure. Soupionis et al. [START_REF] Soupionis | [END_REF] have proposed a combinatorial method for automatic detection and classification of faults and cyber-attacks occurring on the power grid system when there is limited data from the power grid nodes due to cyber implications. In addition, Tao et al. [START_REF] Zhu | Detecting node replication attacks in wireless sensor networks: a survey[END_REF] have described the network attack knowledge, based on the theory of the factor expression of knowledge. They have studied the formal knowledge theory of SCADA network from the factor state space and equivalence partitioning. This approach utilizes the Factor Neural Network (FNN) theory which contains high-level knowledge and quantitative reasoning, used to establish a predictive model including analytic FNN and analogous FNN. This model abstracts and builds an equivalent and corresponding network attack and a defense knowledge factors system. Finally [START_REF] Cruz | Improving Cyber-Security Awareness on Industrial Control Systems: The CockpitCI Approach[END_REF] has introduced a new European Framework-7 project Cockpit CI (Critical Infrastructure) and roles of intelligent machine learning methods to prevent SCADA systems from cyber-attacks. Also, existing multi-task learning methods can be categorized into two main categories: learning with feature covariance and learning with task relations [START_REF] Zhang | A convex formulation for learning task relationships in multi-task learning[END_REF]. Different from prior solutions to distributed multi-task learning, which are focused on the former category [START_REF] Wang | Distributed multi-task learning[END_REF], our proposed multi-task learning falls into the latter category. On the other hand, distributed machine learning has attracted more and more interests recently [START_REF] Xing | Strategies and Principles of Distributed Machine Learning on Big Data[END_REF].

There have been tremendous efforts done on different machine learning problems. Also, online multi-task learning assumes instances from different tasks arrive in a sequence and adversarial chooses task to learn. Cavallanti et al. [START_REF] Cavallanti | Linear algorithms for online multitask classification[END_REF] exploited online multi-task learning with a given task relationship encoded in a matrix, which is known beforehand. Parallel multi-task learning aims to develop parallel computing algorithms for multi-task learning in a shared-memory computing environment.

Recently, Zhang [START_REF] Zhang | Parallel multi-task learning[END_REF] proposed a parallel multi-task learning algorithm named PMTL. In this work, dual forms of three losses are presented and accelerated proximal gradient method is applied to make the problem decomposable, and thus possible to be solved in parallel. Finally, distributed multi-task learning is an area that has not been much exploited.

Wang et al. [START_REF] Wang | Distributed multi-task learning[END_REF] proposed a distributed algorithm for multi-task learning by assuming that different tasks are related through shared sparsity. In another work [START_REF] Baytas | Asynchronous multi-task learning ICDM[END_REF], asynchronous distributed multi-task learning method is proposed for multi-task learning with shared subspace learning or shared feature subset learning. Different from the above-mentioned approaches, our method aims at solving multi-task learning by learning task relationships from data, which can be positive, negative, or unrelated, via a task-covariance matrix.

The proposed Cyber-Typhon Framework

This research proposes the Cyber-Typhon model, which combines the algorithmic power of OS-ELM and RBM in a hybrid mode. It is an innovative Multi-Task Learning approach [START_REF] Demertzis | MOLESTRA: A Multi-Task Learning Approach for Real-Time Big Data Analytics[END_REF], that performs control of network traffic in Critical Infrastructures [START_REF] Demertzis | A Spiking One-Class Anomaly Detection Framework for Cyber-Security on Industrial Control Systems[END_REF][18] [START_REF] Demertzis | A Computational Intelligence System Identifying Cyber-Attacks on Smart Energy Grids[END_REF]. The final target is the detection of vulnerabilities that are usually due to APT attacks [START_REF] Demertzis | The Next Generation Cognitive Security Operations Center: Network Flow Forensics Using Cybersecurity Intelligence[END_REF][21] [START_REF]Cyber-Security and Information Warfare[END_REF]. More specifically, the Cyber-Typhon initially exports features related to network traffic, which are used as input to an OS-ELM neural network. The OS-ELM has been trained with proper data, in order to be able either to classify traffic as normal or (in the opposite case) to identify the threat or the attack type. Obviously, it performs Multiclass classification in order to identify one of the following eight [START_REF] Chen | A model-based approach to self-protection in computing system[END_REF] If the network traffic is normal further communication is allowed. In the opposite case, the type of anomaly is determined and the data flow is redirected to a proper absolutely specialized and dedicated Restricted Boltzmann Machine. If the first RBM does not recognize the specific anomaly for which it is specialized, the data is redirected to the next RBM responsible for the detection of another anomaly and so on till the successful identification is achieved. If detection cannot be done by any of the trained RBM (which are as many as the types of the known anomalies) the network flow data return to the initial OS-ELM, which can perform online sequential learning. Thus, the classification effort can be re-adjusted. The following figure 1 presents the architecture of the Cyber-Typhon. The proposed ELM is a Single-Hidden Layer Feed Forward Neural Network (SHLF2N2) [START_REF] Huang | Extreme learning machine: theory and applications[END_REF] with N hidden neurons, randomly selected input weights and random values of bias in the hidden layer. The output weights are calculated with a single vector matrix multiplication [START_REF] Wang | Distributed multi-task learning[END_REF]. Hidden nodes or hidden level parameters can be randomly created before seeing the training data, while it is remarkable that non-differential activation equations can be handled, and known Neural Network problems such as stopping criterion, learning rate and learning epochs are not addressed. Specifically, the input data is mapped to a random L-dimensional space with a discrete training set N, where (𝑥 𝑖 , 𝑡 𝑖 ), 𝑖 ∈ ⟦1, 𝑁⟧ with 𝑥 𝑖 ∈ 𝑅 𝑑 and 𝑡 𝑖 ∈ 𝑅 𝑐 . The specification output of the network is given by the following equation 6 [START_REF] Huang | Extreme learning machine: theory and applications[END_REF]:

𝑓 𝐿 (𝑥) = ∑ 𝛽 𝑖 ℎ 𝑖 (𝑥) = ℎ(𝑥)𝛽 𝐿 𝑖=1 𝑖 ∈ ⟦1, 𝑁⟧ (6)
Vector matrix 𝛽 = [𝛽 1 , … , 𝛽 𝐿 ] 𝑇 is the output of the weight vector matrix connecting hidden and output nodes. On the other hand, ℎ(𝑥) = [𝑔 1 (𝑥), … , 𝑔 𝐿 (𝑥)] is the output of the hidden nodes for input x, and 𝑔 1 (𝑥) is the output of the ith neuron. Based on a training set {(𝑥 𝑖 , 𝑡 𝑖 )} 𝑖=1 𝑁 , the ELM can solve the Learning Problem 𝛨𝛽 = 𝛵, where 𝛵 = [𝑡 1 , … , 𝑡 𝑁 ] 𝑇 are the target labels and the output vector matrix of the Hidden Layer 𝐻 is given by the following equation 7 [START_REF] Huang | Extreme learning machine: theory and applications[END_REF]:

𝛨(𝜔 𝑗 , 𝑏 𝑗 , 𝑥 𝑖 ) = [ 𝑔(𝜔 1 𝑥 1 + 𝑏 1 ) ⋯ 𝑔(𝜔 𝑙 𝑥 1 + 𝑏 𝑙 ) ⋮ ⋱ ⋮ 𝑔(𝜔 1 𝑥 𝑁 + 𝑏 1 ) ⋯ 𝑔(𝜔 𝑙 𝑥 𝑁 + 𝑏 𝑙 )
] 𝑁×𝑙 [START_REF] Zhang | Parallel multi-task learning[END_REF] The input weight vector matrix of the hidden layer ω (before training) and the bias vectors b are created randomly in the interval [-1,1], by employing equations 8a and 8b.

𝜔 𝑗 = [𝜔 𝑗1 , 𝜔 𝑗2 , … , 𝜔 𝑗𝑚 ] 𝛵 (8𝑎) and 𝛽 𝑗 = [𝛽 𝑗1 , 𝛽 𝑗2 , … , 𝛽 𝑗𝑚 ] 𝛵 (8𝑏)
The output weight vector matrix of the hidden layer H is calculated by the use of the activation equation in the training dataset, based on equation 9a and the output weights β are calculated based on equation 9b.

𝛨 = 𝑔(𝜔𝑥 + 𝑏) (9a) 𝛽 = ( 𝛪 𝐶 + 𝐻 𝑇 𝐻) -1 𝐻 𝑇 𝑋 (9b)
where 𝐻 = [ℎ 1 , … , ℎ 𝑁 ] is the output vector matrix of the hidden layer and 𝛸 = [𝑥 1 , … , 𝑥 𝑁 ] is the input vector matrix of the hidden layer. Indeed, β can be calculated by the following equation 10:

𝛽 = 𝛨 + 𝑇 ( 10
)
where + is the generalized inverse Moore-Penrose vector for matrix H [START_REF] Huang | Extreme learning machine: theory and applications[END_REF]. The Cyber-Tyfon employs the ELM algorithm, which employs the kernel of the following Gaussian Radial Basis Equation K(u,v)=exp(-γ||u-v|| 2 ) (5).

The number of the hidden neurons k is equal to 20 and the assigned random input weights are wi where bi, i=1,…,N are the bias. The calculation of the hidden layer output matrix H has been done by employing equation 11 below [START_REF] Huang | Extreme learning machine: theory and applications[END_REF].

𝐻(ℎ) = [ ℎ(𝑥 1 ) ⋮ ℎ(𝑥 𝑁 ) ] = [ ℎ 1 (𝑥 1 ) ⋯ ℎ 𝐿 (𝑥 1 ) ⋮ ⋮ ℎ 1 (𝑥 𝑁 ) ⋯ ℎ 𝐿 (𝑥 𝑁 ) ] (11)
where h(x) = [h1(x), . . . , hL(x)] is the output (row) vector of the hidden layer with respect to the input x. Also, h(x) maps the data from the D-dimensional input space to the L-dimensional hidden-layer feature space (ELM feature space) H. Thus, h(x) is indeed a feature mapping. ELM aims to minimize the training error ||Hβ -T|| 2 as well as the norm ||β|| of the output weights: where H is the hidden-layer output matrix of the equation 11.

Minimization of the norm of the output weights ||β|| is actually achieved by maximizing the distance of the separating margins of the two different classes in the ELM feature space 2/||β||.

To calculate the output weights β we used equation 12 [START_REF] Huang | Extreme learning machine: theory and applications[END_REF].

β=( I C + H T H) -1 H T T (12)
where the value of C (a positive constant) and the value of T are obtained from the Equation Approximation of the SHLF2N2 with additive neurons (see equation 13 below).

ti=[ti1, ti2,…,tim] T ∈ R m (13), T = [ t 1 T ⋮ t N T ]
The OS-ELM [START_REF] Demertzis | MOLESTRA: A Multi-Task Learning Approach for Real-Time Big Data Analytics[END_REF] is an alternative technique for large-scale computing and machine learning approaches that used when data becomes available in a sequential order to determine a mapping from data set corresponding labels.

The main difference between online learning and batch learning techniques is that in online learning the mapping is updated after the arrival of every new data point in a scale fashion, whereas batch techniques are used when one has access to the entire training data set at once. It is a versatile sequential learning algorithm because of the training observations are sequentially (one-by-one or chunk-by-chunk with varying or fixed chunk length) presented to the learning algorithm. At any time, only the newly arrived single or chunk of observations (instead of the entire past data) are seen and learned. A single or a chunk of training observations is discarded as soon as the learning procedure for that particular (single or chunk of) observation(s) is completed. The learning algorithm has no prior knowledge as to how many training observations will be presented. Unlike other sequential learning algorithms which have many control parameters to be tuned, OS-ELM only requires the number of hidden nodes to be specified [START_REF] Huang | Extreme learning machine: theory and applications[END_REF].

The proposed OS-ELM algorithm consists of two main phases namely:

The Boosting (BPh) and the Sequential Learning (SLPh). Obviously, if the network flow is classified as normal, it is allowed to continue. In the opposite case, there are 7 RΒΜs, as many as the abnormal classes, where each one of them has been trained to perform One-Class Classification (OCC) in order to exclusively recognize one specific network attack [START_REF] El-Yaniv | Optimal single-class classification strategies -google scholar[END_REF]. The OCC also known as Unary Classification Method (UCΜ), implements intelligent categorization of cases belonging to a specific class, among an existing set of records. OCC is learning from a training set that contains only records of one specific class. Usually in this method the negative class is absent because it is not sampled. Thus, the division boundaries are determined effectively only with the knowledge of the positive class. This process is much more difficult than a traditional binary or multiclass classifier [START_REF] Munroe | Multi-class and single-class classification approaches to vehicle model recognition from images[END_REF], as it is trained to accept target objects and to reject the ones that have significant deviation. Minimizing the errors is also a difficult process, because in this type of categorization, cross validation is unavailable since there is no data from the other classes. Finally, it should be stressed that one class problem-solving technique is inverse to the generalization approaches that are pursued in other machine learning problems, as it tends to provide a fully defined configuration of parameters. This can exponentially increase the complexity of the classifier, trying to correctly classify target data. The more complex the model, the smaller the rank range in the target data range, and the less likely it is for the outliers to be categorized correctly. In practice, one can create a complex model by setting all its possible parameters without being at risk from overfitting.

In this sense, the ΟCC is the most appropriate approach for detecting abnormalities and identifying patterns or trends, in a set of data that displays divergent behaviors than expected. ΟCC achieves high levels of successful detection, while maintaining low false error rates (false alarm) [24] [25].

The RBM [START_REF] Zhang | An overview on Restricted Boltzmann Machines[END_REF] belong to the family of energy-based models, where each configuration of the variables of interest corresponds to a finite scalar energy value used for training. The learning process is performed by modifying the energy equation (ENF), so that its shape has desirable. Specifically, RBM is a symmetric graphical model. The units of one layer are connected (and thus dependent) only with units of the next one. The proposed RBM ENF with V visible units and H hidden ones is defined as follows equation 14 [START_REF] Zhang | An overview on Restricted Boltzmann Machines[END_REF]: [START_REF] Xing | Strategies and Principles of Distributed Machine Learning on Big Data[END_REF] where v and h are binary vectors related to the state of visible units, and to the state of hidden units respectively. Moreover, 𝑣 𝑖 and ℎ 𝑗 correspond to the individual state of each visible unit (VU) i, and each hidden unit (HU) 𝑗 respectively, and 𝑤 𝑖𝑗 is the weight assigned to their connection. Finally, 𝑏 𝑖 𝑣 and 𝑏 𝑗 ℎ is the bias of the VU i and the HU j. The reserved probability 𝑝(𝑣|ℎ) is given by the following equation 15:

𝛦(𝑣, ℎ) = -∑ ∑ 𝑣 𝑖 𝐻 𝑗=1 𝑉 𝑖=1 ℎ 𝑗 𝑤 𝑖𝑗 -∑ 𝑣 𝑖 𝑏 𝑖 𝑣 𝑉 𝑖=1 -∑ ℎ 𝑗 𝑏 𝑗 ℎ 𝑉 𝑗=1
𝑝(𝑣|ℎ) = 𝑒 -𝐸(𝑣,ℎ) ∑ 𝑒 -𝐸(𝑔,ℎ) 𝑔 (15)
In the specific case of examining a unit i of the visible level, the assigned reserved probability distribution (if we know the status of the hidden layer h), is calculated by equation 16:

𝑝(𝑣 𝑘 = 1|ℎ) = 1 1+𝑒 -∑ ℎ 𝑗 𝑤 𝑘𝑗 +𝑏 𝑣 𝑘 𝐻 𝑗=1 (16)
The reserved probabilities 𝑝(ℎ|𝑣) and 𝑝(ℎ 𝑘 = 1|𝑣) are defined by equations 17 and 18:

𝑝(ℎ|𝑣) = 𝑒 -𝐸(𝑣,ℎ) ∑ 𝑒 -𝐸(𝑣,𝑔) 𝑔 (17) and 𝑝(ℎ 𝑘 = 1|𝑣) = 1 1+𝑒 -∑
𝑣 𝑖 𝑤 𝑘𝑗 +𝑏 ℎ 𝑘 𝑉 𝑖=1 [START_REF] Demertzis | An innovative soft computing system for smart energy grids cybersecurity[END_REF] These relations express the independence of the units of the two levels from the units of the same level. The training process of the RBM is the process of finding values for its parameters, that maximize the mean logarithmic probability of the occurrence of set C (equation 19).

∑ log 𝑝

∑ 𝑒 -𝐸(𝑣 𝑐 ,𝑔) 𝑔 ∑ ∑ 𝑒 -𝐸(𝑢,𝑔) 𝑔 𝑢 𝐶 𝑐=1 [START_REF] Demertzis | A Computational Intelligence System Identifying Cyber-Attacks on Smart Energy Grids[END_REF] The following equation 20 estimates the cost equation to wij that represents the renewal of the weights: The first term calculates the average values of 𝑣 𝑖 𝑐 , 𝑔 𝑗 when the visible level of RBM is leaded by the data (𝑣 𝑐 ), whereas the second term corresponds to the values of 𝑣 𝑖 , 𝑔 𝑗 when the data are "produced" by the model. An equivalent way of formulating, would suggest that every weight 𝑤 𝑖𝑗 should change to become equal to 𝛥𝑤 𝑖𝑗 (equation 21):

𝛥𝑤 𝑖𝑗 = 𝜀 𝑤 (𝐸 𝑑𝑎𝑡𝑎 [𝑣 𝑖 ℎ 𝑗 ] -𝐸 𝑚𝑜𝑑𝑒𝑙 [𝑣 𝑖 ℎ 𝑗 ]) (21
) Τhe RBM model in this case starts approaching the actual values of the data. The first term 𝐸 𝑑𝑎𝑡𝑎 [𝑣 𝑖 ℎ 𝑗 ] can be calculated easily, since knowing the values of the units at the visible level, we can, by means of the above equations, calculate the reserved probability for each unit at the hidden level. Calculation of the second term, which presupposes the existence of samples from the model itself, is more difficult. In this paper we use optimization through Contrastive Divergence (CD) [START_REF] Ma | Convergence Analysis of Contrastive Divergence Algorithm Based on Gradient Method with Errors[END_REF]. The steps that we used to procedure a single sample can be summarized as follows: 1. Take a training sample v, compute the probabilities of the hidden units and sample a hidden activation vector h from this probability distribution. 2. Compute the outer product of v and h and call this the positive gradient. 3. Use h, to sample a reconstruction v' of the visible units, and then resample the hidden activations h'. (Gibbs sampling step) 4. Compute the outer product of v' and h' and call this the negative gradient.

5. Let the update to the weight matrix W be the positive gradient minus the negative gradient, times some learning rate: 𝛥𝑊 = 𝜀(𝑢ℎ 𝑇 -𝑢 ′ ℎ ′𝑇 ). 6. Update the biases a and b analogously: 𝛥𝑎 = 𝜀(𝑢 -𝑢 ′ ), 𝛥𝑏 = 𝜀(ℎ -ℎ ′ ). The CD method gives lower energy to the actual data and much higher energy to the "reconstructions" resulting from them. This helps the model approach the actual data distribution. Also, the Typhon employs the Multi-Task Grouping and Overlap Learning approach, combined with the optimal use of the OS-ELM and RBMs methodologies. The sliding windows (SLIW) [START_REF] Dietterich | Machine Learning for Sequential Data: A Review[END_REF] are used to partition the data stream. Using the SLIW technique, the system estimates a table of indices, after accepting the data flow vectors as input. The data are divided in partitions of 1,000 samples with 400 of them overlapping between adjoining windows. This enables continuous and unintentional scanning of the data, which results in faster and more accurate control results. This happens because a small SLIW is much more likely to be more uniform than a larger one, and therefore it is more predictable. Each new sample is checked by the OS-ELM algorithm and force it in the appropriate RBM. Once the optimal model has been created, it is then applied to the control window of the SLIW. This process is followed until the input data of each window is completed. The total knowledge of the entire window is stored in MM1 (the first window). The MM1 is transferred in the window 2 as MLT, which means that the total knowledge of the first window is transferred to the second window. The process continues for all the data in window 2, and so on for all other windows, until the full analysis of the data flow is done. The full representation of the proposed process is presented in figure 2 below. When the OS-ELM classification detects an attack type, the network flow data are directed to the corresponding RBM. If it is a false positive case, the data are directed to the next RBM and so on. If the output is a false alarm, then the traffic is redirected to the initial OS-ELM, which has the online learning potential and it re-examines the dataflow from the beginning as if it was a new one. This architecture has been chosen due to the fact that in multi parametric and high complexity problems related to big data (like the one examined herein) the classification results are unstable especially regarding the analysis and the incorporation of the data flows. The introduced architecture is a serious method of a hybrid combined resolution of the APT attacks. The system is flexible and effective and it not only offers a robust approach but it also offers a faster convergence of the overall model. Finally, the proposed architecture can exploit the multi-tasking learning in order to enhance generalization. Each one of the developed RBMs is totally dedicated to the specific problem, which can result in bias and variance reduction. This can result in overfitting elimination and it offers a robust framework capable of coping with high complexity problems.

Datasets

Appropriate datasets were chosen that closely simulate Industrial Control System (ICS) communication and transaction data. They were used in the development and evaluation of the proposed model. Contained preprocessed network transaction data and preprocessed to strip lower layer transmission data, were used (e.g. TCP, MAC) [START_REF] Morris | Industrial Control System Simulation and Data Logging for Intrusion Detection System Research[END_REF]. Specifically, the Gas_Dataset that was chosen for the purpose of this research includes [START_REF] Zhang | An overview on Restricted Boltzmann Machines[END_REF] The dataset is determined and normalized in the interval [-1,1] in order to phase the problem of prevalence of features with wider range over the ones with a narrower range, without being more important. Also, the outliers and the extreme values spotted were removed based on the Inter Quartile Range technique [START_REF] Zwillinger | CRC Standard Probability and Statistics Tables and Formulae[END_REF]. The reader can find details regarding the dataset and the data collection and assessment methodology in [START_REF] Morris | Industrial Control System Simulation and Data Logging for Intrusion Detection System Research[END_REF].

Training the RBMs

The RBMs were trained based on an innovative application of the One Class Classification methodology. More specifically, according to this approach the system is exclusively trained with data related to a vulnerability of the network, in order to be able to identify the specific behavior -anomaly that has a potential relation with APT attacks. This is achieved by estimating the RBMs Reconstruction Error (RER), which is a basic criterion for the safe determination of the classes. It is calculated by using a threshold which is unique for each class and which emerged after several trial and error attempts, aiming in the optimal output of the system. If the error is higher than this predefined limit (which is characteristic for each class) then it is rejected as an unknown class and the data are redirected to the next RBM.

Results

In data cases using a machine learning classifier, for estimating the real error during training, the full probability density of both categories should be known. The classification performance is estimated by the development of a Confusion Matrix (CM), where the main diagonal values (top left corner to bottom right) correspond to correct classifications and the rest of the numbers correspond to very few cases that were misclassified. The following table1 presents the CM results of the OS-ELM:

The following table 2 present the analytical results of the proposed method.

Discussion and Conclusions

This research introduces a highly innovative, reliable and effective anomaly detection system, which is based on advanced computational intelligence approaches [START_REF] Demertzis | A Bio-Inspired Hybrid Artificial Intelligence Framework for Cyber Security[END_REF][33] [START_REF] Demertzis | Evolving Smart URL Filter in a Zone-Based Policy Firewall for Detecting Algorithmically Generated Malicious Domains[END_REF]. The Cyber-Typhon performs a series of sophisticated anomaly detection equations, by using Multi-Task Learning and by effectively combining On-Line Sequential Extreme Learning Machines with Restricted Boltzmann Machines. It ensures (in the most effective and intelligent way) the safe network communication among the interacting devices in a critical infrastructure environment. The proposed system, significantly enhances the security mechanisms of the Critical Infrastructures, which are a constant target due to their high importance. Also, this architecture offers the potential of developing a safe platform that can control and integrate network transactions. This can be done without the need for human intelligence or for a central authority. Also, it has been proven that collective intelligence technologies offer a smart solution for the determination of digital security and for monitoring of assets related to critical infrastructure.

The development of the model was based on the extremely effective OS-ELM algorithm, and on the multi-task learning method, which can perform transfer of knowledge between relative processes that are executed in parallel.

The training process was based on the employment of the RBMs that were trained based on the Unary Classification method. This was done by using specific datasets, each one of them corresponding to the behavior of a certain attack, in order to ensure the absolute reliability of the classifier.

The performance of the proposed system was tested on a multidimensional dataset of high complexity. This has resulted after an extensive research into the operation of critical infrastructure control systems and after comparisons, audits and tests. The target was the identification of the most appropriate limits, which express and realistically represent the classes they represent. The high accuracy of the results has greatly enhanced the validity of the general methodology.

It is important to mention that this particular model is presented for the first time in the literature.

It constitutes an important guideline in further exploitation of intelligent technologies in the automations that compose industrial networks. Proposals for the development and future improvements of this system, should focus on further optimizing the parameters of the RBMs used in order to achieve an even more efficient, accurate and quicker classification, capable of dividing even more precisely the boundaries between the situations of systems.

It would be important to study the equation-extension of the proposed algorithm with metalearning methods. This could further improve the anomaly detection process. Finally, the introduced model can employ adaptive learning in order to gain self-improvement potentials. This would automate 100% the whole process.

Fig 1 .

 1 Fig 1. Graphic representation of the Cyber-Typhon

  The BPh is used to train the SLFNs by applying basic ELM, with a batch of training data in the initialization stage. The boosting training data are discarded as soon as BPh is completed. The volume of the required training data vectors is very small, and it can be equal to the number of hidden neurons. A detailed description of the OS-ELM classifier is done below, where equations 9a and 9b are employed.Phase 1 (BPh)The BPh for a small initial training setN={(xi, ti)|xi ∈ R n , ti ∈ R m , i=1, ••• , 𝑁 ̃} isthe following: a) Assign arbitrary input weights wi and biases bi or centers µi and their corresponding impact widths σi, i=1, ••• 𝑁 ̃, where 𝑁 ̃ is the number of hidden neurons used by the RBF kernel for a specific application. b) Calculate the initial hidden layer output matrix 𝐻 0 = [ℎ 1 ,••• , ℎ 𝑁 ̃ ] T where ℎ 𝑖 = [𝑔(𝑤 1 • 𝑥 𝑖 + 𝑏 1 ), ••• , 𝑔(𝑤 𝑁 ̃ • 𝑥 𝑖 + 𝑏 𝑁 ̃ )] T , 𝑖 = 1, ••• , 𝑁 ̃, and 𝑔 is the activation equation or the RBF kernel. c) Estimate the initial output weight 𝛽 (0) = 𝑀 0 𝐻 0 𝑇 𝑇 0 , where 𝑀 0 = (𝐻 0 𝑇 𝐻 0 ) -1 and 𝑇 0 = [𝑡 1 , … , 𝑡 𝑁 ̃]𝑇 . d) Set 𝑘 = 0 Phase 2 (SLPh) In the second SLPh, the OS-ELM algorithm learns the train data one-by-one, or chunk-by-chunk, and all of the training data are discarded once the learning procedure on them is completed. The actual steps of this phase (for each incoming observation (𝑥 𝑖 , 𝑡 1 )), are described below. a) Calculate the hidden layer output vector ℎ (𝑘+1) = [𝑔(𝑤 1 • 𝑥 𝑖 + 𝑏 1 ), ••• , 𝑔(𝑤 𝑁 ̃ • 𝑥 𝑖 + 𝑏 𝑁 ̃ )] T . b) Calculate the latest output weight 𝛽 (𝑘+1) by employing the Recursive Least-Squares (RLS) algorithm where 𝛽 ̂= (𝐻 𝑇 𝐻) -1 𝐻 𝑇 𝑇. c) Set 𝑘 = 𝑘 + 1, where xi ∈ R n , ti ∈ R m and 𝑖 = 𝑁 ̃+ 1, 𝑁 ̃+ 2, 𝑁 ̃+ 3: xi ∈ R n , ti ∈ R m and 𝑖 = 𝑁 ̃+ 1, 𝑁 ̃+ 2, 𝑁 ̃+ 3.

Fig 2 .

 2 Fig 2. The proposed Multi-Task Grouping and Overlap Learning approach

  independent parameters and 210,770 instances, from which 61,156 normal and 149,614 abnormal (7 categories of attacks: Naïve Malicious Response Injection (ΝΜRI) 16,578, Complex Malicious Response Injection (CMRI) 15,466, Malicious State Command Injection (MSCI) 28,152, Malicious Parameter Command Injection (ΜPCI) 30,548, Malicious Equation Code Injection (ΜFCI) 20,628, Denial of Service (DoS) 11,022 and Reconnaissance (Recon) 27,220). The following figure 3 offers a graphical representation of the abnormal class distribution in the final Gas_Dataset.

Fig. 3

 3 Fig.3 Abnormal class distribution

  The Precision (PRE) the Recall (REC), the F-Score and the Total Accuracy (TA) indices are defined as in equations 25, 26, 27 and 28 respectively[START_REF] Fawcett | An introduction to ROC analysis[END_REF]:

Table 1 .

 1 Confusion Matrix of the OS-ELM All of the performance metrics in this testing, were estimated based on the One Versus ALL approach, because it is a Multi-Class Classification case. The numbers of misclassifications are related to the False Positive (FP) and False Negative (FN) indices appearing in the confusion Matrix. The True Positive rate (TPR) also known as Sensitivity the True Negative rate also known as Specificity (TNR) and the False Positive Rate (FPR) are defined by using equations 22, 23, and 24 respectively[START_REF] Fawcett | An introduction to ROC analysis[END_REF].

	Normal	ΝΜRI	CMRI	MSCI	ΜPCI	ΜFCI	DoS	Recon
	59,826	428	93	289	453	2	65	0
	632	15,944	0	2	0	0	0	0
	40	0	15,426	0	0	0	0	0
	264	0	0	27,888	0	0	0	0
	503	0	0	0	29,900	125	20	0
	2	0	0	0	157	20,469	0	0
	139	0	0	1	24	0	10,858	0
	0	0	0	0	0	0	0	2,220

Table 2 .

 2 Classification Accuracy and Performance MetricsThe 10-fold cross validation (10_FCV) is employed to obtain performance indices. Crossvalidation is a technique to evaluate predictive models by partitioning the original sample into a training set to train the model, and a test set to evaluate it.

	Classifier	Fold	TA	RMSE	Precision	Recall	F-Score	AUC
		1 st	98.51%	0.0548	0.980	0.980	0.9800	0.998
		2 nd	98.63%	0.0541	0.990	0.990	0.9900	0.999
		3 rd	97.96%	0.0482	0.976	0.976	0.9760	0.989
		4 th	98.63%	0.0543	0.990	0.990	0.9900	0.996
		5 th	98.98%	0.0578	0.989	0.989	0.9890	0.997
	OS-ELM	6 th	98.00%	0.0490	0.981	0.981	0.9810	0.995
		7 th	98.60%	0.0549	0.986	0.986	0.9860	0.999
		8 th	98.75%	0.0560	0.987	0.987	0.9870	0.999
		9 th	98.28%	0.0567	0.986	0.986	0.9860	0.999
		10 th	98.30%	0.0536	0.985	0.985	0.9850	0.999
		Avg	98.46%	0.0539	0.985	0.985	0.985	0.997


George Kyriakides
email: ge.kyriakides@uom.edu.gr

Konstantinos Margaritis

Comparison of Neural Network Optimizers for Relative Ranking Retention Between Neural Architectures

Keywords: Deep Learning, Neural Architecture Search, Ranking

Autonomous design and optimization of neural networks is gaining increasingly more attention from the research community. The main barrier is the computational resources required to conduct experimental and production project. Although most researchers focus on new design methodologies, the main computational cost remains the evaluation of candidate architectures. In this paper we investigate the feasibility of using reduced epoch training, by measuring the rank correlation coefficients between sets of optimizers, given a fixed number of training epochs. We discover ranking correlations of more than 0.75 and up to 0.964 between Adam with 50 training epochs, stochastic gradient descent with nesterov momentum with 10 training epochs and Adam with 20 training epochs. Moreover, we show the ability of genetic algorithms to find high-quality solutions of a function, by searching in a perturbed search space, given that certain correlation criteria are met.

Introduction

Optimization and autonomous design of neural architectures have been increasingly investigated in the past years [START_REF] Miikkulainen | Evolving Deep Neural Networks[END_REF][START_REF] Zoph | Neural Architecture Search with Reinforcement Learning[END_REF][START_REF] Zoph | Learning Transferable Architectures for Scalable Image Recognition[END_REF][START_REF] Kyriakides | Neural Architecture Search with Synchronous Advantage Actor-Critic Methods and Partial Training[END_REF][START_REF] Liu | Progressive Neural Architecture Search[END_REF][START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF][START_REF] Cai | ProxylessNAS: Direct Neural Architecture Search on Target Task and Hardware[END_REF][START_REF] Cortes | AdaNet: Adaptive Structural Learning of Artificial Neural Networks[END_REF]. Partly due to the popularity of deep learning [START_REF] Dong | Convolutional Neural Networks for Clothing Image Style Recognition[END_REF][START_REF] Lee | Deep Reinforcement Learning Agent for Playing 2D Shooting Games[END_REF][START_REF] Gatys | Image style transfer using convolutional neural networks[END_REF] and partly due to the cumbersome process of designing and optimizing the networks. It is only natural that we wish to free ourselves of the repetitively unique process of designing and optimizing. Many methods have been proposed.

Although most studies propose methods that produce better networks or converge faster, most of them are concerned with the way that the design or optimization algorithm acts. Nonetheless, the majority of computational resources are consumed during the candidate architectures' evaluation. Thus, many researchers distribute the workload among many machines in order to speed up the process. Although a possible and quite effective solution, not all algorithms favor distributed setups. Furthermore, it requires availability of additional hardware.

On the other hand, the networks are evaluated only to further the search or optimization method, until the final architecture is produced. As such, the actual performance of the network does not matter as an absolute value. It matters in relation to the other alternatives that the design algorithm has at each point in time. As such, methods that enable a quick relative evaluation of network architectures can greatly reduce the computational workload of designing and optimizing architectures.

In this paper, we study the use of various neural network optimizers under a small number of training epochs as a proxy of relative network performance. Furthermore, we conduct Monte-Carlo simulations of searching for the best parameters in 3-d spaces with similar correlation coefficients. The motivation behind this work is to examine the feasibility of using faster evaluation schemes (reduced number of epochs) when evaluating architectures in an autonomous design framework, where absolute performance has a smaller impact than relative performance between candidate architectures.

Related Work

As stated previously, most research focuses on faster and more efficient methodologies. Nonetheless, much of this work has greatly contributed to the field. A well-known approach entails the use of genetic algorithms in order to evolve networks of fixed submodule design [START_REF] Miikkulainen | Evolving Deep Neural Networks[END_REF]. Two versions are proposed, DeepNEAT which is an extension of the older NEAT [12] method for fully-connected networks and CoDeepNEAT which utilizes the co-evolution [START_REF] Moriarty | Forming Neural Networks Through Efficient and Adaptive Coevolution[END_REF] of graph structure and sub-module populations. Another approach, utilizes reinforcement learning in order to design CNN architectures [START_REF] Zoph | Neural Architecture Search with Reinforcement Learning[END_REF] and transferable architectures [START_REF] Zoph | Learning Transferable Architectures for Scalable Image Recognition[END_REF]. In [START_REF] Kyriakides | Neural Architecture Search with Synchronous Advantage Actor-Critic Methods and Partial Training[END_REF] an attempt to reduce the number of trainable parameters by freezing all network layers minus the output layer, results in a positive relative ranking. This indicates that freezing intermediate layers is a viable approach for neural architecture search. Another study attempts to design architectures directly on the target dataset, without the use of proxy tasks [START_REF] Cai | ProxylessNAS: Direct Neural Architecture Search on Target Task and Hardware[END_REF]. Although the authors argue that a proxyless search produces better results, it also demands greater computational resources. It is only logical that a search conducted on the target hardware and training conditions will produce the best results. Furthermore, in [START_REF] Miikkulainen | Evolving Deep Neural Networks[END_REF] researchers noticed that by restricting the available training epochs during evolution, the final produced architecture required less epochs to achieve performance comparable to state of the art.

Other approaches employ random search [START_REF] Bergstra | Random search for hyper-parameter optimization[END_REF], sequential model-based optimization [START_REF] Liu | Progressive Neural Architecture Search[END_REF] and boosting [START_REF] Cortes | AdaNet: Adaptive Structural Learning of Artificial Neural Networks[END_REF].

Methodology

Evolving Architectures

In this paper we use the methodology proposed in [START_REF] Miikkulainen | Evolving Deep Neural Networks[END_REF], called DeepNEAT to evolve neural architectures. Although the original method uses speciated crossover to produce offspring from a population as well as mutations, we use only the mutation operator, as our goal is to produce an array of different architectures of increasing complexity. The methodology initializes an individual architecture comprised only of an input and output node. At each mutation, there is a certain probability that a new connection will be added or an existing connection will be deactivated and replaced by a new module. Each module consists of a convolutional layer, a dropout layer, a weight scaling factor and a max pooling layer with probability of 0.5. Each gene in the chromosome contains information about the module's parameters: the number of filters in the convolutional layer, in the space [32, 256], the dropout rate in [0, 0.7], the weight scaling factor in [0, 2.0], the kernel size in {1, 3} and if the max pooling layer is present or not. In our experiment the probability of adding a new module node was set at 0.05 and the probability of adding a new connection at 0.1. Some of the architectures generated are depicted in Fig. 1. Number 1 is the starting architecture, consisting only of an input and an output layer. Each node denotes a whole module. Due to the fact that many node modules must be scaled in order for the network to be functional, Number 20 has 15 modules that translate to 86 PyTorch layers.

Experimental Setup

For our experiments, we generated several distinct neural networks of increasing complexity by mutating a simple network using DeepNeat for 20 generations, with the goal of classifying the CIFAR10 image recognition dataset. We trained each network using 1, 5, 10, 20, and 50 epochs respectively. We compared a total of seven optimizers; Adadelta, Adam, Adamax, RMSprop, Stochastic Gradient Descent (SGD), SGD with momentum (SGD-M), and SGD with nesterov momentum (SGD-NM). Following, we calculated each network's rank for its specific optimizer and epochs combination. Our goal was to compare the relative rankings of low training cost combinations to the relative rankings of fully trained networks (50 epochs).

In order to compare the networks' rankings, we use Kendall's rank correlation coefficient (τ) [START_REF] Kendall | A New Measure of Rank Correlation[END_REF], which is a measure of similarity between ranks (correlation of ranks). High positive or negative values imply correlation (positive and negative respectively).

To further investigate the ability of optimizers to discern high-performing architectures, we perform the analysis on all architectures as a whole, as well as the last ten generated architectures. The motivation behind this decision is that simple architectures will exhibit similar behavior and minor differences in performance can be attributed to other factors, such as initial conditions. Complex architectures on the other hand, demand more computational resources to train. Thus, design methods will benefit more from a faster but relatively consistent evaluation. Following, we conducted Monte-Carlo simulations of 3-d Rastrigin functions with similar correlation coefficients to those observed in our experiments in order to demonstrate the feasibility of using reduced training epochs as proxies for neural architecture search. We employed a simple genetic algorithm in order to find the local minima, both in the original function as well as the perturbed function. We initialized a population of 10 individuals, with crossover rate of 0.9 and mutation rate of 0.02. The population was selected for crossover by tournament selection. We experimented with 10, 100, and 1000 generations. Each experiment was repeated 1000 times, in order to obtain a reasonable distribution sample. Given the multi-modal nature of Rastrigin function, its search space resembles those in neural architectures.

Implementation

For our experimental implementation we used the NORD benchmarking framework. DeepNEAT was implemented in Python but only the mutation operator was used. Furthermore, we used ESA's Pygmo implementation of a simple genetic algorithm for our Monte-Carlo simulations [START_REF] Biscani | esa/pagmo2: pagmo 2.10[END_REF] The experiments were run on 7 Tesla K40 GPUs.

Results

Neural Architectures

In this study our focus is the retention of relative rankings while using reduced training epochs for candidate architectures. Nonetheless, the achieved accuracy of each architecture, optimizer and epoch combination provides us with an intuitive way to assess the quality of each network. Fig. 2 depicts the accuracy of each optimizer for 1 and 50 epochs of training. It is evident that with 50 epochs the stability is greatly increased. Nevertheless, it seems that there exists a pattern, that both groups follow. A sharp increase in accuracy at generation 10, a relatively stable performance for a small number of generations and then a sharp decrease. This pattern seems to be followed by all optimizers on both groups, except for the 'Ada' family (Adamax, Adadelta, Adam) when given 50 epochs of training. The training time required for each group is depicted in Table 1. As it is evident, the number of training epochs seems to be the most influential factor. Furthermore, required time seems to scale linearly with the number of epochs, as it is expected. Although useful, it is not easy to judge relative performance from accuracy plots and training time tables. In order to better compare the relative performance of the various optimizer and training epoch combinations, we ranked the generated architectures, depending on the accuracy that they achieved on the specific epoch-optimizer combination. We then computed the rank correlation coefficient for each list with the list generated by the same optimizer when the architectures are trained for 50 epochs before evaluation. The results show high correlation coefficients (above 0.65) for most combinations, except for RMSprop_01 and SGD_01 as depicted in Fig. 3 .

Fig. 3. Rank correlation coefficients heatmap.

The results seem promising, although they do not show a clear or extremely strong correlation between the groups. This can be partly attributed to the fact that simple architectures will perform relatively close, independently of the optimizer and the number of epochs trained. If the number of parameters is small and the problem domain large, there are simply too little degrees of freedom to adequately fit the model. Thus, we repeat the analysis, by discarding the first half of the generated architectures. The goal is to test if rank correlations are retained on more complex architectures. As it is evident in Fig. 4, there are now strong relationships between some groups. Adam_01, Adam_05, Adam_20, Adamax_10, Adamax_20, RMSprop_20, SGD-M_10, and SGD-NM_10 have high correlations with Adam (τ > 0.75, p < 0.05). Furthermore, SGD-NM_10 and Adam_20 have τ = 0.964 with p < 0.01, indicating a strong correlation in relative rankings. This confirms that as architecture complexity increases, evaluation schemes come to a more robust agreement about their relative performance.

Monte-Carlo Simulations

We saw the need to further study the ability to optimize through a proxy task with specific rank correlation coefficient, compared to the original search space. We also saw the need to study the behavior of optimization algorithms when they are given more freedom to explore both search spaces (i.e. they are not restricted by computational resources). In order to do so, a controlled simulation environment had to be created. The environment should have many local optima, similar to the search space of neural architectures. Furthermore, we wanted to search on a proxy space, were the proxy and original spaces had similar correlation coefficient to the one found in our original experiment. This would enable to study the feasibility of using any proxy task that is known to produce correlated results with the original search space. Thus, we generated a 3-d Rastrigin function and added uniform noise, until a fixed rank correlation coefficient was achieved, compared to the original function. We followed the same procedure as in our original experiment. The genetic algorithm searched the space of the perturbed function (proxy task) and the final solution was tested on the original. A 2-d example of the two functions can be seen in Fig. 5 showing how the addition of noise does change the function's surface, but leaves major features and patterns intact. Table 2 depicts the average (E) and Table 3 the standard deviation(σ) of the produced solutions, when applied to the original function. We conducted the experiment for three levels of correlation (τ=0.75, τ=0.85, and τ=0.95). In the original space, the genetic algorithm is able to produce solutions close to the top 1% by searching the space of the unperturbed Rastrigin function. For τ=0.75, by conducting the search on the perturbed spaces and then applying the final solutions on the normal function, the genetic algorithm is able to produce solutions near the top 4%. For τ=0.85 and τ=0.95 it is able to produce solutions near the top 3% and 2% respectively. By constructing a 95% confidence interval, the original function's space produces a solution that belongs in the top [0, 0.033]. For the worst-case scenario (τ=0.75), the perturbed function's space produces solutions in the top [0.016, 0.08]. For τ=0.85, the solutions lie in [0.006, 0.07] and for τ=0.95 they exist in [0.003, 0.039]. Furthermore, we can see that the solutions' standard deviation reduces as τ increases, given enough generations to search the given space. As it is evident, the most important factor concerning the final solution's quality is the number of generations that the algorithm was allowed to run.

Limitations

This study provides some interesting results, both in the real-world experiment, as well as the Monte-Carlo simulation. Nonetheless, more thorough research must be conducted in order to assess the feasibility of using reduced epoch training as a proxy for evaluating relative rankings in neural network architectures. The candidate architectures in this study are admittedly few, even though they are diverse. Moreover, we tested our results in a single image recognition task, although it utilized a popular and well-studied dataset. We hope that in our future work we will overcome some of these limitations, by running longer experiments on more recent hardware, thus being able to employ more architectures and datasets.

Conclusions

In this paper we studied the relative ranking correlations of neural architectures by evolving an initial architecture with DeepNEAT's mutation operator and training them for a set number of epochs, utilizing select optimizers. By comparing the 1, 5, 10, and 20 epoch groups to the 50-epoch group, we observed positive and relatively high correlations for most optimizers, except for the RMSprop_01 and SGD_01 optimizerepochs combination. In order to test the most complex architectures of our candidates, we split the candidates in two and tested the relative rankings of the second half. This revealed a strong correlation between Adam and the two optimizer-epochs groups, SGD-NM_10 and Adam_20. The correlation coefficient was very high for both groups (τ = 0.964) and showed high statistical significance (p < 0.01). In order to test the ability of optimization algorithms to produce optimal solutions from correlated proxy tasks, we employed a genetic algorithm on a perturbed Rastrigin function, with a rank correlation coefficient of 0.75 to the original. The genetic algorithm produced solutions that belong on average, to the top 3,5% of the original function's solutions.

From the above, we conclude that there is a positive probability that reduced epoch training can be used in order to evaluate the relative rankings of neural architectures for design and optimization purposes. It is evident that the design/optimization process does not care for the performance of intermediate solutions, only for their relative performance. Furthermore, genetic algorithms perform relatively well when searching in perturbed spaces, given that the correlation with the original space is at least 0.75. Nonetheless, more thorough research must be conducted in order to establish optimal proxies for neural architecture design and optimization in a variety of domains.

Finally, the most decisive factor concerning the algorithm's solution quality is the number of generations that it is allowed to run. Thus, it can be advantageous to use proxy tasks with lower correlation to the original search space, if the reduction in computational resources required to evaluate a single network can be dedicated to running a more exhaustive search.

Future Work

In light of our recent findings, we hope to continue the research into feasible solutions for the relative comparison of neural network architectures. A higher number of candidate architectures and a broader range of datasets will greatly enhance our confidence in the results, given that they agree with the current. Moreover, comparisons with even higher number of training epochs will further enhance our confidence in the results. Finally, a study that employs well-established design methodologies will contribute to the final verdict on the feasibility of using reduced epoch training, by collecting data from the actual problem.

Fig. 1 .

 1 Fig. 1. Sample architectures evolved.

Fig. 2 .

 2 Fig. 2. Comparison of achieved accuracy in 1 epoch of training to 50 epochs.

Fig. 4 .

 4 Fig. 4. Rank correlation coefficients heatmap for the second half of generated architectures.

Fig. 5 .

 5 Fig. 5. 2-D Rastrigin functions, original and perturbed with τ = 0.738

Table 1 .

 1 Average training times (in seconds per architecture).

	Optimizer	1 Epoch	5 Epochs	10 Epochs	20 Epochs	50 Epochs
	Adadelta	36.1	158.6	312.3	622.1	1701.5
	Adam	35.9	157.5	311.7	617.4	1724.8
	Adamax	35.8	157.1	310.3	617.8	1718.6
	RMSprop	35.4	156.5	310.2	614.6	1687.6
	SGD	36.0	157.1	308.3	613.2	1687.3
	SGD-M	35.7	157.2	309.3	613.5	1709.5
	SGD-NM	35.5	156.6	311.2	617.7	1707.4

Table 2 .

 2 Mean of achieved solution's top percentage for the original and perturbed search space, when applied to the original function.

	Generations	E(original)	E(τ=0.95)	E(τ=0.85)	E(τ=0.75)
	10	0.160	0.162	0.184	0.195
	100	0.050	0.058	0.086	0.103
	1000	0.013	0.022	0.032	0.040

Table 3 .

 3 Standard deviation of achieved solution's top percentage for the original and perturbed search space, when applied to the original function.

	Generations	σ(original)	σ(τ=0.95)	σ(τ=0.85)	σ(τ=0.75)
	10	0.075	0.066	0.075	0.084
	100	0.031	0.034	0.045	0.052
	1000	0.010	0.012	0.019	0.021

Acknowledgements

This work was supported by computational time granted from the Greek Research & Technology Network (GRNET) in the National HPC facility -ARIS -under project ID DNAD. Furthermore, this research is funded by the University of Macedonia Research Committee as part of the "Principal Research 2019" funding program.