A clustering algorithm for multivariate data streams with correlated components

Common clustering algorithms require multiple scans of all the data to achieve convergence, and this is prohibitive when large databases, with data arriving in streams, must be processed. Some algorithms to extend the popular K-means method to the analysis of streaming data are present in literature since 1998 (Bradley et al. in Scaling clustering algorithms to large databases. In: KDD. p. 9–15, 1998; O’Callaghan et al. in Streaming-data algorithms for high-quality clustering. In: Proceedings of IEEE international conference on data engineering. p. 685, 2001), based on the memorization and recursive update of a small number of summary statistics, but they either don’t take into account the specific variability of the clusters, or assume that the random vectors which are processed and grouped have uncorrelated components. Unfortunately this is not the case in many practical situations. We here propose a new algorithm to process data streams, with data having correlated components and coming from clusters with different covariance matrices. Such covariance matrices are estimated via an optimal double shrinkage method, which provides positive definite estimates even in presence of a few data points, or of data having components with small variance. This is needed to invert the matrices and compute the Mahalanobis distances that we use for the data assignment to the clusters. We also estimate the total number of clusters from the data.


Introduction
Clustering is the (unsupervised) division of a collection of data into groups, or clusters, such that points in the same cluster are similar, while points in different clusters are different.When a large volume of (not very high dimensional) data is arriving continuously, it is impossible and sometimes unnecessary to store all the data in memory, in particular if we are interested to provide real time statistical analyses.In such cases we speak about data streams, and specific algorithms are needed to analyze progressively the data, store in memory only a small number of summary statistics, and then discard the already processed data and free the memory (Garofalakis et al, 2016).Data streams are for example collected and analyzed by telecommunication companies, banks, financial analysts, companies for online marketing, private or public groups managing networks of sensors to monitor climate or environment, technological companies working in IoT, etc.In this framework, there are many situations in which clustering plays a fundamental role, like customer segmentation in big e-commerce web sites, for personalized marketing solutions, image analysis of video frames for objects recognition, recognition of human movements from data provided by sensors placed on the body or on a smartwatch, monitoring of hacker attacks to a telecommunication system, etc.

Related literature
The methods for cluster analysis present in literature can be roughly classified into two main families: probability-based methods (see e.g.Aggarwal and Reddy (2013)), which are based on the assumption that clusters come from a mixture of distributions, from a given family.In such case the clustering problem is reduced to the parameter estimation.These algorithms are well suited to detect the presence of non-spherical or nested clusters, but are based on specific assumptions on the data distribution, the number K of clusters is fixed at the very beginning, and, more important, they require multiple scans of the dataset to estimate the parameters of the model.Thus they can not be applied to massive datasets or data streams.
The second family of clustering algorithms is composed by distance-based approaches.Given a dataset of size n, grouped into K clusters, such methods have usually the goal to find the K centers of the clusters which minimize the mean squared distance between the data and their closest centers.These methods usually take different names depending on the type of considered distance.If the Euclidean distance is used, the corresponding method is the classical and very popular K-means method (see e.g.Jain (2010)), which is probably the most diffused clustering algorithm, because of its simplicity.Anyway the exact solution of the minimization problem connected with K-means is NP-hard, and only local search approximations are implemented.The method is sensitive to the presence of outliers and to the initial guess for the centers, but improvements both in terms of speed and accuracy of the algorithm have been implemented in K-means++ (Arthur and Vassilvitskii, 2007), which exploits a randomized seeding technique.Unfortunately both the classical K-means and the K-means++ algorithms require multiple scans of the dataset or a random selection from the entire dataset, in order to solve the minimization problem.Since data streams can not be scanned several times and we can not (randomly) access to the entire dataset all together, also these methods are not suitable for clustering data streams.
When the elements to be clustered are not points in R d but more complex objects, like functions or polygons, other clustering algorithms are used, like PAM, CLARA, CLARANS, (Kaufman and Rousseeuw (1990); Ng and Han (2002)), which are based on non-euclidean distances defined on suitable spaces.These methods are looking for medoids, instead of means, which are the "most central elements" of each cluster and are selected from the points in the dataset.Also these algorithms can not be efficiently applied to analyse data streams, since they either require multiple scans of the sample, or the extraction of a subsample to identify the centroids or medoids, then all data are scanned according to such identification and the medoids are not any more updated with the information coming from the whole dataset.Actually such popular methods are suited for data which are very high dimensional (e.g.functions) or for geometrical or spatial random objects, but not for datasets with an high number of (rather small dimensional) data.
The key element in smart algorithms to treat data streams is to find methods to represent the data with summary statistics which are retained in memory, while the single data are discarded.Such summary statistics must be updated when each new observation, or group (chunk) of observations, is processed, since a second scan of the data is not allowed.This strategy to analyse data streams is followed in O'Callaghan et al (2001), where the STREAM algorithm is proposed as an extension of BIRCH (Zhang et al, 1996).The STREAM method solves a so called K-Median problem, which is a generalization of K-means where the Euclidean distance is replaced by a general distance.The performance of the STREAM method with respect to computational costs and quality of the clustering, measured in terms of sum of squared distances (SSQ) of data points from the assigned clusters centers, is also studied, in particular in comparison with K-means, providing good theoretical and experimental results.In the STREAM algorithm the number K of clusters is not specified in advance, and is evaluated by an iterative combination between SSQ and the number of used centers.The main defect of the STREAM algorithm is that it uses a global metric D on the space of the data points, and thus does not take into account that different clusters may have different specific variability.Further the metric D is supposed completely known and is not estimated from the data.
In many situations the quality of the clustering is improved if a local metric is used.A local metric is a distance which takes into account the shape of the "cloud" of data points in each cluster to assign the new points (see Figure 1).A first attempt to use a local distance is given by the Bradley-Fayyad-Reina (BFR) algorithm (Bradley et al (1998); Leskovec et al ( 2014)), which solves the K-means problem by using a distance based on the variance of each component of the random vectors belonging to the different clusters.The BFR algorithm is based on the assumption that the clusters' distribution results from a mixture of multivariate normal distributions, whose parameters are estimated from the data streams.The BFR Algorithm for clustering is based on the definition of three different sets of data: a) the retained set (RS): the set of data points which are not recognized to belong to any cluster, and need to be retained in the buffer; b) the discard set (DS): the set of data points which can be discarded after updating the summary statistics; c) the compression set (CS): the set of summary statistics which are representative of each cluster.
Each data point is assigned to one of these sets on the basis of its local distance from the center of each cluster.Here the Mahalanobis distance is used, computed with respect to the sample covariance matrix of each cluster.
The main weakness of the BFR Algorithm resides in the assumption that the covariance matrix of each cluster is diagonal, which means that the components of the analyzed multivariate data should be uncorrelated.With such assumption, at each step of the algorithm only the means and variances of each component of the clusters centers must be retained, reducing thus the computational costs.Further, in this setting the estimated covariance matrices are invertible even in presence of clusters composed just by two p-dimensional gaussian data points.Anyway such assumptions geometrically imply that the level surfaces (ellipsoids) of the gaussians including the data points in each cluster should be oriented with main axes parallel to the reference system.

Aims and overview of the paper
We here propose a method to clusterize data streams, using a local metric which is estimated in real time from the data.Such metric is based on the Mahalanobis distance of the data points from each cluster center c i , computed using an estimator of the covariance matrix of the corresponding i − th cluster.In the following we will always represent vectors as column vectors and we will assume that our data are vectors in R p .
Definition.Let x be a data point and c i be the center of the i − th cluster.Assume that the elements of the i − th cluster come from a population having covariance matrix Σ i .Then the Mahalanobis distance of x from c i is given by We assume that the data points are vectors in R d with correlated components and we thus estimate all the terms of the covariance matrix of each cluster, including the off diagonal ones.We use the sample mean of each cluster as centers c i .
We divide the data in the same three sets defined in the BFR algorithm, we don't fix a priori the number K of clusters, and we evaluate and update such number using a density condition.Thus in our procedure from time to time new clusters will be formed, composed only by a few data points, not sufficient to obtain a positive definite estimate of the corresponding covariance matrix using the classical sample covariance estimator.We thus use an optimal double shrinkage estimator of the covariance matrix, which provides always positive definite matrices, that are then inverted to compute the Mahalanobis distance.
In our setting we will relax a little bit the assumption of gaussianity stated in the BFR algorithm, assuming that the data come from a mixture of "bell shaped" distributions, but possibly having a bigger multivariate kurtosis (i.e.fatter queues) than a gaussian.
Our algorithm is thus an improvement of the BFR algorithm, relaxing some of its assumptions.Since with our method also the covariance terms of the clusters must be retained, there is an increase in the computational costs with respect to BFR, but such increase can be easily controlled and is affordable if the processed data are not extremely high dimensional.Therefore our algorithm is targeted to problems with data streams composed by data points of "medium" dimension, i.e. a dimension not so small to apply visualization techniques to identify the clusters (2D or 3D problems), which usually work better, but much smaller than the number of available data.
The paper is then structured as follows: in Section 2 we face the problem of the estimate of the covariance matrix of each cluster.We modify a Steinian linear shrinkage estimator in order to obtain a positive definite estimator of the covariance matrix, which can be applied also to non-gaussian cases, and which can be incrementally updated during the data processing.In Section 3 we introduce the summary statistics that will be retained in memory for each cluster, and we show that they can easily be updated when new data streams are processed.We then describe the way by which the data points are assigned to the three sets RS, CS, DS.In Section 4 we describe the secondary compression, that is the way by which the points in RS and CS can be merged to pre-existing clusters or are put together to form new clusters.In Section 5 we apply our method first to synthetic data, and we compare heuristically its performances with the case in which the data points are assumed to have uncorrelated components, like in the BFR algorithm.We then apply our method to cluster the real dataset KDD-CUP'99 ( http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html), a network intrusion detection dataset, that was used also to test the STREAM algorithm.We apply our algorithm to all the variables in the dataset which are declared continuous.Actually some of such variables have a very small variance; anyway our optimal double shrinkage estimator of the covariance matrices of the clusters guarantees positive definite estimates also in this situation, stabilizing thus the local Mahalanobis distances that we use in our procedure.The results are coherent with the structure of the dataset, whose data should be divided into 5 clusters, as we obtain.
In this paper we don't study the asymptotic properties of our algorithm, but we limit ourselves to show heuristically that our algorithm provides better results of other methods to cluster data streams present in literature, just with a small increase in the computational costs.

The covariance matrices of the clusters
Our algorithm is based on the Mahalanobis distance, it is hence crucial to estimate the covariance matrices of the clusters in an optimal way.Let us first observe that when a new cluster is formed, it contains too few data points to obtain a positive definite estimate of the covariance matrix, using the sample covariance matrix, at least until N ≤ p, where N is the number of data in a cluster and p the data points dimension.
To solve this problem in an optimal way, we exploit the optimal double shrinkage estimator given in (Ikeda et al, 2016, Equation (3.6)) by where S is the sample covariance matrix, D S is its diagonal matrix, I p is the identity matrix of order p, and 0 ≤ λI + λD ≤ 1 are weighting the convex combination of the three matrices.This estimator is optimal in terms of quadratic loss (Himeno and Yamada, 2014;Touloumis, 2015;Ikeda et al, 2016), and it leads to covariance matrix estimators that are non-singular, well-conditioned, expressed in closed form and computationally cheap regardless of p.Therefore, in these terms, it is the optimal choice among the possible alternatives, where the first term (1 − λI − λD ) should be initially settled close to 0, and then its value is increasing to 1 when N → ∞.We note that when λI = 0 and λD = 1 we obtain the local distance used in BFR.In (Ikeda et al, 2016), λI , λD are given as functions of the quantities where tr[Σ 2 ] and tr[Σ(Σ − D Σ )] are unbiased estimators of the corresponding quantities tr Σ is the true covariance matrix of the considered cluster, and D Σ its diagonal matrix.Unfortunately (see Touloumis (2015); Ikeda et al ( 2016)), both these estimators are based on the scalar statistics proposed by (Himeno and Yamada, 2014), where is the centroid of the considered cluster, composed by N data points.In data stream framework, we note that is not a function of few summary statistics, which can be updated when a new data point is added to the cluster.In fact, xN is changing with N and Q (N ) must be then recomputed, due to the quadratic term in its definition, using all the data in the cluster when a new point is added.To overcome this problem, we prove in the following section the existence of two unbiased estimators for tr[Σ 2 ] and tr[Σ(Σ − D Σ )] based on the following statistics Q N : The key point is that Q N is defined recursively, and it is a function of Q N −1 , the new added point, and the centroid of the cluster at the time of the update.
In the following we will describe the details of our method and the assumptions that must be satisfied to apply it.

A model for the estimate of the covariance matrices
Our dataset is given by a sequence of p-dimensional vectors x 1 , x 2 , . ... Each observation x n is independent on the others and, if belonging to the cluster k, it is generated as where µ k is the mean vector and Σ k ) is strictly positive definite.The following hypothesis of uncorrelation is assumed on the first four moments: for any integers γ 1 , . . ., γ q satisfying 0 ≤ q 1 γ i ≤ 4, and where z n,i is the i-th component of the vector Assume that the sequence x 1 , x 2 , . . .belongs to the same cluster with Σ is formed by independent vectors with null expectation.Then, as a consequence of ( 5), we have that Moreover, E[y i y j y k y l ] = 0 only in the following situation: where κ 11 is defined in Himeno and Yamada (2014) as Note that κ 11 = 0 for gaussian data, thus it is an indicator of deviation from gaussianity in terms of kurtosis.In case of gaussian data its estimation can be neglected (Fisher and Sun, 2011).
when (i = k) = (j = l) the same as above, since y k y l = y l y k , hence Lemma 1 As a consequence of (6d), Lemma 2 As a consequence of all the relations (6), Lemma 3 As a consequence of (6b), y N y N y i y j = (N − 1)(trΣ) 2 .

Optimal shrinkage estimation
We now use the previous results to solve the problem of finding the optimal estimates of λI , λD in (1), as a function of the statistics S (sample covariance matrix of the data in the same cluster), of Q N given in (4), and of two quantities S N and T N that can be updated inductively.As can be seen in ( 2), the problem here is the unbiased estimation of the terms tr[Σ 2 ] and tr(Σ 2 ) − tr(D 2 Σ ).The derivation of this estimate is given in the next section, after a technical result given hereafter.
We may use the following additional relations in our estimates (Himeno and Yamada, 2014;Ikeda et al, 2016) once we have recalled that the quantity R N is negligible (see, again, Touloumis (2015); Ikeda et al (2016)).
When the data are distributed as gaussians, a direct estimation without κ 11 based on (7a-c) may be done (see Fisher and Sun (2011)), since κ 11 = 0.When this is not the case, we may use the statistics Q N already introduced in (4) and we will prove in Lemma 4 that where S N and T N are two quantities that may be simply calculated inductively as: if a cluster is made by merging two clusters of N 1 and N 2 points; (9) if a new point is added to a cluster of N − 1 points; T N1 + T N2 if a cluster is made by merging two clusters of N 1 and N 2 points; (10) Lemma 4 With the notations of (3), ( 4), ( 9) and (10) we have if a new point is added if a cluster is made by merging two clusters of N 1 and N 2 points; and hence See the Appendix for the proof.
The system composed by (7a-c) and ( 8) may be read as Now, the matrix A may be shown to be invertible, and hence Ẑ = BA −1 X is a linear (in X) unbiased estimator for Z, since For sake of completeness, we give here the elements of the matrix

Summary statistics and primary data compression
In this section we define the summary statistics that will be retained in memory for each cluster and we describe the first phase of our clustering procedure.As in the BFR algorithm, we first perform the primary data compression, that is the identification of items which can be assigned to a cluster, and then discarded (Discard Set, DS), after updating the corresponding summary statistics contained in the Compression Set CS.Data compression refers thus to representing groups of points by their summary statistics and purging these points from RAM.In our algorithm, like in BFR, primary data compression will be followed by a secondary data-compression, which takes place over data points in the Retained Set (RS), not compressed in the primary phase.Assume that data points x 1 , . . ., x N ∈ R p , N ≥ 2 must be compressed in the same cluster.We will retain only the following summary statistics and the statistics Q N , S N , T N defined in ( 4),( 9),(10), respectively.
In particular the statistics s N are needed to compute the sample means xN = 1 N N i=1 x i , that are used as clusters centers, while the matrices Σ N are used to compute the unbiased sample covariance matrices of the clusters S = 1 , which are needed, together with Q N , S N , T N , to compute the optimal double shrinkage estimators described in the previous section.
The summary statistics (11) can also be easily updated when a new data point x N +1 must be added to the cluster, without processing again the already compressed points.In fact while the other summary statistics have already been defined recursively.
Note that the matrix Σ N is symmetric, thus at each step of the algorithm we have to retain in memory only p(p+1) 2 + p + 4 = p 2 2 + 3 2 p + 4 summary statistics for each cluster, where p is the dimension of the data points.Thus, in case of K clusters, our computational costs are of the order of Kp 2 .In addition, note that we should simply sum the corresponding statistics if we want to merge two clusters.
Similarly to the BFR algorithm, in order to assign a point to a cluster we use the squared Mahalanobis distance from its center (sample mean), i.e. we assign a new data point x to cluster h with center xh and estimated covariance matrix Ŝh , if h is the index which minimizes Differently from the BFR algorithm, here we estimate the covariance matrices of the clusters with the optimal double shrinkage estimators described in the previous section.In order to avoid the inversion of a matrix and thus to reduce the computational costs, we observe that the Mahalanobis distance between two points x, y, computed with respect to a covariance matrix S, can be rewritten as follows (see e.g.(Rencher, 1998, Expression A.7.10)): In our algorithm we will actually use expression ( 12) for the computation of all the Mahalanobis distances.We also compare x with each point x o in the retained set RS, if any, by computing where ŜP matrix is the pooled covariance matrix based on Ŝh of all the K clusters: and where n h is the number of points in cluster h.With ŜP , we emphasize the weighted importance of directions that are more significant for the clusters when we compute the distance between two "isolated" points.Since the retained set contains the points which do not belong clearly to one specific cluster, with this comparison we check if they can be aggregated with the new incoming data, to form new clusters.We then approximate locally the distribution of the clusters with a p−variate Gaussian and we build confidence regions around the centers of the clusters (see Hotelling (1931)).Following the approach stated in Bradley et al (1998), which is motivated by the assumption that the mean is unlikely to move outside of the computed confidence interval, we perturb xh by moving it in the farthest position from x in its confidence region, while we perturb the centers of the other clusters by moving them in the closest positions with respect to x and we check if the cluster center closer to x is still xh .If yes, we assign x to cluster h, we update the corresponding summary statistics and we put x in the discard set; otherwise, we put x in the retained set (RS) (see Figure 2).If in the first comparisons the point x is closer to a point x o of the retained set than to any cluster, we form a new secondary cluster with the two points if x o remains the closest to x after the centers' perturbation.In this case we add the corresponding summary statistics to the compressed set CS, and we put x and x o in the discard set.Otherwise we put x and x o in RS (see Figure 3).
Let us see the procedure of centers' perturbation in deeper detail.

Confidence regions
It is well-known (Hotelling (1931)) that a confidence region for the mean µ based on x and Ŝ may be based on the Hotelling's T -squared distribution where F p,n−p is the F-distribution with parameters p and n − p.
Then, if we denote by CI k the confidence region for the mean of cluster k, i.e.
then the perturbation p k (x) for the data point x is Denoting by t α = T 2 p,n−1 (1 − α), if we introduce a Lagrange multiplier λ * , the problems of minimization or maximization stated in the definition of p k (x) can be solved by differentiating the following lagrangian form The resolution ∇ µ L = 0 gives µ = x−λx k 1−λ , where λ = nλ * .In particular, the optimal µ is the linear combination of x and x k in CI k which is farther from x or closer to x, when k = j or k = j, respectively.The constrain reads Denoting by Summarizing we obtain the following perturbations of the clusters centers, referred to the data point x,

Secondary data compression
The purpose of secondary data compression is to identify "tight" sub-clusters of points among the data that we can not discard in the primary phase.In Bradley et al (1998) this is made using the euclidean metric.We adopt a similar idea, but we use a local metric, based on the Mahalanobis distance.We exploit a technique based on hierarchical clustering, mimicking the Ward's method (Gan et al, 2007;Legendre and Legendre, 2012).
Given two clusters h 1 and h 2 with n h1 ≥ 2, n h2 ≥ 2 points, and centroids xh1 and xh2 , respectively, then the squared Mahalanobis distance of one centroid to the other cluster may be measured as ∆ 2 Ŝh i (x h1 , xh2 ), i = 1, 2. Accordingly, to decide whether two clusters are close or not, we compare the weighted combination of those distances with the the squared Mahalanobis distance ∆ 2 Ŝh 1 h 2 (x h1 , xh2 ) of the two centroids, evaluated with the pooled covariance matrix of the two clusters Ŝh1h2 = . The distance between a single retained point x and a cluster h is computed by the squared Mahalanobis distance ∆ 2 Ŝh (x, xh ) between the point and the cluster centroid, based on the estimated covariance matrix of the cluster, while the distance between two Table 1 Results of the application of our proposed algorithm (PA) and of the BFR algorithm to synthetic data.We call chunk the number of processed data out of which we apply secondary compression.
retained points x 1 , x 2 is computed by their squared Mahalanobis distance ∆ 2 ŜP (x 1 , x 2 ) based on the pooled covariance matrix (13) of all the clusters.
Based on the hierarchical tree built with such distances, we sequentially merge two clusters or points only if a suitable density condition is fulfilled.This condition is different for the different types of merging that we can perform: we merge two clusters h 1 and we merge a retained point x and a cluster h if ∆ 2 Ŝh (x, xh ) < θ 1 (tr(S h )); we merge two retained points x 1 and Here θ i , i = 0, 1, 2, are thresholds, chosen by the user.For what concerns θ 2 , we suggest to use a significant quantile of the χ-square distribution that arises under the null hypothesis H 0 : the retained points come from a gaussian distribution with covariance matrix given by the pooled covariance matrix (13) of all the clusters.
5 Results on simulated and real data and discussion

Results on synthetic data
Synthetic data were created for the cases of 5 and 20 clusters.Data were sampled from 5 or 20 independent pvariate Gaussians, with elements of their mean vectors (the true means) uniformly distributed on [−5, 5].The covariance matrices were generated by computing products of the type Σ = U HU T , where H is a diagonal matrix with elements on the diagonal distributed as a Beta(0.5, 0.5) rescaled to the interval [0.5, 2.5], and U is the orthonormal matrix obtained by the singular value decomposition of a symmetric matrix M M T , where the elements of the p × p matrix M are uniformly distributed on [−2, 2].In either cases of 5 or 20 clusters, we generated 10.000 vectors for each cluster, having dimensions p = 5, 10, 20.
This procedure guarantees that these clusters are rather well-separated Gaussians, in particular for higher vector dimensions.
We applied both our procedure and the BFR algorithm to these synthetic data, to compare the performance of the two methods.In both cases, we computed the secondary data compression once out of 25, or out of 50 data points.In the tests on data from 20 clusters we started from a lower number of initial clusters (equal to 10), in order to check the ability of our algorithm to detect the correct number of clusters.The results are reported in Table 1.
We note that the number of clusters is sometimes underestimated by our method, in particular in the case of 20 clusters.In such cases, if the point clouds in different clusters are gathered in rather close ellipsoids, then the correct detection of the clusters may be more difficult.Anyway in all cases the estimates provided by our algorithm are equal or better than those obtained with the BFR algorithm.
We also point out that in the case of 20 clusters with p = 10, and secondary compression performed once out of 50 processed data, the overestimation of the number of clusters obtained with our algorithm is compensated by the presence of two small clusters, composed by a few hundreds of data points, which can then be revisited as groups of outliers.Anyway also in this case our results are better than those obtained with BFR.
The method seems to be sensitive to the frequency of the secondary compression only in presence of many clusters.
Note that our method gives always a correct estimation of the number of clusters in all cases with 5 true clusters, while the BFR method overestimates the correct number in particular when the data dimension is small (p = 5).This is reasonable since in lower dimensional spaces the shape and orientation of the point clouds must be correctly estimated and taken into account to identify the clusters in a proper way (see Figure 4).We tested also cases with bigger values of p, but in such cases both algorithms are able to detect the correct number of clusters, in an equivalent way, since a few clusters in high dimensional spaces are almost always well separated, because of "curse of dimensionality" reasons.

Results on a real dataset
We applied our algorithm to a real dataset to detect network intrusions.Detecting intrusions is a typical data streaming problem, since it is essential to identify the event while it is happening.In our experiments we used the KDD-CUP'99 (http://kdd.ics.uci.edu/databases/kddcup99/kddcup99.html)intrusion detection dataset which consists of two weeks of raw TCP dump data.This dataset is related to a local area network simulating a true Air Force environment with occasional attacks.Variables collected for each connection include the duration of the connection, the number of bytes transmitted from source to destination (and viceversa), the number of failed login attempts, etc.We applied our algorithm to the 34 variables that are declared to be continuous.Some of these variables actually are almost constant, giving thus an estimated zero sample variance in many clusters.In such situation, if the BFR algorithm is applied, singular covariance matrices are estimated for some clusters.Consequently the Mahalanobis distance becomes unstable.Our optimal double shrinkage estimators are thus necessary to overcome this instability, and as a byproduct, they can take into account the deviation of the kurtosis from the Gaussian case.
The same dataset was analysed in O'Callaghan et al ( 2001), via the STREAM algorithm, but they used the Euclidean distance, which is a global distance that gives the same importance to all the variables.
We obtained stable results.We applied the secondary compression every 100 data, starting from 4 clusters composed by less than 20 points.We observed the presence of 6-8 big clusters starting from about 100000 processed data.We processed about 646000 data, ending with 5 big clusters, composed by the following number of points: 133028; 121661; 242206; 53235; 95977.Note that we detected the final correct number of clusters, since in this dataset there are four possible types of attacks, plus no attacks.The four types of attacks are denial-of-service; unauthorized access from a remote machine (e.g.guessing password); unauthorized access to local superuser (root) privileges; surveillance and other probing (e.g., port scanning).
In Figure 5 we show the effectiveness of secondary compression on the stabilization of the number of clusters.Actually when the number of identified clusters is bigger than 8 after secondary compression, the exceeding ones are formed just by a few points, and can then be reinterpreted as groups of outliers.For example when 300113 data have been processed, we find 15 clusters composed respectively by the following number of points: 95141; 22451; 50098; 79943; 30683; 11834; 7762; 1228; 712; 118; 100; 33; 4; 4; 2. Note that 7 out of 15 clusters are quite small, containing less than 1000 data points.

Conclusion
We have introduced a new algorithm to cluster data streams with correlated components.Our algorithm in some parts imitates the BFR algorithm, since, like BFR, it uses a local distance approach, based on the computation of the Mahalanobis distance.In order to compute such distance, positive definite estimators of the covariance matrices of the clusters are needed, also when the clusters contain just a few data points.We obtained such estimators by considering a Steinian double shrinkage method, which leads to covariance matrix estimators that are non-singular, well-conditioned, expressed in a recursive way and thus computable on data streams.Further such estimators provide positive definite estimates also when some components of the data points have a small variance, or the data distribution has a kurtosis different from the Gaussian case.
We applied both our proposed method and the BFR algorithm to synthetic gaussian data, and we compared their performance.From the numerical results we conclude that our method provides rather good clustering on synthetic data, and performs better than the BFR algorithm in particular in presence of few clusters in spaces of rather low dimension.This is reasonable since the BFR algorithm approximates the "clouds" of data with ellipsoids having axes parallel to the reference system, and this leads to a wrong classification when the clusters are elongated, not much separated, and with axes rotated with respect to the reference system.In such situations our algorithm is able to capture in a more proper way the geometry of the clusters, and thus improves the classification.
Anyway the secondary compression could be possibly improved by applying some incremental modelbased technique (see Fraley et al (2005)), but modified in such a way to avoid multiple scans of the sample.
We also applied our algorithm to a real dataset, obtaining good results in terms of correct identification of the number of clusters, and stability of our algorithm.
The advantage of our algorithm with respect to other methods present in literature, like BFR or STREAM, is that it relaxes the assumptions on the processed data streams, and can thus be effectively applied to a wider class of cases, on which it performs better.In the cases where the assumptions of the other methods are satisfied, our algorithm provides equivalent results.It can then be systematically substituted to other methods to analyze data streams, in all cases in which the data points are not too much high dimensional.

Proof of Lemma 4
For N = 2, as a consequence of the model: Since E[(y 1 y 1 )] = trΣ, by (6a) and (6d), we obtain the first part of the thesis.
Let us add a point to a cluster of N − 1 points.We obtain

Fig. 1
Fig.1A typical situation in which the shape of the point clouds must be taken into account for the points assigment: two clusters of gaussian points are represented, both composed by 100 data points.The new data point P is closer to the cluster center A if we use a global metric, like e.g. the Euclidean distance, but the point is more likely to belong to the cluster centered at B

Fig. 2
Fig.2The 2 possible situations after the perturbation of the clusters centers.Left: the point x is assigned to cluster h; right:x is moved to RS

Fig. 4
Fig.4Example of a typical situation in 2D: the BFR algorithm is approximating the shape of the clusters using ellipsoids parallel to the main axes.Our proposed algorithm (PA) is using ellipsoids with the correct orientation, thus the uncertainty region (overlapping of the ellipsoids) is reduced.

Fig. 5
Fig. 5 Effect of the secondary compression.Red line: number of clusters obtained before secondary compression.Blue line: number of clusters obtained after secondary compression.Secondary compression is applied once out of 100 iterations.