Sheng-Hui Wang

Huai-Ting Li

En-Jui Chang

Andy Wu
email: andywu@access.ee.ntu.edu.tw

Entropy-Assisted Emotion Recognition of Valence and Arousal Using XGBoost Classifier

Keywords:

Emotion recognition is an essential function to realize humanmachine interaction devices. Physiological signals which can be collected easily and continuously by wearable sensors are good inputs for emotion analysis. How to effectively process physiological signals, extract critical features, and choose machine learning model for emotion classification has been a big challenge. In this paper, an entropy-based processing scheme for emotion recognition framework is proposed, which includes entropy domain feature extraction and prediction by XGBoost classifier. We experiment on AMIGOS database and the experimental results show that the proposed scheme for multi-modal analysis outperforms conventional processing approaches. It achieves approximately 80% and 68% accuracy of prediction for two affect dimensions, valence and arousal. For one modality case, we found that galvanic skin response (GSR) channel is the most potential modality for prediction, which leads to best performances.

Introduction

Affective computing is a key technology for human-computer interaction (HCI) in the future era of Internet of thing (IoT) [START_REF] Picard | Affective Computing for HCI[END_REF], which makes it possible for machines and computers to realize human's emotion and mentality in real time. Furthermore, it can give appropriate responses and services based-on human's current mental status. Related databases have been developed for related researches on affective computing. Among the open-source databases available on the Internet, FEEDTUM [START_REF] Wallhoff | Facial expressions and emotion database[END_REF] and Berlin Database of Emotional Speech (EmoDB) [START_REF] Burkhardt | A Database of German Emotional Speech[END_REF], released in 2005, used facial video and speech for emotion recognition, respectively. They predicted basic emotion of people, such as anger and happy, by single type of signal.

As intelligent IoT develops, more and more wearable devices are equipped with different kinds of sensors. We can effectively get various signals sensed from subjects and have the access for continuously monitoring. For this scenario, physiological signals can be good inputs for affective computing framework, since we do not express our emotions on faces and make sounds every time. Integrating multiple-channel of physiological signals for emotion recognition and aggregating machine learning classifiers for ensemble learning have gain attention in recent years. DEAP database [START_REF] Koelstra | DEAP: A Database for Emotion Analysis; Using Physiological Signals[END_REF] proposed in 2012 has been experimented in many works [START_REF] Chen | Identifying valence and arousal levels via connectivity between EEG channels[END_REF]- [START_REF] Wu | Estimation of valence of emotion using two frontal EEG channels[END_REF], it contains several types of physiological signals which were recorded when subjects were watching music videos. In DEAP, emotions perceived by subjects are self-assessed on two affect dimensions, valence and arousal. Valence tells how positive or negative the emotion is and arousal indicates the intensity of emotion. In 2016, ACERTAIN database [START_REF] Subramanian | ASCERTAIN: Emotion and Personality Recognition using Commercial Sensors[END_REF] was proposed that first used commercial wearable sensor to collect data and added personality into experiments. AMIGOS database [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF] proposed in 2017 contains stable data, several types of videos, and elements for other related research, such as social context and mood. These all indicate that multi-modality analysis of physiological signals is a trend of affective computing. Therefore, how to effectively process different physiological signals and further extract critical features from different channels for emotion classification have been a big challenge.

The aim of this paper is to establish reliable emotion recognition framework based on physiological signals. To improve the accuracy of emotion recognition, we firstly extract the entropy domain features to quantify the regularity and randomness of signal, which has high potential to represent the different levels of emotion. Next, we apply the XGBoost [START_REF] Chen | XGBoost[END_REF] to enhance the performance of classification which has high scalability and efficiency for training and has the ability to learn from high dimension data without using high complexity feature selection algorithms.

The rest of the paper is organized as follows. Section 2 describes common methods for emotion recognition task. Section 3 illustrates the enhanced methods for feature extraction and machine learning engine. Section 4 shows the experimental results and Section 5 concludes.

Common Approaches of Emotion Recognition Framework

The flow chart of the emotion recognition framework in this study is shown in Fig. 1. Three kinds of physiological signals are inputs and the output is the prediction for high or low of affect dimensions. Three other blocks including pre-processing, feature extraction, and machine learning engine will be introduced as below.

Pre-Processing

We use the processed data provided by the official website of AMIGOS dataset. In AMIGOS database, each subject watched both long and short videos as stimulus for emotion. Because long-length video data contain multiple values for valence In this work, we do binary classification for valence and arousal. Selfassessments of each subject on two affect dimensions provided in AMIGOS database is used for labels. They are originally values between [START_REF] Picard | Affective Computing for HCI[END_REF][START_REF] Chen | XGBoost[END_REF] and needed to be transformed to either positive class or negative class label, which represents high and low, respectively. Fig. 2 shows two different ways to process. Directly cut by threshold 5 is intuitive, but hard to distinguish one's relative high and low emotion. In this way, values may highly accord to personal tendency to rate high or low scores, which results in imbalance of labels and decrease in performances. Our method is to define labels by subject-dependent mean value among all 16 videos of each person. Using mean value instead of median value can avoid ambiguity whether the median is labeled positive class or negative class.

Feature Extraction

Feature extraction is implemented on all three kinds of signals with different ways to reduce the input dimension before using machine learning methods. Features illustrated in this Section are slightly different from the feature set in AMIGOS database. Features extracted among physiological signals are shown in Table 1. Features from raw ECG (60): spectral power 0-6 Hz (per 0.1Hz) HRV features (10): mean, standard deviation, skewness, kurtosis, root mean square, root mean square of the successive differences, percentage of points higher than mean + standard deviation, percentage of points lower than mean -standard deviation, low frequency spectral power (0.01-0.08 Hz), middle frequency spectral power (0.08-0.15 Hz), high frequency spectral power (0.15-0.50 Hz). HR features (6): mean, standard deviation, skewness, kurtosis, percentage of points higher than mean + standard deviation, percentage of points lower than mean -standard deviation EEG signal (105 features) PSD features (70): average 5 bands (theta, slow alpha, alpha, beta, and gamma) PSD in 14 channels.

Asymmetry features (35): difference of 5 bands PSD between 7 pairs of channel.

GSR signal: Partly follow methods in [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF] and [START_REF] Kim | Emotion recognition based on physiological changes in music listening[END_REF], skin response (SR) related features [START_REF] Tsai | Empirical mode decomposition based detrended sample entropy in electroencephalography for Alzheimer's disease[END_REF] are first extracted and skin conductance (SC) signal is got by computing SR signal's reciprocal. SC signal is then normalized and used to extract features (4). Next, skin conductance slow response (SCSR) and skin conductance very slow response (SCVSR) is got by low-pass filtering normalized SC signal at 0.2 and 0.08 Hz, respectively. Afterwards, some SCSR related features are extracted (4). Last, after de-trending SCSR and SCVSR by using empirical mode decomposition (EMD) methods [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] to remove the last half intrinsic mode functions (IMFs), some related features are extracted.

ECG signal: Partly follow methods in [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF] and [START_REF] Kim | Emotion recognition based on physiological changes in music listening[END_REF], 0-6 Hz spectral power features are first extracted (60, per 0.1). Then, R-R interval (RRI) series is calculated after detecting R peaks in ECG. Heart rate variability (HRV) and heart rare (HR) time series can be computed by using RRI. Finally, HRV related features [START_REF] Huang | The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] and HR related features (6) are extracted.

EEG signal: Follow the methods in [START_REF] Koelstra | DEAP: A Database for Emotion Analysis; Using Physiological Signals[END_REF], average power spectral density (PSD) of theta band (4-7 Hz), slow alpha band (8-10 Hz), alpha band (8-13 Hz), beta band (14-30 Hz), gamma band (31-47 Hz) of each EEG channel are extracted (5 × 14 = 70). Also, asymmetry of PSD of 5 bands between 7 pairs of EEG channel are extracted (5 × 7 = 35).

Machine Learning Engine

This block consists of a feature selection method and a machine learning classifier, as shown in Fig. 3a. Feature space which contains 214 features is comparatively big while we only have 528 (16 subjects * 33 videos) data in AMIGOS database. Thus, feature selection is applied to eliminate redundant features and reduce the model complexity of machine learning algorithm in next stage. This is helpful for enhancing overall performance and sparing lots of computation resources. In Fig. 3b, overview of machine engine block in AMIGOS database is shown. For feature selection method in original paper of Amigos database [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF], Fisher's linear discriminant (FLD) J [START_REF] Song | Feature Selection Based on Linear Discriminant Analysis[END_REF] is calculated over each feature, which is defined as (

Afterward, user can decide how much discriminant features they would like to select based on validation sets and pick up features with the highest J values. For machine learning classifier, Gaussian Naïve Bayes (GaussianNB) classifier is used. Assuming features are independent and Gaussian distributed, GaussianNB is given as

1 1 ( ,... , ) = ( ) ( ) n n c i i i G f f argmax p C c p F f | C c      , (2)
and calculate the probability of each sample for each class, where F is feature set and C is the classes set. These two methods are not powerful and lead to relatively low f1-score shown in [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF].

In this work, we first implement support vector machine (SVM) [START_REF] Cristianini | An Introduction to Support Vector Machines and Other Kernel-based Learning Methods[END_REF], [START_REF] Luts | A tutorial on support vector machine-based methods for classification problems in chemometrics[END_REF] based approach as the benchmark, as shown in Fig. 3c. SVM is a famous tool in the field of machine leaning and it has good performance on classification and prediction over many applications. Compared with other classifier, SVM can deal with small-size datasets well, since it uses only support vectors to construct hyperplane. When it comes to feature selection, recursive feature elimination (RFE) algorithm [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF] can efficiently remove irrelevant features and be suit for SVM. RFE-SVM eliminates features with smallest weight in SVM model on a sequential backward selection (SBS) based process. The progress continues until one feature is left and the features combined with SVM model with the highest performance is outputted.

Proposed Framework with Entropy Domain Features and XGBoost Classifier

Refine Composite Multiscale Entropy (RCMSE) [16]:

Multiscale entropy (MSE) [START_REF] Costa | Multiscale entropy analysis of biological signals[END_REF] has been used widely to evaluate physiological control mechanisms, such as atrial fibrillation and Alzheimer's disease [START_REF] Tsai | Empirical mode decomposition based detrended sample entropy in electroencephalography for Alzheimer's disease[END_REF]. RCMSE is an improved version of MSE. It reduces the possibility of undefined value problem when the signal length is short and has a better accuracy for entropy estimation. The concept of RCMSE is using different scales of local matching pattern to compute the regularity of the signal. There are two steps of RCMSE. First step: time series of signal is coarse-grained into multiscale series. For each scale factor  ,  series are generated by average  points in non-overlapping windows and each of series overlap  -1 points with neighbor series. The j-th point of k-th coarse-grained series,

  ( ) ( ) ( ) ( ) ( ) 1 2 3 ... k k , k , k , k , p y y y y       y
,with scale factor  of signal x is defined:

1 () ( 1) 1 1 1 jk k , j i i j k N x , j , k                y . (3)
In the conventional MSE algorithm, output of coarse graining is only first series of each scale factor. Second step: averaged sample entropy of each is calculated as below:

1 1 1 ( ( ) m k, k m k, k ) n RCMSE , ,m,r ln n            x , ( 4
)
where m is the matching pattern length and r is the similarity criterion. In (4), m k, n  is the number of two sets of simultaneous data points of length m in kth series have the difference < r. When the ratio of

1 1 m k, k n      and 1 m k, k n    
is small, sample entropy would be large, which means high complexity of signal. Value of RCMSE is undefined only when

1 1 m k, k n      = 0, which means m k,
n  of k series are all zero. By summation of matching pattern of multiple series, RCMSE has less probability to be undefined than conventional MSE.

In this work, RCMSE is applied to RRI series of ECG signals. Because the length of ECG signals are relatively short (each video 55s-155s), the scale factor is set up to be 3. In addition, we set the matching pattern m up to be three and the similarity criterion r is set to be 0.2 of standard deviation.

Turning Points Ratio (TPR):

TPR is proposed on the basis of nonparametric "Runs Test" to evaluate the randomness in a time-series [START_REF] Wallis | A Significance Test for Time Series Analysis[END_REF] and the idea was used in RRI of ECG signal [START_REF] Dash | A statistical approach for accurate detection of atrial fibrillation and flutter[END_REF]. The concept of TPR is to measure the complexity of the signal by number of turning points compared to total points. Turning point is found by comparing each point with left and right neighbor points and TPR is calculated as follows:

  1 11 2 ( - )( - ) > 0 2 N i i i i i x x x x TPR N       , ( 5
)
where N is the length of the signal x. Besides original TPR, we extracted modified TPR (MTPR) on signals. The procedure of MTPR is to use EMD methods to extract trend of the signal first, in order to remove trivial peaks of the signal. Next, compute TPR as (3) on the extracted trend. In this work, TPR and MTPR are calculated on RRI of ECG signals and GSR signals. EEG signal has less information on time series and is skipped. Shannon Entropy: Shannon entropy is commonly calculated in the domain of information theory. It is defined for a given discrete probability distribution, using probability of each symbol to measure the uncertainty or randomness of the data.

However, almost every point has different values. We can't see each of different points as a new symbol, otherwise almost the same value of Shannon entropy we would get from different signals. In other way, we classify every points into one group of the group set. First, outliers that have larger differences than three standard deviations with mean are removed. Second, we sort rest data points into N groups, which are equally divided between max and min value. Last, i p of each ith group is calculated and Shannon entropy we can get by: 1 ( )

N ii i Shannon Entropy p log p    . ( 6
)
The concept of Shannon entropy is to observe the complexity of the signal by overall distribution. Each group is seen as a symbol.

In this work, Shannon entropy is calculated on RRI of ECG signals and GSR signals. Since the optimal group number is data dependent, we apply total 4, 8, 16, 32, 64 groups for simulation.

For these three enhanced feature extraction methods, total 26 features are added into original feature space. The details are listed in Table 2. 12): m = 0:  = {1, 2, 3}, sum of  = {1, 2, 3} m = 1:  = {1, 2, 3}, sum of  = {1, 2, 3} m = 2:  = {1, 2, 3}, sum of  = {1, 2, 3} TPR (2): TPR, MTPR. Shannon Entropy (5): number of groups: 4, 8, 16, 32, 64

Extreme Gradient Boosting (XGBoost)

In this part, we would like to change the content of the machine learning blocks mentioned before in Section 2.3, where RFE based SVM approach was used. Among all the machine learning algorithms, gradient boosting tree based model [START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF] has shown in many applications in different domains. XGBoost [START_REF] Chen | XGBoost[END_REF] is an efficient and scalable gradient boosting machine, which has won lots of machine competitions in recent years [START_REF] Adam-Bourdarios | The Higgs Machine Learning Challenge[END_REF], [START_REF] Phoboo | Machine Learning wins the Higgs Challenge[END_REF]. It is an ensemble model consisting of sets of classification and regression tree (CART). While XGB is used for supervised learning problems and we use training data i x to predict a target variable i y , the model can be described in the form:

1 ( ) , K i k i k k ŷ f x f F    , ( 7
)
where K is the total number of trees, fk for kth tree is a function in the functional space F, and F is the set of all possible CARTs. In the training, each of newtrained CART will try to complement the so-far residual. Objective function been optimized at (t+1)th CART is described:

() 11
( ) + ( ) nt t i i i ii ôbj l y , y f    , (8)
Where l( ) denotes the training loss function, i y the is ground truth, and () t i ŷ is the prediction value at step t.

( )

 given by 2 1 1 ( ) = 2 T j j f T w     (9)
is the regularization term, where T are the number of leaves and j w is the score on jth leaf. When (9) is optimized, Taylor's expansion is used so that gradient descent can be used for different loss functions. Furthermore, feature selection is no need when we use XGBoost approach. During training period of XGBoost, good features would be chosen as node in trees, which means features not used are abandoned.

In this work, we use the scikit-learn API for XGBoost classification. The inputs of XGBoost are total 240 features (214 traditional features + 26 entropy domain feature) and the outputs are prediction results for valence or arousal. For loss function in [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF], logistic loss function is set. The details of used and performances are shown in next Section.

Experimental Settings and Results

Single trial classification for two affect dimensions, arousal and valence, is experimented and the flow is shown in Fig. 1. After signal is pre-processed and transformed into features, all features are normalized to [1 -1]. Leave-one-subjectone approach is conducted to evaluate performance. That is, every time one of the 33 subjects is leaved as test set and machine learning engine is trained using remaining 32 subjects' features. When the procedure is repeated 33 times, final performance is calculated by averaging 33 values in each procedure.

There are 3 processing schemes that are compared.  Scheme_1: Using the feature set and the method described in original AMIGOS database [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF], where FLD and GaussianNB were used for machine learning engine block.



Scheme_2: Using the feature set illustrated in Section 2.2 and RFE-SVM for machine learning engine block.  Scheme_3: Using new feature set with entropy domain features and XGBoost for machine learning engine block. Table 3 shows the comparison of f1-score and accuracy over three schemes, F1-score is the harmonic mean of precision and recall and here we average f1scores of positive and negative class as final value. Four scenarios were experimented: using only GSR, ECG, EEG, and using all modalities. Table 4 shows the features used in trees of XGBoost classifier for classification of valence and arousal. , HRV_ root mean square, HRV_ low frequency spectral power, HRV_ high frequency spectral power, HRV_percentage of points lower than mean -standard deviation, TPR, Shannon group4, RCMSE_m=2 sum EEG: T8_gamma, O2_gamma, T7_gamma, O2_beta, T7-T8_slow alpha, T7-T8_beta, , T7-T8_theta, O1-O2_beta Arousal GSR: spectral power 0-0.2 1.6-1.8, SC_mean of derivative, SC_mean of second derivative, SCSR_number of local minimum ECG: spectral power 1.8-1.9 1.9-2.0 3.0-3.1, HRV_percentage of points lower than mean -standard deviation, HRV_ low frequency spectral power, RCMSE_m=1, τ=1 EEG: FC6_alpha F1-score for Scheme_1 is directly obtained from [START_REF] Miranda-Correa | AMIGOS: A Dataset for Mood, Personality and Affect Research on Individuals and Groups[END_REF]. For dimension valence and arousal, best accuracy is about 80% and 68%, respectively. Scheme_3 outperforms Scheme_1 and Scheme_2 in almost every scenario on both affect dimensions, especially in valence dimension, outperforms by more than 10% of accuracy. Using all-modalities for prediction can have highest performances except for f1-score in arousal. In our experiments (Scheme_2 and Scheme_3), using GSR channel to predict can have the best performance for both valence and arousal dimensions, when only one channel is used. However, the features extracted in GSR channel is the least among three modalities, which indicates that the quality but not the quantity of features lead to better prediction ability.

Conclusion

In this paper, we proposed the scheme including both entropy domain features and XGBoost. We enhance feature extraction methods which are entropy domain and helpful for evaluating the complexity of physiological signals. On the other hand, XGBoost classifier which gains popularity in recent years is used for learning and prediction. The proposed scheme can reach the performance of approximately 80% and 68% accuracy on valence and arousal dimension, respectively. It outperforms the processing scheme in original paper of AMIGOS database and the scheme that contains common features and traditional SVM model.

Fig. 1 .

 1 Fig. 1. Flow chart of general emotion recognition framework. and arousal, only short-length video data are used in this work to reduce the uncertainty. There are total 40 subjects and 16 short-length videos in database. Each video of each subject is seen as one data. Each of data has 14 channels for Electroencephalography (EEG), 2 channels for Electrocardiography (ECG), and 1 channel for galvanic skin response (GSR). 7 subjects whose physiological signals of videos contain some missing values are removed (ID number: 9, 12, 21, 22, 23, 24, 33). Therefore, there are totally (40-7) (subjects)× 16 (videos)= 528 data. We further checked the data stability. First Channel of ECG signals suffer from fewer noises compared to second channel and are used for further analyzing. In addition, GSR signals are filtered by the low-pass filter with cut-off frequency 3Hz to remove abnormal high frequency noises.In this work, we do binary classification for valence and arousal. Selfassessments of each subject on two affect dimensions provided in AMIGOS database is used for labels. They are originally values between[START_REF] Picard | Affective Computing for HCI[END_REF][START_REF] Chen | XGBoost[END_REF] and needed to be transformed to either positive class or negative class label, which represents high and low, respectively. Fig.2shows two different ways to process. Directly cut by threshold 5 is intuitive, but hard to distinguish one's relative high and low emotion. In this way, values may highly accord to personal tendency to rate high or low scores, which results in imbalance of labels and decrease in performances. Our method is to define labels by subject-dependent mean value among all 16 videos of each person. Using mean value instead of median value can avoid ambiguity whether the median is labeled positive class or negative class.

Fig. 2 .

 2 Fig. 2. (a) Define labels by comparing with threshold score 5. (b) Define labels by comparing with subject-dependent mean among all videos.

Fig. 3 .

 3 Fig. 3. (a) Overview of machine learning engine block. (b) Machine learning engine block in AMIGOS paper. (c) Machine learning engine block of simple benchmark.

  entropy domain features, such as sample entropy, permutation entropy  based features, are widely used on physiological signals. Extracted entropy value can help quantify the regularity of signal and thus be applied on medical diagnosis. For this emotion recognition work, three types of entropy domain features, including refined composite multiscale entropy (RCMSE), turning point ratio (TRP), and Shannon entropy, are applied to measure the complexity of physiological signals. The details are shown as follows.

Table 1 .

 1 Features extracted among physiological signals mean, mean of derivative, mean of negative derivative, proportion of negatives in derivative, number of local minimum, mean of each rise time series, spectral power 0-2.4 Hz (12, per 0.2Hz). SC features (4): mean, standard deviation, mean of derivative, mean of second derivative. SCSR features (4): mean, standard deviation, mean of derivative, mean of second derivative. SCSR & SCVSR features (6): zero crossing rate of SCSR, zero crossing rate of SCVSR, average magnitude of SCSR occurrences, average magnitude of SCVSR occurrences, average magnitude of all SCSR and SCVSR occurrences, ratio of number of SCSR occurrences to SCVSR.

	Modality	Features Extracted
	GSR signal (32 features) SR features (18): ECG signal
	(77 features)

Table 2 .

 2 Enhanced features extracted among physiological signals

	Modality	Features Extracted
	GSR signal	TPR (2): TPR, MTPR.
	(7 features)	Shannon Entropy (5): number of groups: 4, 8, 16, 32, 64
	ECG signal	RCMSE (
	(19 features)

Table 3 .

 3 Performance of scheme_1-scheme_3 on emotion recognition framework. (F1score is mean for positive and negative class. Red value indicates the highest accuracy or f1-score in each affect dimension. Scheme_1: Amigos feature set + FLD + GaussianNB; Scheme_2: commom feature set + RFE-SVM; Scheme_3: new feature set + XGBoost)

			Accuracy		F1-score
		GSR ECG EEG	All	GSR ECG EEG	All
	V Scheme_1	-	-	-	-	0.531 0.535 0.576 0.570
	Scheme_2 0.640 0.615 0.589 0.680 0.610 0.588 0.556 0.666
	Scheme_3 0.776 0.634 0.581 0.801 0.774 0.630 0.577 0.800
	A Scheme_1	-	-	-	-	0.548 0.550 0.592 0.585
	Scheme_2 0.644 0.623 0.566 0.663 0.643 0.612 0.557 0.660
	Scheme_3 0.682 0.542 0.579 0.684 0.707 0.566 0.604 0.698

Table 4 .

 4 Features selected in XGBoost classifier for classification of two affect dimensions. (Entropy domain features are marked in red color)

	Dimension	Features
	Valence	GSR: spectral power 0-0.2 0.4-0.6 0.6-0.8 0.8-1.0 1.4-1.6, SC_mean of
		derivative, SCSR_mean, SR_MTPR
		ECG: spectral power 1.2-1.3 2.5-2.6

Acknowledgements. This work was supported by the Ministry of Science and Technology of Taiwan (MOST 106-2221-E-002-205-MY3 and MOST 106-2622-8-002-013-TA), National Taiwan University and Pixart Imaging Inc.


Sotirios-Filippos Tsarouchis

Maria Th Kotouza
email: maria.kotouza@issel.ee.auth.gr

Fotis E Psomo- Poulos

Pericles A Mitkas
email: mitkas@auth.gr

A multi-metric algorithm for hierarchical clustering of same-length protein sequences

Keywords: Hierarchical clustering, amino acid sequences, sequence similarity, sequence identity

The identification of meaningful groups of proteins has always been a major area of interest for structural and functional genomics. Successful protein clustering can lead to significant insight, assisting in both tracing the evolutionary history of the respective molecules as well as in identifying potential functions and interactions of novel sequences. Here we propose a clustering algorithm for same-length sequences, which allows the construction of subset hierarchy and facilitates the identification of the underlying patterns for any given subset. The proposed method utilizes the metrics of sequence identity and amino-acid similarity simultaneously as direct measures. The algorithm was applied on a realworld dataset consisting of clonotypic immunoglobulin (IG) sequences from Chronic lymphocytic leukemia (CLL) patients, showing promising results.

Introduction

One of the main challenges in computational biology concerns the extraction of useful information from biological data. This requires the development of tools and methods that are capable of uncovering trends, identifying patterns, forming models, and obtaining predictions of the system [START_REF] Larranaga | Machine learning in bioinformatics[END_REF]. The majority of such tools and methods exist within the field of data mining. Clustering [START_REF] Berkhin | A survey of clustering data mining techniques[END_REF] is a data mining task that divides data into several groups using similarity measures, such that the objects within a cluster are highly similar to each other and dissimilar to the objects belonging to other clusters based on that metric. From a machine learning perspective, the search for meaningful clusters is defined as unsupervised learning due to the lack of prior knowledge on the number of clusters and their labels. However, clustering is a widely used exploratory tool for analyzing large datasets and has been applied extensively in numerous biological, genomics, proteomics, and various other omics methodologies [START_REF] Belacel | CLUSTERING: UNSUPERVISED LEARNING IN LARGE SCREENING BIOLOGICAL DATA[END_REF]. Genomics is one of the most important domains in bioinformatics, whereas the number of sequences available is increasing exponentially [START_REF] Larranaga | Machine learning in bioinformatics[END_REF]. Often the first step in sequence analysis, clustering can help organize sequences into homologous and functionally similar groups, can improve the speed of data processing and analysis, and can assist the prediction process.

There have been several approaches in the past, attempting to address the issue of identifying meaningful groupings of sequences of identical lengths. While sequence clustering has a long history in the field of bioinformatics ( [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences[END_REF], [START_REF] Edgar | Search and clustering orders of magnitude faster than BLAST[END_REF]), there are few attempts in literature that can be successfully applied to sequences of the same length. One of the most notable approaches is the Teiresias algorithm [START_REF] Isodore Rigoutsos | Combinatorial pattern discovery in biological sequences: The TEIRESIAS algorithm[END_REF], [START_REF] Darzentas | A dif-ferent ontogenesis for chronic lymphocytic leukemia cases carrying stereotyped antigen re-ceptors: molecular and computational evidence[END_REF], that discovers rigid patterns (motifs) in biological sequences based on the observation that if a pattern spans many positions and appears exactly k times in the input, then all fragments (sub patterns) of the pattern have to appear at least k times in the input. The main drawback of this algorithm is that pairwise comparisons are employed between all the sequences of the dataset, leading to an exponential increase in execution time and memory requirements for large-scale datasets.

In this paper, we introduce a method for clustering amino acid sequences of identical length, using an approach that does not demand pairwise comparisons between the sequences, but it is instead based on the usage of a matrix that contains the amino acid frequencies for each position of the target sequences.

Methodology

The proposed clustering method uses both sequence identity and amino-acid similarity as similarity measures to form the clusters (both concepts are further defined below).

Ultimately, a binary top-down tree is constructed by consecutively dividing the frequency amino acid matrix of a given cluster into two sub-matrices, until only two sequences remain at each cluster at the leaf-level.

Binary Tree construction

The first phase consists of a top down hierarchical clustering method. Hierarchical clustering is one of the most commonly used approaches for sequence clustering [START_REF] Cai | ESPRIT-Forest: Paral-lel clustering of massive amplicon sequence data in subquadratic time[END_REF].

At the beginning of the process, it is assumed that all Ν sequences belong to a single cluster, which is consequently split recursively while moving along the different levels of the tree. Ultimately, the constructed output of the clustering process is presented as a binary tree. The right side of the tree is expected to be much longer than the left side due to the constraints posed by the split process; the sequences with the highest similarity percentage at a specific sequence position are assigned to the right side, whereas the remaining sequences are assigned to the left. The process of this phase (Algorithm 1, Fig. 2) can be formally described in the following steps and further detailed below:

1. Create frequency and frequency-similarity based matrix (FM, FSM) 2. Compute average identity of the matrices (𝑖𝑑 ̅̅̅̅ , 𝑖𝑑𝑆 ̅̅̅̅̅ ) 3. Split each frequency matrix into two sub matrices 4. Update the Level matrix and the Identity matrices (Y, I, IS) 5. Check for branch break

Step 1: Frequency amino acid and frequency-similarity based amino acid matrix. The first step of the proposed method aims to construct a frequency amino acid matrix. This is defined as a 2-dimensional matrix, with number of rows equal to the number of the different amino acids (i.e. 20 rows) and number of columns equal to the length (L) of the sequences provided as input. Each element (i,j) of the matrix corresponds to the number of times amino acid i is present in position j for all sequences. The count matrix (CM) contains the absolute values, whereas the frequency matrix (FM) contains the corresponding frequencies (Eq. 1). In addition to CM, a second frequency matrix is constructed using the same approach, but instead of the 20 amino acids, groups of similar amino acids are used under given schemes. As a use case in this paper, the 11 IMGT physicochemical classes are taken into consideration, as shown in Fig. 1, and an 11 x L frequency matrix is constructed. Step 2: Compute Identity. The identity is a similarity metric that is computed for each cluster based on the corresponding amino-acid frequency matrix. This metric, calculated as a percentage, indicates how compact the cluster is. Its maximum value (100%) corresponds to the case when the number of unique sequences that belong to a cluster is equal to one. The overall identity is equal to the average identity of each position in the given sequence (Eq. 2), and it is produced based on the CThr matrix, that contains only the elements of the amino acid matrix (CM) that correspond to amino acids that appear more than once in the corresponding column (Eq. 3).

FM = CM / Ν (1)
𝑖𝑑 ̅̅̅̅ = (∑ 𝑚𝑎𝑥(𝐶𝑇ℎ𝑟[, 𝑗]) /𝛮) 𝐿 𝑗=1 /L (2) CThr[i,j] = { 𝐶𝑀[𝑖, 𝑗], 𝑖𝑓 𝐶𝑀[𝑖, 𝑗] > 1 0, 𝑖𝑓 𝐶𝑀[𝑖, 𝑗] ≤ 1
, i=1,..,20, j=1,..,L

Steps 3-4: Split of frequency matrices. Each cluster is divided into two distinct subsets according to the following criteria and in the order that they are listed:

1. Select the element of the frequency amino acid matrix with the highest percentage. An example of the division of cluster0 using this criterion is shown in Fig. 3. 2. If there exist more than one columns of the frequency amino acid matrix that contain the same highest percentage value, the selection is applied using the entropy criterion, defined below. 3. In the case where more than one columns exhibit the exact same entropy value, criterion 1 is applied to the frequency similarity amino acid matrix. 4. In the case of non-unique columns, criterion 2 is applied to the frequency similarity amino acid matrix. 5. If the number of columns is still more than one, one column from the above sub group of columns is randomly selected.

Entropy criterion. The entropy is computed for each column of the frequency matrix and represents the diversity of the column. A lower entropy value indicates a more homogeneous column, therefore the column with the lowest entropy value is selected during the splitting process. Step 5: Branch breaking. A cluster is no further divided into subsets when the number of sequences that belong to the cluster is less than 3.

Fig. 3. The division of Cluster 0 using the first criterion.

Software Implementation

The algorithm outlined in the paper is implemented in R, a programming language that is widely used for statistical computing, graphics and data analytics. In order to produce an interactive and user-friendly tool, thus making it easier for users to interact with the data, the analysis, and the visualization of the results, an R Shiny application was built using the shiny package. In practical terms, a Shiny App is a web page / UI connected to a computer running a live R session (Server). The users can select personalized parameters via the UI. These parameters are passed on to the Server, where the actual calculations are performed and the UI's display is updated according to the produced results.

R Shiny applications consist of at least two R scripts; the first one implements the User Interface (ui.R) by controlling the layout of the page using html commands and other nested R functions, and handles the input parameters provided by the users. The second script implements the Server (server.R) and contains essential commands and instructions on how to build the application and process the data. Except from those two essential scripts, a helper script is defined (helpers.R), which includes all the functions needed for further processing of the data and achieving the desirable plot formats.

Apart from the shiny package, several other packages have been utilized in order to add further functionality to the application and visualize the produced results. Indicatively, the DiagrammeR, data.tree and collapsibleTree libraries were very useful towards the visualization of the constructed tree. The latter is more interactive and gives the user the opportunity of collapsing branches and focusing on the branch or level of interest. Another special graph of our application is the logo graph, which contains the common letters of a specific cluster or level and can be produced through the use of the ggseqlogo library. Finally, our R Shiny Application is publicly available from the following URL: https://github.com/mariakotouza/H-CDR3-Clustering

Results

Our method was applied on a real-world dataset comprising 123 clonotypic immunoglobulin (IG) amino acid sequences (deduced from the corresponding IG gene rearrangement sequences) from patients with chronic lymphocytic leukemia (CLL). All sequences utilized the IGHV4-34 gene with 111 (90.2%) being assigned to 6 distinct biologically relevant groups (subsets), and had an identical length of 20 amino acids. These subsets are characterized by the presence of common amino acid sequence patterns within the VH CDR3 of the clonotypic IG. Subset #4 is a major subset and patients belonging to this subset display an indolent clinical course, while the other subsets are minor. In detail, 101 sequences were assigned to subset #4, 2 to subset #207, 2 to subset #4-34/20-1, 4 to subset #4-34-16, 2 to subset #4-34-18 whereas 12 sequences carried heterogeneous receptors and, thus, were not assigned to any subset. Through the application of our tool, a binary tree with 19 levels was constructed. The average value and the standard deviation of each level's identity using the 20 amino-acid and the 11 IMGT Physicochemical classes are summarized in Table 1. The table shows that the identity value increased towards the leaves of the tree. Notably, when the 11 classes were used instead of the individual amino-acids, the total identity value was a little higher as expected.

In more detail, both Average Identity and Average Similarity values are equal to zero at the root of the tree and they almost consistently increase at each level, reaching the value of 98 percent on level 19. At any given point, the Average Similarity value is equal or greater than the Average Identity value, which is reasonable assuming that every similarity group contains one or more individual letters. Regarding the Identity Standard Deviation and Similarity Standard Deviation, it can be observed that both have an initial value of NA at the root, they are slightly increased for levels 1 and 2 and then their values decrease for every level until reaching the leaves. This is an expected outcome of the process, since the standard deviation is a measure of how spread out numbers are. Finally, the Identity Standard Deviation is equal to or greater than Similarity Standard Deviation, because the number of individual letters is greater than the number of similarity groups and therefore the amount of spread is higher. Subset #4 is the largest subset in the present data series comprising 101 clonotypic IG sequences. Most of them (93/101, 92%) formed cluster 13, at level 4 of the clustering process, with level 0 being the root of the tree. Cluster 13 also included 3 non-subset IG sequences. The identity and similarity rates of this cluster were 20% and 30%, respectively. The CDR3 of cluster 13 is 20 amino acids long and consisted of 4 conserved positions and 16 positions that are characterized by variability. Subset #4-34/20-1 consisted of 2 clonotypic IG sequences. These sequences were grouped together at level 9 (cluster 57) with high identity and similarity rates (80% and 95%, respectively). Only 4 positions of the CDR3 were encoded by different amino acids:

A R Y _ _ V V T A _ _ N Y Y Y Y G M D V.
More uniform groups and less variability at a position level is noticed at higher levels of clustering. Consequently, identity and similarity rates increase from one level to another. Subset #4-34-16 consisted of 4 clonotypic IG sequences. Three out of 4 IG sequences (75%) were grouped together at level 5 forming cluster 31 with 45% identity and 50% similarity. Cluster 31 also included a non-subset IG sequence. At the next level of clustering (level 6) the sequences assigned to subset #4-34-16 (3/4, 75%) formed cluster 39 with higher identity and similarity rates (70%|80%). Six positions were characterized by variability.

Discussion

In this study, we applied a new clustering technique to 123 immunoglobulin (IG) sequences of the CDR3 region from patients with chronic lymphocytic leukemia (CLL), all of which express the IGHV4-34 gene. This is a hierarchical method that constructs a top-down binary tree by iteratively breaking the amino-acid frequency matrix into two sub-matrices and using identity and amino-acid similarity as measures for the clustering. The results of this analysis showed that the identity value increases as one transitions from the root towards the leaves of the tree. In the case when the 11 IMGT Physicochemical classes were used instead of the 20 individual amino-acids, the identity value was a little higher as expected.

The proposed method can extract patterns from the clusters at the different tree levels. For example, the pattern that corresponds to cluster 13 (as shown in Fig. 4a) is _ R _ _ _ _ _ _ _ _ _ _ Y _ Y _ G _ _ _, whereas the consensus pattern of cluster 39 (Fig 4d

) is _ _ _ V V P A A _ _ P _ Y Y Y Y G M D V.
These patterns can be directly used by biomedical experts in order to evaluate the cluster and infer potential meaningful interaction at the amino acid level. In our case, clustering is performed to discover previously unknown patterns (unsupervised learning) and then, as a next step, sequence labeling can be performed to assign new sequences to the existing classes.

Future steps involve the application of graph theory on the produced tree, in conjunction with string distance metrics, in order to further refine the clustering process and uncover potential connections between different clusters / nodes of the tree. Moreover, additional testing will be required, especially focusing on large-scale datasets and optimizing the code for speed and memory usage.

Fig. 1 .

 1 Fig. 1. The 11 IMGT Physicochemical classes for the 20 amino acids [9].

Ε 2 )Fig. 2 .

 22 Fig. 2. Block diagram of the tree construction.

Fig. 4 Fig. 4 .

 44 Fig. 4. CDR3 region logos of four clusters of different levels of the tree.

Fig. 5 .

 5 Fig. 5. The binary tree with 12 levels and values for each cluster with the format [cluster id, number of cluster sequences, identity, similarity]

Table 1 .

 1 Average and Standard deviation identity and similarity value of the clusters of the 19 tree levels.

		Identity	Identity	Similarity	Similarity
		Mean	Standard	Mean	Standard
		Value	Deviation	Value	Deviation
	level.0	0.000	NA	0.000	NA
	level.1	5.000	0.000	7.500	3.535
	level.2	38.750	41.708	41.250	39.660
	level.3	34.166	33.078	37.500	31.741
	level.4	37.500	25.083	41.500	25.500
	level.5	55.500	31.837	61.000	28.848
	level.6	71.875	29.753	79.375	23.969
	level.7	65.833	31.211	71.666	26.204
	level.8	66.666	26.394	72.500	22.079
	level.9	68.125	23.594	78.125	19.628
	level.10	70.833	24.579	77.500	18.907
	level.11	73.750	20.310	80.000	16.903
	level.12	75.000	15.811	80.500	14.615
	level.13	79.583	15.877	84.166	13.953
	level.14	80.416	13.221	83.333	12.851
	level.15	85.000	12.792	86.250	13.164
	level.16	86.250	111.505	87.083	11.957
	level.17	87.500	9.414	88.750	8.822
	level.18	93.750	7.723	94.166	7.334
	level.19	98.000	4.216	98.000	4.216

Acknowledgments

The authors would like to thank Katerina Gemenetzi, Andreas Agathangelidis and Kostas Stamatopoulos from the Institute of Applied Biosciences (INAB|CERTH) for providing the dataset as well as for their critical review of the paper.


Shah Nawaz Khan
email: s.khan@fbk.eu

Roberto Riggio
email: r.riggio@fbk.eu

On Edge Cloud Architecture and Joint Physical Virtual Resource Orchestration for SDN/NFV

Keywords: Edge Cloud, SDN/NFV, Orchestration, Kubernetes

5G networks will incorporate virtualized network infrastructures and new technologies such as Software Defined Networking and Network Function Virtualization. The developments have predominantly been made in the core and back-haul network segments. However, the ambitious Key Performance Indicators (KPIs) set for 5G networks will necessitate a renewed focus on the network edge in view of the virtualized infrastructure and SDN/NFV adoption. Several challenges are anticipated to be addressed for cloud at the edge, prime among which is the heterogeneity that spans, among others, hardware, software, radio, networking and virtualization domains. In this paper, we investigate cloud at the edge architecture, its unique challenges and the technology enablers. We present a prototype heterogeneous edge cloud implementation with focus on joint physical and virtual resource management and orchestration while supporting multi-tenancy for edge services. A scheduling & orchestration module is presented that interfaces with Kubernetes cloud management system to support service deployment per two scheduling polices of load balancing and energy saving. We present details of the considered edge cloud platform, the scheduling & orchestration module and its functions inside the edge cloud. Finally, we present some preliminary results and comparisons of the implemented orchestration policies in the context of heterogeneous edge services.

Introduction

Wireless communication is an integral part of modern day lives across the world. In the past decade, requirements on wireless networks to support increasing number of new services beyond the simple voice calls and connecting larger number of end-users have been consistently intensifying. Until recently, mobile communication networks have been addressing these requirements and challenges with continuous evolution of the radio access technologies and by incrementally increasing the overall network capacity. However, with the recent focus on and proliferation of social media, networking, augmented and virtual reality and Internet of Things applications, the requirements and challenges have intensified such that mere evolution of 4G networks will not fully address them. These requirements and challenges are more stringent and diverse ranging from ultra low latency, to massive machine type communications to enhanced multimedia and broadband applications. Moreover, the pace at which the new network services and applications are developing, it is difficult for network operators to keep pace with by upgrading the core infrastructure, back-haul and radio access capacities. It is anticipated that in the next few years, the annual global IP traffic volume will reach an unprecedented 3.3 ZettaBytes (1 ZB = 1000 Exabytes) [START_REF]CISCO: Cisco Visual Networking Index: Forecast and Methodology[END_REF]. These numbers and the associated network requirements have made 5G, the next big step in the mobile wireless communication networks, one of the most important areas of research and development.

Several key technologies are anticipated to play a major role in defining the 5G networks' characteristics including Cloud based virtualized network infrastructure, Software Defined Networking (SDN) and Network Function Virtualization (NFV). These technologies and their practical realization has been an area of immense interest among researchers from both academia and industry. A large number of initiatives have been taken under the 5G networks R&D umbrella developing virtualization, SDN/NFV, Management and Network Orchestration (MANO) solutions [START_REF]OpenStack: Open source software for creating private and public clouds 2018[END_REF], [START_REF]OSM: OpenSource Management & Network Orchestration[END_REF] [START_REF]CORD: Central Office Re-architected as a Datacenter[END_REF]. Most of these initiatives have targeted the core and back-haul network segments with centralized cloud infrastructures. The network edge, that is, the radio, networking and communication infrastructure closest to the end-users (Access Points, Base Stations) has not received the attention it warrants. Some of this can be attributed to the issues such as complexity, lack of standardized platforms, node heterogeneity, distributed nature of the edge infrastructure and the large number of Points of Presence (PoPs). These factors make it challenging to develop and manage cloud infrastructure for the edge and orchestrate the resources for multiple tenants as anticipated for 5G networks. Providing cloud architectural blueprints for resource constrained edge devices, technology components for virtualization and resource orchestration (both physical and virtual) will help bring the benefits of virtualization and SDN/NFV closer to the network edge. In this paper, we focus on lightweight virtualization and resource management at the network edge. We characterize the unique features of network edge, the limitations of state of the art solutions and present a prototype heterogeneous edge cloud with a joint physical and virtual resource orchestrator module. The orchestrator can pro-actively scale the physical and virtual resources available to multi-tenant services per a given policy of load-balancing or energy saving. Our implementation is based on Kubernetes [START_REF]Kubernetes: Production-Grade Container Orchestration System[END_REF] container orchestration system and integrates with its main control elements for placing VNFs at the edge. The prototype edge cloud consists of nodes having different hardware architectures and computing resources integrated into a unified cloud fabric supporting service deployment, multi-tenancy, resource isolation and scaling. We present details of the hardware platform, the cloud virtualization stack and the orchestrator module and evaluate its performance. The rest of the paper is organized as follows. Section 2 discusses state of the art on cloud virtualization and SDN/NFV related technologies in view of their limitations for the network edge. Section 3 details the main contribution of this paper and presents an edge cloud architecture and virtualization stack addressing key constraints and function requirements. Preliminary results from the evaluation platform are presented in section 4. The paper concludes with a summary and an outlook for future work in section 5.

State of the Art

Recently, the network edge has attracted significant interest from researchers as evident from the abundant literature available under synonymous concepts but different terminologies. Mobile Edge Computing, FoG Computing, Multi-Access Edge Computing, Cloudlets etc., are a few well-known umbrella terms which in essence, aim to realize computing and network services closer to the end-users. In spite of these parallel developments in nomenclature and predominantly theoretic research works, there is a general dearth of clarity on the platform for edge computing such as hardware, software and architecture. A common consensus from the SDN/NFV community is that the benefits of cloud and virtualization must be realized at the network edge to support 5G use-cases and features such as multi-tenancy, resource slicing and orchestration etc. However, the available virtualization and resource orchestration technologies tailored towards data-center oriented centralized clouds and Virtual Machine (VM) based VNF abstraction do not suit the constraints of the network edge. For creating the cloud abstraction, the existing commercial and open source hypervisor technologies for SDN/NFV such as VMware ESX and ESXi, Microsoft Hyper-V, Xen, KVM etc., have a large hardware resource requirements footprint [START_REF] Blenk | Survey on Network Virtualization Hypervisors for Software Defined Networking[END_REF]. Usually in order of multiple high power processing cores and Gigabytes of memory, these hypervisors are virtually impossible to deploy on dispersed and resource constrained nodes. From the 5G networks and SDN/NFV perspective, the additional management and control elements required in the ETSI MANO architecture [START_REF]ETSI: Network Functions Virtualisation[END_REF] such as Virtual Infrastructure Managers (VIM), VNF Manager (VNFM) and NFV Orchestrators need to be deployed. Many realizations of these components such as OpenDaylight, OSM, OpenO, OpenBaton etc., are tied to the centralized data-center oriented cloud model. The hardware resources (storage, compute, networking) required for realizing the same cloud architecture at the network edge are literally unavailable as edge nodes are resource constrained devices and are spread out over larger geographical areas. The level of hardware abstraction realized for centralized clouds where VNFs are completely devoid of information about the underlying hardware features may also prove limiting for the edge. For certain services, it may make more sense to expose the features available at the nodes such as hardware architecture, computing resources, dedicated hardware etc., to tailor their scheduling towards those nodes. Another limiting factor for cloud at the edge using existing virtualization solutions is the VM based VNF realization. A VM adds many layers of software abstraction (e.g., guest OS, unnecessary libraries) between the physical processor and the actual software that carries out a virtual network function. This makes VM based VNF realization unsuited to the limited resources available at the edge cloud. Moreover, with large sizes of VM based VNFs and the limited memory and storage resources of individual edge nodes, the scheduling and VNF placement problem becomes extremely complicated. Another constraint that proves a hindrance is the dispersed nature of edge nodes and the variance of networking links among them (throughput, firewalls etc.). This necessitates a more distributed realization of edge cloud rather than the centralized solutions available for core SDN/NFV implementations. A federation umbrella could tie together dispersed edge clouds and provide an interface to service providers through which network services could be deployed across distributed edge cloud. Finally, a number of advances have been made in the cloud resource management and cloud-native software design domains which have not been considered in the cloud for SDN/NFV context. These developments include the micro-services based software design, DevOps, Serverless Functions etc., which place a higher emphasis on resource utilization efficiency, robust cloud-native designs and rapid deployment. These developments are more suited to the edge cloud context and call for a design where VNFs are composed of small services interfaced together to create the end-to-end network service.

Edge Cloud Architecture and Resource Orchestration

Considering the challenging KPIs for 5G networks addressed in section 1 and the limitations of tailoring the existing virtualization solutions to the network edge addressed in section 2, it can be argued that cloud at the edge must be based on lightweight virtualization technologies and must integrate heterogeneous hardware while supporting the 5G service requirements. Furthermore, the ETSI MANO stack where specific control elements such as VIM, VNFM and NFVO are responsible specific SDN/NFV management tasks should be realized, potentially in a converged form and with small resource requirements footprints. Considering these and the recent developments in the cloud and virtualization domains such as Micro-Services Architectures, DevOps and Function as a Service (FaaS) etc., where resources are utilized to their fullest, we condense the requirements of edge clouds into following list.

-Heterogeneity: Cloud at the edge must integrate hardware nodes of different architecture and resources (compute, storage, memory, networking) into a unified abstraction providing similar interfaces to the service providers as in centralized cloud infrastructures. Moreover, cloud at the edge should hide the complexity of these underlying hardware differences from the deployed services i.e., a service provider should not be forced to design the service template according to the hardware architecture of the host infrastructure. Any mapping between the requirements of a VNF and the capabilities of the underlying hardware should be done by the cloud control and management functions. -Resource Awareness: Cloud at the edge should be several order of magnitude lighter in terms of hardware requirements compared with centralized cloud technologies. This implies that the edge cloud and MANO features should be realized using virtualization technologies that can run efficiently on hardware resource constrained devices deployed in end-user premises. -MANO Features: MANO is an essential requirement from the SDN/NFV management and control point of view. Service providers must have a simple interface through which they can deploy an edge resident network service and manage it in isolation from other services. Therefore, the essential features of VIM, VNFM, and NFVO must be realized albeit considering the discussed resource constraints and heterogeneity of the network edge. Moreover, management polices for cloud resources (physical and virtual) and tenant-specific service orchestration must be supported. -Services: Cloud at the edge must support the distinctive features of centralized SDN/NFV clouds such as resource slicing, multi-tenancy, service isolation, auto-scaling etc. Moreover, support for new cloud-native software design principles such as micro-services architecture should be supported to utilize the dispersed computing resources more effectively.

Considering these requirements and constraints, we focus on container based virtualization where VNFs can be realized in container images instead of VMs.

Containers have a significantly low overhead not only in terms of storage, memory and processing requirements compared with VMs but they also eliminate the need for the hypervisor layer. This resource efficiency usually comes at the cost of security and service isolation but several measures can be taken to address these concerns including container specific Linux kernels and more sophisticated access policies. Beyond the container based VNF abstraction, we use Kubernetes to emulate the role of VIM and NFVM by extending it to support (a) heterogeneous hardware nodes, (b) custom scheduling and orchestration module and (c) allowing for mixed service deployment in architecture agnostic manner. The remainder of this section details the considered virtualization stack and MANO feature realization. The weave element is deployed alongside Kubernetes core control elements to provide overlay connectivity among service VNFs. The details of these control elements is beyond the scope of this work and can be found in [START_REF]Kubernetes: Production-Grade Container Orchestration System[END_REF]. Kubernetes uses a Master-Worker based architecture where the master node hosts all the main control plane elements while the worker nodes host a subset which enable communication with the master node. We host the master node functions on a Raspberry Pi and integrate the rest of ARM and x64 nodes as workers. However, by default, Kubernetes does not integrate nodes of different hardware architecture into a single cloud abstraction i.e., the worker nodes must be of the same architecture as the master node. This limitation arises during the Kubernetes installation where the master node downloads control plane elements container images on the worker nodes according to its own architecture. Therefore, if the master node elements are instantiated on an ARM node as in our case, Kubernetes will download ARM based images for the rest of cluster members including on the x64 nodes. To overcome this problem, we modified the configuration for master node allowing it to download architecture specific images on worker nodes. To this end, we modified the control daemon configuration for kube-proxy by duplicating it for x64 architecture and specifying the correct architecture images.. This problem is expected to be natively addressed in future Kubernetes versions where control plane images are expected to be architecture specific using the Docker image manifests. Most of the ETSI MANO control functions can be realized inside the Kubernetes cluster including VIM, VNFM and NFVO by either using the default Kubernetes control elements or by extending the cluster with new control elements. We utilize both approaches in our prototype edge cloud setup. For the VIM and NFVO functions, we have developed a dedicated scheduling and orchestration module (detailed in the next subsection). This module, depicted in yellow color in figure 1, takes the role of VIM by elastically scaling the infrastructure resources to the requirements of the deployed services following a scheduling policy. To this end, the available physical resources of the nodes and requirements of the deployed services are actively monitored. The same module is responsible for network service orchestration including deployment, scheduling VNFs, chaining and scaling. The VNFM role is shared between the Kubernetes control elements that actively probe the state of deployed VNFs and restart a service VNF if detected to be dysfunctional. Finally, to ensure service isolation and that tenants are allocated a specific set of resources, both physical and virtual, we use the Kubernetes namespaces feature together with labeling of the nodes. A Kubernetes namespace is a virtual container for deployed services in the edge cloud which can enforce resource restrictions on the services falling in a particular namespace. With namespaces and node labels, a service can not only be confined to specific nodes e.g., deploy on ARM or x64 nodes but also be restricted to a subset of the resources on those nodes.

Scheduler & Orchestrator Architecture

Figure 2 shows the internal architecture and functions of the scheduler and physical & virtual resource orchestrator component integrated with the kubernetes control plane elements. The orchestrator receives a policy parameter as input at start-up time and spins up two concurrent processes namely the Policy function and the Resource Monitor function. The policy parameter determines the scheduling behavior and the placement of service VNFs on the worker nodes. Currently, two options are available for policy parameter targeting energy saving and load balancing. In the energy saving mode, the orchestrator targets running the least number of worker nodes in running state in order to conserve energy expended by the cluster. Practically, this implies a single worker node is kept running for as long as the deployed service VNFs are providing the required quality of service. A new worker node is boot up only when the running services are scaled beyond the resources of the running worker nodes or when new deployments are requested that either require more physical resources or are targeting a particular node/set of nodes in the cluster with specialized hardware. A service deployment may target a specific set of nodes for several reasons such as due to assigned quota, resource constraints or performance reasons. For example, a service might request a VNF to be deployment on a node that has SSD storage or a GPU for graphic processing. The master node keeps a record of the hardware attributes of the worker nodes in the cluster. The orchestrator can map these attributes to the requirements of a service deployment. The other option for policy parameter is load balancing which keeps all the worker nodes in running state and priorities minimal delay in service deployment and scaling latency. This is achieved by compromising on the overall energy consumption of the cluster. The whole service deployment and scheduling process is executed as follows:

-Service deployment is requested by providing a service template file to the Kubernetes API Server component. To do this, a service provider needs to authenticate its access with the API Server. The service template is provided as Yaml [START_REF]YAML: Yaml ain't markup language[END_REF] file describing the service at a high abstraction level such as the required VNFs (Docker containers to run), hardware resources required, service and target node attributes. -The API server receives the deployment requests and pushes a scheduling request trigger to the orchestrator component. At this point, the scheduling and orchestration component has to determine the criteria of VNF placement resource scaling. These include (a) find the worker nodes in the cluster that have the required resources available to host the VNFs requested by the deployment (b) determine, for each VNF, a single node out of the possibly multiple qualified nodes as host (c) determine the state of the node which can be either ready or switched off (d) bind the VNF to the selected node if the node is in ready state (in load balancing policy) or first boot the node (in energy saving policy) and then bind the VNF to the selected node once it is in ready state.

-For task (a), a Resource Mapper function takes the requested resources of a given service deployment request and compares them to the active resource map maintained by the Resource Monitor process instantiated by the scheduler. It should be noted that the Resource Monitor keeps a real-time load of running nodes only and an indication of the physical capacity of the nodes that are offline but part of the cluster. This way, the resource mapping process considers all nodes in the cluster instead of the online ones if the scheduling policy is load balancing or the VNFs image is for a particular hardware architecture. The resource mapper will return the list of all qualified nodes that can host the service VNFs.

-For task (b), the policy function is responsible for reducing the possible list of several qualified worker nodes to a single node based on the policy used by the scheduler. It is worth pointing out here that this process is sequential i.e., a deployment may contain several VNFs (containers) and the deployment of those are carried out sequentially in order to fulfill the load balancing or energy saving policy requirements.

-For task (c), the node activator module takes care of booting up a node based on the node ID in case it is offline. It uses the Wake On Lan feature of NICs on the worker nodes. -Finally for task (d), the binder function does the actual scheduling of the VNF on the selected node. For this, the binder interacts with the kubernetes control plane elements to download the VNF container image (if not present locally), and then boot it up with the required resources provisioned. Lastly the Binder function logs the success and failure events of the binding process.

Evaluation

We have carried out preliminary evaluation of the resource consumption (Storage, Memory, CPU) of the control elements in our prototype edge cloud environment and its performance in terms of service deployment and deletion latencies.

Figure 3 shows the image sizes of Kubernetes core control elements and our addon purpose-built scheduler and orchestrator (labeled as Sch+Orch). The add-on component has been implemented in Golang and compiled into a static binary.

As evident from figure 3, the container image sizes of the core Kubernetes elements are small enough to be deployed on tiny storage resource constrained devices. Second, it is very easy to create add-on elements for specific functionality such as Sch+Orch that has our statically linked binary embedded inside a tiny sized Docker image (5.3MB only). Figure 4 shows the CPU and Memory utilization on the Master node (Raspberry Pi 3) of the three main control elements that are used during deployment of a service containing two VNFs. In general, all elements have a very small resource utilization footprint during an actual service deployment. However, since our custom Sch+Orch element takes away the scheduling and orchestration duties from the Kubernetes default scheduler, there is no observable variance in its CPU and Memory consumption during the deployment. In fact, the only significant variance is in the memory footprint of the API server and the CPU utilization of our scheduler VNF. This is expected since a service deployment is an active interaction process between the Kubernetes API Server and the scheduler VNF. To evaluate the scheduling and service deployment capabilities of our platform and to analyze the effects of scheduling policies on the deployment latencies, we deployed a couple of prototype web and video streaming edge services in our evaluation platform. Figure 5 shows the two prototype services and their VNF composition. On the left side, a service that combines web services and access services is depicted presented. This service targets the ARM nodes in the cluster as the Docker images for Nginx VNF and WiFi access point VNFs have been developed for that platform. The association between VNF architecture and the nodes having that physical architecture is done by label matching between the service template (yaml file) and the labeled nodes. On the right side in figure 5, a mixed architecture streaming service is depicted. This service template exemplifies a video streaming scenario where users of different privileges are provisioned video streaming service of different quality. The VNFs in this service includes an access VNF (WiFi access) and two streaming VNFs that stream video files in different quality. The video files are provided to both streaming services by a back-end video store VNF. To eliminate some of the latency fac- tors of service deployment, all the Docker images for required service VNFs were downloaded prior to the actual deployment. Figure 6 shows the service deployment and deletion time for the two services per the load-balancing and energy saving scheduling policies. The web-service is entirely hosted on ARM nodes and uses the load-balancing scheduling policy. As evident, the service takes approximately 26 seconds to be active in the edge cloud including providing the WiFi Access and the Nginx web-server. The video streaming service is deployed on different architecture nodes where most of the streaming VNFs are deployed on x64 machines while the WiFi access VNF is hosted on ARM node. This service uses the energy efficient scheduling policy where the x64 machine is first booted up from an OFF state before the streaming VNFs can be deployed. The x64 node takes around 50 seconds to boot up and associate with the Master node. The load balancing and energy saving policies can be mixed depending on the service requirements. Figure 6 also shows the respective deletion times for the two example edge services.

S c h + O r c h K -

Conclusion

In this paper, we characterized the core attributes of network edge and its implications for brining cloud abstraction and SDN/NFV benefits to the end-users. We also presented a candidate architecture for heterogeneous edge clouds supporting multi-tenant edge services. Moreover, a prototype VNF scheduler and joint physical virtual resource orchestrator was presented and evaluated for prototype edge service deployment. In future we intend to extend the orchestration umbrella across multiple distributed edge-clouds by tying together a federated

Fig. 2 .

 2 Fig. 2. Software architecture and functional diagram

Fig. 3 .Fig. 4 .

 34 Fig. 3. Container image sizes for core control elements

Fig. 5 .

 5 Fig. 5. Prototype edge services for web hosting and video streaming

  3.1 Physical Architecture & Virtualization StackFigure1depicts a subset of the physical architecture of edge cloud infrastructure and the virtualization stack used in this work. The hardware platform comprises five ARM architecture based Raspberry Pi 3 nodes and two Intel x64 based Laptops assembled in a cluster and local network. All nodes run platform specific Linux distribution and a container runtime daemon (Docker in our case). The selection of Docker in general and container based virtualization in particular was made considering the resource limitations and hardware heterogeneity of the nodes. To integrate the nodes into a cluster where services can be deployed at cloud abstraction level, we use Kubernetes[START_REF]Kubernetes: Production-Grade Container Orchestration System[END_REF]. Kubernetes is an open source container management and orchestration platform that can run on several types of physical nodes and has a considerably small resource requirements footprint compared with hypervisor and VM based virtualization platforms. The core Kubernetes control elements such as Etcd, Controller Manager, API Server, DNS Fig. 1. Physical architecture and virtualization stack of the edge cloud and default Scheduler run in Docker containers on a single Raspberry Pi node.

	MASTER: ARM-RPi-3		Worker: ARM-RPi-3	Worker: Intel x64
	Scheduler + Phy/Virt Resource Orchestrator	VNF 1	VNF 2	VNF 1	VNF 2
	Etcd	API Server	Cotroller Manager	DNS Scheduler Weave	Proxy Kubelet Weave	Proxy Kubelet Weave
			Container Runtime	Container Runtime	Container Runtime
			Linux		Linux	Linux
				Hardware Switch

cloud environment. We also intend to investigate the interplay between selforganizing edge services and resource orchestration in host infrastructure using Artificial Intelligence algorithms.