We present a real-time multivariate anomaly detection algorithm for data streams based on the Probabilistic Exponentially Weighted Moving Average (PEWMA). Our formulation is resilient to (abrupt transient, abrupt distributional, and gradual distributional) shifts in the data. This novel anomaly detection routines utilize an incremental online algorithm to handle streams. Furthermore, our proposed anomaly detection algorithm works in an unsupervised manner, eliminating the need for labelled examples. Our algorithm performs well and is resilient to concept drift.
Anomaly detection is the task of classifying patterns that depict abnormal behaviour. Outliers can arise due to (human/equipment) errors, faulty systems, and others.
Anomaly detection is well-suited to unbalanced data scenarios, where the ideal scenario is to predict minority class behaviour. There are numerous applications for detecting loan default, fraud detection, and network intrusion detection. An anomaly detection algorithm can work in Unsupervised, Supervised, or hybrid modes.
There are diferent types of anomalies described as follows. • Static: These algorithms work with static data.
Every item is loaded into memory at the same time in order to perform computations. • Online: These algorithms work in real-time data streams. Items are incrementally loaded into memory and processed in chunks. • Static + Online: The model may operate in two stages, as initial parameters get estimated in the static setting. The parameters are incrementally updated as more data arrives. Our work is of this type.
PEWMA was introduced in the work [1] for online anomaly detection on univariate time series. Drawing inspiration from their work, we have provided extensions to support real-time anomaly detection of a multivariate data stream.
1 2. Background • Point Anomaly: the algorithm identifies a single instance as an anomaly concerning the entire data set. • Contextual Anomaly: a data instance can be anomalous based on the context (attributes and position in the stream) and proximity of the chosen anomaly. This anomaly type is ideal for multivariate data, e.g. in the snapshot reading of a machine, an attribute of a single data point may seem abnormal but can be normal behaviour based on consideration of the entire data. • Collective Anomaly: the algorithm decides a set of data points that are anomalies as a group, but individually these data points exhibit normal behaviours.
EWMA limitations motivated the discovery of the Probabilistic Exponentially Weighted Moving Average (PEWMA). PEWMA [1] improves on EMWA by adding a parameter that includes the probability of the data in the model, as shown in Equation 2. PEWMA works for every shift including abrupt transient shift, abrupt distributional shift, and gradual distributional shift respectively.
PEWMA and EWMA use a damped window [8].
← (1 − ) end if //moving average 1 ← 1 + (1 − ) 2 ← 2 + (1 − ) 2 //incremental mean +̂ 1 ← 1 //incremental standard deviation +̂1 ← √ 2 − 12
The parameters of the anomaly detection algorithm +1 = (1 − ) −1 + (1 − (1 − )) (2) consist of the current data, the mean of the data, ̂ is the mean of the data, ̂ the current standard deviation, 3. Method the probability density function, +̂ 1 the mean of the next data (incremental aggregate), +̂1 the next standard Our formulation provides an implementation of the on- deviation (incremental aggregates), the data size, and line covariance matrix in Subsection 3.2, alongside an on- a point in . Initialize the process by setting the initial = ∗
= ∗ 2. Perform Cholesky factorization on the initial covariance matrix (positive-definite), as shown in
Equation 4. 3. The updated covariance in the presence of new data is equivalent to the weighted average of the past covariance without the updated data and the covariance matrix of the transformed input, as shown in Equation 5.
+1 = + ∗ Where = ∗ and ∈ is understood in our implementation as the current data. and are positive scalar values. 4. Increment the Cholesky factor of the covariance matrix as shown in Equation 6. +1 = √ + ⎛ ⎝ √ ‖ ‖2 ⎜√ ⎜ 1 +
2 ‖ ‖ ⎞ ⎟ ⎠ − 1⎟ ∗ (6) 5. There are dificulties with setting /ℎ
and respectively. + = 1
as an explicit form of exponential moving average coeficients. The author chose to set the values of , using the data stream statistics, as shown in Equation 7.
The parameters are set as = 2, = 1 − and is the size of the original data in the static
2 settings.
Where = √1 −
and +1 = + ⎛ ⎝ ⎜ 1 + ‖ ‖2 ⎜√ (3) (4) (5) = 2+6
.
Our work used the following parameters , , and as seen in Subsection 3.2. These anomaly thresholds are chosen based on the criteria that outliers are ≥ 3 times the standard deviation in normally distributed data.
1. Estimate covariance matrix for initial data, ∈ × .
Initial covariance matrix, where ∈ × , is the number of samples, is the number of dimensions as shown in Equation 3. (8) (9) (10) (11) (12) (13) 1. Estimate covariance matrix for initial data, ∈ × .
The initial covariance matrix, where ∈ × , is the number of samples, is the number of dimensions as shown in Equation 8.
= ∗ −1 = ( ∗ )−1 = ∗ Inverse the covariance matrix, −1 as shown in Equation 9. 2. Perform Cholesky factorization on the initial covariance matrix, as shown in Equation 10. 3. Increment the Cholesky factor of the covariance matrix
+−11 = ( + ∗ )−1 +−11 = −1( + ∗ )−1 Let us fix, ̂ = expression to
using Sherman−Morrison Formula reduces the ∗ . The resulting simplification 1 +−11 = ( −1 − −1 ̂ −1 1 + ( ̂ −1 ) )
The probability density function utilizes ideas from hypothesis testing for deciding on a threshold to set the confidence level. This threshold is used for determining the acceptance and rejection regions of the Gaussian distribution curve.
variance matrix, +1
−1. 1. Use the covariance matrix, +1 and inverse co2. We increment the mean vector, as new data Equation 14. arrives. It is possible to simplify the Covariance matrix, , which captures a number of system dynamics. Let represent the current data count before updated data arrives. Also, ̂: is the most recent data, +1 : moving average as shown in +1 = ( ∗ ) + ̂ + 1 (14) 3. Set a threshold to determine the acceptance and rejection regions. Items in the acceptance region are considered normal behaviour as shown in Equation 15. () =
1 √(2 ) || exp (− 1 ( − ) −1( − ) ) 2 (15) Where is mean vector, is the covariance matrix, || is the determinant of matrix, ∈ is data vector, and is the dimension of respectively.
Furthermore, we have provided detailed experiments on the proposed algorithms in diferent realistic scenarios. However, we maintain the statistical framework provided by the work [1] with theoretical guarantees.
We have experimented to determine the usefulness of our algorithm by creating a simulation with 10000000 random vectors with dimensions of 15. The repeated trial shows that our algorithm is not sensitive to initialization seeds and matrix dimensions. This requirement was a deciding factor in the choice of the evaluation metric, as shown in Equation 16. We have provided more information on the metric in Section 5.
Our experiment will check the efect of varying the size of the initial static window versus the update window. This is shown in Subsection 4.1 and Subsection 4.2.
We evaluated the trade-of between the static window and the update window. The experiment setup is as follows: • Split the data into 5 segments train on 1st segment(static), update covariance on 2nd (online), compare with static covariance, and calculate the error. • Train on 1, 2 segments (static), update covariance on 3rd (online), compare with static covariance and calculate the error. • Train on 1, 2, 3 segments (static), update covariance on 4th (online), compare with static covariance and calculate the error. • Train on 1, 2, 3, 4 segments (static), update covariance on 5th (online), compare with static covariance and calculate the error. with static covariance and calculate error on segments (5)
Our random matrices get flattened into vectors and used as input. The length of the flattened vector is applied as a normalization factor to make the loss metric agnostic of the matrix dimension. The loss function used in the evaluation is Absolute Average Deviation (AAD) because it gives a tighter bound on the error than mean squared error (MSE) or mean absolute deviation (MAD) as shown in Equation 16. We take the average of the residuals divided by the ground truth for every sample in our evaluation set. If the residual is close to zero, we contribute almost nothing to the measure. On the contrary, if the residual is high, we want to know the diference from the ground truth. 377–380. 25 (2013) 2283–2301. and is the length of the flattened matrix.
We can observe that building your model with more data in the init (static) phase leads to lower errors than having fewer data in the init phase and using more data for an update. The observation matches our intuition because when you operate online, you tend to use smaller storage space. However, there is still a performance tradeof compared to batch mode.
The error at the beginning of our training is significant in both charts. This insight shows that rather than performing the expensive operation of converting a covariance matrix to have positive definiteness, it is better to use random matrices that are positive definite. More data would help us achieve convergence as more data arrives.
The success of the experiments has given us confidence that our multivariate anomaly detection algorithm would have similar characteristics to the univariate case described in the work [1].
There is no generic anomaly detection that works for every task. The underlying assumption in this work is that the features in use capture relevant information about the underlying dynamics of the system. Our proposed implementation is an anomaly detection algorithm for handling multivariate streams, even with challenging shifts. In future work, we will extend support for nonstationary distributions in multivariate data streams. (16)