webtrafficHMM: a State Model for Mixed-Type Variables

Models web traffic (or similar data) using Gaussian Mixture Hidden Markov Models. Such models are notoriously difficult to fit. What follows is a proposal to obtain good fits for the Gaussian Mixture Hidden Markov Model (GMMHMMs) for a class of mixed-type features. The GMMHMM is adaptive, setting the number of states and mixtures at a level appropriate for the complexity exhibited by the input data.

1. Introduction

HMMs are powerful models that attempt to infer hidden states based on changes in observed behaviour. For example, suppose we wish to infer a person's state from a set of acceleratorometer measurements (source: Dai et al., 2017. Recurrent Hidden Semi-Markov Model. ICLR). The HMM does a remarkable job of isolating the various states based on the disparate behaviour that occurs in each state.

In much the same way, we seek to infer modes of network behaviour using either Proxy or Netflow logs. The ability to infer states enables detection of a subtle form on anomaly known as "collective anomalies". There are three types of anomalies: Collective Anomalies are notoriously difficult to detect. One way to detect such anomalies in time series data is to look for anomalous inferred states derived from an HMM.

Network Features are often of mixed type. That is, features exhibit both continuous-like and discrete-like behavior and their cumulative distribution functions (cdfs) consist of regions where the cdfs are strictly monotonically increasing whilst also being occasionally punctuated with discontinuities. For the mixed-type features prevalent in Network Data Logs, Gaussian Mixture Model HMMs (GMMHMMs) provide a promising modelling framework: the traditional Gaussian Mixture Models (GMM) can naturally capture the continuous-like behavior of the feature-in-question whilst any discrete-like behavior can be captured by setting one of the mixture components of the GMM to a Dirac Delta function centered around the value in question. The diagram below shows a GMM fit to observations from a mixed-type feature.

Furthermore, a GMMHMM can allow the process modeled to evolve over time, exhibiting dramatically different behavior at different points in time. The ability of an HMM to capture the behavior of the various states is one of its key selling points as it enables one to infer the modes-of-behaviour that the system exhibits. This can be valuable as, for example, examining such behavior modes could potentially enable one to pick up entrenched malicious activity in a network – that is, malicious network activity that has occurred for long enough duration so that any anomaly detection algorithm would simply view the behavior as normal activity. A protracted period of such malicious activity could in principle be captured by the HMM within one of its states. A subsequent closer examination of the states could then show that unusual combinations of network features are characteristic of the "malicious state" (e.g, unusual flow of bytes on port 3389, which is meant primarily for remote desktop connections or very small inter-arrival times between traffic events).

Unfortunately, the high level of flexibility of the Gaussian Mixture Model HMM (GMMHMMs) comes at the expense of a high sensitivity to the model initialization. This means that any traditional HMM estimation procedure is virtually guaranteed to yield an ill-fitted GMMHMM and consequently a model which gives a highly distorted representation of the traffic dynamics. For such a model, each state could well be a mixture of disparate modes of behavior that have no real connection to each other. We look to initialize the GMMHMM in a way so that the final model is likely to yield states that correspond to some concrete behavior of the network.

2. Method

Four guiding principles inform the way we go about initializing the GMMHMMs:

1. The states of the HMM are persistent. That is, we expect the network to remain in all but the rare states for protracted periods of time.
2. Each discrete-like behaviour by the network will manifest itself as a High Density Region (HDR) within the network feature space and each HDR is likely to arise from a single state of the HMM. At the same time, it is possible for multiple HDRs to arise from the same state if the HDRs tend to occur simultaneously. Enabling multiple HDRs to arise from the same state can make the states of the GMMHMM more persistent.
3. Each of the mixture components in any HMM state should have its relative levels of activity remain roughly constant over the time points that its state is active. Let . The GMMHMM implicitly assumes that is constant through time for all states and all mixture components. However, due to the sensitivity of the GMMHMMs to their initialization, a fitted GMMHMM can sometimes have sample estimates for 's that vary dramatically over time for a subset of state-mixtures. We need to explicitly guard against time-inhomogeneous 's for the GMMHMM.
4. States should not occur simultaneously across the entire observation period. Any errant states should be aggregated into a single state.

The first principle causes us to consistently initialize a GMMHMM with a strongly diagonal transition matrix. The second principle is foundational in the search for an appropriate initialization of the GMMHMMs. This leads us to search for the "discrete-like" behaviour in the traffic data in order to use them to seed some of the mixture components for the HMM states. We expect such discrete-like behaviour to present itself as high density regions (HDRs) in the feature space. We have a initialization pipeline for a GMMHMM that is composed of three, interdependent initialization procedures, each of which has differing levels of emphasis to the four guiding principles. These are:

• HDR-seeded Initialization & Fitting (Initialization Stage): This is the first stage of the GMMHMM initialization. Particular attention is paid in this stage to finding High Density Regions (HDRs) and using them to initialize the appropriate states within the GMMHMM. A strongly persistent transition matrix is used in the initialization and the GMMHMM is subsequently fitted to the data and HDRs that occur simultenously over the entire observation period are placed as separate mixture components in the same state. As such, this stage pays particular attention to adhering to the first, second and fourth principles shown above.
• Iterative, Missing-HDR-seeded Initialization & Fitting (Iteration Stage): This stage takes the last fitted GMMHMM and partitions the data according to the GMMHMM state that each observation is predicted to lie in. We then search for undiscovered HDRs in each of these partitions. As the data partitions should have more homogeneous behaviour than the pooled data, this stage will often find HDRs that have been missed in earlier stages. The resulting expanded set of HDRs are subsequently used in conjunction with a strongly persistent transition matrix to initialize the corresponding states of a GMMHMM and the model is fitted to the data. This stage assumes that any newly-found HDRs lie in the same state as the partition which they form a part of. This stage attempts to adhere to the first and second principles, though it is possible for a missing HDR to lie in multiple partitions and hence to manifest itself in multiple states of the GMMHMM.
• Enforcing Time-Homogeneity in Mixture Initialization & Fitting (Enforcement Stage): This stage takes the state-mixture combinations from an existing GMMHMM and checks that the mixture components tend to occur at the same times. The mixture components are reorganized so as to ensure that state-mixture components that occur roughly simultaneously are placed in the same state. Once the mixtures components are reallocated between the states, mixtures with very similar distributions and which are allocated to the same state are amalgamated into a single mixture. As with the other two stages, the resulting set of mixtures are used to initialize a GMMHMM and fitted to the data. This stage pays particular attention to the third principle described above.

It was noted that a pipeline using the aforementioned three stages could occasionally yield a GMMHMM with a subset of states that occur virtually simultaneously. This is due to some of the state-mixture components in the Enforcement stage having changed considerably in the final model fit. Such state-mixtures that originate from different simultaneous states were adjudged in the Enforcement Stage to not be simultaneous but were since changed in the terminal model fit so that in the final model these same state-mixtures occur simultaneously. Such state-mixture components give rise to simultaneous states, which is undesirable as it makes the states less persistent. Thus, in order to more closely adhere to the first principle, we choose to aggregate such simultaneous states into a single state whose mixture components are determined by the dominant mixture components across the member states. We call this the "Aggregation Stage". The Aggregation Stage pays particular attention to the fourth principle.

The diagram below shows hows the aforementioned stages may be arranged in order to ensure a good fit to the GMMHMM.

3. Pipeline Details

We describe the technical details of each of the four stages as well as how we go about obtaining a regularized fit for the GMMHMM.

3.1. Initialization Stage

We present an overview of the Initialization Stage:

We describe each of these stages in more detail below.

Algorithm

1. Model Initialization. We start by setting the means, the covariances and the model weights using the hmmlearn package's initialization function. This effectively sets the covariances equal to the sample covariance of the entire data set, the weights as uniformly equal and the means as the centroids returned by the k-Means algorithm with k equal to the number of states. k-Means is known to be very sensitive to the initialization condition and hence this sensitivity affects the initialization stage.
2. Fit GMM to find HDRs. A Gaussian Mixture Model (GMM) is fitted to the time series. In the process, we "flatten" the time series and assume the observation at each time point is independent and identically distributed. The HDRs are defined to be those Gaussian components that have covariances below some pre-determined threshold (thus exhibiting discrete-like behaviour for the feature(s) in question) and where the number of observations predicted to lie in the Gaussian component exceeds some pre-determined frequency threshold. The means, covariances and weights of the HDRS are returned. As with the k-Means algorithm, the GMM is very sensitive to the initialization. Hence the HDRs found can vary dramatically with each run of this step.
3. Perform Spectral Clustering on the HDRs. A similarity matrix is calculated based on the relative simultaneity of each pair of HDRs. The observation window is split into time bins whose width is determined by both the width of the observation period and the standard deviation of all the HDR's recorded time indices (as suggested by Scott, 1979). For each HDR, we calculate a vector showing the relative activity of the HDR in each time bin. A cosine similarity can be calculated between each pair of HDRs using their activity vectors and this can be used as the basis to cluster the HDRs. All HDRs placed in the same cluster will be part of the same state in the GMMHMM. Each of the HDRs in the same cluster will initialize a separate mixture component of the same state. As each HDR cluster will be used to initialize a separate state of the GMMHMM, we check if the number of HDR clusters exceed the maximum number of states allowed for the GMMHMM. In such a case, we redo the spectral clustering, making it easier for previous clusters to merge together. In much the same vein, we also need to check that both the number of clusters and the maximum cluster size are a reasonable size compared to the maximum number of states and mixtures allowed for the GMMHMM. If not, we adjust the thresholds that define the HDRs and redo the exercise of finding the HDRs by re-fitting a GMM to the flattened time series. Note that the clustering solution yielded by this stage allows us to adaptively reduce the complexity of the GMMHMM (i.e. the number of state and mixture components) if the dynamics of the time series are relatively simple. This is made possible as the majority of the model complexity of the GMMHMM is intended to account for the discrete-like behaviour in the network traffic. Thus, if there are not many HDRs discovered, we can simplify the GMMHMM.
4. Initialize states of GMMHMM. The HDR clustering solution is used to change a subset of the means, covariances and weights initialized by the model initialization procedure of hmmlearn. We make each state of the GMMHMM correspond to one of the HDR clusters, with the state's mixtures components initialized by the means, covariances and weights of the HDRs within the cluster. By design, at least one mixture component in each state as well as one state will not be initialized using the characteristics of the extracted HDRs. This gives the GMMHMM some flexibility to learn time series behaviour that does not arise from any HDR. However, this explicit allowance for continuous behaviour may make this stage run into numerical instabilities when analyzing a time series with very little continuous-like behaviour.

3.2. Iteration Stage

As the initialization stage uses both the k-Means algorithm as well as fits a GMM to the time series data, the output from the Initialization Stage is quite sensitive to the starting conditions. In particular the HDRs used to initialize the GMMHMM can change with each run. The Iteration Stage is meant to mitigate some of the vagaries of the Initialization Stage by conducting a second search through the time series data. The time series is partitioned according to the predicted states for the time series according to the fitted GMMHMM from the previous stage and we search for missing HDRs in each partition. We present an overview of the processing performed by the iteration stage for each state (that is, partition) of a fitted GMMHMM.

We describe each of the stages in more detail below.

Algorithm

1. Data Partition Predicted to lie in State k. We find the predicted states for the observed time series based on the previously-fitted GMMHMM and extract all observations predicted to lie in State k.
2. Fit GMM to find HDRs. A Gaussian Mixture Model (GMM) is fitted to the time series partition. The HDRs are defined to be those Gaussian components that have covariances below some pre-determined threshold and where the number of observations predicted to lie in the Gaussian component exceeds some pre-determined frequency threshold. The means, covariances and weights of the HDRS are returned. As with the k-Means algorithm, the GMM is very sensitive to the initialization. Hence the HDRs found can vary dramatically with each run of the system. However, the partitioning of the data into (hopefully) more homogeneous subsets should make the HDRs found for each partition less sensitive to the starting conditions.
3. Match HDRs to Mixtures of GMMHMM. The HDRs were found for the subset of the time series predicted to lie in State k. We wish to use the HDRs to initialize the relevant mixture components of the k-th state of the GMMHMM. We rank the HDRs in descending order of its frequency and subsequently match each HDR to the mixture component that has the smallest Kullback-Leilbler distance to the HDR. Once a mixture component has been matched to an HDR, it is removed from the pool of mixture components that can be matched to the remaining HDRs.
4. Initialize states of GMMHMM using found HDRs. We initialize the matching mixtures of the matching states of the GMMHMM with the corresponding HDR's mean, covariance and weight.

3.3. Enforcement Stage

We present an overview of the Enforcement Stage:

Algorithm

1. Predict State-Mixture for each data point. We predict the most likely mixture that gave rise to each observation by finding the Mixture that maximizes where is the most likely state for the i-th observation according to the Viterbi algorithm and is the Gaussian component corresponding to the i-th mixture of the most likely state for the i-th observation . Hence, we assume that the most likely State-Mixture for each observation must be one of the mixtures from the most likely State for that same observation. Though this assumption may not always hold, it does greatly speed-up the calculations. Even still, this step is often the main computational bottleneck for the Enforcement Stage.
2. Get dominant State-Mixture & their indices. A dominant state-mixture is defined to be one that has at least 200 observations predicted to lie within it. We check if the dominant state-mixtures cumulatively account for at least 95% of observations. If not, we drop the frequency threshold that defines the dominant state-mixtures.
3. Choose appropriate width for Time Bins. We use the standard deviation of the time indices together with the sample size of the dominant state-mixtures to calculate the appropriate width for the time bins (see Scott, 1979) to use to represent the activity profile through time for each of the dominant State-Mixtures.
4. Calculate Mean & Covariance for Rare State-Mixtures. Any state-mixture not meeting the frequency threshold to be considered a "dominant" state-mixture is rather considered to be a rare state-mixture. We pool all such state-mixtures together and calculate the overall sample mean and covariance for such observations. These shall be used to initialize GMMHMM components that don't have a corresponding dominant state-mixture.
5. Calculate relative simultaneity between pairs of State-Mixtures. The activity vector for each of the dominant state-mixtures is calculated. This shows the proportion of the state-mixture's activation that occurs in each time bin. A cosine similarity matrix is subsequently calculated so that each matrix element compares the activity vectors of a pair of state-mixtures. This forms a measure of the relative simultaneity of the state-mixtures.
6. Cluster dominant State-Mixtures based on simultaneity. The cosine similarity matrix is used to calculate an adjacency matrix. Two dominant state-mixtures are connected if the cosine distance between them is less than some threshold and at least one of the two points has the other point as one of its two closest points. Spectral clustering is applied to the adjacency matrix to yield clusters of dominant state-mixtures that all have similar activity profiles over time. The goal is that each of such clusters should initialize a state in the GMMHMM. Hence, the number of clusters must be less than the number of states for the GMMHMM. If this is not the case, we increase the threshold for the cosine similarity thus making it easier for the state-mixture clusters to merge.
7. For each State-Mixture cluster, cluster the Mixtures based on Distribution Similarity. For each of the state-mixture clusters from the previous step, we calculate the Jensen-Shannon divergence between each pair of state-mixtures that lies in the cluster (provided the cluster size is at least four). If the cluster has less than four elements within it, we just assume each element (i.e. each state-mixture) forms its own cluster. The Jensen-Shannon divergence is closely related to the Kullback-Leibler divergence. However, unlike the KL divergence, the Jensen-Shannon divergence is a valid distance metric. We subsequently calculate an adjacency matrix where a pair of mixture components are connected if their Jensen-Shannon divergence is below some threshold and at least one of the two mixtures has the other as one of its three closest mixtures. The adjacency matrix is used to cluster the mixtures based on their distributional similarity. Thus, at the end of this step, we have a two-level nested clustering of the state-mixtures. The top level is derived in step 6 and clusters the state-mixture components if they have similar time-activity profiles (that is, the events from the state-mixture components tend to occur simultaneously). The next level goes through the mixtures in each of the top-level clusters and groups them together based on their distributional similarity. Each of these nested clusters will be used to initialize one of the mixtures in one of the GMMHMM's states. As such the size of the largest top-level cluster cannot be larger than the number of mixtures allowed for the GMMHMM. If this happens to be the case, we increase the Jensen-Shannon threshold, thus making it easier for the nested clusters to merge and thus reduce the size of the top-level cluster. Note that the clustering solution yielded by this stage allows us to adaptively reduce the complexity of the GMMHMM (i.e. the number of state and mixture components) if the dynamics of the time series are relatively simple. We make the number of states at most one more than the number of top-level clusters found and the number of mixtures at most one more than the size of the largest top-level cluster found. As such, the model complexity of the GMMHMM is set by an analysis of the underlying network traffic dynamics of the user in question.
8. Set the weights for Rare State-Mixtures. A separate state is set aside to account for behaviour outside the dominant state-mixtures. The weights for the mixture components of this "rare" state will be determined by the relative proportion of observations predicted to lie in the most common rare state-mixtures (i.e. the top rare state-mixtures where is the smaller of the number of rare state-mixtures & the allowed number of mixtures in the GMMHMM). If the number of rare state-mixtures is less than the number of mixtures allowed by the GMMHMM, we allocate 90% of the weight to the mixture components corresponding to the top rare state-mixtures (split according to their relative proportions) and the remaining 10% is split evenly among the remaining mixture components. Otherwise, all the weight is split amongst the rare state-mixtures.
9. Initialize the State-Mixture of GMMHMM using the Nested Clustering solution. As alluded to in step 7, the top-level of the clustering solution has the clusters split so that each cluster's elements tend to occur simultaneously. Each of these top-level clusters correspond to a particular state of the GMMHMM. The second level clustering solution takes the state-mixtures that tend to occur simultaneously and clusters them according to their distributional similarity. Each of these nested clusters thus correspond to one of the mixtures in one of the states in the GMMHMM. We choose to initialize the GMMHMM's mean and covariances using a weighted average of the state-mixtures in the corresponding nested cluster. The weight is proportional to the total number of times that the GMMHMM predicts that an observation will lie in one of these corresponding state-mixtures. For all GMMHMM components that don't have a corresponding dominant state-mixture, we initialize them using the behaviour of the pooled, rare state-mixtures.

3.4. Aggregation Stage

We present an overview of this Stage:

Algorithm

1. Choose appropriate width for Time Bins. We use the standard deviation of the time indices for each state together with the sample size of the states to calculate the appropriate width for the time bins (see Scott, 1979) to use to represent the activity profile through time for each of the states.
2. Calculate relative simultaneity between pairs of States. The activity vector for each of the states is calculated. This shows the proportion of the state's activation that occurs in each time bin. A cosine similarity matrix is subsequently calculated so that each matrix element compares the activity vectors of a pair of states. This forms a measure of the relative simultaneity of the states.
3. Cluster States based on simultaneity. The cosine similarity matrix is used to calculate an adjacency matrix. Two states are connected if the cosine distance between them is less than some threshold and at least one of the two points has the other point as one of its two closest points. Spectral clustering is applied to the adjacency matrix to yield clusters of states that all have similar activity profiles over time. The goal is that each of such clusters should initialize a state in the GMMHMM and most (if not all) of the clusters should be of size one with the larger clusters composed of states that occur simultaneously and hence should be combined.
4. Initialize the State-Mixtures of GMMHMM using the Clustering solution. Each of the clusters will correspond to a particular state of the GMMHMM. For those clusters of size one, we set a state in the GMMHMM with the corresponding weights, means and covariances of the single element (i.e. state) in the matching cluster. For states of the GMMHMM corresponding to clusters with more than one element, we order the mixture components of the member states in descending order of likelihood and use the top mixture components across these member states to initial the mixture components of the single, aggregated state.

3.5. Regularized Fit

At the end of each of the above stages, the initialization values are given for the GMMHMM's means, covariances and weights. We however need to do some extra work before fitting the GMMHMM.

1. Initializing the Transition Matrix. The transition matrix is also initialized. In keeping with the first principle from section 2, we initialize the transition matrix to be strongly diagonal (i.e. have 90% recurrence probability for each of the states) whilst splitting the remaining 10% probability evenly amongst the remaining states.
2. Regularize the Covariance Matrix. We take the covariance matrices suggested by the preceding stage (i.e. either the Initialization, Iteration, Enforcement or Aggregation stage) and ensure that the variance is never less than some minimum threshold. We also add a small number along the diagonal of the covariance matrix as this is known to make analysis on the covariance matrix more numerically stable. The correlations implied by the covariance matrix is also not allowed to be too strong (i.e. more than 0.9 or less than -0.9).
3. Fit the GMMHMM. We take the resulting means, covariances, weights and transition matrix to initialize the GMMHMM. Subsequently, the GMMHMM is fit to the data.
4. Regularize the Covariance Matrix. As the fitting procedure can result in a numerically unstable covariance matrix, we again choose to regularize the estimated covariance matrices using the same procedure as in step 2.
5. Regularize the Starting Probabilities. It is common for the fitted GMMHMM to yield a highly concentrated starting probability (i.e. one state which has probability almost one and all other states with probability almost zero). We choose to regularize the starting probabilities to reflect the overall relative occurrences of the various states.
6. Score the GMMHMM model. We can take the fitted, regularized GMMHMM and calculate the log-likelihood it yields on the observed data. These scores are used to decide whether the Iteration or Enforcement stages have succeeded. That is, these stages are only judged to have succeeded if they have improved on the log-likelihood score from the previous stage.

References

• Chandola, V., Banerjee, A., Kumar, V. 2009. Anomaly Detection: a Survey. ACM Computing Surveys.
• Dai, H., Dai, B., Zhang, Y., Li, S., Song, L. 2017. Recurrent Hidden Semi-Markov Model. ICLR.
• Scott. D.W. 1979. On Optimal and Data-Based Histograms. Biometrika 66, 605-610.