Nageswara S V Rao

Neena Imam
email: nimam@ornl.gov

Zhengchun Liu
email: zhengchun.liu@anl.gov

Rajkumar Kettimuthu
email: kettimut@anl.gov

Ian Foster
email: foster@anl.gov

Machine Learning Methods for Connection RTT and Loss Rate Estimation Using MPI Measurements Under Random Losses *

Keywords: Round Trip Time, Loss Rate, Message Passing Interface, Machine Learning, Generalization Bounds, Regression, Information Fusion

Scientific computations are expected to be increasingly distributed across wide-area networks, and Message Passing Interface (MPI) has been shown to scale to support their communications over long distances. Application-level measurements of MPI operations reflect the connection Round-Trip Time (RTT) and loss rate, and machine learning methods have been previously developed to estimate them under deterministic periodic losses. In this paper, we consider more complex, random losses with unform, Poisson and Gaussian distributions. We study five disparate machine leaning methods, with linear and non-linear, and smooth and non-smooth properties, to estimate RTT and loss rate over 10Gbps connections with 0-366ms RTT. The diversity and complexity of these estimators combined with the randomness of losses and TCP's non-linear response together rule out the selection of a single best among them; instead, we fuse them to retain their design diversity. Overall, the results show that accurate estimates can be generated at low loss rates but become inaccurate at loss rates 10% and higher, thereby illustrating both their strengths and limitations.

Introduction

Computations distributed across geographically dispersed facilities, such as multiple supercomputer sites connected over a wide-area network, are of increasing interest in science applications. Their execution times are effected by network latencies and loss processes, often in a complex way, due to the close coupling between computations and communications in these applications. Recently, Message Passing Interface (MPI) has been shown to be effective in supporting communications over wide-area connections, including ones long enough to span the globe, under external packet loss rates up to 20% Table 1: RMS errors of RTT estimation by individual and fuser methods. [START_REF] Rao | A case study of MPI over long distance connections[END_REF]. In contrast with executions at a single facility, these distributed computations need to account for the longer and more varied executions times of MPI operations to avoid inefficiencies due to unbalanced computing and networking operations; for example, MPI join operation over connections with wide ranging latencies will be delayed by the longest. Motivated by such considerations, Round Trip Times (RTT) and loss rates of wide-area connections are estimated using execution time measurements of MPI primitives in distributed computations [START_REF] Rao | Estimation of rtt and loss rate of wide-area connections using mpi measurements[END_REF]. A main contributor to these execution times is the Transmission Control Protocol (TCP) which is a dominant underlying transport mechanism of MPI for wide-area connections. In particular, at increased loss rates and randomness, the execution time variations are dominated by TCP's highly non-linear response dynamics [START_REF] Hassan | High Performance TCP/IP Networking: Concepts, Issues, and Solutions[END_REF][START_REF] Lakshman | TCP/IP performance with random loss and bidirectional congestion[END_REF].

Machine Learning (ML) methods have been developed for a number of networking tasks for science data flows, for example, detecting flow anomalies [START_REF] Giannakou | Flowzilla: A methodology for detecting data transfer anomalies in research networks[END_REF] and classifying elephant and mice flows [START_REF] Chhabra | Classifying elephent and mice flows in high-speed scientific networks[END_REF]. In particular, ML methods are developed to estimate the connection RTT and loss rate under deterministic periodic losses in [START_REF] Rao | Estimation of rtt and loss rate of wide-area connections using mpi measurements[END_REF] for 10Gbps emulated connections with 0-366ms RTT; these connections represent local, cross country, continental and round the earth distances. In this paper, we consider more realistic, complex scenarios with random losses, in particular under unform, Poisson and Gaussian distributions up to 20% loss rates, to study the strengths and limitations of ML methods. We study five disparate ML methods, with linear and non-linear, and smooth and non-smooth properties, to estimate connection RTT and loss rate. They include four non-linear estimators, namely, smooth Support Vector Machine (SVM) and Gaussian Process Regression (GPR), and non-smooth Ensemble of Trees (EOT) and Regression Trees (RT), in addition to the baseline Linear Regression (LR) method. The diversity and complexity of these estimators combined with the randomness of losses and TCP's non-linear response rule out the identification of a single best among the estimators. Analytical results establish the finitesample limits in asserting the performance superiority of any such method based on samples [START_REF] Devroye | A Probabilistic Theory of Pattern Recognition[END_REF]. In particular, the training error is an insufficient indicator of estimator's performance due to potential over-fitting that leads to poor generalization performance on future datasets.

Over-fitting is often specific to an estimator method and is less likely to occur across estimators of radically different designs. In several cases, by fusing diverse estimators both the performance and diversity of design are preserved [START_REF] Rao | On fusers that perform better than best sensor[END_REF]. However, the fused estimators are also subject to finite sample limits since they are also estimators. We study linear regression fusion (LR-F) and GPR fusion (GPR-F) methods, and our results show that the latter achieves lowest Root Mean Square Error (RMSE) among all estimators for RTT in three out of four scenarios. We develop analytical characterization of the performance improvements of fused estimates over individual RTT estimates under finite sample, distribution-free framework [START_REF] Vapnik | Statistical Learning Theory[END_REF].

By using MPI execution times as the independent variable, we formulate the problem of estimating RTT and loss rate as a regression estimation problem. The overall results are illustrated using RTT estimates with smallest RMSE among individual methods and fusers in Figure 1 for four loss rates. In each plot, datasets are concatenated at four loss rates in increasing order and at each loss rate RTT is increased left to right, and measurements are repeated 10 times at each RTT value as shown in Figure 2. Among individual RTT estimates, GPR has the lowest RMSE in all four scenarios, and GPR-F fuser achieved even lower RMSE in three out of four loss scenarios, as shown in Table 1 while encompassing the design diversity of individual methods. Overall, our results show that accurate estimates can be generated at low loss rates but become inaccurate at loss rates 10% and higher, wherein the datasets appear much too complex for these methods (as in the case of deterministic periodic losses [START_REF] Rao | Estimation of rtt and loss rate of wide-area connections using mpi measurements[END_REF]). In addition, they reveal some subtle performance effects including over-smoothing by some estimators in achieving lower RMSE, and bleeding effects of RTT in loss rate estimates.

The organization of this paper is as follows. The testbed used in collecting MPI execution time measurements is described in Section 2. An analytical formulation of the underlying regression problem is presented in Section 3. Various datasets of execution time measurements are qualitatively described in Section 4. RTT estimators are described in Section 5, wherein five different ML methods are described in Section 5.1 and two fusers are described in Section 5.2. Generalization equations of the fusers for RTT estimation are described in Section 6. Loss rate estimators are described in Section 7. The performance of the estimators is qualitatively interpreted in the context of datasets at lower and higher loss rates in Section 8. A summary of results and directions for future work are described in Section 9. A computing cluster with InifiniBand (IB) interconnect is expanded to constitute a testbed to run MPI codes across the wide-area Ethernet connections. Additional Ethernet Network Interface Cards (NIC) are installed in two cluster computing nodes (tait1 and tait2), which are connected to Ethernet switches and a hardware-based Ethernet emulator in the configuration shown in Figure 3. The IB connections of the cluster are subject to 2.5 ms latency limit, and hence MPI measurements over IB are not indicative of the performance over long distance connections. Specifically, the shorter distances combined with credit-based IB protocol flow control do not adequately reflect the complex variations of TCP over wide-area connections, particularly under packet losses. Furthermore, due to their latencies, wide-area networks are more prone to more losses compared to IB networks.

Test Configuration

Typical wide-area connections consist of a number of switches and routers whereas IB connections have fewer IB switches. This testbed connection consists of two Ethernet switches between the source computing node and a port of the emulator, which reflects a site connection. Similarly at the other end, the connection consists of two Ethernet switches between the second port of the emulator and destination computing node. Thus, this symmetric end-to-end connection consists of four Ethernet internal cross-connections, six short Ethernet segments and one emulated long distance Ethernet connection with variable latency, loss rate and type of loss distribution.

ANUE/Ixia hardware-based emulator is used to collect MPI measurements over Ethernet connections with 11 Round Trip Times (RTT) in 0-366 ms range. These RTT values are strategically chosen to represent three ranges: (a) smaller values represent cross-country connections, for example, computing facilities distributed across the US, (b) 93-183 ms represent inter-continental connections, and (c) 366 ms represents a connection spanning the globe, which is mainly used as a limiting case. External periodic and random packets losses are introduced using ANUE/Ixia devices at four different loss rates. These emulators delay the packets as per the specified RTT value, and thus closely emulate the physical long distance path. Equally importantly, these emulations closely match TCP dynamics of physical connections with corresponding RTTs, which is a critical factor in assessing MPI performance over long distance connections. In particular, these emulations lead to different TCP dynamics and responses under deterministic (periodic) and random losses of Ethernet segments [START_REF] Rao | A case study of MPI over long distance connections[END_REF], which result in a wider spread of the execution times at high loss rates under random losses (Section 8).

Analytical Formulation

We now provide a formal description of the underlying estimation problems to support subsequent analytical treatment of RTT and loss rate estimation methods. Let E be a random variable representing the execution time of MPI Send Receive primitive; it is distributed according to the joint probability distribution P E,R,L , where R and L are the random variables representing RTT and loss rate, respectively. In general, the distribution P E,R,L is quite complex since it depends on the properties of the network connection and host systems, and also the software stack consisting of the operating system, networking and MPI modules. Given an execution time measurement E = e, the conditional distribution P R,L|e = P R,L|E=e , characterizes the distribution of RTT and loss rate at this value e. Then, RTT-regression function f R is defined as the expected value of RTT at E = e given by f R (e) = RdP R,L|e , which is averaged over both R and L at each e. The loss-regression function f L is the expected value of the loss rate at E = e, which is similarly given by f L (e) = LdP R,L|e .

In general, these regressions cannot be obtained even in theory since the underlying distribution P E,R,L is unknown. In stead, ML methods are employed to estimate their approximations using a training sample (E i , R i , L i ), i = 1, 2, . . . , l, wherein E i is the execution time measured over a connection with RTT R i and loss rate L i . The distributions of the connection parameters R and L are determined by the design of connection configurations, and are fixed while the measurements of E are repeated. Thus, the distribution of E encompasses factors due to the properties of physical connection parameters as well as operating system, TCP and MPI modules.

Then, RTT and loss rate estimation problems can be cast as estimating the regression functions f R and f L , respectively, using measurements. We consider that RTTregression estimate f R A is obtained by method A ∈ A = {EOT,GPR,LR,RT,SVM,LR-F,GPR-F} using the measurement pairs (E i , R i ), i = 1, 2, . . . , l. Similarly, the loss-regression estimate f L A is obtained by method A ∈ A using the measurement pairs (E i , L i ), i = 1, 2, . . . , l. At a given execution time E = e, f R A (e) and f L A (e) are the estimates of RTT and loss rate, respectively, provided by method A.

Execution Time Measurements

The execution times of MPI Send Receive operations collected at the application-level are shown as a function of RTT in Figure 4 for loss rates, 0.1, 1, 10 and 20%, of externally induced losses under four loss scenarios, one deterministic periodic and three random, namely uniform, Poisson and Gaussian. Their increasing trend as a function of RTT is evident at lower loss rates, 0.1% and 1%, but it becomes less prominent at 10% loss rate, and essentially disappears at 20% loss rate as outliers dominate. Overall, the execution times as well as their variations increase as loss rate is increased, which is an indication of the increased complexity of their estimation at higher loss rates.

In terms of losses, the execution times are shown as function of loss rates 0.1, 1, 10 and 20% in Figure 5 under the four loss scenarios. The measurements at any loss rate encompass all 11 RTT values, and 10 repeated measurements at each RTT value. Their increasing trend as a function of loss rate is evident overall but is sharper for periodic losses and is more diffused with overlaps across loss rates in all three random loss scenarios. The ranges of execution times are much wider for random losses compared to periodic losses. Also, the execution times as well as their variations increase overall as loss rate is increased for deterministic periodic loss scenario. But, in random loss scenarios the variations are more subtle: their spread is similar at all loss rates except at 20%, wherein a few measurements are very large, which indicate the complexity of loss rate estimation in these scenarios.

RTT Estimators

We present RTT estimates in the form of traces that are indexed by groups of 440 measurements that correspong to increasing loss rates, and within each group we have 11 sub-groups that correspond to increasing RTT values as shown in Figure 2; each sub-group corresponds to 10 repeated measurements at a fixed pair of loss rate and RTT. The 440 measurements for each loss scenario shown in Figure 6 will be used to compare qualitatively with the corresponding RTT estimate traces. We utilize the regression estimation codes from matlab statistics toolbox.

Five Estimators

The five estimation methods are chosen to reflect different characteristics of the underlying regressions, namely, smooth and non-smooth functions, respectively. GPR and SVM with Gaussian kernels [START_REF] Scholkopf | Advances in Kernel Methods[END_REF] are based on non-linearly transforming the feature space X into regression space of Y . They both provide smooth regression functions f GPR and f SVM , and their respective function classes F GPR and F SVM consist of smooth functions as a result of Gaussian kernels. EOT [START_REF] Breiman | Random forests[END_REF][START_REF] Hastie | The Elements of Statistical Learning: Data Mining, Inference, and Prediction[END_REF] method is based on boosting of a collection of classification trees that are customized to fit the training data using the AdaBoost method. RT [START_REF] Breiman | Classification and Regression Trees[END_REF] methods is also based on trees that are customized to fit the training data. They both lead to a highly non-smooth regression functions f EOT and f RT , and their function classes F EOT and F RT consists of a collection of decision tree. LR is a smooth and linear method and leads to f LR from the function classes F LR , which is effective in RTT estimation under no losses [START_REF] Rao | Estimation of rtt and loss rate of wide-area connections using mpi measurements[END_REF] but is quite limited under losses as indicated by its RMSE in Tables 1 and2.

The estimators under periodic, uniform, Poisson and Gaussian loss scenarios are shown in Figures 7,8, 9 and 10, respectively. Under periodic and Gaussian losses, all estimates are more accurate at 0.1% and 1% loss rate but are inaccurate at higher loss rates; in particular, they capture the increasing trends in RTT at low loss rate but exhibit high variation at 10% and 20% loss rate. Under unform and Poisson losses, GPR method does not capture the increasing RTT trend at any loss rate, while other non-linear estimates captured it. Interestingly, GPR achieved lowest RMSE among individual estimators which is due to the inclusion of measurements at higher loss rate that resulted in "averaging" across all loss rates. This undesirable artifact of low RMSE but less accurate estimate at low loss rates is illustrated in Figures 8 and9. LR and SVM methods have highest and second highest RMSE among the twenty cases in Table 1.

Estimator Fusers

The estimators from five individual methods are used as 5-dimensional input to a fuser which produces RTT as its output. The linear regression fuser (LR-F) is a linear combination of the individual estimators, and the GPR fuser (GPR-F) is obtained using GPR method based on the outputs of five estimators corresponding to the training sample. GPR-F achieved lower RMSE than best individual estimator GPR in all except under Gaussian losses, whereas LR-F has lower RMSE for Poisson and uniform losses. As shown in Figure 1, the fused estimates were able to capture the increasing RTT trend at lower loss rates while achieving lower RMSE error than GPR under uniform and Poisson losses, unlike GPR estimator with lowest RMSE among individual estimators.

Generalization Equations for Fused Estimates

We consider five individual estimates, indexed by A ∈ A I = {EOT,GPR,LR,RT,SVM}, such that the fuser input vector X consists of five real-valued components, X A , A ∈ A I , and output Y is a real-valued estimate of RTT. RMSE values in Table 1 have been used to compare the performance of fusers and individual estimators in previous sections, which are subject to statistical variations since they are computed based on a sample. We now derive confidence bounds for these RMSE values which provide analytical justification for their use. For simplicity of presentation of analytical results, we use MSE in place of RMSE following the common practice in finite sample theory [START_REF] Vapnik | Statistical Learning Theory[END_REF].

Regression Problem: Finite Sample Generalization

In a generic regression estimation problem the feature vector X ∈ ℜ d and the output vector Y ∈ ℜ are distributed jointly accordingly to an unknown distribution P X,Y . The expected error of a regression function f is

I( f ) = ( f (X) -Y ) 2 dP X,Y .
The expected best regression estimator f * minimizes I(.) over F , i.e., I( f * ) = min

f ∈F I( f ).
The empirical error Î( f ) based on training data (X i ,Y i ), i = 1, 2, . . . , l, is defined as

Î( f ) = 1 l l ∑ i=1 ( f (X i ) -Y i ) 2
It is an approximation of I( f ) computed based on the training data. The empirical best regression estimator f minimizes Î(.) over F , i.e., Î( f

) = min f ∈F Î( f ).
The joint distribution P X,Y of data is complex, domain specific, and is only partially known. In our context, it depends on the finer details of the underlying software and hardware components, which may manifest as additional random variables. For an individual estimator A ∈ A I , X and Y correspond to execution time E and RTT R, respectively, and P X,Y corresponds to P E,R,L which involves additional random variable of the loss rate L. For fusers, X and Y correspond to 5-dimensional vector consisting of outputs of estimators and RTT R, respectively. In general, an optimal f * cannot be computed precisely with probability one even in principle, since P X,Y is either unknown or not computationally conducive. Under certain conditions, Vapnik's generalization theory [START_REF] Vapnik | Statistical Learning Theory[END_REF] establishes that there exists a confidence function δ (.) such that for a "suitable" estimator f obtained from training data we have

P l X,Y I( f ) -I( f * ) > ε < δ (ε, ε, l) (1)
where ε, ε > 0, 0 < δ < 1, and Î( f ) = min f ∈F Î( f ) + ε. This condition ensures that "error" of f is within ε of optimal error (of f * ) with probability 1δ , irrespective of the underlying measured or computed data distribution P l X,Y . Furthermore, under these conditions, the confidence parameter δ (ε, ε, l) approaches 1 as the sample size l approaches infinity.

Consider the fuser class F F used in fusing the estimators f A ∈ F A , A ∈ A I . Let f F denote the regression function obtained by composing f A 's with the fuser function from F F . The error reduction ∆ F of the fused estimate over the best individual classifier is defined as

∆ F = min A∈A I I( f A ) -I( f F ).
Then, if F F has the isolation property [START_REF] Rao | Measurement-based statistical fusion methods for distributed sensor networks[END_REF], then ∆ F ≥ 0. The best error reduction is given by

∆ * F = min A∈A I I( f * A ) -I( f * F ).
and its estimate based on a sample is given by ∆F = min

A∈A I Î( f A ) -Î( f F ).
The error reduction values ∆ F based on measurements are shown in Table 1 for the two fusers LR-F and GPR-F. GPR-F has positive ∆F values in three scenarios indicating that the fused estimate has lower RMSE than the lowest of its constituent estimators (namely, GPR). LR-F das positive ∆F values in two scenarios, which might be attributed to the lack of the required statistical independence in estimator outputs. We show in the next section that the estimate ∆F reflects the optimal improvement ∆ * F achievable by the fuser within a formal framework.

Estimator Fusers: Generalization Equations

The generalization bound δ (ε, ε, l) applicable to five individual estimators can be derived using various properties of the corresponding estimator classes [START_REF] Rao | Finite-sample generalization theory for machine learning practice for science[END_REF]. In particular, these bounds for GPR and SVM with Gaussian kernels could be based on fat-shattering index [START_REF] Scholkopf | Advances in Kernel Methods[END_REF], and for EOT and RT they may be based on bounded total variation [START_REF] Anthony | Neural Network Learning: Theoretical Foundations[END_REF]. In Theorem 1, we assume that these generalization bounds are available from existing works, and their detailed derivations are beyond the scope of this paper.

We now show that the estimate ∆F is within ε of the optimal ∆ * F with a probability that improves with the training data size l independent of the underlying distribution P Y,X .

Theorem 1. Consider that there exists δ B (ε, εB , l) such that based on i.i.d. l-sample, we have

P l X,Y I( fB ) -I( f * B ) > ε < δ B (ε, εB , l) . (2)
for all individual estimators B ∈ A I , N A I = |A I |, and both fusers B = LR-F, GPR-F such that δ B (ε, εB , l) → 0 as l → ∞. Then, the probability that the closeness between ∆F and ∆ * F is within ε is bounded as

P l X,Y | ∆F -∆ * F | < ε > 1 -δ D (ε/2, εD , l) -∑ A∈A I δ A ε/(2N A I ), εA , l ,
for both fusers D = LR-F, GPR-F.

Proof. We first note that for

D = LR-F, GPR-F | ∆F -∆ * F | ≤ | Î( fD ) -I( f * D )| + min A∈A I Î( fA ) -min A∈A I I( f * A ) ,
which establishes that the condition | ∆F -∆ * F | > ε implies that at least one term on the right hand side is greater than ε/2. We now have

| Î( fD ) -I( f * D )| ≤ | Î( fD ) -I( fD )| + |I( fD ) -I( f * D )|,
which in turn establishes that the condition | Î( fD

) -I( f * D )| > ε/2
implies that at least one term on the right hand side is greater than ε/4. Then, by hypothesis in Eq (2), both conditions are simultaneously satisfied with probability at most δ D (ε/4, εd , l). Similarly, we have In each case, the condition that left hand side is larger than ε/2 implies at least one of the terms under the summation is greater ε/(2N A I ). Under the hypothesis of this theorem in Eq (2), both conditions are satisfied with probability at most

∑ A∈A I δ A (ε/(2N A I ), εa , l).
By combining the above terms together, we have

P l X,Y | ∆F -∆ * F | > ε < δ D (ε/2, εD , l) + ∑ A∈A I δ A ε/(2N A I ), εA , l ,
which proves the theorem.

The confidence bound in this theorem is distribution-free in that it does not depend on P X,Y . It is expressed in terms of the precision parameter ε and the confidence

parameter 1 -δ D (ε/2, εD , l) -∑ A∈A I δ A ε/(2N A I ), εA , l ,

Loss Rate Estimators

The loss estimates of four non-linear estimators are shown in top left plot of Figure 11 for periodic losses. For random losses, SVM estimates have extreme variations and hence are omitted in the plots. Also, LR estimator is omitted in all plots due to its extremely large variation under all loss scenarios. Qualitatively, smooth SVM and nonsmooth RT methods both exhibit large variations, which indicate the underlying properties of the data rather than these methods; indeed GPR is the only method that did not produce large variations. RMSE of five estimators and linear fuser LR-F are shown in Table 2. Methods with lowest RMSE for loss rate estimation are shown in Figure 12, which are EOT for periodic losses and GPR for all others. Both of them exhibited lower variations at low loss rate and an increasing trend as RTT is increased at a fixed loss rate. As in the case of RTT estimation, the "averaging" by GPR resulted in lower RMSE but less accurate estimates at low loss rates; but, interestingly, this effect is more dominant for Gaussian errors unlike for RTT estimation. In almost all cases, at fixed loss rate, the estimators showed an increasing trend as RTT is increased.

The large variations shown in the scatter plots in Figure 5 indicate high RMSE by any estimate since its output is a function and data dispersed around it contributes to RMSE. In particular, a smooth estimate will not be able to capture these variations as applicable to GPR and SVM methods. While tree-based methods in principles can capture such variations, they require a large number of leaf nodes, and those with smaller number will result in large RMSE. In summary, the results indicate the challenging nature of the underlying datasets, which in some sense expose the limitations of the conventional ML approaches for loss rate estimation.

Execution Time Measurements: Data Regressions

Qualitative insights into the performance of regression estimators can be gained by examining the scatter plots of measurements separately at low and high loss rates. Overall increasing trend of RTT when plotted as a function of execution time is evident under 0.1% and 1% loss rates, but not under and 10% and 20% loss rates in all four scenarios, as shown in Figure 13 for periodic and unform losses, and in Figure 14 for Poisson and Gaussian losses. Even GPR-F estimate with the lowest overall RMSE captures the increasing trend only in the former case but not in the latter case. Qualitatively, the wide spread of measurements at high loss rates indicates the lack of information needed to estimate RTT by any method that uses regression function, smooth or non-smooth. For loss rates, the scatter plots are shown in Figure 15, which have significant variations in execution times at fixed values of loss rate. The non-smooth EOT method with lowest RMSE under periodic losses captures several loss rates at 10 and 20% loss rates as shown in top left plot. The smooth GPR estimator is plotted along with data in all three random loss scenarios in which it has lowest RMSE; several estimated points are in between the loss rate values, and the estimator shows a continuous trend in the mapping from measurements to loss rate. Similar to RTT estimates, an overall increasing trend of loss rate when plotted as a function of execution time is evident under 0.1% and 1% loss rate; but, the loss rate is fixed at two values 0.1 and 1% as shown in Figure 16 for periodic and unform losses, and in Figure 17 for Poisson and Gaussian losses. This is a "bleed over" artifact due to increasing RTT at both 0.1 and 1%loss rates. Under 10% and 20% loss rates in all four scenarios there is a significant scatter in loss estimates, indicating the underlying complexity of regression estimation.

Conclusions

Rich datasets of MPI measurements are becoming increasingly available as more and more computations are distributed over wide-area networks. These measurements exhibit certain characteristics, such as longer executions times and large variations, that are atypical of conventional MPI applications executed on single computing systems with Inifiniband or custom interconnects. Losses are integral to wide-area networks as TCP that supports MPI utilizes self-induced losses to pace its flows. Consequently, these distributed computations need to mitigate the inefficiencies due to network de- lays and their variations. These computations may be distributed across geographically dispersed nodes that are dynamically identified; consequently, the RTTs and loss rates of the underlying connections may not be a priori known. The MPI measurements collected at the application-level reflect the connection length and losses, and have been shown to be useful in estimating RTT and loss rate using ML methods, albeit accurately only at low loss rates.

Complementing previous works under deterministic periodic loss scenarios, we studied five ML methods to estimate the connection RTT and loss rates under random losses, which are more reflective of practical scenarios. As in previous works [START_REF] Rao | Estimation of rtt and loss rate of wide-area connections using mpi measurements[END_REF], the results show that accurate estimates can be generated at low loss rates but they become inaccurate at loss rates 10% and higher. However, this randomness manifests in subtle ways, resulting in different performances of non-linear estimators; in particular, GPR that achieves low RMSE does not provide accurate RTT estimates at low loss levels, unlike others with higher RMSE. These effects are mainly due to the highly non-linear response of the underlying TCP dynamics that "amplify" the randomness of losses. Furthermore, it is equally complex to assess the performance of ML methods due to their non-linear nature, and their fusers are only effective in some scenarios for RTT estimation. In another direction, these results highlight the strengths and limitations of ML methods for network-level estimation problems using application-level measurements.

Fig. 1 :

 1 Fig. 1: RTT estimates with lowest RMS error of individual and fuser methods. loss type EOT GPR LR RT SVM LR-F GPR-F ∆LR-F ∆GPR-F periodic 109.20 87.78 102.27 95.31 103.49 92.94 85.50 -1.63 2.28 poisson 101.89 91.90 104.32 104.31 120.35 89.85 85.69 2.05 6.21 gaussian 84.33 73.13 2020.55 73.28 137.29 88.22 87.06 -15.06 -13.92 uniform 109.11 99.59 106.47 107.62 121.81 90.12 84.90 9.46 14.69

Fig. 2 :

 2 Fig. 2: Index representing increase of loss rate and RTT from left to right.

Fig. 3 :

 3 Fig. 3: Configuration for long Ethernet connection between compute nodes of IB cluster.

Fig. 4 :

 4 Fig. 4: Execution times of MPI Sendrecv operations as function of RTT under four external loss scenarios, namely, periodic, uniform, Poisson and Gaussian.

Fig. 5 :

 5 Fig. 5: Execution times of MPI Sendrecv operations as function of loss rate under four external loss scenarios, namely, periodic, uniform, Poisson and Gaussian.

Fig. 6 :

 6 Fig. 6: Traces of execution times in seconds of MPI Sendrecv operations under four external loss cases.

Fig. 7 :

 7 Fig. 7: Periodic losses.

Fig. 8 :

 8 Fig. 8: Uniform losses.

Fig. 9 :

 9 Fig. 9: Poisson losses.

Fig. 10 :

 10 Fig. 10: Gaussian losses.

  fA )min A∈A I( f * A ) , which in turn establishes that the condition min A∈A I Î( fA )min A∈A I I( f * A ) > ε/2 implies that at least one term on the right hand side is greater than ε/4. Then, we consider the two upper bounds min A∈A I Î( fA )min A∈A I( fA ) ≤ ∑ A∈A I Î( fA ) -I( fA ) fA ) -I( f * A ) .

Fig. 11 :

 11 Fig. 11: Loss rate estimates with lowest RMSE of individual methods.

Fig. 12 :

 12 Fig. 12: Loss rate estimates with lowest RMSE of individual methods.

Fig. 13 :

 13 Fig. 13: Data regressions under periodic and uniform losses.

Fig. 14 :

 14 Fig. 14: Data regressions under Poisson and Gaussian losses.

Fig. 15 :

 15 Fig. 15: Data regressions of loss estimators with lowest RMSE.

Fig. 16 :

 16 Fig.16: Loss regression at low and high loss rates for periodic and uniform losses.

Fig. 17

 17 Fig. 17: Loss regression at low and high loss rates for Poisson and Gaussian losses.

:

  Fig. 17: Loss regression at low and high loss rates for Poisson and Gaussian losses.

Table 2 :

 2 which approaches 1 with increasing number of measurements l. Loss rate estimates RMSE of five estimators and linear fuser.

	loss type EOT GPR LR RT SVM LR-F
	gaussian 5.19 6.39 5.34 6.82 7.17 6.42
	periodic 7.26 6.62 6.82 7.00 6.91 6.81
	poisson 6.99 6.61 7.42 6.93 6.72 6.74
	uniform 7.26 6.62 6.82 7.00 6.91 6.81

This work constitutes only initial steps in understanding the complexity of estimating network-level parameters using application-level measurements, and the performance of various ML solutions, including individual and fused estimates. Future work may involve studying the random losses due to external traffic in production networks, which may not follow known random processes. Since there is no universal way to choose among various ML methods from sample performance only, it would be of future interest to investigate into domain specific customizations, hyper-parameter tuning, fusers and other approaches to RTT and loss rate estimation [12,13].

* This work is funded by RAMSES project and Applied Mathematics program, Office of Advanced Computing Research, U.S. Department of Energy, and by Extreme Scale Systems Center, sponsored by U.S. Department of Defense, and performed at Oak Ridge National Laboratory managed by UT-Battelle, LLC for U.S. Department of Energy under Contract No. DE-AC05-00OR22725.


Ola Salman

Imad H Elhajj

Ayman Kayssi

Ali Chehab
email: chehab@aub.edu.lb

Denoising Adversarial Autoencoder for Obfuscated Traffic Detection and Recovery

Keywords: Machine Learning, Network Security, Traffic Classification, Obfuscation, Deep Learning, IoT, Autoencoder, Generative Adversarial Network

Traffic classification is key for managing both QoS and security in the Internet of Things (IoT). However, new traffic obfuscation techniques have been developed to thwart classification. Traffic mutation is one such obfuscation technique, that consists of modifying the flow's statistical characteristics to mislead the traffic classifier. In fact, this same technique can also be used to hide normal traffic characteristics for the sake of privacy. However, the concern is its use by attackers to bypass intrusion detection systems by modifying the attack traffic characteristics. In this paper, we propose an unsupervised Deep Learning (DL)-based model to detect mutated traffic. This model is based on generative DL architectures, namely Autoencoders (AE) and Generative Adversarial Network (GAN). This model consists of a denoising AE to de-anonymize the mutated traffic and a discriminator to detect it. The implementation results show that the traffic can be denoised when different mutation techniques are applied with a reconstruction error less than 10 -1 . In addition, the detection rate of fake traffic reaches 83.7%.

Introduction

In the IoT era, billions of things are being connected to the Internet. This results in an unprecedented growth in the amount of generated traffic. This traffic presents different QoS and security challenges. Traffic classification emerges as a key enabling tool to meeting some of these challenges. For example, IoT devices' vulnerabilities have been exploited to perform critical network attacks (e.g. Mirai) [START_REF] Perrone | The day after mirai: A survey on mqtt security solutions after the largest cyber-attack carried out through an army of iot devices[END_REF]. Tracking these devices and detecting abnormal traffic are key to prevent harmful network attacks.

In this context, Machine Learning (ML) based methods were proposed for traffic classification and intrusion detection. Many supervised and unsupervised methods have been employed accordingly. While supervised methods are devoted to classifying the (attack) traffic based on known class labels, the unsupervised ones allow the detection of unknown traffic. The traditional ML methods rely on well-structured and hand-designed features. These features are extracted using statistical traffic measures (e.g. maximum, minimum, standard deviation, etc. of the packet sizes, and packets interarrival times).

However, mutation techniques alter the traffic statistical features, making it very challenging to know the original traffic type. Moreover, the mutation techniques might change the packets characteristics while maintaining the flow statistical features unchanged. In this case, the detection of abnormal traffic can be evaded.

Recently, Deep Learning (DL) has acquired a lot of attention due to its representation learning capabilities. Using a new data representation in this paper, we propose an unsupervised DL model to detect abnormal traffic and de-anonymize the mutated one. Generative DL architectures, namely Autoencoders (AE) and Generative Adversarial Networks (GAN), have been applied mainly in the computer vision domain to detect abnormality and to denoise images. AE is a DL architecture for extracting data representation, and GAN has the capability to generate fake data samples to enhance the discrimination between real and fake data. In this paper, we combine AE and GAN to detect abnormal traffic and deanonymize the mutated one. The proposed architecture consists of an encoder, a decoder, and a discriminator. The encoder-decoder pair form a denoising AE responsible for learning the original data representation and to denoise the mutated one. In parallel, the discriminator is trained to differentiate between the mutated traffic (abnormal) and the denoised traffic (normal). The training of the proposed model relies on data collected from real IoT devices and IoT attacks. The testing results show the robustness of the proposed method to detect mutated traffic and to recover the original one. Note that the proposed model is not limited to IoT traffic and can be applied to any type of network traffic.

The rest of the paper is organized as follows: in section 2, we present an overview of related concepts. In section 3, we present our proposed model. Section 4 details the implementation and presents the evaluation results. In section 5, we discuss the results and present our future work. Finally, we conclude in section 6.

Background and Related Work

In this section, we present the related work including a review of the generative DL architectures, their application in detecting abnormality, traffic classification and intrusion detection methods, along with the obfuscation techniques that can affect their accuracy.

Unsupervised Deep Learning

In this subsection, we explain the DL architectures that we based our work on. These architectures can learn the input data representation and are therefore called generative DL models.

Generative Adversarial Network: Introduced by Goodfellow et al. in [START_REF] Goodfellow | Generative adversarial nets[END_REF], GAN consists of two parts: the generator (G) and the discriminator (D). From a game theoretic perspective, GAN can be interpreted as a zero-sum or minmax game between the generator and the discriminator. The generator tries to learn the input data representation to generate data samples very similar to the real ones. The discriminator tries to maximize the probability of distinguishing between fake and real input. The GAN objective function can be presented as follows:

V (G, D) = E x∼p data (x) (log(D(x)) + E z∼pz(z) (log(1 -D(G(z)))))
where x is the input data, p data (x) is the data distribution, D(x) is the discriminator output, p z (z) is the fake data distribution, z is a sample from p z (z), and G(z) is the generator output. The generator aims at minimizing the probability of fake data detection by the discriminator, which means that the G objective is to find min

G (E z∼p(z) (log(1 -D(G(z))))).
The discriminator aims at maximizing the probability of detecting real data as real and fake data as fake, which means that the D objective is to find max

D (E x∼p data (x) (log(D(x))) + E z∼p(z) (log(1 - D(G(z))))). Thus, the GAN objective is to find min G max D (V (G, D))
. Primarily, GAN is applied for synthetic data generation. GAN has been applied also for image anomaly detection [START_REF] Zenati | Efficient ganbased anomaly detection[END_REF][START_REF] Li | Anomaly detection with generative adversarial networks for multivariate time series[END_REF][START_REF] Deecke | Image anomaly detection with generative adversarial networks[END_REF]. In fact, the adversarial learning permits the discriminator to detect abnormal input data. Furthermore, the generative learning permits the generator to learn the real data representation, which makes GAN suitable for image denoising [START_REF] Tripathi | Correction by projection: Denoising images with generative adversarial networks[END_REF][START_REF] Bengio | Improving generative adversarial networks with denoising feature matching[END_REF].

Autoencoders: Being a generative model, AE is a type of DL networks that is specialized in extracting the input data representation. The AE consists of two parts: an encoder and a decoder. The encoder extracts a compressed data code by estimating a function f , in such a way z = f (x), where x is the input data, and z is the extracted representation or latent variable. The decoder aims at reconstructing the input data by relying on the extracted representation. In other terms, the decoder tries to infer the inverse function g, in such a way that y = g(f (x)) = g(z). The AE objective function can be represented at the minimization of the difference between g(f (x)) and x. In other terms, the AE aims at minimizing the reconstruction error. A type of AEs is the probabilistic autoencoder, which aims to infer the distribution of x, p θ (x) by means of another distribution q φ (z/x) (probability of the latent variable z knowing the input x). A well-known type of probabilistic AEs is the Variational Autoencoder (VAE), which imposes a prior restriction on p(z) to be a normal distribution. In this case, the problem reduces to maximizing the Evidence Lower Bound (ELBO) or maximizing the KullbackLeibler (KL) divergence between q φ (z) and p θ(z/x) , represented by KL(q φ (z/x)||p(z)). Having the ability to extract the real data distribution, VAEs have been applied for anomaly detection. Indeed, when the reconstruction error is large, anomaly is detected in the input data. Borrowing the adversarial concept from GAN, Adversarial AEs (AAE) were introduced by Makhzani et al. in [START_REF] Makhzani | Adversarial autoencoders[END_REF]. Similar to the GAN, AAE includes a discriminator that tries to differentiate between the data sampled from the latent variable prior p(z) and the real data. In this case, the discriminator aims to minimize

L dis = -1/N N -1 i=0 log(d x (z i )) + 2N j=N log(1 -d x (z j ))
, and the generator tries to minimize

L prior = 1/N N -1 i=0 (log(1 -d x (z i ))
, where N is the number of samples, and d x (z i ) is the discriminator output of the latent space variable. In this case, if the total loss function is optimized, q φ (z/x) will be very similar to p(z), or in other terms, KL(q φ (z/x)||p(z)) will be minimized, and thus the log likelihood of the original data distribution will be maximized. AAEs were applied also for anomaly detection in images [START_REF] Ger | Autoencoders and generative adversarial networks for anomaly detection for sequences[END_REF][START_REF] Wang | Self-adversarial variational autoencoder with gaussian anomaly prior distribution for anomaly detection[END_REF][START_REF] Vu | Anomaly detection with adversarial dual autoencoders[END_REF][START_REF] Beggel | Robust anomaly detection in images using adversarial autoencoders[END_REF]. Furthermore, a recent work has considered to add the denoising function to the AAEs for image denoising [START_REF] Creswell | Denoising adversarial autoencoders[END_REF]. In this case, the corrupted data x is considered as input and two methods were proposed for model representation. The first consists of matching qφ (z/x) to p(z), and the second consists of matching q φ (z/x) to p(z). However, in our work, we use a sparse AE that aims to minimize the Mean Square Error (MSE) between the reconstructed data and original data. In addition, applying the adversarial concept, we choose to train a discriminator to detect abnormal traffic when the reconstruction error is high. Thus, unlike previous work, the generator part of GAN is omitted [START_REF] Chalapathy | Deep learning for anomaly detection: A survey[END_REF].

ML based Traffic Classification and Intrusion Detection

Traffic classification is an essential network function, which is necessary for traffic engineering, QoS management, and security management. Different methods have been proposed for traffic classification using different sets of features [START_REF] Dainotti | Issues and future directions in traffic classification[END_REF][START_REF] Nguyen | A survey of techniques for internet traffic classification using machine learning[END_REF][START_REF] Finsterbusch | A survey of payload-based traffic classification approaches[END_REF][START_REF] Cao | A survey on encrypted traffic classification[END_REF][START_REF] Pacheco | Towards the deployment of machine learning solutions in network traffic classification: a systematic survey[END_REF][START_REF] Callado | A survey on internet traffic identification[END_REF]. More recently, DL has been applied for traffic classification using new features and new data representations [START_REF] Rezaei | Deep learning for encrypted traffic classification: An overview[END_REF][START_REF] Wang | End-to-end encrypted traffic classification with one-dimensional convolution neural networks[END_REF][START_REF] Chen | Seq2img: A sequence-to-image based approach towards ip traffic classification using convolutional neural networks[END_REF][START_REF] Smit | Looking deeper: Using deep learning to identify internet communications traffic[END_REF][START_REF] Fadlullah | State-of-the-art deep learning: Evolving machine intelligence toward tomorrows intelligent network traffic control systems[END_REF]. On the other hand, intrusion detection is an essential network security component. Several approaches have been adopted to detect and prevent network attacks. Traditional Intrusion Detection Systems (IDSes) are rule-based, where the attack signature is known by identifying some patterns in the packets fields. However, this method fails to detect unknown attacks (e.g. zero-day attacks) and requires the inspection of the packet header. In addition, traffic anonymization can be used to thwart detection by traditional IDS. ML methods have been proposed to detect network attacks and traffic abnormality by means of statistical features [START_REF] Buczak | A survey of data mining and machine learning methods for cyber security intrusion detection[END_REF]. However, the traditional methods require the extraction of the features by computing statistical measures that might be data dependent. In addition, the detection of abnormal traffic is not a straightforward task, given that the abnormal traffic might present similar statistical behavior to the normal one [START_REF] Baddar | Behavioral-anomaly detection in forensics analysis[END_REF]. Recently, DL has been applied in the network domain for intrusion detection [START_REF] Alom | Network intrusion detection for cyber security using unsupervised deep learning approaches[END_REF][START_REF] Vinayakumar | Deep learning approach for intelligent intrusion detection system[END_REF]. More specifically, IoT, which presents aggravated security challenges, has acquired a special attention from the intrusion detection perspective [START_REF] Chaabouni | Network intrusion detection for iot security based on learning techniques[END_REF][START_REF] Da Costa | Internet of things: A survey on machine learning-based intrusion detection approaches[END_REF]. Convolutional Neural Network (CNN) and Recurrent neural Network (RNN) have been applied for payload-based attack detection in [START_REF] Liu | Cnn and rnn based payload classification methods for attack detection[END_REF]. While supervised learning was mainly considered for identifying specific attacks [START_REF] Xiao | An intrusion detection model based on feature reduction and convolutional neural networks[END_REF], the detection of unknown and zero-day attacks call for unsupervised-learning-based methods. In fact, the proposed unsupervised [START_REF] Munir | Deepant: A deep learning approach for unsupervised anomaly detection in time series[END_REF] or semi-supervised [START_REF] Umer | A two-stage flow-based intrusion detection model for next-generation networks[END_REF] DL methods for intrusion detection use statistical features. To the best of our knowledge, no previous work considered the recovering (denoising) of the mutated traffic. In addition, this is the first work to consider the question of detection of mutated traffic, in addition to the (unknown) attack traffic.

Traffic Obfuscation Techniques

In the aim of protecting user privacy, a new research direction is considering traffic obfuscation to thwart classification. In this context, many obfuscation techniques have been proposed [START_REF] Cheng | Evasion techniques: Sneaking through your intrusion detection/prevention systems[END_REF][START_REF] Khraisat | Survey of intrusion detection systems: techniques, datasets and challenges[END_REF]. These methods can be classified in seven categories: steganography, tunneling, anonymization, mutation, morphing, and physical layer obfuscation. While steganography and physical layer obfuscation techniques require specific protocols to recover the original data, the remaining techniques can be applied without imposing changes to the current network protocols. Anonymization and tunneling hide some packet-related information by encryption and establishment of virtual connections (port numbers and internet addresses). However, statistical ML methods still have power to classify the anonymized or tunneled traffic. Mutation and morphing are two techniques that consider modifying the statistical traffic characteristics considering the modification of the packet size and the packet Inter-Arrival Times (IAT) [START_REF] Chaddad | App traffic mutation: Toward defending against mobile statistical traffic analysis[END_REF]. This has great impact on the accuracy of ML-based classification and intrusion detection. Even though the obfuscation techniques were intended to protect the user privacy, attackers might use them to perform their attacks without being detected. In this context, padding and traffic shaping are proposed to modify the packet size and IAT respectively [START_REF] Qu | An empirical study of morphing on network traffic classification[END_REF][START_REF] Qu | An empirical study of morphing on behavior-based network traffic classification[END_REF]. Recently, the GAN DL architecture has been employed to adapt malware traffic to the normal one [START_REF] Verma | Network traffic obfuscation: An adversarial machine learning approach[END_REF][START_REF] Zhang | Deep adversarial learning in intrusion detection: A data augmentation enhanced framework[END_REF][START_REF] Lin | IDSGAN: generative adversarial networks for attack generation against intrusion detection[END_REF][START_REF] Rigaki | Bringing a gan to a knife-fight: Adapting malware communication to avoid detection[END_REF]. On the contrary, in this paper, we will consider a discriminative denoising DL model to detect abnormality by relying on the trained discriminator and to de-anonymize mutated traffic by means of a denoising AE. Our contributions in this paper can thus be summarized by the following points:

-The discriminative part of GAN with a denoising AE are combined to allow mutated traffic detection and recovery. -A new data representation is used to permit the de-anonymization of mutated traffic by applying DL-based techniques.

Proposed Abnormal Traffic Detection

In this section, we detail our proposed DL model, the attack model, and the traffic representation method.

Attack Model

Our attack model consists of an attacker trying to modify the packet size (padding) and IAT (shaping), in such a way to hide any information that serves for attack detection or traffic classification. The mutations are therefore of two types, padding and shaping [START_REF] Lin | IDSGAN: generative adversarial networks for attack generation against intrusion detection[END_REF][START_REF] Liu | Cnn and rnn based payload classification methods for attack detection[END_REF], as shown in Fig. 1. In the following, we summarized the packet padding techniques listed in [START_REF] Qu | An empirical study of morphing on network traffic classification[END_REF][START_REF] Qu | An empirical study of morphing on behavior-based network traffic classification[END_REF][START_REF] Chaddad | Mobile traffic anonymization through probabilistic distribution[END_REF], with s being the packet original size and m(s) being the packet size after mutation. The last method, described in [START_REF] Dyer | Peek-a-boo, i still see you: Why efficient traffic analysis countermeasures fail[END_REF], combines shaping and padding.

10. Fixed size and fixed IAT: this method consists of padding the size of all the flow packets to the MTU and sending all the packets at a fixed time interval.

This model considers two types of adversaries, including:

-Malicious adversary: this type of attackers aims at mutating the attack traffic to evade intrusion detection. -Benign adversary: this type of attackers aims at hiding the traffic characteristics to protect the user privacy.

Deep Learning Model

Our proposed model, illustrated in Fig. 2, consists of two parts: a denoising AE and a discriminator. The denoising AE consists of an encoder and a decoder. The encoder aims at estimating the function f , that maps the input space of X to the latent space of the latent variable Z, where Z is a compressed representation of X. The denoising AE is fed with a dataset containing the mutated version x of the initial data X. The AE aims at minimizing the reconstruction error L(X, X). In our case, the loss function is the MSE function, L(X, X) = (Xg(f ( X))) 2 , where X is the original data and g(f ( X)) is the reconstructed data by considering the mutated data as input. The discriminator aims at differentiating between normal and abnormal traffic. In fact, the discriminator will be trained on classifying the mutated data (i.e. the mutated input data) as abnormal and the reconstructed data (i.e. the decoder output) as normal. The output function of the discriminator is a sigmoid function p(y = y (j) /x) = 1/(1 + e -y (j) ), where y (j) is the output class that can take the two values: 0 for j = 0 and 1 for j = 1, and the loss function is a cross entropy function L D = 1/m 1 j=0 y (j) log(y (j) ) + (1 -y (j) )log(1 -y (j) )).

Proposed Model Workflow

After training the model, the proposed scheme workflow, illustrated in Fig. 3, consists of: 1) passing the traffic to the discriminator; 2) if the traffic is normal, it is passed to the classifier of normal traffic; 3) If not, it is passed to the denoiser. 4) After denoising, it is passed again to the discriminator. 5) If a normal traffic is detected, it is passed to the classifier; 6) if not, the traffic is detected as abnormal, so it might be obfuscated or attack traffic. The same workflow is repeated for the abnormal traffic to know the attack type or detect an unknown attack traffic. In fact, two cases are considered. First, when training the model on normal traffic, the aim is to detect attack traffic and the mutated one as abnormal. In this case, the denoising aims at recovering the normal traffic for classification. However, when training the model with attack traffic, at the testing phase the aim is to detect unknown attacks and to denoise the mutated attack traffic to know the exact attack type. Note that the classifier will be trained in a supervised mode to classify the traffic based on the IoT device type in the normal traffic case and based on the attack type in the attack traffic case.

Data Representation

In our case, the data consists of collected network traffic. This traffic is filtered by flows, where a flow is defined as being the set of packets having the same: source IP address, source port number, destination IP address, destination port number, and protocol (TCP or UDP). For each packet, we extract three features: size (s), IAT (t), and direction (d). For each flow, we consider the first 4 * 4 packets in either direction. In fact, as shown in Fig. 3, the extracted features can be visualized in 4 * 4 RGB images, where the i th pixel RGB coordinates are Fig. 3. Proposed scheme workflow represented by the i th packet features [s, t, d]. Thus, every feature represents a color channel. In our case, R is given for the size , G is given for the interarrival time, and B is given for the direction. This data representation is explained in detail in [START_REF] Salman | A multi-level internet traffic classifier using deep learning[END_REF].

Implementation

In this section, we detail the data collection and preprocessing, the model implementation, and the evaluation results.

Data Collection and Preprocessing

For collecting normal IoT traffic, we used a set of IoT devices, including Wi-Fienabled devices and hub-connected devices. The Wi-Fi devices are: D-Link HD 180-Degree Wi-Fi Camera DCS-8200LH, D-Link Wi-Fi Smart Plug DSP-W215, D-Link Wi-Fi Siren DCH-S220, and D-Link Wi-Fi Water Sensor DCH-S160, while the hub-connected devices are components of the Samsung SmartThings Home Monitoring Kit, including a motion sensor, a multi-purpose sensor, and a smart plug. These devices were installed in a home environment and were left to function normally. The Wi-Fi enabled devices are routed through a laptop to the Internet. A bridge is created at this laptop to forward the incoming traffic to the Ethernet interface connected to the home gateway. For the hub-connected devices, the hub is connected to the laptop by Ethernet and this laptop is configured to forward the incoming traffic to the wireless interface connected to the home gateway. In both cases, the traffic is collected at the laptop. One day of traffic for each device was considered for training and one day of traffic was used for testing. For IoT attack traffic training data, we consider the dataset collected in [START_REF] Koroniotis | Towards the development of realistic botnet dataset in the internet of things for network forensic analytics: Bot-iot dataset[END_REF], this dataset consists of multiple PCAP files for each type of attack. We choose one file for each type of attack for training and one file for testing. To preprocess the data and extract the flows, the dpkt Python library was used [1]. In total, for training, we have 10320 flows of normal traffic and 3308 flows of attack traffic. For testing, we have 258 flows of normal traffic, 350 flows of attack traffic, and 2000 flows of (unknown) attack traffic. The normal traffic is categorized into five classes based on the device model while the attack traffic is categorized into three classes: data theft, denial of service and scanning.

Model Implementation and Training

The proposed model was implemented in Python using the tensorflow library [START_REF]tensorflow[END_REF]. The encoder, decoder and the discriminator consist of a 2-layer fully connected neural network with 1000 neurons each. The output layer of the decoder and the discriminator is a sigmoid layer, with the difference that the decoder output is of the same dimension of the input; however, the discriminator output is of dimension one. The model is trained with 100 epochs and a batch size of 100, the learning rate is 10 -3 , and the momentum decay is 0.9 (beta1). The Rectified Linear Unit (ReLU) was used as the activation function for all hidden layers, and the weights are optimized using the Adam optimizer. For generating mutated data, we implement the mutation techniques listed in section III-B ((1) →( 10)).

The training data is randomly mutated, where each data sample is uniformly mutated to one of the 10 mutation techniques. For testing the effectiveness of the denoising process on the classification accuracy, we implement a CNN classifier with three convolutional layers, two max pooling layers, and two fully connected layers with one dropout layer with 50% dropout probability. Similarly, ReLU is applied for activation in the hidden layers, and the Adam algorithm is used for optimization with a learning rate of 10 -3 . In addition, we apply cross validation for the classifier training with 4 folds. The performance metric used to evaluate the classifier is the accuracy, which is the ratio of the correctly classified samples to the total number of samples.

Evaluation and Results

For each of the experiments (1) →( 10), the corresponding mutation technique is applied to the testing data. First, the traffic is passed to the discriminator. Then, it is passed to the denoising AE. Furthermore, the mutated and denoised traffic are passed to a classifier, that classifies the flow based on the device model for normal traffic and based on the attack type for attack traffic. Samples of the mutated traffic and resulted images after denoising are included in Table 5. It is visually clear from the included images that the denoiser succeeds in learning the original data representation disregarding the mutation level. The difference between the denoised images and the original ones shows that the model does not overfit the data, however it learns the representation and the noise pattern.

For the normal traffic, the denoiser succeeds in recovering flows very similar to the original ones for most of the mutation techniques, except for the mutations ( 8) and [START_REF] Cao | A survey on encrypted traffic classification[END_REF]. Similarly, for the attack traffic, the denoiser recovers the main flow representation, except for the mutation (2). In Table 1, the MSE between the original data and the mutated one (mutation degradation), the reconstruction loss, and the abnormality detection rate are reported for the normal traffic and the attack traffic. The results present a high detection rate for the different mutation techniques, except for the mutation techniques (4), [START_REF] Buczak | A survey of data mining and machine learning methods for cyber security intrusion detection[END_REF], and ( 9). This can be explained by the fact that ( 6) and ( 9) are normal distribution based mutations. In these cases, the main traffic characteristics will remain unchanged and this will harden the differentiation between the original and mutated traffic. ( 4) is a mice-elephant mutation of the packet sizes. However, in our case (i.e. IoT traffic), most of the packets are of small size and therefore the mutation (4) will have limited affect on the traffic characteristics. Moreover, it can be noticed that the MSE between the recovered data and the original one (reconstruction loss) Table 2 presents the accuracy of the classification based on the traffic label: device model for normal traffic and attack type for attack traffic. The mutation techniques (1), ( 4), [START_REF] Qu | An empirical study of morphing on network traffic classification[END_REF], and (10) affect the accuracy and misleads the classification noticeably in the normal traffic case; however, after denoising, the accuracy increases. However, for the techniques (5), ( 6), [START_REF] Callado | A survey on internet traffic identification[END_REF], and (9), the denoiser fails to reconstruct the original traffic. This is due to the fact that the randomness will create a denoised traffic of random type. Moreover, for the mutation technique (2), the mutation is linear and this does not affect the CNN classifier accuracy, being immune to the linear degradation of the image. However, for statistical based machine learning methods, the mutation techniques (2), [START_REF] Beggel | Robust anomaly detection in images using adversarial autoencoders[END_REF], and (6) affect noticeably the classification accuracy. To see the effect of the mutation on the statistical based classification, we include the results of the mutated and denoised traffic statistical based classification in Table2. The statistical classification uses a subset of the Moore features (see Table3) and Random Forest (RF) as classifier. It can be noticed that overall our representation method with CNN classifier outperforms the RF method before and after denoising in the normal traffic case. However, in the attack traffic case, our method gives better results after denoising.

Discussion and Future Work

The results show that unsupervised DL architectures are very powerful in learning the data representation. AE is a well-known DL generative model, that is highly effective in extracting compressed representation from image-type data. Consequently, after representing network traffic as images, we applied AE to extract the representation patterns from IoT traffic. Moreover, the capability of AE to denoise images is used to overcome the mutation technique challenges.

The results show the effectiveness of the proposed method to recover the original traffic representation for different levels of mutation, some of which are rather severe (e.g. fixed packet size and IAT). In fact, the considered mutation techniques cause any statistical classifier to fail. However, one limitation of this work is that we assume that we know ahead of time the mutation technique used by the attacker, so we can decide to choose to apply the classifier before or after denoising. While we considered mutation techniques in this paper, we will study as future work additional obfuscation techniques. In addition to the abnormal Feature Description total f wd pkt Total packets in the forward direction total f wd bytes Total bytes in the forward direction total bck pkt Total packets in the backward direction total bck bytes Total bytes in the backward direction min pckt size f wd The min packet size in the forward direction mean pckt size f wd The mean packet size in the forward direction max pckt size f wd The max packet size in the forward direction std pckt size f wd

The standard deviation packet size in the forward direction min pckt size bck

The min packet size in the backward direction mean pckt size bck The mean size of packets in the backward direction max pckt size bck The max packet size in the backward direction std pckt size bck

The standard deviation packet size in the backward direction total size

The total flow size min pckt size

The min packet size in either direction mean pckt size

The mean size of packets in either direction max pckt size

The max packet size in either direction std pckt size

The standard deviation packet size in either direction min iat f wd

The minimum interarrival time in the forward direction mean iat f wd

The mean interarrival time in the forward direction max iat f wd

The maximum interarrival time in the forward direction std iat f wd

The standard deviation interarrival time in the forward direction min iat bck

The minimum interarrival time in the backward direction mean iat bck

The mean interarrival time in the backward direction max iat bck

The maximum interarrival time in the backward direction std iat bck

The standard deviation interarrival time in the backward direction total time

The duration of the flow min iat

The minimum interarrival time in either direction mean iat

The mean interarrival time in either direction max iat

The maximum interarrival time in either direction std iat

The standard deviation interarrival time in either direction avg pckt f wd

The average number of packets in the forward direction avg bytes f wd

The average number of bytes in the forward direction avg pckt bck

The average number of packets in the backward direction avg bytes bck

The average number of bytes in the backward direction avg iat f wd

The proportion of flow time in the forward direction to the total flow time avg iat bck

The proportion of flow time in the backward direction to the total flow time Table 3. Feature Set for statistical classification traffic detection using the discriminative part of the GAN architecture, in the morphing case for example, we need to consider the generative part also. The generator will be trained to generate morphed traffic similar to the original one. In this case, the model convergence will not be a straightforward task.

Conclusion

Applying ML for network management and network security has gained a lot of interest in the last decade. However, new traffic obfuscation techniques have been developed to thwart classification and avoid detection in the sake of user privacy. These techniques have been employed by attackers to mount their attacks without being detected. While the obfuscation techniques might be detected by behavioral or statistical ML techniques, the recovery of the initial traffic is unfeasible. In this paper, inspired from a promising DL application, which is image denoising, we transform the traffic to images then combine two well-established DL architectures (AE and GAN) to reconstruct the original traffic and detect abnormal one. The test results show the effectiveness and robustness of the proposed model to detect abnormal traffic in all its variants.

Fig. 1 . 7 . 8 .

 178 Fig. 1. Attack model

Fig. 2 .

 2 Fig. 2. Proposed model

Fig. 4 .

 4 Fig. 4. Data Representation

Fig. 5 .

 5 Fig. 5. Visualized images representing network traffic

Table 1 .

 1 Testing evaluation results

	Mutation	Autoencoder Loss Mutation Loss Discrimination Rate
	Technique	Normal Attack Normal Attack Normal	Attack
	(1)	0.042	0.0357	0.2713 0.1698 100%	100%
	(2)	0.03218	0.1569	0.0171 0.082	100%	100%
	(3)	0.0149	0.1569	0.0009 0.0065 88.99%	83.13%
	(4)	0.0543	0.0463	0.2704 0.1696	50%	49.94%
	(5)	0.067	0.0519	0.0171 0.082	100%	100%
	(6)	0.0621	0.0703	0.0129 0.0411	49%	49.88%
	(7)	0.0378	0.1569	0.2974 0.3245 100%	100%
	(8)	0.0595	0.0257	0.09878 0.1016 99.02%	100%
	(9)	0.1728	0.0303	0.2825 0.0012	50%	41.59%
	(10)	0.0369	0.0519	0.5687 0.4943 100%	100%

Table 2 .

 2 Classification Results

			Before Denoising			After Denoising
	Mutation Technique	Normal	Attack	Normal	Attack
		CNN	RF	CNN	RF	CNN	RF	CNN	RF
	(1)	30.14% 54.47% 51.55% 80.84% 65.19% 38.78% 53.1% 37.05%
	(2)	93.13% 25% 36.27% 36.24% 79.16% 68.68% 33.41% 33.41%
	(3)	99.01% 87.19% 76.25% 90.33% 87.99% 78.79% 33.41% 33.41%
	(4)	29.9% 55.2% 51.55% 80.1% 67.4% 40.68% 43.79% 42.57%
	(5)	93.13% 25.06% 36.99% 36.21% 20.09% 20.09% 63.12% 38.21%
	(6)	95.09% 25.06% 56.55% 62.94% 56.12% 33.88% 36.99% 42.69%
	(7)	75.49% 85.6% 33.17% 19.63% 77.69% 55.39% 33.41% 33.41%
	(8)	78.18% 81.92% 33.15% 15.96% 62.99% 41.85% 62.05% 38.48%
	(9)	80.14% 82.35% 34.24% 39.64% 20.09% 20.09% 62.76% 30.31%
	(10)	19.6% 28.18% 33.17% 25.65% 78.43% 49.14% 60.6% 43.22%

Acknowledgments

Research funded by the AUB University Research Board, the Lebanese National Council for Scientific Research, and TELUS Corp., Canada.


Jaekyung Yang
email: jkyang@jbnu.ac.kr

Jong-Yeong Lee

Myoungjin Choi

A new approach to determine the optimal number of clusters based on the Gap statistic

Keywords: Clustering, Number of Clusters, Data Mining. 1

Data clustering is one of the most important unsupervised classification method. It aims at organizing objects into groups (or clusters), in such a way that members in the same cluster are similar in some way and members belonging to different cluster are distinctive. Among other general clustering method, kmeans is arguably the most popular one. However, it still has some inherent weaknesses. One of the biggest challenges when using k-means is to determine the optimal number of clusters, k. Although many approaches have been suggested in the literature, this is still considered as an unsolved problem. In this study, we propose a new technique to improve the gap statistic approach for selecting k. It has been tested on different datasets, on which it yields superior results compared to the original gap statistic. We expect our new method to also work well on other clustering algorithms where the number k is required. This is because our new approach, like the gap statistic, can work with any clustering method.

Introduction

There are still many open challenges in the clustering task. Those challenges are getting even worse in the current big data era, where data is collected from many sources at high speed. This paper focuses on answering the question: how to decide on the number of clusters k? Being one of the oldest question in the clustering literature, the question has been tackled by hundreds of researchers with many solutions that have been proposed. Among these solutions, the gap statistic is one of the most modern approaches. It is backed by the rigorous theoretical foundation and has been shown to outperform many other heuristic-based approaches such as elbow or silhouette. However, there are still several drawbacks to the original design of the gap statistic, which limits its applicability in real applications. This paper introduces a new technique to mitigate those limitations. The technique can improve the effectiveness of the gap statistic in multiple dimensions. The gap statistic that uses the newly proposed technique is called the "new gap" for short. The following few subsections describe literature reviews

The Elbow Approach

The oldest method called 'elbow' has been proposed to determine the number of clusters for k-mean clustering algorithm [START_REF] Trupti | Review on determining number of Cluster in k-means Clustering[END_REF]. This is a visual method. The idea of the elbow method is to run clustering method on the dataset for a range of values of k (for example from 1 to 10), and for each value of k calculate clusters and internal index (it could be the sum of squared error (SSE), the percentage of variance, etc.). Then plot a line chart of the internal index for each value of k. At some value of k the value of internal index drops dramatically, and after that, it reaches a plateau when k is increased further. This is the best k value we can expect. Fig 1 illustrates how the elbow method work. In Fig 1, the line chart goes down rapidly with k increasing from 1 to 2, and from 2 to 3, and reaches an elbow at k = 3. After that, it decreases very slowly. Looking at the chart, it looks like maybe the right number of cluster is three because that is the elbow of this curve.

Fig. 1. Identification of Elbow point

However, the elbow method does not always work well. Sometimes, there are more than one elbow, or no elbow at all.

Average Silhouette Approach

Average silhouette method computes the average silhouette of observations for different values of k [2] [START_REF] Rousseeuw | Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[END_REF]. The optimal number of clusters k is the one that maximizes the average silhouette over a range of possible values for k [START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF]. Given a clustering result with k clusters (k > 1), we can estimate how well an observation i is clustered by calculating its silhouette statistic 𝑠𝑠 𝑘𝑘 (𝑖𝑖). Let a(i) be the average distance from observation i k (the number of clusters) Internal Index "Elbow" to other points in its cluster, and b(i) be the average distance from observation i to points in its the nearest cluster, then the silhouette statistic 𝑠𝑠 𝑘𝑘 (𝑖𝑖) is calculated by: 𝑠𝑠 𝑘𝑘 (𝑖𝑖) = 𝑏𝑏(𝑖𝑖) -𝑎𝑎(𝑖𝑖) max{𝑎𝑎(𝑖𝑖), 𝑏𝑏(𝑖𝑖)} A point is well clustered if 𝑠𝑠 𝑘𝑘 (𝑖𝑖) is large. The average silhouette score 𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎(𝑘𝑘) gives an estimation of the overall clustering quality when clustering the dataset into k clusters:

𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎(𝑘𝑘) = 1 𝑛𝑛 � 𝑠𝑠 𝑘𝑘 (𝑖𝑖) 𝑛𝑛 𝑖𝑖=1
, where n is the number of data points. Therefore, we select k so that it maximizes the average silhouette score. However, this average silhouette is only a heuristic metric, which can be shown to perform poorly in many cases. Note that avgS(k) is not defined at k = 1.

Hartigan Statistic

Hartigan (1975) proposed the statistic [START_REF] John | Clustering algorithms[END_REF]:

𝐻𝐻(𝑘𝑘) = 𝑊𝑊 𝑘𝑘 𝑊𝑊 𝑘𝑘+1 - 1
𝑛𝑛 -𝑘𝑘 -1 , where 𝑊𝑊 𝑘𝑘 is the average within-cluster sum of squares around the cluster means. The formula to calculate 𝑊𝑊 𝑘𝑘 is given in the next section about the gap statistic.

The idea is to start with k = 1 and keep adding a cluster until H(k) is sufficiently large. Hartigan suggested the "sufficiently large" cut-off is 10. Hence the estimated number of clusters is the smallest 𝑘𝑘 ≥ 1 such that 𝐻𝐻(𝑘𝑘) ≤ 10.

Gap Statistic

Gap statistic was introduced in 2001 by Tibshirani et.al. [START_REF] Tibshirani | Estimating the number of data clusters via the Gap statistic[END_REF] and is still a state-ofthe-art method for estimating k. It has been shown to outperform the elbow, average silhouette, and Hartigan methods in both synthesized and real datasets [START_REF] Tibshirani | Estimating the number of data clusters via the Gap statistic[END_REF] [START_REF] Chiang | Intelligent Choice of the Number of Clusters in K-Means Clustering: An Experimental Study with Different Cluster Spreads[END_REF]. The method works by assuming a null reference distribution. It then compares the change in within-cluster dispersion with the expected change if the null distribution is true. If when k = K and the within-cluster dispersion starts decreasing slower than the expected rate of the reference distribution, the gap statistic returns k as the expected number of clusters. The formal definition of the gap statistic is given as follows:

Let 𝑑𝑑 𝑖𝑖𝑖𝑖 = �𝑥𝑥 𝑖𝑖 -𝑥𝑥 𝑖𝑖 � 2 denotes the Euclidean distance between observation i and j, 𝐷𝐷 𝑟𝑟 is the sum of the pairwise distance for all points in a given cluster 𝐶𝐶 𝑟𝑟 containing 𝑛𝑛 𝑟𝑟 points.

𝐷𝐷 𝑟𝑟 = � � 𝑑𝑑 𝑖𝑖𝑖𝑖 𝑖𝑖∈𝐶𝐶 𝑟𝑟 𝑖𝑖∈𝐶𝐶 𝑟𝑟
Then measure of compactness of clusters 𝑊𝑊 𝑘𝑘 is the average within -cluster sum of squares around the cluster means:

𝑊𝑊 𝑘𝑘 = � 1 2 𝑛𝑛 𝑟𝑟 𝐷𝐷 𝑟𝑟 𝑘𝑘 𝑟𝑟=1
The purpose of clustering is with a given K finding the optimal 𝑊𝑊 𝑘𝑘 , when k increases, 𝑊𝑊 𝑘𝑘 decreases. But the speed reduction of 𝑊𝑊 𝑘𝑘 also decreases. The idea of elbow method is to choose the k corresponding to the "elbow" (finding k that point has the most significant increase in goodness-of-fit). The problems when using elbow method is no reference clustering to compare, and the differences 𝑊𝑊 𝑘𝑘 -𝑊𝑊 𝑘𝑘-1 's are not normalized for comparison.

The main idea of the gap statistic is to standardize the graph of log(𝑊𝑊 𝑘𝑘 ) by comparing it with its expectation under an appropriate null reference distribution of the data. Estimate of the optimal number of clusters is then the value of k for which log(𝑊𝑊 𝑘𝑘 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 ) falls the farthest below this reference curve log�𝑊𝑊 𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 �:

𝐺𝐺𝑎𝑎𝐺𝐺 𝑛𝑛 (𝑘𝑘) = 𝐸𝐸 𝑛𝑛 * {log�𝑊𝑊 𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 � -log (𝑊𝑊 𝑘𝑘 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 )} With 𝐸𝐸 𝑛𝑛 * is
, b = 1, 2 …, B, k = 1, 2 …, K. Compute the (estimated) gap statistic: 𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘) = 1 𝐵𝐵 � 𝑙𝑙𝑙𝑙𝑎𝑎𝑊𝑊 𝑘𝑘𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 𝐵𝐵 𝑘𝑘=1 -𝑙𝑙og (𝑊𝑊 𝑘𝑘 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 )
Those 𝑙𝑙𝑙𝑙𝑎𝑎𝑊𝑊 𝑘𝑘𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 from the B Monte Carlo replicates exhibit a standard error 𝑠𝑠𝑑𝑑(𝑘𝑘) which, accounting for the simulation error, is turned into the quantity

𝑠𝑠 𝑘𝑘 = � 1 + 1 𝐵𝐵 . 𝑠𝑠𝑑𝑑(𝑘𝑘)
Finally, the optimal number of cluster K is the smallest k such that 𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘) ≥ 𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘 + 1) -𝑠𝑠 𝑘𝑘+1 The above rule to select k is presented in the original gap statistic paper and called the "Tibs2001SEmax" rule in the R clustering implementation of the gap statistic. Since 2001, several other alternatives to this rule have been proposed, such as the "firstSEmax" rule [START_REF] Maechler | firstSEMax rule in R document[END_REF] or the "globalSEmax" rule [START_REF] Dudoit | A prediction-based resampling method for estimating the number of clusters in a dataset[END_REF]. In this study, the Tibs2001SEmax rule in all experiments was used as the baseline approach. In this paper, the term "gap statistic" refers to the function Gap(k) with the Tibs2001SEmax is used as the k-selecting rule.

Fig. 2 provides an example of how the gap statistic works. Fig. 2a plots the example dataset with two well-separated clusters. Fig. 2b shows the line representing the within sum of squares 𝑊𝑊 𝑘𝑘 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 , which is a downward trend in number of cluster k. Fig. 2c shows the log of the expected rate log(𝑊𝑊 𝑘𝑘𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 ) using an assumed null distribution (uniform distribution in this case). Fig. 2d shows the gap statistic, which is calculated by subtracting the log expected rate log(𝑊𝑊 𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 ) for the log(𝑊𝑊 𝑘𝑘𝑘𝑘 𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 ). The optimal number of k is the smallest k such that there is a significant chance that Gap(k) is higher than Gap(k+1), which is k = 2 in this case. Tibshirani used one standard deviation 𝑠𝑠 𝑘𝑘+1 to determine when the chance is significant.

Methodology

Although being backed by a rigorous theoretical foundation (unlike other heuristicbased methods like elbow or silhouette), the Gap statistic still has several drawbacks that limit its applicability to practical applications. In this section, we conduct several experiments with synthesized datasets to demonstrate those limitations. Based on the insights learned from those experiments, we then introduced a new technique to improve the gap statistic.

The Gap Statistic Limitations

By design, the gap statistics can only work well when all the clusters in the dataset are well-separated from each other. However, this is rarely the case in practice, where clusters usually overlap up to a certain degree. This "non-overlapping" assumption is one of the main reason that limits the gap statistics effectiveness in real applications. Fig. 3 shows how the gap statistics fail to identify the correct K in simple synthesized datasets, that the clusters only barely overlap each other. However, clusters should not overlap with each other too much. Otherwise, the notion of "cluster" will become very fuzzy. This is because the data density in the overlapping area is the sum of the data density of the two clusters in that area. This can potentially make the overlapping area become another cluster. In some applications, we indeed want to recognize that overlapping space as a cluster, while that behavior is unexpected in other applications. Fig. 4 illustrates this confusion in the case of two strongly overlapping clusters. Besides the non-overlapping assumption, the gap statistic also assumes that there is no hierarchical clustering structure in the dataset. This means in the dataset; there is no cluster that consists of many smaller clusters. In addition, the gap statistics require a lot of computing power to compute the expected 𝑊𝑊 𝑘𝑘 under the null reference distribution 𝐸𝐸 𝑛𝑛 * {log(𝑊𝑊 𝑘𝑘 )}. It has to sample the null reference distribution B times (𝐵𝐵 ≥ 50), for each sample b, we run the clustering algorithm. In this case, the clustering algorithm is PAM, which takes 𝑂𝑂(𝑛𝑛 2 ) with n is the number of data points. In total, the complexity of the algorithm to estimate the 𝐸𝐸 𝑛𝑛 * {log(𝑊𝑊 𝑘𝑘 )} is 𝑂𝑂(𝐵𝐵𝑛𝑛 2 ). This would make it impossible to apply gap statistic on dataset with more than several thousands of data points.

The new gap

As described in the previous section, the gap statistic method has largely three limitations. However, we only focus on the overlapping issue to produce a new gap. The other limitation issues will be covered in the further research.

The 1stDaccSEmax Rule for Overlapping Clusters. The Tibs2001SEmax rule returns the smallest k such that the gap at that point has a significant chance (one standard error) to be higher than the next gap. As shown in the previous section, this rule is very sensitive to overlapped clusters. In fact, when there are overlapping clusters in the dataset, the gap does not decrease but slightly increase after 𝑘𝑘 = 𝐾𝐾 (where K is the real number of clusters in the dataset). This results in over-estimation of K. Therefore, instead of using the gap statistic directly, we propose to use the deceleration of the gap statistic (Dacc statistic for short). The Dacc is calculated as follows:

𝐷𝐷𝑎𝑎𝐷𝐷𝐷𝐷(𝑘𝑘) = [𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘) -𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘 -1)] -[𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘 + 1) -𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘)] = 2𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘) -𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘 -1) -𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘 + 1) Fig. 5 shows how the Dacc(k) statistic can be computed from the Gap(k) statistic.

We designed this statistic based on the insight that when k is going from 1 to K, the Gap(k) increases with constant or accelerated speed, up to the point where k = K. At that point, the Gap(k) will suddenly slow down its speed of increasing or start to decreasing (negative speed). Fig. 6 illustrates how the Dacc(k) looks like in different scenarios.

Conclusion

This study focuses on improving the gap statistic for the task of predicting the number of clusters k of a dataset. It identifies and demonstrates three main limitations of the gap statistic, including the overlapping clusters problem, the hierarchical clustering structure problem, and the big dataset problem. Based on these insights, we proposed the new technique to tackle the overlapping problem: the 1stDaccSEmax rule. The performance of the new method is evaluated with several synthetic datasets. It is believed that the performance of the new gap method would be shown to be better than all other traditional approaches. The further numerical experiments will be done on several real datasets with some other new techniques to overcome the other gap limitations.

Fig. 2 .

 2 Fig. 2. How the gap statistic works on a dataset with two well-separated clusters

Fig. 3 .

 3 Fig. 3. Overlapping clusters problem with gap statistic (a) the ovl2Gauss dataset: 400 data points in 2 dimensions that sampled equally from the two 2D Gaussian distributions: 𝒩𝒩 �� 0 0 � , � 1 0.7 0.7 1 �� and 𝒩𝒩 �� 4 0 � , � 1 -0.7 -0.7 1 ��. (b) gap statistic with Tibs2001SE rule suggests k = 3 instead of 2 for the ovl2Gauss.

Fig. 4 .

 4 Fig.4. Two strongly overlapping clusters can be correctly seen as one, two, or 3 clusters.

Fig. 5 .Fig. 6 .

 56 Fig.5. How to compute Dacc(k) from Gap(k) k k-1 k+1

Fig. 6f )Fig. 8 .Fig. 9 .

 6f89 Fig. 6f) In dataset with K clusters where some clusters strongly overlap each other: the definition between clusters becomes very fuzzy. Two strongly overlapping clusters can

  the expectation under a sample size of n from reference distribution, we estimate 𝐸𝐸 𝑛𝑛 * �log�𝑊𝑊 𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 �� by an average of B copies log (𝑊𝑊 𝑘𝑘 𝑛𝑛𝑛𝑛𝑛𝑛𝑛𝑛 ), each of which is computed from a Monte Carlo sample from reference distribution. Cluster the Monte Carlo samples into k groups and compute 𝑙𝑙𝑙𝑙𝑎𝑎𝑊𝑊 𝑘𝑘𝑘𝑘

ACKNOWLEDGMENT • This work was supported by the Industry-Academic Cooperation R&D grant (Dec. 24 2018 ~ Dec. 23 2019) funded by the LX Spatial Information Research Institute (SIRI). • This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (2017R1D1A1B03034475).

be correctly considered as one, two, or three clusters. Therefore, both Dacc and Gap statistic behave unpredictably in this case.

To take into account the sampling error occurring when estimating the expected 𝑊𝑊 𝑘𝑘 under the null distribution, I incorporate the standard error 𝑠𝑠 𝑘𝑘 to the 𝐷𝐷𝑎𝑎𝐷𝐷𝐷𝐷(𝑘𝑘) to get the 𝐷𝐷𝑎𝑎𝐷𝐷𝐷𝐷𝑎𝑎𝐸𝐸(𝑘𝑘) as follows:

𝐷𝐷𝑎𝑎𝐷𝐷𝐷𝐷𝑎𝑎𝐸𝐸 The Gap(k) chart can have multiple peaks, especially when the dataset has a hierarchical clustering structure. Therefore, instead of selecting k where the DaccSE(k) reaches its global maxima, we select the k where DaccSE(k) reaches its first local maxima. This is similar to the idea of searching for the first local maxima of the Tibs2001SEmax rule introduced in the original gap paper. This new rule is called the 1stDaccSEmax rule. Generally, the 1stDaccSEmax rule keeps looking for k with the highest positive DaccSE, with k sequentially running from k = 2 to k = kmax and stop at the point where 𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘) higher than 𝐺𝐺𝑎𝑎𝐺𝐺(𝑘𝑘) -𝑠𝑠 𝑘𝑘+1 . Fig. 8 shows how the 1stDac-cSEmax rule works in different situations.

Note that although the DaccSE(k) statistic does not define when k = 1, the 1stDac-cSEmax rule can still detect if there is no cluster in the dataset. This can happen in two situations, which are illustrated in Fig. 8. In Fig. 8b, Gap(1) > Gap(2) by a margin