Mohammed Ismail Smahi
email: mohamedismail.smahi@univ-tlemcen.dz

Fethallah Hadjila
email: fhadjila@univ-tlemcen.dz

Chouki Tibermacine
email: chouki.tibermacine@lirmm.fr

Abdelkrim Benamar
email: abenamar@univ-tlemcen.dz

A deep learning approach for collaborative prediction of Web Service QoS

Keywords: Web Service, QoS Prediction, Deep Autoencoder, Self-Organizing Map

Web services is the corner stone of many crucial domains, such as cloud computing and the Internet of things. In this context, QoS prediction for Web services is a highly important and challenging issue. In fact, it allows for building value-added processes as compositions and workflows of services. Current QoS prediction approaches, such as collaborative filtering methods, mainly suffer from the problems of data sparsity and cold-start obstacles. In addition, previous studies have not explored in depth the impact of geographical characteristics of services/users and QoS rating on the prediction problem. To address these difficulties, we propose a deep learning-based approach for QoS prediction. The main idea consists of combining a matrix factorization model based on a deep auto-encoder (DAE) and a clustering technique based on the geographical characteristics to improve the prediction effectiveness. The overall method proceeds as follows: First, we cluster the QoS data using a self-organizing map that incorporates the knowledge of geographical neighborhoods; by doing so, we allow for the reduction of the data sparsity while preserving the topology of input data. Besides that, the clustering step effectively handles the cold start problem. Second, for each cluster, we train a DAE that minimizes the squared loss between the ground truth QoS and the predicted one. Third, the missing

Introduction

Quality of service (QoS) is an important factor for performing Web service selection and recommendation over the internet [START_REF] Huang | An optimal qos-based web service selection scheme[END_REF][START_REF] Yu | Efficient algorithms for web services selection with end-to-end qos constraints[END_REF][START_REF] Zheng | Collaborative web service qos prediction via neighborhood integrated matrix factorization[END_REF]. In this context, QoS prediction constitutes a crucial step in building a recommendation system. In practice, QoS values of services constantly change due to environmental constraints (e.g., IT infrastructure and network load), and this fact makes QoS prediction a challenging task [START_REF] Chen | Web service recommendation via exploiting location and qos information[END_REF]. In addition, we observe that end users may have invoked only a small number of web services; consequently, the user-service QoS data are likely to be sparse [START_REF] Huang | Applying associative retrieval techniques to alleviate the sparsity problem in collaborative filtering[END_REF], which may have a significant impact on the accuracy of QoS prediction.

To address the QoS prediction issue, many existing works have leveraged Collaborative Filtering (CF) [START_REF] Adomavicius | Toward the next generation of recommender systems: A survey of the state-of-the-art and possible extensions[END_REF] to infer the missing data. These approaches are based on the exploration of the historical QoS data of Web services recorded from previous interactions. They typically utilize a user-service QoS matrix to model all historical QoS records. The current CF methods can be divided into two classes: memorybased [START_REF] Chen | Collaborative filtering for orkut communities: discovery of user latent behavior[END_REF][START_REF] Deshpande | Item-based top-n recommendation algorithms[END_REF][START_REF] Zheng | Qos-aware web service recommendation by collaborative filtering[END_REF] and model-based [START_REF] Zheng | Collaborative web service qos prediction via neighborhood integrated matrix factorization[END_REF][START_REF] Koren | Matrix factorization techniques for recommender systems[END_REF]. Memory-based approaches consist of two stages: the first one computes the similarity between services (or users) through the use of QoS data, while the second stage computes the missing QoS using average weighting of the known values of similar services (or users). Model-based approaches learn a set of latent factors from the QoS matrix and make predictions. Matrix Factorization (MF) techniques are the most representative approaches of this class. Despite the advantages and the ease of use of CF methods, most of them do not handle contextual attributes (e.g., geographical information of users/services) to enhance the accuracy of prediction; additionally, they suffer from the cold-start problem [START_REF] Wei | Collaborative filtering and deep learning based recommendation system for cold start items[END_REF]. The cold-start problem rises when we want to predict the QoS for new services that have no past QoS records or only few QoS data. The sparsity of the QoS matrix also limits the effectiveness of CF methods and results in a poor accuracy performance. To address the above mentioned issues, we propose an extended CF-based approach that uses both Matrix Factorization (MF) and clustering to handle data sparsity and QoS fluctuations.

This paper is an extension and enhancement of a work originally presented in the European Conference on Service-Oriented and Cloud Computing [START_REF] Smahi | An encoder-decoder architecture for the prediction of web service qos[END_REF]. The proposed approach extends our previous work through the use of both a DAE and clustering. The contributions of this paper are summarized as follows:

(A) To consolidate the related works, we added new approaches that involve deep learning methods, matrix factorization techniques, and contextual recommendation. (B) To face the QoS sparsity issue, we perform a clustering of the original QoS matrix and yield a set of submatrices that share a subset of contextual characteristics (such as geographic information or service providers).

After that, each DAE is trained on a single sub-matrix. The clustering step is realized using a Self-Organizing Map (SOM) [START_REF] Kohonen | Self-organization and associative memory[END_REF]. This choice is motivated by the ability of SOM to preserve the topological properties of QoS data (see the motivation Section for more details). The learned neurons of the map are grouped using K-means to create more representative clusters. We notice that the use of SOM allows us to address the cold-start problem (i.e., estimation of QoS for new services or new users).

More precisely, we initialize the new service (or new user) with the QoS data of the closest cluster-heads (that share the same contextual attributes or the nearest ones). (C) To improve the optimization results, the different hyperparameters of DAEs are fine-tuned through the use of k-fold cross-validation [START_REF] Refaeilzadeh | Cross-validation[END_REF]. (D) To reduce the model sensitivity to over-fitting which can harm the final performances, we add a random noise to the input data during the training phase (i.e., we use a denoising autoencoder [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF] which is a variant of the standard autoencoder). (E) To evaluate the effectiveness and the accuracy of the proposed DAE, we perform several experiments by adopting different densities and different hyper-parameters.

The remainder of the paper is structured as follows. We first introduce in Section 2 some background material on both autoencoders and Self Organizing Maps. In Section 3, we present the proposed method for Web service QoS prediction. After that, we show in Section 4 the experiment design, the results and the analysis. In Section 5, we present the related works. Finally in Section 6, we summarize the contribution of the paper and put forward our future work.

Background

We introduce the basic background material on both autoencoders and Self Organizing Maps.

Autoencoders

The autoencoder [START_REF] Rumelhart | Learning internal representations by error propagation[END_REF][START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF] is an unsupervised neural network. It allows to reproduce the original data from the input layer to the output layer through one or several hidden layers. This unsupervised neural network learns a representation of the inputs that yields the least deformation. Typically, the mechanism consists of the encoding data from the input layer until a central hidden layer, called latent factors layer, through one or more hidden layers. After that, the decoding of the central layer data is performed by a set of decoding layers (having the same size as in the encoding layers). As can be seen in Figure 1, the autoencoder's architecture consists of three parts: the encoder, code (or latent factors) and decoder. The latent factors part represents a single layer of artificial neural network (ANN) where the size should be less than the input layer (which is denoted by K). This hidden layer must be compact and meaningful. The encoder part is represented by a fully connected feed-forward non recurrent ANN. It can be viewed as an encoding function F e that takes a vector x ∈ R D as input data and maps it to a hidden representation z ∈ R K (z = F e (x)). Finally, the decoder part is the reverse mapping of the encoder operation. It is a fully connected feed-forward non recurrent ANN. It is considered as a decoding function F d that takes the resulting latent representation z then maps it back to a reconstructed vector x ∈ R D ( x = F d (z)). Thus, the output will be:

x = F d (F e (x)) (1)
If we assume a basic version of the autoencoder (which uses one hidden layer), the encoding decoding functions will be defined as:

x = σ (W z + b ) with z = σ (W x + b) (2)
such that: σ and σ are transfer functions which can be linear or non-linear 1 . W and W are two weight matrices having the dimensions (K, D) and (D, K) respectively. b and b are the offset vectors with K and D dimensions respectively.

In terms of learning, the autoencoder is trained to minimize the dissimilarity function (or the reconstruction error)

argmin W,W ,b,b n ∑ i=1 δ (x i , F d (F e (x i ))) (3)
where δ is a dissimilarity function, such as the square loss or the cross entropy loss, x i is an example that belongs to the learning dataset and n is the size of the dataset.

Denoising Autoencoders

The denoising autoencoder [START_REF] Vincent | Stacked denoising autoencoders: Learning useful representations in a deep network with a local denoising criterion[END_REF][START_REF] Goodfellow | Deep Learning[END_REF] is a variant version of the basic autoencoder presented above. This kind of autoencoder aims to tackle the noise problem. The goal is to reconstruct the original input data using a corrupted version of it. Note that several corruption types can be applied on the input data, such as: (a) Gaussian noise: where an isotropic Gaussian noise is added to a subset of the input x. (b) Masking noise: A fraction v of the original input x is randomly forced to zero. (c) Salt-and-pepper noise: A fraction v of x is randomly set to their minimum or maximum possible value.

In terms of learning, the denoising autoencoder is still trained to minimize the same reconstruction error between the cleaned input and its reconstruction:

argmin W,W ,b,b n ∑ i=1 δ (x i , F d (F e (x i ))) (4)
where x is a copy of x that has been corrupted by some kind of noise cited above. In this case, the denoising training forces the encoding function F e and decoding function F d to implicitly learn the structure of the original input data [START_REF] Alain | What regularized auto-encoders learn from the data-generating distribution[END_REF].

Self-Organizing Map

Self-Organizing Map (SOM) or Kohonen feature map [START_REF] Kohonen | Self-organization and associative memory[END_REF] is a kind of unsupervised neural network based on a competitive learning paradigm. This method is usually used in performing tasks such as clustering, outlier detection and dimensionality reduction [START_REF] Kiang | Self-organizing map network as an interactive clustering tool-an application to group technology[END_REF]. The objective of this network is to reduce a high-dimensional input data space into a lowdimensional map of neurons (usually two dimensions) through the use of self-organizing neural networks while preserving their topological structure.

Self-Organizing Map principles

In terms of architecture (see Figure 2), the network is composed of two layers: the input layer and the map layer; the layers are fully connected. The map layer comprises a set of neurons M = {n 1 , . . . , n k } arranged in a two l-dimensional map of neurons (the map is usually represented graphically as a rectangle with k neurons (k = x dim ×y dim )). Each neuron has an associated weight vector W known as the code-book. The code-book describes the typical data profile at a given training step t, where T is the number of training time step instances 2 :

W (t) = { - → w 1 (t), . . . , - → w k (t)} t ∈ {1, . . . , T }.
The competitive learning process attempts to construct a nonlinear topology that ensures the mapping of the input vectors X = { - → x (t)|t ∈ {1, . . . , T }}, into a set of neurons M on the map. In terms of learning, the SOM Algorithm consists of two steps, a competitive and an update step. In the first one, each data vector -→ x (t), for a particular training instance t, is mapped to its best matching neuron on the map : bm( -→ x (t)) = n i ∈ M. The winning neuron is denoted as:

n i ( - → x (t)) = argmin j∈1,...,k - → x (t) -- → w j (t) (5)
where each neuron n i codes the subset of input data space, whose elements are closer to n i .

In the update step, the weight vectors of both the best match neuron n i and its whole neighborhood neurons Γ (n i ) are adjusted toward the presented input vector -→ x (t) using the following equation:

w j (t + 1) = w j (t) + η(t) • h i, j (t) • (x(t) -w j (t))
where j ∈ Γ (n i ), 0 < η < 1 is the learning rate, and h i, j is the neighborhood function. It is important to mention that the neighborhood of the wining neuron (n i ) is progressively reduced until it reaches a set of size one. A frequently used neighborhood function is the Gaussian:

h i, j (t) = exp - d(i, j) 2 2σ 2 (t) ) (6)
2 SOM training stops when the weight vectors are stabilized or the maximum number of iterations is reached.

where d(i, j) is the distance between the neuron n i and neuron n j on the map and the radius σ (t) = σ 0 • exp t T max decreases after each iteration to restrict the area of the neighborhood.

Batch learning

A batch algorithm is a useful variant of SOM learning process presented above, which is widely used in parallel implementations [START_REF] Wittek | somoclu: An efficient parallel library for self-organizing maps[END_REF]. If the whole training data is available beforehand, SOM has a batch formulation of update of all weights:

w j (T ) = ∑ T t =1 h i, j (t )x(t ) ∑ T t =1 h i, j (t ) (7)
It means that the best matching neurons for all input vectors are selected and hereafter their weights and the weights of their respective neighbors are updated (according to Equation 6).

Initialization Methods

The results and the performances of SOM method are strongly influenced by the initialization of neurons' weights. To this end, two types of initializations can be used. The first one consists of a random initialization of the code-book. In this method, random weight vectors are selected from the data points. Technically this method is easy to implement, but its main weakness lies in its high running time and its non deterministic side-effects on results. To tackle the previous problem, a linear initialization technique can be used. The initial values of code-book are selected from the subspace spanned by the Principal Component Analysis (PCA) algorithm [START_REF] Jolliffe | Principal component analysis[END_REF]. The behavior of SOM method using PCA initialization converges more rapidly towards a better clustering and ensures a deterministic clustering result [START_REF] Akinduko | Som: Stochastic initialization versus principal components[END_REF].

Proposed Method

First, we introduce the architecture model of the proposed approach. Then, we motivate the choices made in the design of this architecture model. At last, we present in depth how QoS prediction is made thanks to this model.

Architecture model

Figure 3 illustrates the architecture model of the proposed approach. First, we assume that there are m users and n Web services in our system. QoS values are collected from user's invocation of Web services. We use various sources to get this data (such as social networks, a monitoring system, or direct feedbacks). The collected values constitute a m × n Fig. 3 Overview of the prediction system architecture rating matrix M = {r u,s } m×n . We assume that QoS matrix contains missing values because it is impossible to invoke all services by all users. In order to predict the missing values (first challenge), we propose an approach that is composed of the following three stages:

1. First, we cluster the initial QoS rating matrix according to the geographical characteristics. More specifically, we will perform a clustering based on SOM Algorithm. As mentioned previously, this method allows us to create clusters whose members are closer in terms of contextual attributes. This step aims to produce a set of clusters that share the same contextual characteristics. Additionally, it allows to decrease the sparsity of the QoS matrix (i.e., the QoS matrix of each cluster is less sparse in comparison to the initial QoS matrix) by considering the top-k clusters. 2. Second, we perform the learning of latent factors for each cluster by leveraging a deep denoising autoencoder.

During our experiments, we used a deep autoencoder with three hidden layers (this choice is explained in the experimental study). We assume that all the clusters are trained with the same number of hidden layers. 3. Third, the learned deep autoencoder predicts missing QoS values and stores them in the initial dataset.

The second challenge that interests us is the cold-start problem. Based on the results of the Data Clustering and Autoencoder Building phases, this task describes the arrangements to ensure the best solution for this problem (later in this section, we will describe the detailed mechanism).

Motivations behind architecture design choices

Autoencoder

Our proposition uses a DAE that allows for predicting QoS values of services. This choice is mainly motivated by the high accuracy achieved by matrix factorization methods, such as autoencoders [START_REF] Li | Deep collaborative filtering via marginalized denoising auto-encoder[END_REF]. According to [START_REF] Koren | Matrix factorization techniques for recommender systems[END_REF][START_REF] Xu | Context-aware qos prediction for web service recommendation and selection[END_REF], matrix factorization techniques are more accurate (in terms of prediction performance) than the remaining CF methods. In addition, the flexibility of DAE (e.g., the number of layers and the activation functions) is leveraged to ensure the best performance. Indeed, it is possible to vary the activation functions of each layer (e.g, relu, so f t plus, leaky-relu, ELU) in order to maximize the prediction accuracy.

SOM

To boost the accuracy of our prediction model, we have to address the sparsity issue of the QoS matrix. A higher sparsity results in a lower quality of prediction. To handle this situation, we cluster the QoS data into groups using the SOM method.

We have chosen this method for different reasons, including the following:

-The SOM method has the ability to preserve the density and the topological structure (the form) of the original data; -SOM-based clustering allows us to address the cold-start problem. More specifically, if we handle a new service (or user) (without QoS values), we first compare its contextual attributes with the context of the cluster-heads; then we initialize the QoS of the new service (or user) with the QoS values of the closest cluster-head; -SOM network ensures a more faithful representation (in terms of contextual attributes such as country, provider and autonomous systems) of the services (or users), that is, the services (or users) having similar contextual attributes are more likely to belong to the same cluster; -Their ability to preserve the neighborhood relationships (i.e, two adjacent services/users are represented by two adjacent cluster-heads, and two neighboring cluster-heads share a part of their geographical characteristics while the other properties are different).

According to [START_REF] Chen | Service recommendation based on separated time-aware collaborative poisson factorization[END_REF][START_REF] Chen | Exploiting web service geographical neighborhood for collaborative qos prediction[END_REF][START_REF] Tang | Collaborative web service quality prediction via exploiting matrix factorization and network map[END_REF], QoS values are directly influenced by users and services geographical characteristics. The more the services share the same geographical characteristics the more the QoS values are similar.

In practice, the geographical characteristics of Web services can be categorized into three different groups: Provider group (P), AS3 group (A ) and Country group (C ). One group is a subset of another (i.e, P(i) ⊆ A (i) ⊆ (i), where i ∈ {service, user}). This implies that the more the group is restricted, the more the QoS values are likely to be similar.

As we will see in the experimental evaluation Section, an empirical analysis is conducted to demonstrate the influence of geographical characteristics on the clustering results of SOM.

Collaborative QoS prediction

The first algorithm (Algorithm 1) enables us to cluster the input data (services S or users U) according to the batch version of SOM principle. Firstly (line 1), we initialize the weight vector for all neurons on the map (code-book initialization). To this end, we use a personalized initialization mechanism to ensure a deterministic behavior and to preserve the neighborhood relationships (more details are given in Section 4.3.3). After that (from line 3 to 5), we determine for each input its winner neuron (BMU) on the map. In next steps (from line 6 to 12), we adjust the weight vectors of all BMUs and their neighborhoods according to Equation 7. The process is repeated as many times as necessary (epoch times as specified in line 2). At the end, we create the final clusters by processing the code-book with a k-means Algorithm [START_REF] Hartigan | Algorithm as 136: A k-means clustering algorithm[END_REF].

max = 8 14 < C 1 , • • • ,C max >= KmeansCluster(W, max) postprocess the code-books with k-means algorithm 15 return < C 1 , • • • ,C max > list of services clusters
In the second algorithm (Algorithm 2), we learn a model that achieves the highest QoS prediction performance. The goal is to infer the best learning model for our deep autoencoder. This algorithm builds a deep autoencoder for each cluster C i . In line 3, we split the current cluster C i data into k-f old folders. Each folder represents a given percentage (density%) of the available cluster data. This parameter represents the ratio between the training dataset and the entire dataset. For example, if density = 20% then the training set represents 20% of all cluster data and the validation set represents the remaining 80% of QoS data. After that in line 7, we learn a deep autoencoder for each folder and with respect to all training sets (∀S ∈ T ). The goal of this step is to minimize the squared error between the model prediction and the desired values. It is important to mention that the training is controlled by the noise parameter. If this parameter is null (i.e, no noise), we use the deep autoencoder architecture. Otherwise, the selected architecture of the model is a denoising autoencoder. In line 8, we calculate the deep autoencoder error achieved on the validation dataset. In line 11, the k-f older results are averaged to produce a single validation error (mve). Finally, we return the best training model as well as the mean validation error.

Algorithm 2: Autoencoder k-fold cross-validation

Inputs : C = {C 1 , • • • ,C max } , density, noise, k-f old 1 for i → max do 2 T = C i training set 3 < f older 1 , ..., f older k-f old >= Partition(C i , density, k-f old) 4 for k : 1 → k-f old do 5 T = T -f older k 6 V = f older k validation set 7 Model k = argmin ∀Sm∈T   1 |T | |T | ∑ m=1 |DAE(S m , noise, |C i |) -S m | 2   8 err k = 1 |V | |V | ∑ m=1 |Model k (S m ) -S m | 2 9 end 10 Model k = argmin ∀k∈k-f old (err k ) get the k-th best model in term of err k 11 mve = 1 k-f old k-f old ∑ k=1 err k 12 Θ [i] =< Model k , mve > 13 end 14 return Θ
In the last algorithm (Algorithm 3), we address the coldstart problem. To solve that, and taking into consideration that QoS values are highly influenced by geographical characteristics (as discussed in the end of Section 3.2.2), we compute the initial QoS values V (s) = {qos 1 , • • • , qos m } for each new service according to its geographical characteristics. Technically, we check if the actual provider of the new service (P(s)) belongs to the initial set of providers (of our dataset) (line 1). If so, in line 2, we compute the representativeness rate of the actual provider for each service cluster (α c,s , ∀c ∈ C S ). Afterward, we average that with respect to all cluster-heads (H (c), ∀c ∈ C S ). If not, we check if the AS number of the new service (A (s)) is in our initial set of AS numbers (line 3). If this is the case, we calculate the representativeness rate of the actual AS number for each services cluster, and we average the initial QoS values V (s) with respect to all cluster-heads. Otherwise, we proceed to compute the representativeness rate of the service country (C (s)) and we average its initial QoS values. In case where the geographical characteristics of the new service are not available, we search the provider having the nearest service (as mentioned in line 8), where D denotes a geodesic distance function, and we compute its QoS values as mentioned in line 9. This vector (V (s)) is then provided as an input to all the trained DAEs. Finally, we average the prediction results of the initial QoS value for the new service (line 11). Note that the same algorithm is used for a new user (V (u) = {qos 1 , • • • , qos n }), with the appropriate changes as well.

Algorithm 3: cold-start Algorithm

Inputs : s new Service C = {C 1 , • • • ,C max } list of services clusters H (C) = {H (C 1 ), • • • , H (C max )} list of cluster-heads Θ
best training models Data : P, A , C List of Providers, ASNs and Countries from dataset

1 if (P(s) ∈ P) then 2 V (s) = 1 ∑ c∈C α c,s ∑ c∈C α c,s H (c) with α c,s = 1 |c| |c| ∑ j=1 1, if P(s) = P(s j ) 0, otherwise 3 else if (A (s) ∈ A ) then 4 V (s) = 1 ∑ c∈C α c,s ∑ c∈C α c,s H (c) with α c,s = 1 |c| |c| ∑ j=1 1, if A (s) = A (s j ) 0, otherwise 5 else if (C (s) ∈ C ) then 6 V (s) = 1 ∑ c∈C α c,s ∑ c∈C α c,s H (c) with α c,s = 1 |c| |c| ∑ j=1 1, if C (s) = C (s j ) 0, otherwise
∀s i ∈S (D(s, s i ))) compute the provider of nearest service 9 V (s) = 1 ∑ c∈C α c,s ∑ c∈C α c,s H (c) with α c,s = 1 |c| |c| ∑ j=1 1, if P(s) = P(s j ) 0, otherwise 10 end 11 V (s) = 1 |Θ | ∑ Model ∈Θ
Model (V (s)) average prediction of the initial QoS values 12 return V (s) service initial vector Figure 4 shows an example that explains the principles of the proposed approach: The subplot (A) shows a l by k grid that covers the different clusters of the SOM map; Each line of the map represents a country and involves k neurons, and each neuron is a subset of services. The subplot (B) shows the initialization process of the neurons' weights of a given country, the neurons are assigned in a self-organized mode. The subplot (C) shows the training phase of the different DAEs that are related to the learned clusters. The subplot (D) details the cold start scenario.

Experimental Evaluation

In this section, we conduct experiments to evaluate the performance of the proposed approach and its different variants (deep autoencoder, clustered deep autoencoder and Top-k

- → n 1 • • • • • • • • • - → n i • • • • • • • • • - → n k
Compute the maximal geodesic distance:

D(S2, S4) = 25 km 1
Update the 1 st and k th neurons:

- → n 1 = - → S 2 - → n k = - → S 4 2 Compute the index of S5: i = D(S2,S5)×k D(S2,S5)+D(S4,S5) 3 Update the i th neuron: - → n i = - → S 5 4 (A) Grid of l × k neurons each line represents a country C 1 . . . . . . . . . C i . . . C l - → n 1 • • • • • • • • • - → n i • • • • • • • • • - → n k CLUSTER1 CLUSTER2 • • • • • • • • • CLUSTERn (C)
CLUSTER1 S P1 1 S P1 2 S P2 3 S P2 4 H (cl1)
. . .

S P1 new =? CLUSTERn S P2 5 S P1 6 S P2 7
H (cl1)
Compute the clusters heads:

H (cl1) = ( - → S P1 1 + - → S P1 2 + - → S P2 3 + - → S P2 4 )/4 . . . H (cln) = ( - → S P2 5 + - → S P2 6 + - → S P1 7 )/3 1
The representativeness rate of the ptovider P1:

αcl 1 = 1 + 1 + 0 + 0 = 2 . . . αcln = 0 + 1 + 0 = 1 2
Average initial vector for the new service:

V (Snew) = 1 (αcl 1 +•••+αcln ) αcl 1 H (cl1) + • • • + αcln H (cln) 3
Average initial vector for the new service:

V (Snew) = 1 n (Model1(V (Snew)) + • • • + Modeln(V (Snew))) 4
In this example, we assumed that the provider of S P1 new belongs to the initial set of our providers. Otherwise, we will follow the other steps mentionned in algorithm 3.

NB Fig. 4 Illustrative example describing the principles of the proposed approach clustered deep autoencoder). The objective is to answer the following questions: (1) How sensible are our proposed models to the hyper-parameters? (2) What is the impact of the geographical characteristics on the SOM clustering? (3) How do our approach compare to the state-of-the-art methods under different scenarios?

Experimental setup

The experiments were conducted on the MESO@LR-Platform of the University of Montpelier, France. To this end, we used 4 nodes (14 cores) with 128 GB of RAM. All the learning programs were implemented in Tensorflow Python. Note that, for the training of the SOM algorithm, we turn to a high-performance implementation called Somoclu 4 , which is considered as a massively parallel tool for training a batch formulation of self-organizing maps on large datasets [START_REF] Wittek | somoclu: An efficient parallel library for self-organizing maps[END_REF].

Data collection

To evaluate the proposed approach, we conducted experiments on a large-scale real-world Web service QoS repository named WS-Dream, released by [START_REF] Zheng | Investigating qos of real-world web services[END_REF]. This repository contains two datasets. Note that the dataset contains about 26% of missing values, which is a quite high date sparsity rate. This problem specific to this dataset may pose the risk to having an invalid clustering (a cluster with services with null values). In order to reduce this risk, we considered the following two rules for a given service invocation QoS matrix M c,t = {QoS u,s } m×n (for a selected criterion c and at given time slot t):

1. For each QoS c,t s,u value for a given invocation matrix M c,t , if this value is invalid, then we replace it by the average of the valid QoS c,t s,u values on all the previous matrices time slot: QoS c,t s,u =

1 t -1 t-1 ∑ t =1
QoS c,t s,u 2. For each QoS c,t s,u value for a given matrix invocation M c,t , if this value is invalid and all its previous values (regarding time slots) are invalid, then QoS c,t s,u = 0. Other changes are made on the first dataset, in order to reduce the missing values on the services attributes in terms of geolocational information (AS number, IP address, latitude and longitude). To do this, we used the geolocational information derived from GeoIP25 , and IP2Location6 databases.

Table 2 shows the different improvements operated on the initial dataset.

Evaluation metrics

To evaluate the performance of the method, we use two wellknown metrics frequently used in collaborative filtering: mean absolute error:

MAE V = 1 |V | ∑ (u,s)∈V
|r u,s -ru,s | and root mean squared

error: RMSE V = 1 |V | ∑ (u,s)∈V |r u,s -ru,s | 2
, where V represents the validation dataset. r u,s is the actual QoS score for service s given by user u, and ru,s the predicted one.

Since we used k-fold cross-validation, we assess the prediction on the average of MAE and RMSE too. The average MAE is defined as follows:

MAE AV G = 1 allV k-f old ∑ allV k-f old (|V k-f old | × MAE V k-f old ) (8)
The average RMSE (which is also denoted as mve in step 11 of Algorithm 2) is defined as:

RMSE AV G = 1 allV k-f old ∑ allV k-f old (|V k-f old | × RMSE V k-f old ) (9)
where V k-f old represents a validation set in both previous equations.

Effects of training parameters

Autoencoder models are full of hyper-parameters and setting a proper parameter initialization is a great challenge. This requires expertise and extensive trial and error. In our work, we focused on two important hyper-parameters: the size of each layer in the autoencoder and the choice of the activation function on these layers. In the following, we present a series of comparative experiments conducted in order to obtain the optimal parameters of our model.

Impact of autoencoder layer sizes

The size of the common layer (latent factor layer) is a hyperparameter that we must set before training the autoencoder. Like any feed-forward networks, it has been proven that deep autoencoders yield much better compression than corresponding shallow or linear autoencoders [START_REF] Hinton | Reducing the dimensionality of data with neural networks[END_REF]. So, additional layers can learn complex representations by approximating any mapping from input to code arbitrarily well. For this reason, and aside from the latent factor layer, we used a three layer net for both encoding and decoding networks. This choice is based on the power of three-layer neural network experimentally proven in [START_REF] Hecht-Nielsen | Theory of the backpropagation neural network[END_REF]. The final architecture of our autoencoder is composed of four fully connected layers (three layers + latent factor layer) for both encoding and decoding networks.

In order to determine the ideal number of neurons for each layer, many rules-of-thumb have been taken into account [START_REF] Panchal | Behaviour analysis of multilayer perceptrons with multiple hidden neurons and hidden layers[END_REF].

The code size should be: (1) between the size of the input layer and the size of the output layer; (2) two-thirds the size of the input layer, plus the size of the output layer;

(3) less than twice the size of the input layer. We considered these three rules as a starting point for determining the upper bound on this parameter. Figure 5(A) shows the results of four different values that operate on the layers sizes of the autoencoder. As we can see, the configuration that performs best on our 4-layer-model is 1024, 512, 256 and 128. Note that this configuration depends on the dimension of our input (with 4500 input services) for the autoencoder model. For the clustered autoencoder, this configuration is updated by the rule of three on the input size (which represents the cluster size).

Impact of activation functions

The second important hyper-parameter used in the training of our autoencoder is the target activation. To examine the influence of this parameter, we empirically used some of the most popular choices in deep learning. So, we used a sigmoid such as standard activation function and exponential linear units [START_REF] Clevert | Fast and accurate deep network learning by exponential linear units (elus)[END_REF] (ELU), standard rectified linear unit (ReLU) [START_REF] Sun | Deeply learned face representations are sparse, selective, and robust[END_REF] and finally, a leaky rectified linear unit (Leaky ReLU) [START_REF] Xu | Empirical evaluation of rectified activations in convolutional network[END_REF] known as non-standard activation functions. The effect of the different activation functions on RMSE training metric is shown in Figure 5(B). In order to accelerate the convergence speed we opted for a non-saturated activation function: Leaky ReLU.

Results and discussion

In deep learning, it has been experimentally proven that the classical model of autoencoders degenerate into identity networks and they fail to learn the latent relationship between data [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. To tackle this issue, we used a variant model of deep autoencoder by corrupting the inputs, training the model to denoise the final outputs.

To enhance the robustness of the model, we adopted a masking noise technique to set a random fraction of the input to zero. For this reason, we explore four rates of noise: 0% (i.e., no noise), 20%, 50% and 80%, applied on these models with the whole dataset in its input. Note that the same training process is applied on experiments for each clustered deep autoencoder after SOM clustering phase on all input data.

Deep autoencoder

To extract the latent factor from the dataset, we used a 4layer encoder network with the following configuration on the layers sizes: 1024-512-256-128. The same symmetric configuration is used for the decoder network to reconstruct the input data.

To ensure more robust prediction results, we lead our experiments with a multiple folders according to the k-fold cross-validation principle. Indeed, the experiments are performed using five partitions (folders) with different density values. Each partition represents a given percentage (density%) of the available dataset. The density represents the ratio between the training and the validation datasets. To this end, we explore three possibilities for density ratio: 20%, 50% and 80%. Finally, the validation results are averaged over the 5-folders through the application of equations 8 and 9. In Figure 6 (A), we illustrate the average RMSE/MAE after five runs operated on three different data densities. For each density value, we applied four noise variations on the input data. Note that the QoS prediction accuracy increases continuously when the matrix density increases (due to the increase of the noise rate).

In Figure 6 (B), we can see that the training converged for all the possibilities, after fewer than 100 epochs. Also, we can easily observe that there is a stability (no increase after converging), and no overshoot (no increase before converging).

Clustered deep autoencoder

In order to provide yet a better prediction accuracy, we used a Clustered Autoencoder. This autoencoder variant is designed after the partition of all services into eight clusters that are homogeneous with regard to their geographical characteristics. As mentioned above, we used the same architecture of our deep autoencoder (4-layer-model configuration: 4-layer sizes, and the activation functions are Leaky ReLU).

The main difference is that the size of the four layers that make up each clustered autoencoder is depending on the size of the inputs (number of services in the cluster). The layer sizes are updated by the rule of three on the initial configuration (1024-512-256-128) according to the cluster size (input size). To do so, we performed eight trainings for each cluster. We applied the same process on each cluster under various input data densities (20%, 50% and 80%). For each density case, we applied a set of noise variations (0%, 20%, 50% and 80%). Thereafter, we average the validation over those clusters. As shown in the results from Figure 7(A average RMSE/MAE for different densities increases continuously when the rate of corrupted input data increases (except for the case when density is 20%). In Figure 7 (B), we notice that the training converges a little slower than the deep autoencoder model.

SOM Initialization

This subsection demonstrates the performance of SOM as a method for geographical clustering instead of using other alternative methods like k-means with Euclidean distance, Hierarchical clustering, etc. As already noted, the final results of clustered autoencoder (see Figure 7 (A)), are the average result for the training of eight different autoencoders. Each autoencoder is built according to its input vector size (number of services). The histogram depicted in Figure 8 As we can see, the best results are for the clusters with the highest number of services (clusters with 1803 and 1749 services). Their results are far less than the average result. However, the third and the fourth ranked autoencoders are the clusters with the smallest number of services (60 and 73 respectively). For the remaining clusters, despite the fact that they contain a considerable number of services, we note that their results are significantly higher then the average result.

The major significance of these results is on the manner that we operate for the initialization of the code-book. In fact, we opted for a Kohonen feature map network with 80 × 80 neurons to perform the clustering of 4500 services. Each line on this map contains 80 neurons and represents one country from our dataset. Note that, in ws-dream dataset we detected 81 different countries (we merged two countries). The size of each neuron vector is equal to the number of users available in our dataset (i.e., 142 users). So, each line for a selected country i on the map represents semantically the maximal geodesic distance (according to their longitude/latitude values) between all services pairs which compose this country i (i.e., argmax(D(s 1 , s 2 )), ∀s 1 , s 2 ∈ Country i ), where D denotes a geodesic distance function. This means that the first neuron is initialized with the vector of s 1 divided by the size of all services which compose the selected country. The same thing is operated for the last neuron with the vector s 2 . The rest services are distributed on the rest of neurons according to their geodesic distance. Note that, if there are several services with the same geographic coordinates, we put on the corresponding neuron the average vector of those services. Finally, the services which do not have the geographical information (i.e., latitude/longitude) are randomly dispatched throughout the whole neurons of their corresponding country line.

In support of our initialization outcomes (Figure 8 (A)), other experiments were performed with a random initialization of the code-book in SOM algorithm. As we can see in Figure 8 (B)), the average results in terms of RMSE and MAE are less than the results of experiments operated on SOM clustering results based on Geographical characteristics initialization.

Other interesting observations on our proposed initialization method are: -Firstly, in terms of RMSE and MAE, we observe that about 80% of all services are clustered in the top-3 clusters for clustering based on Geographical characteristics initialization (figure 8 (A)), and only 69% for a random initialization (figure 8 (B)). -Secondly, we observe that there is a directly proportional relationship between the training results and the sparsity on data of each cluster. The more the density of input data is the better the results are. Table 3 summarizes the sparsity (column two) ratio on each cluster (column one), knowing that the sparsity of all data is about 23%. -The third important observation is about the percentage of missing geolocational information. We observed that most of services without geographical characteristics are clustered in the latest clusters. The third column in Ta-ble 3 shows the missing rate according to the clusters size. -Finally, with our proposed initialization mechanism, we ensured a deterministic behavior of clustering result.

Based on this analysis, and in order to preserve the topological properties of the dataset and reduce the effect of QoS sparsity, we take the advantage of SOM algorithm to create the appropriate clusters with our proposed initialization mechanism. Table 4 summarizes the results in terms of RMSE, MAE and TIME for the three architectures: A deep autoencoder, a clustered autoencoder and clustered* autoencoder with only top-3 clusters. As already mentioned, the results represent the comparison between the three different approaches that we have implemented, so, the results are grouped according to those methods. The Deep AE group concerns the results on test data after training the deep autoencoder without using any clustering. The Clustered AE group concerns the average results for eight clusters of the deep autoencoder. The third group (Clustered* AE) presents the average testing results on the top-3 best clusters (about 80% of all services are clustered in the top-3 clusters). All the experiments are conducted with three different data density values, and by adding a random noise to the input data for each of them. Note that the technique of adding the noise, as specified in [START_REF] Neelakantan | Adding gradient noise improves learning for very deep networks[END_REF], not only helps us to avoid over-fitting, but also can result in lower training loss (especially when density is about 20%). This clearly demonstrates the interest of using the clustering phase before. Indeed, the results show that the method with SOM clustering behaves better than the first method where the prediction of QoS values is performed on the whole dataset. Consequently, using deep-autoencoder on the top-3 clusters outperforms the other methods. We also note that along with the data densities increasing, the prediction metrics become smaller.

A detailed version of experiments operated on the whole dataset according to k-fold cross-validation principle can be found at the following links:

-Deep Autoencoder: http://bit.ly/2JrIGSs -Clustered deep autoencoder: http://bit.ly/2XKRMxO

We additionally observe, that not only the prediction accuracy decreases on the clustered AE, but the training time improves significantly from the deep autoencoder method to the clusterd AE method.

Performance Comparison

With the aim of evaluating the advantages of our method and its variants, we compare them with the following stateof-art baseline QoS prediction methods. UPCC [START_REF] Zhao | User-based collaborative-filtering recommendation algorithms on hadoop[END_REF]: Userbased CF using PCC. This method is a user-based model using Pearson Correlation Coefficient for recommendation and prediction of Web services. IPCC [START_REF] Sarwar | Itembased collaborative filtering recommendation algorithms[END_REF]: Item-based CF using PCC. This method uses similar services for the QoS prediction using Pearson Correlation Coefficient method. ARIMA [START_REF] Godse | Automating qos based service selection[END_REF]: this method is often considered as the baseline method. It is a statistical method adapted to QoS web service prediction. WSRec [START_REF] Zheng | Wsrec: A collaborative filtering based web service recommender system[END_REF]: It is a hybrid CF algorithm that combines a user-based prediction [START_REF] Zhao | User-based collaborative-filtering recommendation algorithms on hadoop[END_REF] model with item-based prediction model [START_REF] Sarwar | Itembased collaborative filtering recommendation algorithms[END_REF]. Lasso [START_REF] Wang | A spatial-temporal qos prediction approach for time-aware web service recommendation[END_REF]: this approach optimizes the recommendation problem by adapting the lasso penalty function. Country-clustered Autoencoder [START_REF] Smahi | An encoder-decoder architecture for the prediction of web service qos[END_REF]: this is our first approach to predict QoS of Web services based on historical data. We used a simple (not deep) autoencoder architecture (with only one hidden layer) to predict the QoS values on sets of clusters based on the country ID. For that, all the experiments for our previous approach have been reconducted in line with the actual execution parameters (number of clusters 7 and number of folders).

We note that our work is not compared with QoS Prediction methods that use a Topk based model (since they do not consider the entire dataset).

Table 5 presents the MAE and RMSE results of different prediction methods on response-time criterion when the training set densities take two different values: 80% and 50%. From these results, we notice the following observations:

1. The prediction accuracy of Lasso is better then Country AE, WSRec and ARIMA methods; 2. The Deep AE method is slightly better then Lasso method for the two density values. However, when the rate of corruption data increases the performance of our Deep AE decreases, compared to Lasso method; 3. Compared with Lasso method, the Clustred AE (when no noise) can obtain as high as 9.6% and 11% improvements in prediction accuracy when data density is 80% and 50% respectively. Furthermore, the performances are almost identical when the input noise is about 50%. 4. We clearly remark that the top-3 clustering deep autoencoders outperform the others methods in all cases (density and noise variations) even when the noise rate is around 80%. 5. Globally, when the noise is lower or equal to 20%, we observed that the results of the three variants of our proposed method outperform other prediction approaches in prediction accuracy in all cases in terms of RMSE metric (except for Deep AE when density is about 80% and noise is equal to 20%), and in the most cases from MAE metric. First, we use the geographical characteristics of a service to determine the most representative cluster (in terms of country, autonomous system, or provider); this step is performed by using the first part of Algorithm 3 (from line 1 to line 9). Then, we use the new QoS values (given by the first step) as an input of all trained DAEs and take the average of the returned results (line 11). Note that, in the learning phase, we use about 60% of entire dataset as examples to fit the parameters of the SOM classifier. In the testing phase, we use about 40% of examples to assess the performance of our fully-specified classifier.

As depicted in Table 6, the RMSE, MAE, and testing TIME, show the postive impact of the contextual characteristics on the resolution of the cold-start problem. For that, we are not only focused on the checking of the most represen- tative cluster in terms of geographical characteristics (step 1), but we also use already the trained DAE (after using Algorithm 2) to average the returned results. We show that the application of DAE (step 2), can significantly improve the results.

Threats to validity

Regarding the recommendations in [START_REF] Wohlin | Experimentation in software engineering[END_REF], there are several elements to discuss about the validity of the experiment. For that, we consider the threats to internal and external validity of our study:

-Internal Validity threats are intrinsically tied to experimental realism. To select the different parameters for our models, we have tried to be as objective as possible. However, despite the fact that we conduct several pre-trainings for selecting the best configuration of those parameters, other different values may conduct to different results. For the deep autoencoder model, since we split the dataset into training and testing sets with varying proportions and use the k-fold cross-validation, this can alleviate the actual threat. For the clustered deep autoencoder model, and in addition to previous remarks, internal validity threats are mitigated since the number of neurons is proportional with the size of clusters. -External Validity threats are about the ability to generalize our proper findings and conclusions to other contexts. In this study, we performed our experiments on a large dataset (ws-dream), which enables us to train correctly our prediction method. To the best of our knowledge, this is the only large Web service QoS dataset publicly available. Although we used only one dataset, we are quite confident that the impact of using another dataset is reduced with the use of k-fold cross-validation which operates on multiple data densities. On the other side, the fact of using a denoising autoencoders by corrupting the data by randomly forcing some of the input values to zero, decreases the risk of external validity.

Related Works

In [START_REF] Zhang | Deep learning based recommender system: A survey and new perspectives[END_REF][START_REF] Singhal | Use of deep learning in modern recommendation system: A summary of recent works[END_REF][START_REF] Jalili | Evaluating collaborative filtering recommender algorithms: A survey[END_REF], the authors draw up a detailed survey and new perspectives on "deep learning"-based recommendation systems. For the two first articles, the authors treat all categories of recommendation systems. However in the third article, the authors focus only on the Collaborative Filtering recommendation algorithms.

In this paper, we focused on the Collaborative Filtering methods using a deep learning mechanism and/or exploiting the context information (and more particularly geographical information) for clustering Web services. For this reason, we present recent related works according to these two different aspects:

Clustering-based Methods in Collaborative Filtering

In the literature, several clustering-based works have been carried out [START_REF] Tang | Collaborative web service quality prediction via exploiting matrix factorization and network map[END_REF][START_REF] Chen | Your neighbors alleviate coldstart: On geographical neighborhood influence to collaborative web service qos prediction[END_REF][START_REF] Chen | Exploiting web service geographical neighborhood for collaborative qos prediction[END_REF][START_REF] Chen | Your neighbors are misunderstood: On modeling accurate similarity driven by data range to collaborative web service qos prediction[END_REF][START_REF] Ma | Variation-aware cloud service selection via collaborative qos prediction[END_REF][START_REF] Liu | A personalized clustering-based and reliable trust-aware qos prediction approach for cloud service recommendation in cloud manufacturing[END_REF]. Those studies are trying to discover a set of clusters based on the neighborhood characteristics of data.

In [START_REF] Ma | Variation-aware cloud service selection via collaborative qos prediction[END_REF], the authors focus on Cloud Services Selection according to users' non-functional requirements. Their method is based on time series QoS data. In order to identify the user clusters, they use the double Mahalanobis distances to improve the similarity measurement of QoS cloud models during multiple periods. They assume two challenges: Challenge in exactly identifying the neighboring users for a current user and Challenge in selecting the appropriate cloud service with optimal QoS meeting user's period preferences.

In the work presented in [START_REF] Liu | A personalized clustering-based and reliable trust-aware qos prediction approach for cloud service recommendation in cloud manufacturing[END_REF], the authors assume that most existing CF approaches ignore the influence of task similarity among different users on QoS prediction results and vice versa (i.e. the influence of similar users on different tasks). To address this problem, the authors proposed a novel clustering-based and trust-aware method for personalized and reliable QoS values prediction. For that, they combine two contributions to make a more personalized QoS prediction. For the first one, they develop a clustering-based algorithm to identify a set of similar users from the point of view of task similarity. The tasks similarities are computed by incorporating both explicit textual information and rating information as well as implicit context information. For the second one, they made assumptions that the QoS values may be contributed by unreliable users. For that, they design a trust-aware CF approach by merging local and global trust values, to reconstruct trust network of the clustered users. A series of experiments on two real-world datasets were conducted to evaluate the proposed approach.

The work presented in [START_REF] Chen | Exploiting web service geographical neighborhood for collaborative qos prediction[END_REF] aims to improve QoS prediction accuracy. To this end, the authors have considered the factor of geographical neighborhood to perform both matrix factorization and service clustering. The authors built a set of service clusters based on geo-neighbors similarity. The clusters are identified using both a bottom-up hierarchical neighbor-hood and the contextual attributes (like the country and service provider). Accordingly, the authors leveraged a neighborhood-based term to decompose the QoS matrix and predict the missing values.

Deep Learning Methods in Collaborative Filtering

Recently, deep learning-based models are rapidly developed in the era of recommendation. A wide range of Collaborative Filtering approaches have been proposed like in [START_REF] Kanagawa | Cross-domain recommendation via deep domain adaptation[END_REF][START_REF] Xiong | Deep hybrid collaborative filtering for web service recommendation[END_REF][START_REF] Bai | Dltsr: A deep learning framework for recommendations of long-tail web services[END_REF][START_REF] Wang | Online reliability prediction via motifs-based dynamic bayesian networks for service-oriented systems[END_REF]. We focus next on two works that we consider as very related to our work.

In [START_REF] Xiong | Deep hybrid collaborative filtering for web service recommendation[END_REF], the authors proposed a hybrid deep learning approach by combining Matrix Factorization with Content filtering to improve Web service recommendation. To predict the probability of a particular service to be invoked by a particular web application or mashup (user in our case), they used two feed-forward neural networks for both Collaborative filtering and Content filtering components. After that, these two components are combined by concatenating their last hidden layers to compose a third multi-layer perceptron to learn the interaction between them as well as their functionalities (it is well known that the interaction treated with Collaborative filtering and Content filtering is about multiple factors such as invocation history and functionalities). In order to evaluate their method, the authors conducted several experiments using a real-world Web service dataset. However, despite the fact that for those experiments the training data is randomly selected under different percentages of training data (20% to 80% with a step of size 10%), they kept only 20% in test data for each experiment round. The second remark is about the lack of information on their choices of the number of layers and the size of each of them and how those choices can affect the final results. A final remark is about the learning time, which is relatively high due to the fact that they try to learn three different neural networks.

In [START_REF] Bai | Dltsr: A deep learning framework for recommendations of long-tail web services[END_REF], the authors used an autoencoder architecture to learn low-dimensional representations and perform feature extraction. To this end, they build a deep learning framework to recommend long-tail Web services. To tackle the problem of the unsatisfactory quality of the low-quality of the content description, they used a stacked denoising autoencoder to perform feature extraction for long-tail Web services with severe sparsity on both their content description and their historical usage data. The proposed approach is tested on a real-world dataset and compared with several state-of-art baselines. In summary, this framework can be considered as a content-based recommendation system. However, from the viewpoint of matrix factorization and collaborative filtering paradigm, the interaction between the applications (mashups) and services is not addressed here.

Positioning of the approach

In contrast to the aforementioned works, our contribution leverages both contextual attributes and QoS data to perform clustering. By doing so, we limit the chances of having sparse or small-sized clusters. More specifically, the initialization of the cluster-heads (or weights of neurons) is performed using contextual attributes; in addition, this initialization is accomplished in a self-organized manner (i.e., the neighboring neurons of the map are also neighbors in the physical network). Besides that, the training of the SOM model exploits the QoS data (instead of contextual attributes) to update the weights of the cluster-heads.

To the best of our knowledge, our contribution is the first work that designs a deep autoencoder architecture (and its variant denoising autoencoder) for QoS prediction Web services.

Conclusion

QoS is a key factor for building successful service-oriented applications. Most existing QoS-predicting works do not handle data sparsity, the cold-start problems, and the contextual information of services and thus are likely to perform worse in real scenarios.

In this paper, we tackled the aforementioned issues by adopting a DAE architecture. To optimize prediction capabilities, the hyper-parameters of the DAE are tuned through cross-validation. We also introduced a SOM-based clustering step as a pre-processing to the DAE learning. This step aims to reduce QoS matrix sparsity and improve prediction accuracy. We conducted a set of experiments to evaluate the performance of this prediction method. We explored many variants of this method and compared their performance with the state-of-the-art methods. These experiments showed that our method outperforms the existing methods in terms of accuracy.

As perspectives to this work, we would like to investigate first the use of alternative architectures of deep learning models like stacked autoencoders or deep belief networks. These architectures may ensure more accurate prediction results than the autoencoders used in this work. Second, we plan to use deep learning in predicting Web service Quality of Experience (QoE). We intend to use some powerful text mining tools to analyse the users' feedbacks in order to estimate QoE scores for Web services. We argue that the combination of QoS and QoE will contribute to improve Web service recommendation systems.

Fig. 1

 1 Fig. 1 Deep Autoencoder architecture.

Fig. 2

 2 Fig. 2 Kohonen feature map network. (The connections between the input and output layers indicate the relationships between the input and output vectors.)

Algorithm 1 : 3 foreach s ∈ S do 4 compute

 134 Data Clustering According to Batch SOM Algorithm Inputs : S = {s 1 , • • • , s n } , set of n service vectors Data : W = {w 1 , • • • , w k } set of services weight vectors 1 Initialize all w ∈ W 2 for epoch : 1 → N epochs do

7 else 8 P

 78 (s) = P(argmin

Fig. 5

 5 Fig. 5 Training RMSE per minibatch on the second dataset, with learningrate = 0.05, epoch = 20, for an optimal parameters selection. (A) Influence of layers size, (B) influence of activation functions.

Fig. 6

 6 Fig.6with a training/testing set of (density%)/(100 -density%) for a deep autoencoder prediction (on various noises and densities). (A) The average RMSE and MAE results, and (B) The training performance graph.

Fig. 8

 8 Fig. 8 Detailed results of clustered autoencoder with density = 80% and noise = 0%. clustering results: (A) based on Geographical characteristics initialization. (B) with random initialization

  QoS prediction performance

					(D) Cold-start process
		2 Training the best model for each cluster	Compute the initial QoS values for the new service S P1 new =?
	build model for cluster1 1	x 1 ENCODER1	DECODER1	x1	for model1 3 Optimize the squared error loss
	build model for cluster2 1	x 2 ENCODER2	DECODER2	x2	for model2 3 Optimize the squared error loss
		. . .
		. . .
	1 build model for clustern	xn ENCODERn	DECODERn	xn	3 Optimize the for modeln squared error loss

Table 1

 1 Table 1 summarizes some characteristics of WS-Dream repository data. Information details of ws-dream repository (RT : response Time, TP: Throughput Time )

	The first one contains 1 974 675 in-
	vocations from 339 users on 5825 Web services. The sec-
	ond one contains about 30 million Web service invocation
	records collected from invocation of 4500 services by 142

users at 64-time intervals. Each time interval t takes 15 minutes. The innovation, on both datasets, considers response 4 https://github.com/peterwittek/somoclu time and throughput criteria (c ∈ {rt,th}).

Table 2

 2 The results of different improvements on dataset

	Attributes	% of missing values (before)	% of messing value (after)
	QoS values	26%	23%
	AS number	21%	15%
	IP Address	24%	13%
	Latitude/Longitude	22%	15%

  (A) shows the detailed results of training autoencoders for each cluster separately. In fact, this graph represents the results in terms of RMSE and MAE metrics of eight autoencoders. The average RMSE and average MAE are plotted with blue and red lines respectively. Note that the x-axis corresponds to the size of input data for those autoencoders.

	2.33 3 4 5	1.577	1.819	3.697	2.060	3.781	2.665	(A)	4.409	3.056	4.567	2.752	5.101	3.619	5.762	4.303	u	2.217	2.543	3.697	2.115	3.881	2.635	(B)	4.500	3.120	4.803	2.903	6.204	4.492	6.500	4.866	3.23 4 5 u
	1.19	0.590	0.756															0.913	1.104														1.78
		1803	1749	60	73		196	356	99	164		1557	1382	144	420		50	552	91	304
						RMSE		MAE												RMSE		MAE

Table 3

 3 Information details of ws-dream repository

	Cluster size % of sparsity % missing geoloc. info.
	1803	22,34%	13,81%
	1749	23,10%	9,43%
	60	19,30%	20,00%
	73	24,44%	13,70%
	196	28,06%	32,14%
	356	26,63%	33,99%
	99	22,53%	15,15%
	164	26,69%	24,39%

Table 4

 4 Users accuracy comparison methods on RT criterion (training time is specified in minutes)

	Method	Noise	Density=20% RMSE MAE TIME RMSE MAE TIME RMSE MAE TIME Density=50% Density=80%
		0%	3,318	1,819	121,2	2,761	1,407	144,5	2,562	1,270	194,2
	Deep AE	20% 50%	3,229 3,554	1,717 1,972	97,98 92,56	2,819 3,026	1,501 1,681	128,9 115,8	2,616 2,807	1,367 1,511	167,6 143,4
		80%	4,090	2,337	94,92	3,575	2,057	108,4	3,310	1,894	131,1
		0%	2,852	1,571	53,74	2,558	1,383	69,29	2,325	1,188	85,34
	Clustered AE	20% 50%	2,765 2,977	1,514 1,647	54,16 54,34	2,581 2,743	1,398 1,501	70,20 70,30	2,379 2,577	1,236 1,380	86,26 86,12
		80%	3,173	1,724	53,97	2,981	1,612	69,68	2,862	1,527	85,40
		0%	1,908	0,840	42,20	1,762	0,771	54,70	1,700	0,748	66,90
	Clustered* AE	20% 50%	1,997 2,087	0,907 0,936	42,25 42,32	1,806 1,917	0,802 0,858	55,10 55,20	1,731 1,829	0,771 0,822	67,90 68,00
		80%	2,224	0,956	42,31	2,128	0,945	55,10	2,071	0,930	67,50

Table 5

 5 Performance Comparisons of Prediction methods on RT with two data density values

	Method	Density=50% RMSE MAE RMSE MAE Density=80%
	UPCC		3,034	1,470	3,032	1,467
	IPCC		2,951	1,396	2.925	1.372
	ARIMA		3,401	1,236	2,986	1,028
	WSRec		2,945	1,380	2,925	1,372
	Lasso		2,872	1,021	2,572	0,893
	Country AE		3,825	1.892	2,803	1,250
		noise = 0%	2,761	1,407	2,562	1,270
	Deep AE	noise = 20% noise = 50%	2,819 3,026	1,501 1,681	2,616 2,807	1,367 1,511
		noise = 80%	3,575	2,057	3,310	1,894
		noise = 0%	2,558	1,383	2,325	1,188
	Cluster AE	noise = 20% noise = 50%	2,581 2,743	1,398 1,501	2,379 2,577	1,236 1,380
		noise = 80%	2,981	1,612	2,862	1,527
		noise = 0%	1,762	0,771	1,700	0,748
	Cluster* AE	noise = 20% noise = 50%	1,806 1,917	0,802 0,858	1,731 1,829	0,771 0.822
		noise = 80%	2,128	0,945	2,071	0,930
	4.4 Cold-start situation
	To handle the cold-start problem, we adopt the strategy de-
	scribed in algorithm 3. More specifically, we estimate the
	initial QoS values of a new service (or new user) as follows:

Table 6

 6 Cold-start problem (time is specified in seconds)

	Cold-Start Steps RMSE MAE TIME
	Step 1	3.877	1.951	1.201
	Step 2	3.010	1.420	2.194

The autoencoder is called linear if the transfer functions are linear, otherwise it is non-linear. x dim y dim . . .neuron i i n p u t l a y e r

An AS is either a single network or a group of networks that is controlled by a common network administrator (or group of administrators) on behalf of a single administrative entity (university, a business enterprise, etc.)

Maxmind GeoIP2 Geolocational Databases. Retrieved on May 2019 from http://dev.maxmind.com/geoip/geoip2/geolite2/

IP2Location LITE Databases. Retrieved on May 2019 from http://lite.ip2location.com

For this work, the clustering is performed on country ID. In order to have only eight clusters we grouped some of countries in the same cluster

Acknowledgment

This work has been performed with the support of the High Performance Computing Platform MESO@LR, financed by the Occitanie / Pyrénées-Méditerranée Region, Montpellier Mediterranean Metropole and the University of Montpellier, France.


Aymen Laadhari
email: aymen.laadhari@ku.ac.ae

Ahmad Deeb

Numerical Approach Based on the Composition of One-Step Time-Integration Schemes For Highly Deformable Interfaces

In this work, we propose a numerical approach for simulations of large deformations of interfaces in a level set framework. To obtain a fast and viable numerical solution in both time and space, temporal discretization is based on the composition of one-step methods exhibiting higher orders and stability, especially in the case of stiff problems with strongly oscillatory solutions. Numerical results are provided in the case of ordinary and partial differential equations to show the main features and demonstrate the performance of the method. Convergence properties and efficiency in terms of computational cost are also investigated.

INTRODUCTION

We are interested in numerical solutions of strongly coupled systems of PDEs involving highly nonlinear membrane forces, see e.g. [START_REF] Laadhari | An operator splitting strategy for fluid-structure interaction problems with thin elastic structures in an incompressible Newtonian flow[END_REF][START_REF] Laadhari | Exact Newton method with third-order convergence to model the dynamics of bubbles in incompressible flow[END_REF]. For such problems with high nonlinearities, fine meshes and spatial discretizations by highorder finite elements are required [START_REF] Doyeux | Simulation of two-fluid flows using a finite element/level set method. Application to bubbles and vesicle dynamics[END_REF]. Therefore, numerical time integrations with higher orders are also necessary to allow viable numerical approximations. In this preliminary work, we present the formalism in the case of the level set problem, where we seek to produce solutions with a higher-order time approximations.

Different temporal integration techniques with high order approximations exist to study dynamical systems. Popular numerical methods are those associated with discrete flow, where solutions are sought over discrete times using a numerical flow composition [START_REF] Ernst Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF]. Fixed or adapted time steps were used. In the existing literature, several studies on continuous flows have investigated the Borel-Padé-Laplace integrator for stiff and non-stiff problems. See, for example, [START_REF] Ernst Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF][START_REF] Deeb | Comparison between Borel-Padé summation and factorial series, as time integration methods[END_REF]. Indeed, an integrator calculates an approximation over a continuous time interval and has the ability to increase the order of approximation as much as necessary by changing a single parameter in the integrator [START_REF] Deeb | Performance of Borel-Padé-Laplace integrator for the solution of stiff and non-stiff problems[END_REF].

To allow computational savings and improve precision while maintaining numerical accuracy, adaptive timestepping strategies are commonly employed where the level of adaptability may depend on the rate of change between consecutive solutions [START_REF] Guo | A novel adaptive crank-nicolson-type scheme for the time fractional allen-cahn model[END_REF] or on some error estimates [START_REF] Ernst Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF][START_REF] Ernst Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF]. However, the computational efficiency depends on the order of the scheme and how the error estimate is computed. Hence the need for high-order discrete flow to improve numerical efficiency, especially for stiff problems. For a given numerical approximation, the convergence depends on the error in time and in space. Although the adaptive meshing techniques allow a better precision of the numerical solution, the global error estimate will be limited by the order of the numerical scheme in time. Thus, it is important to increase the order of approximation in time to allow high-order of precision.

One can distinguish one-step and multi-step methods based on the number of previous numerical solutions used in the temporal discretization of differential equation. For a detailed discussion of the classical theory of multi-step methods as well as the properties of order, convergence, stability and symmetry, we refer the reader to [START_REF] Ernst Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF][START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF][START_REF] Hairer | Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]. Among the most used schemes for stiff problems with oscillatory terms, we can mention the implicit linear multi-step methods based on the backward difference formulas (BDF) presenting a stability limited to lower orders [START_REF] Iserles | A First Course in the Numerical Analysis of Differential Equations[END_REF][START_REF] Curtiss | Integration of stiff equations[END_REF][START_REF] Butcher | Coefficients for the study of Runge-Kutta integration processes[END_REF][START_REF] Butcher | A history of Runge-Kutta methods[END_REF]. Indeed, only first-order (implicit Euler) and second-order schemes are stable. The BDF schemes of order 3 to 6 exhibit a weaker stability property, called A-stability, and formulas of order greater than 6 are unstable. A generalization of BDF using a second derivative allows to obtain implicitly stable schemes up to the order 10, see e.g. [START_REF] Ernst Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF]. Runge-Kutta methods, referred to as RK, are one-step schemes [START_REF] Runge | Ueber die numerische Auflösung von differentialgleichungen[END_REF]. To cope with the lack of A stability of explicit RK methods, which are not suitable for solving stiff equations, implicit and stable RK methods are also developed. We refer the reader to [START_REF] Ernst Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems[END_REF][START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF] for a note on A-stability and the Gauss, Radau IIA and Lobatto IIIA, IIIB and IIIC schemes. Composition methods have emerged in the literature, with several variants, as a powerful numerical tool based on the composition of lower-order basic time-integration methods with the aim of increasing the order of approximation.

The aim of this work is to provide a higher-order approximation of the level set solution. Various higher-order finite element approximations are used and allow to increase the spatial precision. Since the error in time takes over in case of low-order schemes, composition methods are applied on simple and low-order discrete flows to increase the order of approximation in time, while keeping a reasonable computational cost. Numerical examples are presented to validate and show the main characteristics of the method in the cases of the Lotka-Volterra differential system and the level set advection problem. This paper is organized as follows. Section 1 presents the preliminary concepts and mathematical setting. In Section 2, we present the composition methods for the implicit first-order accurate backward Euler and second-order accurate Heun's method. In Section 3, we provide an assessment of the composition technique through numerical examples.

MATHEMATICAL FORMULATION

Level set problem

Let T > 0 and d = 2 be respectively the simulation period and the spatial dimension. For any time t ∈ (0, T ), the interface is denoted by Γ(t). Its deformations are described implicitly as the iso-value zero of a level set function ϕ. The interface is assumed smooth enough and is enclosed in a larger domain Λ so that

Γ(t) = (t, x) ∈ (0, T ) × Λ : ϕ(t, x) = 0 .
Given a velocity field u, the evolution of the level set function satisfies the advection equation:

∂ t ϕ + u • ∇ϕ = 0, in (0, T ) × Λ. (1)
Appropriate initial and boundary conditions are considered: ϕ = ϕ b on (0, T ) × Σ -and ϕ(0) = ϕ 0 in Λ, where

Σ -= {x ∈ ∂ Λ : u • ν(x) < 0}
is the inflow boundary and ν denotes the unit normal exterior vector to Λ. For an initial internal domain Ω 0 , ϕ is initialized as a signed distance function to Γ(t = 0) so that |∇ϕ 0 | = 1:

ϕ 0 (x) = inf{|y -x|; y ∈ ∂ Ω 0 }, if x / ∈ Ω 0 , -inf{|y -x|; y ∈ ∂ Ω 0 }, if not.
However, the signed distance property is not preserved when solving (1), thus leading to numerical instabilities if the gradient of the level set function becomes either very small, or very large, in particular in the vicinity of the interface or on ∂ Λ [START_REF] Osher | Level Set methods: An overview and some recent results[END_REF]. To restore the signed distance property, an auxiliary redistancing problem, is commonly solved while maintaining the zero-level set position to avoid the loss of mass characterizing Eulerian methods. Higher-order methods and the use of fine meshes in the vicinity of the interface can help enforcing the local mass conservation, see for example [START_REF] Doyeux | Simulation of two-fluid flows using a finite element/level set method. Application to bubbles and vesicle dynamics[END_REF][START_REF] Laadhari | Implicit finite element methodology for the numerical modeling of incompressible two-fluid flows with moving hyperelastic interface[END_REF][START_REF] Laadhari | Fully implicit finite element method for the modeling of free surface flows with surface tension effect[END_REF].

Time integration by composition method

Let y(t) be the solution of the following initial value problem:

y = f (t, y (t)) , for all t ∈ R, with the initial condition y(0) = y 0 . (2)
For t > 0, the exact flow of ( 2) is defined by the map Y t as follows:

Y t : R → R y 0 → Y t (y 0 ) = y(t). (3)
For any given numerical scheme of order p, we can associate a numerical flow, denoted

Φ [p]
∆t , such that a sequence of numerical values y n is constructed on a discrete set of points t n = n∆t as the approximations of the solution at time t n :

y n+1 = Φ [p]
∆t y n . We say that the numerical flow is of order p if the following equality holds:

Y ∆t (y 0 ) -Φ [p] ∆t (y 0 ) = O(∆t p+1 ). The superscript [p] in Φ [p]
∆t refers to the order p of the numerical flow. As stated in [START_REF] Hairer | Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] (Theorem 4.1, page 43), we can increase the order of the approximation using a s times composition of the same numerical flow

Φ [p] ∆t . The resulting discrete flow Φ [p+1] ∆t = Φ [p] a 1 ∆t • Φ [p] a 2 ∆t • . . . • Φ [p]
a s ∆t is at least or order p + 1 if the conditions hold:

a 1 + a 2 + • • • + a s = 1 and a p+1 1 + a p+1 2 + . . . a p+1 s = 0. ( 4
)
There are no real solutions for system (4) when p is odd. However, if p is even, the smallest composition that increase the order by 1 is when s = 3. This composition is called the the triple jump, see [START_REF] Hairer | Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] (page 44, section II.4). One can also continue to increase the order of construction of the discrete flow by continuing to compose discrete flows, see [START_REF] Muruq | Order conditions for numerical integrators obtained by composing simpler integrators[END_REF]. However, a problem for compositions with real coefficients comes from negative coefficients which result in a negative time step. Castella et al. consider solutions for the system with complex coefficients in the case of parabolic equations [START_REF] Castella | Splitting methods with complex times for parabolic equations[END_REF]. In the following, we will use the composition technique in the case s = 2.

Let us consider any discrete flow Φ

[p] ∆t of order p and introduce the following composition:

Φ [p+1] ∆t = Φ [p] a [p] 2 ∆t •Φ [p] a [p] 1 ∆t
.

We recall a useful theorem which states that (see e.g. [START_REF] Casas | Compositions of pseudo-symmetric integrators with complex coefficients for the numerical integration of differential equations[END_REF]):

Theorem 1 Let us consider a differential system of type (2) and any numerical flow

Φ [p] ∆t of order p. If a [p] 2 = a [p] 1 = 1 2 + i 2 sin 2l + 1 p + 1 π 1 + cos 2l + 1 p + 1 π , with    - p 2 l p 2 -1 if p is even, - p + 1 2 l p -1 2 if p is odd.
(

) 5
then the approximation defined by the real part of the output, denoted by Re(.) and obtained by the composition of flow, is an approximation of the solution of order p + 1, i.e.:

Y ∆t (y 0 ) -Re Φ [p] a [p] 2 ∆t • Φ [p] a [p] 1 ∆t (y 0 ) = O ∆t p+2 . (6)
First, we consider the numerical flow Φ BE1 of the backward Euler scheme which is of order one. We define the two time composition as follows:

Φ BE2 ∆t := Φ BE1 a [1] 2 ∆t • Φ BE1 a [1] 1 ∆t . (7)
Here, the coefficients a [START_REF] Laadhari | An operator splitting strategy for fluid-structure interaction problems with thin elastic structures in an incompressible Newtonian flow[END_REF] i are given by formula (5) for l = 0: a

[1] 1 = a [1] 2 = 1 2 + i 2 .
Hence, the approximation given by the real part of the output obtained by the composition ( 7) is an approximation of the solution of order 2, i.e.:

Y ∆t (y 0 ) -Re Φ BE2 ∆t (y 0 ) = O ∆t 3 .
Thus, the double composition of the discrete scheme BE1 of order 1 gives a scheme of order 2. Second, we consider a numerical flow of order 2 associated to the Heun's method Φ HM1 ∆t (second-order RK):

Y ∆t (y 0 ) -Φ HM1 ∆t (y 0 ) = O ∆t 3 ,
thus the composition method:

Φ HM2 ∆t := Φ HM1 a [2] 2 ∆t • Φ HM1 a [2] 1 ∆t
, with 6a

[2] 1 = 6a [2]
2 = 3 + i √ 3 leads to a numerical flow of order 3 and we have:

Y ∆t (y 0 ) -Re Φ HM2 ∆t (y 0 ) = O ∆t 4 .
We can continue increasing the order of the approximation by constructing new numerical flows using the above composition technique. This was thoroughly analysed in [START_REF] Casas | Compositions of pseudo-symmetric integrators with complex coefficients for the numerical integration of differential equations[END_REF] and presented with a comprehensive study. In addition, it is demonstrated that the numerical flow 1 2 Re Φ

[p]

a [p] 1 ∆t • Φ [p] a [p] 2 ∆t + Φ [p] a [p] 2 ∆t • Φ [p] a [p]
1 ∆t is of order p + 2. When composing two discrete flows with permuted coefficients, the imaginary parts are conjugated and are therefore canceled by the summation, resulting in only a real part. Remark that the A-stability of the resulting composition scheme is not studied here but we refer the interested reader to the reference [START_REF] Fedoseev | Preference and stability regions for semi-implicit composition schemes[END_REF]. Roughly speaking, the stability function S r (z) of the resulting composition scheme is the product of two functions: S r (z) = S(a 2 z) S(a 1 z), where S(z) represents the stability function of the initial diagram. Therefore, the region of stability, defined by {z ∈ C, such that |S r (z)| 1}, can be represented as the following intersection:

{z ∈ C : |S(a 2 z)| 1} ∩ {z ∈ C : |S(a 1 z)| 1}.
Note that each set in this intersection has a wider region than the original scheme with a slight translation in the complex plane.

NUMERICAL APPROXIMATION

Space discretization by finite elements. Let T h be a partition of Λ consisting of geometrically conforming open simplicial elements K (triangles for d = 2), such that Λ = ∪ K∈T h . Lagrange finite element polynomials are considered for the space discretization of the level set function on a mesh of size h = max K∈T h diam(K). Let us denote by ϕ h an approximation of ϕ, while X h = ψ h ∈ C 0 (Λ) : ψ h|K ∈ P κ , ∀K ∈ T h , with κ 1 represents the finite dimensional space of admissible level set. The Galerkin scheme associated with (1) consists in finding ϕ h ∈ X h such that ϕ h (t, x) = ∑ i y i (t) ψ h,i (x), for a certain finite element basis.

A semi discrete equation is then obtained by multiplying by a test function and integrating aver Λ. We employ the Stabilization Upwind Petrov Galerkin, referred to as SUPG [START_REF] Loch | The level set method for capturing interfaces with applications in two-phase flow problems[END_REF], method where the SUPG test function is defined as v h,i = ψ h,i +τ K u•∇ψ h,i in order to add a diffusion in the streamline direction. We set a streamline diffusion parameter proportional to the local mesh size so that

τ K = C h K max {|u| 0,∞,K , tol/h K } ,
where C represents a scaling constant and the term tol/h K avoids the division by zero. The semi-discrete system arising from the discretization with finite elements writes:

M(t) y + K(t) y = 0, with y = (y 1 , . . . , y i , . . .) T . (8)
Here, M(t) is the stabilized mass matrix, while K(t) represents the transport matrix, both time-dependent; y is the unknown degrees of freedom for the level set field. In the sequel, we omit the subscript h referring to the spatial approximation to alleviate the notations.

Time advancing scheme. Let us divide [0, T ] into N subintervals t n ,t n+1 of constant time step ∆t with n = 0, ..., N -1. For n > 0, the unknown ϕ n h approximates ϕ at t n . Considering the differential system (8) and for given degrees of freedom, y n represents the vector field which approaches ϕ at time t n . To increase the accuracy in time, we focus on composition methods for certain numerical methods: the backward Euler and Heun's methods, because of their simplicity and accuracy, designated respectively by the acronyms BE and HM. On the one hand, we approximate the solution at time t n using the backward Euler scheme as follows:

Φ BE1 ∆t y n-1 := y n = M(t n ) + ∆tK(t n ) -1 M(t n-1 )y n-1 .
On the other hand, we use the one-step Heun's method as stated by the following numerical flow:

Φ HM1 ∆t y n-1 := y n = I - ∆t 2 [M (t n )] -1 K t n-1 + [M(t n )] -1 K (t n ) I -∆t [M(t n )] -1 K t n-1 y n-1 . Φ BE4 Φ BE2 Φ BE1
u(t) v(t) 0 1 2 u(t) 0 1 2 v(t)

NUMERICAL EXAMPLES Example 1: Validation in the case of Lotka-Volterra equations

We apply the above numerical flow compositions to the Lotka-Volterra system of first order nonlinear differential equations describing the evolution over time of two species: predators and prey. Dynamically speaking, if we assume that predators and prey are not connected, then the number u of prey increases exponentially, with a factor α, and v for predators decreases exponentially, with a factor δ . The Lotka-Volterra system describes the predator-prey interactions, when connected, as follows:

u (t) = α u(t) -β u(t)v(t) v (t) = -δ v(t) + γ u(t)v(t)
, with initial conditions u(0) = u 0 , v(0) = v 0 .

Here, β represents the rate of prey mortality by predators and γ is the growth rate of predators while eating prey. We present in Fig. 2 the evolution of the numerical solution of the Cauchy problem (u(0) = 2, v(0) = 1) in the interval

[0, 55] with α = δ = 2 3 , β = 4 3
and γ = 1. This problem has an invariant prime integral F t , which is therefore identified using the initial conditions, i.e. F t u(t), v(t) = F(u 0 , v 0 ), for any time t > 0. The invariant is given by:

F t u(t), v(t) = β v(t) + γu(t) -α log(v(t)) -δ log(u(t)).
To study the convergence properties of the numerical strategy, we calculate the error between the approximate invariant and the reference invariant at the initial time F 0,ref (u 0 , v 0 ). The error quantification corresponds to the time evolution of the quantity F t (u(t), v(t)) given by:

e ∆t l = ∑ N n=1 F 0,ref (u 0 , v 0 ) -F t (u(t n ), v(t n )) 2 ∑ N n=1 |F 0,ref (u 0 , v 0 )| 2
The parameter ROC represents the convergence rates necessary to set a computational effort to establish a certain precision, in which ∆t is the time step in a given time discretization level l. To investigate the accuracy of different composition methods, we run a series of simulations using the methods of Backward Euler and Heun and their recursive composition for several time steps. The table [START_REF] Laadhari | An operator splitting strategy for fluid-structure interaction problems with thin elastic structures in an incompressible Newtonian flow[END_REF] and the figure 1 report the errors (9) over a complete period for different compositions of numerical flows. The estimated convergence orders are reported in the ROC column, showing that we retrieve the expected theoretical convergence orders of the basic schemes. Moreover, the increase in order by double composition of flows is clearly depicted. Thereafter, we study the computational efficiency by evaluating the CPU times for different time steps and composition methods. Thus, a correlation can be established between the error sought by the user and the CPU time that the numerical flow needs to carry out the simulation, for the two discrete schemes BE and HM and their compositions. Results in Fig. 3 show that the resulting composition scheme allows the user to produce numerical results with a lower computational time than the original scheme, while keeping the same precision. This has been observed for different composition levels, for a given original basic scheme (Heun or Backward Euler method). Note that we do not compare in this work such compositions with other schemes having the same convergence orders.

CPU time [s]

Φ BE1 Φ BE2 = Φ BE1 • Φ BE1 Φ BE4 = Φ BE2 • Φ BE2

CPU time

[s] Φ HM 1 Φ HM 2 = Φ HM 1 • Φ HM 1 Φ HM 4 = Φ HM 2 • Φ HM 2

Example 2: Level set problem by composition methods

In this example, the aforementioned composition of one-step time integration schemes is used for simulations of largely deformable interfaces in a level set framework. The numerical computations are performed using the opensource software FEniCSx [START_REF] Alnaes | The FEniCS project version 1.5[END_REF]. We first consider a reversible vortex test case and evaluate the properties of convergence with respect to the temporal discretization. The computational domain is Λ = [0, 1] 2 . An initially circular interface of radius R = 0.15 is centered at (0.7, 0.7) and is periodically stretched into thin filaments by a vortex flow field. The interface unrolls and reaches its maximum deformations at t = T /2, before resuming its circular shape at time T = 4. The advection velocity is given by: u(t, x) = -2 sin(πx) 2 sin(πy) cos(πy) cos(πt/T ) , 2 sin(πy) 2 sin(πx) cos(πx) cos(πt/T ) T , with x = (x, y) T ∈ Λ. Convergence history for various time steps and time composition schemes.

e ∆t Φ BE1 Φ BE2 = Φ BE1 • Φ BE1 Φ BE4 = Φ BE2 • Φ BE2
Some snapshots showing the interface deformations are reported in Fig. 4. We investigate the order of accuracy with respect to the time discretization using a relatively simple composition method. Let π h be the Lagrange interpolation operator. We compute the errors in the L 2 (Λ) norm for different time step sizes with respect to an exact reference solution π h ϕ at time t = T after one deformation period. The convergence of the composition technique is assessed using the implicit backward Euler scheme, as well as two-and four-time backward Euler compositions, respectively called Φ BE2 and Φ BE4 . Finally, we assess the robustness of our level set solver in the case of the standard Zalesak's rotating disk test, using a one-time composition of the backward Euler scheme. We consider the same setting presented in [START_REF] Ta | An implicit high order discontinuous Galerkin level set method for two-phase flow problems[END_REF], where the slotted disk recovers the initial position after a period T = 4. We choose ∆t = 10 -3 and h = 10

CONCLUSION

We have presented an application of numerical methods of flow composition to the simulation of the level set problem in order to increase the accuracy of the solution and to capture the dynamics of stiff problems. A validation is presented in the case of Lotka-Volterra differential equations. By using simple and lower-order numerical schemes, we CONCLUSION have shown numerically that composition techniques allow to increase the accuracy of the approximation while presenting a lower computational cost. This is part of an ongoing work on modeling the dynamics of highly-deformable biomembranes [START_REF] Laadhari | Fully implicit methodology for the dynamics of biomembranes and capillary interfaces by combining the Level Set and Newton methods[END_REF][START_REF] Gizzi | A three-dimensional continuum model of active contraction in single cardiomyocytes[END_REF][START_REF] Kolahdouz | A numerical model for the trans-membrane voltage of vesicles[END_REF], where temporal integration schemes with high orders are required. It is also planned to implement the method for fluid-structure interaction problems with large structural deformations [START_REF] Laadhari | Eulerian finite element method for the numerical modeling of fluid dynamics of natural and pathological aortic valves[END_REF].

FIGURE 1 :

 1 FIGURE 1: Error versus the time step size for the composition of backward Euler scheme (left) and Heun's scheme (right).

FIGURE 2 :

 2 FIGURE 2: Numerical solution of Lotka-Volterra equations for the first five periods obtained with the BE4 scheme.

FIGURE 3 :

 3 FIGURE 3: CPU time versus error for different compositions of backward Euler (left) and Heun (right) schemes.

FIGURE 4 :

 4 FIGURE 4: (Left) Snapshots showing interface deformations at t ∈ {0, 0.6, 1, 2, 2.6, 3.4, 3.8, 4}, for h = 10 -2 .Convergence history for various time steps and time composition schemes.

  Fig 4 (right) reports the convergence of the computed errors versus the time discretization, showing in particular an increase in the order of convergence by composition of discrete flows.

  -2 . A few snapshots of the rotating disk are shown in Fig 5, showing good preservation of the sharp angles.

FIGURE 5 :

 5 FIGURE 5: Zalesak's rotating disk. Snapshots showing interface deformations at t ∈ {0, 0.6, 1, 2, 2.6, 3.4, 3.8, 4}.

TABLE 1 :

 1 Time convergence history for different time integration composition methods.

	∆t	e [BE1]	ROC e [BE2] ROC e [BE3]	ROC e [HM1] ROC e [HM2] ROC e [HM4] ROC
	4.17E-1 2.04E-1	4.125 1.022	--2.148E-1 --3.899E-3 --2.051E-1 --6.650E-3 --2.888E-4 --1.955 5.231E-2 1.979 3.761E-4 3.277 4.476E-2 2.132 6.021E-4 3.365 1.809E-5 3.881
	1.01E-1 4.156E-1 1.279 1.270E-2 2.013 4.002E-5 3.186 1.046E-2 2.067 6.338E-5 3.201 1.129E-6 3.945
	5.03E-2 1.899E-1 1.122 3.113E-3 2.013 4.595E-6 3.100 2.527E-3 2.034 7.253E-6 3.105 7.056E-8 3.971
	2.51E-2 9.101E-2 1.057 7.698E-4 2.008 5.503E-7 3.051 6.211E-4 2.017 8.667E-7 3.054 4.412E-9 3.985
	1.25E-2 4.457E-2 1.028 1.914E-4 2.005 6.733E-8 3.025 1.540E-4 2.009 1.059E-7 3.027 2.758E-10 3.992
	6.25E-3 2.206E-2 1.014 4.770E-5 2.002 8.328E-9 3.013 3.833E-5 2.004 1.309E-8 3.014 1.725E-11 3.996
	3.13E-3 1.097E-1.007 1.191E-5 2.001 1.035E-9 3.006 9.561E-6 2.002 1.627E-9 3.007 1.084E-12 3.999
	1.56E-3 5.473E-3 1.003 2.975E-6 2.001 1.291E-10 3.003 2.388E-6 2.001 2.027E-10 3.003 6.774E-14 3.999
	7.81E-4 2.733E-3 1.002 7.434E-7 2.000 1.612E-11 3.001 5.966E-7 2.001 2.531E-11 3.002 1.385E-14 --3.91E-4 1.36E-3 1.001 1.858E-7 2.000 2.013E-12 3.002 1.491E-7 2.000 3.156E-12 3.003 2.510E-14 --1.95E-4 6.827E-4 1.000 4.645E-8 2.000 2.231E-13 3.173 3.727E-8 2.000 3.710E-13 3.088 5.621E-14 --9.73E-5 3.400E-4 1.000 1.152E-8 2.000 6.210E-14 1.835 9.245E-09 2.000 4.395E-14 3.060 8.845E-15 --Order of convergence 1 --2 --3 --2 --3 --4

ACKNOWLEDGMENTS

The authors acknowledge financial support from KUST through the grant FSU-2021-027.