Sequential Outlier Detection based on Incremental Decision Trees

TL;DR

This paper presents an online outlier detection algorithm that uses incremental decision trees to model normal data and adaptively identify anomalies in a data stream, improving detection accuracy over existing methods.

Contribution

The paper introduces a novel incremental decision tree-based multi-modal density estimation approach for real-time outlier detection with adaptive thresholding.

Findings

Achieves performance comparable to the optimal fixed threshold in hindsight.

Demonstrates significant improvements over state-of-the-art methods.

Effectively models complex data distributions with growing decision trees.

Abstract

We introduce an online outlier detection algorithm to detect outliers in a sequentially observed data stream. For this purpose, we use a two-stage filtering and hedging approach. In the first stage, we construct a multi-modal probability density function to model the normal samples. In the second stage, given a new observation, we label it as an anomaly if the value of aforementioned density function is below a specified threshold at the newly observed point. In order to construct our multi-modal density function, we use an incremental decision tree to construct a set of subspaces of the observation space. We train a single component density function of the exponential family using the observations, which fall inside each subspace represented on the tree. These single component density functions are then adaptively combined to produce our multi-modal density function, which is shown to achieve the performance of the best convex combination of the density functions defined on the subspaces. As we observe more samples, our tree grows and produces more subspaces. As a result, our modeling power increases in time, while mitigating overfitting issues. In order to choose our threshold level to label the observations, we use an adaptive thresholding scheme. We show that our adaptive threshold level achieves the performance of the optimal pre-fixed threshold level, which knows the observation labels in hindsight. Our algorithm provides significant performance improvements over the state of the art in our wide set of experiments involving both synthetic as well as real data.

Figures

Tables

Table 1. TABLE I : “Log-Loss”, “AUC” and “Running Time” of the algorithms over the datasets described in Section VI-A . The AUC and Running Time values are in the format of “mean value ± plus-or-minus \pm standard deviation.

Algorithm	Log-Loss	Area Under Curve	Running Time (ms)
ITAN	$2.174$	$0.8281\pm 0.1383$	$313.02\pm 5.82$
wGMM	$3.056$	$0.8394\pm 0.0210$	$2275.67\pm 123.23$
wKDE	$3.022$	$0.5678\pm 0.0228$	$273.16\pm 6.92$
ML	$3.087$	$0.2700\pm 0.0943$	$38.56\pm 1.33$

Table 2. TABLE II : “Log-Loss”, “AUC” and “Running Time” of the algorithms over the datasets described in Section VI-B . The AUC and Running Time values are in the format of “mean value ± plus-or-minus \pm standard deviation.

Algorithm	Log-Loss	Area Under Curve	Running Time (ms)
ITAN	$0.833$	$0.7962\pm 0.0757$	$315.85\pm 13.67$
wGMM	$1.303$	$0.7381\pm 0.0285$	$7167.78\pm 140.99$
wKDE	$1.455$	$0.6863\pm 0.0321$	$283.21\pm 1.10$
ML	$1.337$	$0.6859\pm 0.0323$	$34.82\pm 0.43$

Table 3. TABLE III : “Log-Loss”, “AUC” and “Running Time” of the algorithms over the Vehicle Silhouettes dataset. The running time values are in the format of “mean value ± plus-or-minus \pm standard deviation”.

Algorithm	Log-Loss	Area Under Curve	Running Time (ms)
ITAN	$3.507$	$0.7483$	$466.33\pm 13.45$
wGMM	$6.346$	$0.5682$	$4646.87\pm 140.14$
wKDE	$7.684$	$0.4797$	$276.58\pm 6.04$
ML	$4.702$	$0.6806$	$34.71\pm 4.41$

Equations

$$\hat{p}_{t}(\cdot)=f(\mbox{$\mathbf{x}$}_{1},\mbox{$\mathbf{x}$}_{2},...,\mbox{$\mathbf{x}$}_{t-1},d_{1},d_{2},...,d_{t-1}).$$

$$l_{P}(\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}))=-\log(\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}))).$$

$$r_{P,t}(\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),p(\mbox{$\mathbf{x}$}_{t}))=-\log(\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}))+\log(p(\mbox{$\mathbf{x}$}_{t})),$$

$$R_{P,T}=\sum_{t=1}^{T}r_{P,t}(\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),p(\mbox{$\mathbf{x}$}_{t})).$$

$$\hat{d}_{t}=\text{sign}(\tau_{t}-\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t})),$$

$$l_{A}(\tau_{t},\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t})=C_{d_{t}}\mathbf{1}_{\{\text{sign}(\tau_{t}-\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}))\neq d_{t}\}}.$$

$$r_{A,t}(\tau_{t},\tau)=l_{A}(\tau_{t},\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t})-l_{A}(\tau,\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t}).$$

$$R_{A,T}=\sum _{t=1}^T{r}_{A,t}({\tau }_t,\tau ).$$

$$f(\mbox{$\mathbf{x}$}_{t})=h(\mbox{$\mathbf{x}$}_{t})\exp{(\langle\eta,\mbox{$\mathbf{s}$}_{t}\rangle\ -A(\eta))},$$

$$f(\mbox{$\mathbf{x}$}_{t})=\exp{(\langle\eta,\mbox{$\mathbf{s}$}_{t}\rangle\ -A(\eta))}.$$

$$\hat{f}_{t}(\mbox{$\mathbf{x}$}_{t})=\exp{(\langle\hat{\eta}_{t},\mbox{$\mathbf{s}$}_{t}\rangle\ -A(\hat{\eta}_{t}))}.$$

$$\hat{\eta}_{t}=\operatorname*{\arg\!\min}_{\eta}{\sum_{\tau=1}^{t-1}l({\eta},\mbox{$\mathbf{x}$}_{\tau})},$$

$$l({\eta},\mbox{$\mathbf{x}$}_{\tau})=-\langle\eta,\mbox{$\mathbf{s}$}_{\tau}\rangle\ +A(\eta).$$

$$m_{{\hat{\eta }}_t}=\frac{{\sum }_{\tau =1}^{t-1}{s}_{\tau }}{t-1},$$

$$p(\mbox{$\mathbf{x}$}_{t})=\sum_{n=1}^{N}\alpha_{n}f_{n}(\mbox{$\mathbf{x}$}_{t}),$$

$${\tilde{p}}_t(\cdot )=\sum _{i=1}^N{\tilde{\alpha }}_{t,i}{\tilde{f}}_{t,i}(\cdot ),$$

$$\tilde{\alpha}_{t+1,i}=\tilde{\alpha}_{t,i}\exp\left(\theta\frac{\tilde{f}_{t,i}(\mbox{$\mathbf{x}$}_{t})}{\tilde{p}_{t}(\mbox{$\mathbf{x}$}_{t})}\right),$$

$$\sum_{t=1}^{T}l_{P}\left({p}^{*}_{t}(\mbox{$\mathbf{x}$}_{t})\right)\leq AT,$$

$$R_{P,T}({\tilde{p}}_t(\cdot ),p_t^{\ast }(\cdot ))\le \sqrt{2AT\ln N}+\frac{R^2\ln N}{2}.$$

$$D\le \ln N-H(\left({\alpha }_1^{\ast },...,{\alpha }_N^{\ast }\right)),$$

$$\theta =\frac{2\sqrt{\ln N}}{R\sqrt{2AT}+R^2\sqrt{\ln N}},$$

$$\hat{d}_{t}=\begin{cases}+1,&\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t})<\tau_{t}\\ -1,&\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t})\geq\tau_{t}.\end{cases}$$

$$\tilde{l}(\tau_{t},\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t})=C_{d_{t}}\log(\exp((\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t})-\tau_{t})d_{t})+1).$$

$$\tilde{R}_{T_{f}}=\sum_{t\in T_{f}}\tilde{l}(\tau_{t},\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t})-\min\limits_{\tau\in\mathbf{G}}\sum_{t\in T_{f}}\tilde{l}(\tau,\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t}),$$

$$\tau_{t+1}=\begin{cases}\mathbb{P}_{\mathbf{G}}\left(\tau_{t}-\alpha_{t}\nabla_{\tau}\tilde{l}(\tau_{t},\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t})\right),&\text{if $d_{t}$ is known}\\ \tau_{t},&\text{otherwise},\end{cases}$$

$${\mathbb{P}}_{\mathbf{G}}(a)=\underset{b\in \mathbf{G}}{\mathrm{argmin}}\Vert b-a\Vert .$$

$${\alpha }_t=\frac{(1+\exp ({\mathcal{D}}_{\mathbf{G}}){)}^2}{t{C}_{\min }\exp ({\mathcal{D}}_{\mathbf{G}})},$$

$${\tilde{R}}_T\le \frac{\exp ({\mathcal{D}}_{\mathbf{G}})C_{\max }^2}{2C_{\min }}(1+\log T),$$

$$\frac{\partial\tilde{l}(\tau_{t},\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t})}{\partial\tau_{t}}=\frac{-d_{t}C_{d_{t}}}{1+\exp((\tau_{t}-\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}))d_{t})},$$

$$\frac{\partial^{2}\tilde{l}(\tau_{t},\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}),d_{t})}{\partial\tau^{2}_{t}}=\frac{C_{d_{t}}\exp((\tau_{t}-\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}))d_{t})}{(1+\exp((\tau_{t}-\hat{p}_{t}(\mbox{$\mathbf{x}$}_{t}))d_{t}))^{2}}.$$

Full text

EDICS Category: MLR-SLER, MLR-APPL, MLR-LEAR.

I Introduction

I-A Preliminaries

We study sequential outlier or anomaly detection [1], which has been extensively studied due to its applications in a wide set of problems from network anomaly detection [2, 3, 4] and fraud detection [5] to medical anomaly detection [6] and industrial damage detection [7]. In the sequential outlier detection problem, at each round $t$, we observe a sample vector \mbox{\mathbf{x}}_{t}\in\mathbb{X} and label it as “normal” or “anomalous” based on the previously observed sample vectors, i.e., \mbox{\mathbf{x}}_{t-1},...,\mbox{\mathbf{x}}_{1}, and their possibly revealed true labels. After we declare our decision, we may or may not observe the true label of \mbox{\mathbf{x}}_{t}. The objective is to minimize the number of mislabeled samples.

For this purpose, we use a two-stage “filtering” and “hedging” method [8]. In the “filtering” stage, we build in an online manner “a model” for “normal” samples based on the information gained from the previous rounds. Then, in the “hedging” stage, we decide on the label of the new sample based on its conformity to our model of normal samples. A common approach in constructing the aforementioned model is to assume that the normal data is generated from an independent and identically distributed (i.i.d.) random sequence [1]. Hence, in the first stage of our algorithm, we model the normal samples using a probability density function, which can also be considered as a scoring function [8]. However, note that the true underlying model of the normal samples can be arbitrary in general (or may not even exist) [1]. Therefore, we approach the problem in a competitive algorithm framework [9]. In this framework, we define a class of models called “competition class” and aim to achieve the performance of the best model in this class. Selecting a rich class of powerful models as the competition class enables us to perform well in a wide set of scenarios [9]. Hence, as detailed later, we choose a strong set of probability functions to compete against and seek to sequentially learn the best density function which fits to the normal data. Hence, while refraining from making any statistical assumptions on the underlying model of the samples, we guarantee that our performance is (at least) as well as the best density function in our competition class.

We emphasize that there exist nonparametric algorithms for density estimation [10], the parametric approaches have recently gained more interest due to their faster convergence [11, 12]. However, the parametric approaches fail if the assumed model is not capable of modeling the true underlying distribution [9]. In this context, exponential-family distributions [13] have attracted significant attention, since they cover a wide set of parametric distributions [8], and successfully approximate a wide range of nonparametric probabilistic models as well[14]. However, single component density functions are usually inadequate to model the data in highly challenging real life applications [15]. In this paper, in order to effectively model multi-modal distributions, we partition the space of samples into several subspaces using highly effective and efficient hierarchical structures, i.e., decision trees [16]. The observed samples, which fall inside each subspace are fed to a single component exponential-family density estimator. We adaptively combine all these estimators in a mixture-of-experts framework [17] to achieve the performance of their best convex combination.

We emphasize that the main challenge using a partitioning approach for multi-modal density estimation is to define a proper partition of the space of samples [15]. Here, instead of sticking to a pre-fixed partition, we use an incremental decision tree [16] approach to partition the space of samples in a nested structure. Using this method we avoid overtraining, while efficiently modeling complex distributions composed of a large number of components [16]. As the first time in the literature, in order to increase our modeling power with time, we apply a highly powerful incremental decision tree [16]. Using this incremental tree, whenever we believe that the samples inside a subspace cannot belong to a single component distribution, we split the subspace into two disjoint subspaces and start training two new single component density estimators on the recently emerged subspaces. Hence, our modeling power can potentially increase with no limit (and increase if needed), while mitigating the overfitting issues.

In order to decide on the label of a given sample, as widely used in the literature [8], we evaluate the value of our density function in the new data point \mbox{\mathbf{x}}_{t} and compare it against a threshold. If probability density is lower than the threshold, the sample is labeled as anomalous. While this is a shown to be an effective strategy for anomaly detection, setting the threshold is a notoriously difficult problem [8]. Hence, instead of committing to a fixed threshold level, we use an adaptive thresholding scheme and update the threshold level whenever we receive a feedback on the true label of the samples. We show that our thresholding scheme achieves the performance of the best fixed threshold level selected in hindsight.

I-B Prior Art and Comparisons

Various anomaly detection methods have been proposed in the literature that utilize Neural Networks [18], Support Vector Machines [19], Nearest Neighbors [20], clustering [21] and statistical methods including parametric [22] and nonparametric [23] density estimation. In the case when the normal data conform to a probability density function, the anomaly detection algorithms based on the parametric density estimation method are shown to provide superior performance [24]. For this reason, we adopt the parametric probability estimation based approach. In [8], authors have introduced an online algorithm to fit a single component density function of the exponential-family distributions to the stream of data. However, since the real life distributions are best described using multi-modal PDFs rather than single component density functions [25], we seek to fit multi-modal density functions to the observations. There are various multi-modal density estimation methods in the literature. In [15], authors propose a sequential algorithm to learn the multi-modal Gaussian distributions. However, as discussed in their paper, this algorithm provides satisfactory results only if a temporary coherency exists among subsequent observations. In [25], an online variant of the well-known Kernel Density Estimation (KDE) method is proposed. However, no performance guarantees are provided for any of the algorithms. In this paper, we provide a multi-modal density estimation method using an incremental tree with strong performance bounds, which are guaranteed to hold in an individual sequence manner through a regret formulation [8].

Decision trees are widely studied in various applications including coding [26], prediction [27, 28], regression [29] and classification [30]. These structures are shown to provide highly successful results due to their ability to refrain from overtraining while providing significant modeling power. In this paper, we adapt a novel notion of incremental decision trees [31] to the density estimation framework. Using this decision tree, we train a set of single-component density estimators with carefully chosen sets of data samples. We combine these single-component estimators in an ensemble learning [32] framework to approximate the underlying multi-modal density function and show that our algorithm achieves the performance of the best convex combination of the single component density estimators defined on the, possibly infinite depth, decision tree.

Adaptive thresholding schemes are widely used for anomaly detection algorithms based on density estimation [1]. While most of the algorithms in the literature do not provide guarantees for their anomaly detection performance, a surrogate regret bound of $O(\sqrt{t})$ is provided in [8]. However, since in real life applications the labels are revealed in a small portion of rounds [33], stronger performance guarantees are highly desirable. We provide an adaptive thresholding scheme with a surrogate regret bound of $O(\log{t})$. Hence, our algorithm steadily achieves the performance of the best threshold level chosen in hindsight.

I-C Contributions

Our main contributions are as follows:

• •

For the first time in the literature, we adapt the notion of incremental decision trees to the multi-modal density estimation framework. We use this tree, which can grow to an infinite depth, to partition the observations space into disjoint subspaces and train different density functions on each subspace. We adaptively combine these density functions to achieve the performance of the best multi-modal density function defined on the tree.

• •

We provide guaranteed performance bounds for our multi-modal density estimation algorithm. Due to our competitive algorithm framework, our performance bounds are guaranteed to hold in an individual sequence manner.

• •

Due to our individual sequence perspective, our algorithm can be used in unsupervised, semi-supervised and supervised settings.

• •

Our algorithm is truly sequential, such that no a priori information on the time horizon or the number of components in the underlying probability density function is required.

• •

We propose an adaptive thresholding scheme that achieves a regret bound of $O(\log{t})$ against the best fixed threshold level chosen in hindsight. This thresholding scheme improves the state-of-the-art $O(\sqrt{t})$ regret bound provided in [8].

• •

We demonstrate significant performance gains in comparison to the state-of-the-art algorithms through extensive set of experiments involving both synthetic and real data.

I-D Organization

In section II, we formally define the problem setting and our notation. Next, we explain our single-component density estimation methods in section III. In section IV, we introduce our decision tree and explain how we use it to incorporate the single-component density estimators and create our multi-modal density function. Then, we explain the anomaly detection step of our algorithm in section V, which completes our algorithm description. In section VI we demonstrate the performance of our algorithm against the state-of-the-art methods on both synthetic and real data. We finish with concluding remarks in section VII.

II Problem Description

In this paper, all vectors are column vectors and denoted by boldface lower case letters. For a $K$-element vector $\mathbf{u}$, ${u}_{i}$ represents the $i^{\text{th}}$ element and \lVert\mbox{\mathbf{u}}\rVert=\sqrt{\mbox{\mathbf{u}}^{T}\mbox{\mathbf{u}}} is the $l^{2}$-norm, where \mbox{\mathbf{u}}^{T} is the ordinary transpose. For two vectors of the same length $\mathbf{u}$ and $\mathbf{v}$, \langle\mbox{\mathbf{u}},\mbox{\mathbf{v}}\rangle=\mbox{\mathbf{u}}^{T}\mbox{\mathbf{v}} represents the inner product. We show the indicator function by $\mathbf{1}_{\{\text{condition}\}}$, which is equal to $1$ if the condition holds and [math] otherwise.

We study sequential outlier detection problem, where at each round $t\geq 1$, we observe a sample vector \mbox{\mathbf{x}}_{t}\in\mathbb{R}^{m} and seek to determine whether it is anomalous or not. We label the sample vector \mbox{\mathbf{x}}_{t} by $\hat{d}_{t}=-1$ for normal samples and $\hat{d}_{t}=1$ for anomalous ones, where $d_{t}$ corresponds to the true label which may or may not be revealed. In general, the cost of making an error on normal and anomalous data may not be the same. Therefore, we define ${C}_{d_{t}}$ as the cost of making an error while the true label is $d_{t}$. The objective is to minimize the accumulated cost in a series of rounds, i.e., $\sum_{t=1}^{T}{C}_{d_{t}}\mathbf{1}_{\{\hat{d}_{t}\neq d_{t}\}}$.

In our two step approach, we first introduce an algorithm for probability density estimation, which learns a multi-modal density function that fits “best” to the observations. This density function can be seen as a scoring function determining the normality of samples. Due to the online setting of our problem, at each round $t$, our density function estimate, denoted by $\hat{p}_{t}(\cdot)$, is a function of previously observed samples and their possibly revealed labels, i.e.,

[TABLE]

Note that in general, even if the samples are not generated from a density function, e.g., deterministic framework [34], our estimate $\hat{p}_{t}(\cdot)$ can be seen as a scoring function determining the normality of the samples. As widely used in the literature [35], we measure the accuracy of our density function estimate $\hat{p}_{t}$ by the log-loss function

[TABLE]

In order to refrain from any statistical assumptions on the normal data, we work in a competitive framework [9]. In this framework we seek to achieve the performance of the best model in a class of models called the competition class. We use the notion of “regret” as our performance measure in both density estimation and anomaly detection steps. The regret of a density estimator producing the density function $\hat{p}_{t}(\cdot)$ against a density function $p(\cdot)$ at round $t$ is defined as

[TABLE]

where selection of $p(\cdot)$ will be clarified later. We denote the accumulated density estimation regret up to time $T$ by

[TABLE]

In order to produce our decision on the label of observations being “normal” or “anomalous”, at each round $t$, we observe the new sample \mbox{\mathbf{x}}_{t} and declare our decision by thresholding \hat{p}_{t}(\mbox{\mathbf{x}}_{t}) as

[TABLE]

where $\tau_{t}$ is the threshold level. After declaring our decision, we may or may not observe the true label $d_{t}$ as a feedback. We use this information to optimize $\tau_{t}$ whenever we observe the correct decision $d_{t}$. We define the loss of thresholding \hat{p}_{t}(\mbox{\mathbf{x}}_{t}) by $\tau_{t}$ as

[TABLE]

We define the regret of choosing the threshold value $\tau_{t}$ against a specific threshold $\tau$ (which can even be the unknown “best” threshold that minimizes the cumulative error) at round $t$ by

[TABLE]

We denote the accumulated anomaly detection regret up to time $T$ by

[TABLE]

We emphasize that the main challenge in “two-step” approaches for anomaly detection is to construct a density function $\hat{p}_{t}(\cdot)$, which powerfully models the observations distribution. For this purpose, in section III, we first introduce an algorithm, which achieves the performance of a wide range of single component density functions. Based on this algorithm, in section IV, we use a nested tree structure to construct a multi-modal density estimation algorithm. In section V, we introduce our adaptive thresholding scheme, which will be used on the top of the density estimator described in section IV to form our complete anomaly detection algorithm.

III Single Component Density Estimation

In this section we introduce an algorithm, which sequentially achieves the performance of the best single component distribution in the exponential family of distributions [13]. At each round $t$, we observe a sample vector \mbox{\mathbf{x}}_{t}\in\mathbb{R}^{m}, drawn from an exponential-family distribution

[TABLE]

where

• •

$\eta\in\mathbf{F}$ is the unknown “natural parameter” of the exponential-family distribution. Here, $\mathbf{F}\subset\mathbb{R}^{d}$ is a bounded convex set.

• •

h(\mbox{\mathbf{x}}_{t}) is the “base measure function” of the exponential-family distribution.

• •

$A(\eta)$ is the “log-partition function” of the distribution.

• •

\mbox{\mathbf{s}}_{t}\in\mathbb{R}^{d} is the “sufficient statistics vector” of \mbox{\mathbf{x}}_{t}. Given the type of the exponential-family distribution, e.g., Gaussian, Bernoulli, Gamma, etc., \mbox{\mathbf{s}}_{t} is calculated as a function of \mbox{\mathbf{x}}_{t}, i.e., \mbox{\mathbf{s}}_{t}=T(\mbox{\mathbf{x}}_{t}).

With an abuse of notation, we put the “base measure function” h(\mbox{\mathbf{x}}_{t}) inside the exponential part by setting \mbox{\mathbf{s}}_{t}=[\mbox{\mathbf{s}}_{t};\log(h(\mbox{\mathbf{x}}_{t}))] and $\eta=[\eta;1]$. Hence, from now on, we write

[TABLE]

At each round $t$, we estimate the natural parameter $\eta$ based on the previously observed sample vectors, i.e., \{\mbox{\mathbf{x}}_{1},\mbox{\mathbf{x}}_{2},...,\mbox{\mathbf{x}}_{t-1}\}, and denote our estimate by $\hat{\eta_{t}}$. The density estimate at time $t$ is given by

[TABLE]

In order to produce our estimate $\hat{\eta}_{t}$, we seek to minimize the accumulated loss we would suffer following this $\hat{\eta}_{t}$ during all past rounds, i.e.,

[TABLE]

where

[TABLE]

This is a convex optimization problem. Finding the point in which the gradient is zero, it can be seen that it suffices to choose the $\hat{\eta}_{t}$ such that

[TABLE]

where $m_{\hat{\eta}_{t}}$ is the mean of \mbox{\mathbf{s}}_{t} when \mbox{\mathbf{x}}_{t} is distributed with the natural parameter $\hat{\eta}_{t}$.

Note that the memory demand of our single-component density estimator does not increase with time, as is suffices to keep the sample mean of the “sufficient statistic vectors”, i.e., $s_{\tau}$’s, in memory. The complete pseudo code of our single component density estimator is provided in Alg. 1.

IV Multimodal Density Estimation

In this section, we extend our basic density estimation algorithm to model the observation vectors using multi-modal density functions of the form

[TABLE]

where each $f_{n}(\cdot)$ is an exponential-family density function as in (9) and $(\alpha_{1},...,\alpha_{N})$ is a probability simplex, i.e., $\forall n:\alpha_{n}\geq 0$, $\sum_{n=1}^{N}\alpha_{n}=1$.

In order to construct such model, we split the space of sample vectors into several subspaces and run an independent copy of the Alg. 1 in each subspace. Each one of these density estimators observe only the sample vectors, which fall into their corresponding subspace. We adaptively combine the aforementioned single component density estimators to produce our multi-modal density function. In the following, in Section IV-A, we first suppose that a set of subspaces is given and explain how we combine the density estimators running over the subspaces. Then, in Section IV-B, we explain how we construct our set of subspaces using an incremental decision tree.

IV-A Mixture of Single Component Density Estimators

Let $\mathcal{S}=\{S_{1},...,S_{N}\}$ be a given set of $N$ subspaces of the observation space. For instance, in Fig. 1 a set of $11$ subspaces in $\mathbb{R}^{2}$ is shown. We run $N$ independent copies of the Alg. 1 in these subspaces and denote the estimated density function corresponding to $S_{i}$ at round $t$ by $\tilde{f}_{t,i}(\cdot)$. We adaptively combine $\tilde{f}_{t,i}(\cdot)$, $i=1,...N$, in a mixture-of-experts setting using the well known Exponentiated Gradient (EG) algorithm [36]. At each round $t$, we declare our multi-modal density estimation as

[TABLE]

where $\tilde{\alpha}_{1,i}$’s are initialized to be $1/N$ for $i=1,...,N$. After observing \mbox{\mathbf{x}}_{t}, we suffer the loss l(\tilde{p}_{t},\mbox{\mathbf{x}}_{t})=-\log(\tilde{p}_{t}(\mbox{\mathbf{x}}_{t})) and update the mixture coefficients as

[TABLE]

where $\theta$ is the learning rate parameter. The following proposition shows that in a $T$ rounds trial, we achieve a regret bound of $O(\sqrt{T})$ against the multi-modal density estimator with the best $\tilde{\alpha}$ variables (in the log-loss sense), i.e., the best convex combination of our single component density functions.

Theorem 1**.**

For a $T$ round trial, let $R$ be a bound such that $\max_{t,n}\{\tilde{f}_{t,n}\}\leq R,$ for all $t,n$. Let ${p}^{*}_{t}(\cdot)=\sum_{n=1}^{N}{{\alpha}^{*}_{n}}\tilde{f}_{t,n}(\cdot)$ be the optimal (in the accumulated log-loss sense) convex combination of $\tilde{f}_{t,i}$’s with fixed coefficients $\left({\alpha}^{*}_{1},...,{\alpha}^{*}_{N}\right)$ selected in hindsight. If the accumulated log-loss of {p}^{*}_{t}(\mbox{\mathbf{x}}_{t}) is upper bounded as

[TABLE]

we achieve a regret bound as

[TABLE]

Proof.

Denoting the relative entropy distance [37] between the best probability simplex $\left({\alpha}^{*}_{1},...,{\alpha}^{*}_{N}\right)$ and the initial point $\left(\tilde{\alpha}_{1,1},...,\tilde{\alpha}_{1,N}\right)$ by $D$, since $\tilde{\alpha}_{1,n}=1/N,\forall n=1,...,N$, we have

[TABLE]

where $H(\left(\alpha^{*}_{1},...,\alpha^{*}_{N}\right))$ is the entropy of the best probability simplex. Since the entropy is always positive, we have $D\leq\ln{N}$. Using Exponentiated Gradient [36] algorithm with the parameter

[TABLE]

we achieve the regret bound in (19). ∎

Remark 1**.**

We emphasize that one can use any arbitrary density estimator in the subspaces and achieve the performance of their best convex combination using the explained adaptive combination. However, since the exponential family distribution covers a wide set of parametric distributions and closely approximates a wide range of non-parametric real life distributions, we use the density estimator in Alg. 1.

As shown in the theorem, no matter how the set of subspaces $\mathcal{S}$ is constructed, our multi-modal density estimate in (16) is competitive against the best convex combination of the density functions defined over the subspaces in $\mathcal{S}$. However, the subspaces themselves play an important role in building a proper model for arbitrary multi-modal distributions. For instance, suppose that the true underlying model is a multi-modal PDF composed of several components, which are far away from each other. If we carefully construct subspaces, such that each subspace contains only the samples generated from one of the components (or these subspaces are included in $\mathcal{S}$), then the best convex combination of the subspaces will be a good model for the true underlying PDF. This scenario is further explained through an example in Section VI-A.

In the following subsection, we introduce a decision tree approach [16] to construct a growing set of proper subspaces and fit a model of the form (15) to the sample vectors. Using this tree, we start with a model with $N=1$ and increase $N$ as we observe more samples. Hence, while mitigating overfitting issues due to the $\ln{N}$ bound in (19), our modeling power increases with time.

IV-B Incremental Decision Tree

We introduce a decision tree to partition the space of sample vectors into several subspaces. Each node of this tree corresponds to a specific subspace of the observation space. The samples inside each subspace are used to train a single component PDF. These single component probability density functions are then combined to achieve the performance of their best convex combination.

As explained in Section IV-A, our adaptive combination of single component density functions will be competitive against their best convex combination, regardless of how we build the subspaces. However, in order to closely model arbitrary multi-modal density functions of the form (15), we seek to find subspaces that contain only the samples from one of the components. Clearly, this is not always straightforward (or may not be even possible), specially if the centroids of the component densities are close to each other. To this end, we introduce an incremental decision tree [16] which generates a set of subspaces so that as we observe more samples, our tree adaptively grows and produces more subspaces tuned to the underlying data. Hence, using its carefully produced subspaces, we are able to generate a multi-modal PDF that can closely model the normal data even for complex multi-modal densities, which are hard to learn with classical approaches. We next explain how we construct this incremental tree. We emphasize that we use binary trees as an example and our construction can be extended to multi branch trees in a straightforward manner.

We start building our binary tree with a single node corresponding to the whole space of the sample vectors. As an example, consider step 1 in Fig. 2. We say that this node is the $1^{\text{st}}$ node in level [math], and denote this node by a binary index of $(0,1)$, where the first element is the node’s level and the second element is the order of the node among its co-level nodes. We grow the tree by splitting the subspace corresponding to a specific node into two subspaces (corresponding to two new nodes) at rounds $t={\beta}^{k}$ for $k=1,2,...$. Hence, at each round $t$, the tree will have $\lceil{\log_{\beta}t}\rceil$ nodes. We emphasize that, as shown in Theorem 1, selecting the splitting times as the powers of $\beta$, we achieve a regret bound of $O(\sqrt{T\log{\log{T}}})$ against the best convex combination of the single component PDFs (see (19)). Moreover, this selection of splitting times leads to a logarithmic in time computational complexity. However, again, we note that our algorithm is generic so that the splitting times can be selected in any arbitrary manner.

To build subspaces (or sets), we use hyperplanes to avoid overfitting. In order to choose a proper splitting hyperplane, we run a sequential $2$-means algorithm [38] over all the nodes as detailed in Alg. 2. These $2$-means algorithms are also used to select the nodes to split as follows. At each splitting time, we split the node that has the maximum ratio of “distance between $2$ centroids” to “$2^{\{\text{level of the node}\}}$”, where “level of the node” is the number of splits required to build the node’s corresponding subspace as shown in Fig. 2. Note that as this ratio increases, it’s implied that the node does not include samples from a single component PDF, which makes it a good choice to split. This motivation is illustrated using a realistic example in Section VI-A. We split the nodes using the hyperplane, which is perpendicular to the line connecting the two centroids of the 2-means algorithm running over the node and splits this line in half. The splitted node keeps a portion of its $\tilde{\alpha}$ value for itself and splits the remaining among its children. This portion, which is a parameter of the algorithm is denoted by $\xi$. We emphasize that using the described procedure, each node may be splitted several times. Hence, if the splitting hyperplane is not proper due to lack of observations, the problem can be fixed later by splitting the node again with more accurate hyperplanes in the future rounds. As an example, consider Fig. 2. At the last step, node $(2,3)$ is splitted again with a slightly shifted splitting line. This is illustrated in more detail using an example in Section VI-A. The algorithm pseudo code is provided in Alg. 2.

Remark 2**.**

We use linear separation hyperplanes to avoid overtraining while the modeling power is attained by using an incremental tree. However, our method can directly used with different separation hyperplanes.

As detailed in Alg. 2, at each round $t$, the tree nodes declare their single component PDFs, i.e., $\tilde{f}_{t,i}(\cdot),i=1,...,N$. We combine these density functions using (16) to produce our multi-modal density estimate $\tilde{p}_{t}(\cdot)$. Then, the new sample vector \mbox{\mathbf{x}}_{t} is observed and we suffer our loss as (2). Subsequently, we update the combination variables, i.e., $\tilde{\alpha}_{t,i},i=1,...,N$, using (17). The centroids of the 2-means algorithms running over nodes are also updated as detailed in Alg. 2. Finally, the single component density estimates at the nodes are updated as detailed in Alg. 1. At the end of the round, if $t=\beta^{k}$, we update the tree structure and construct new nodes as explained in Section IV-B.

In the following section, we explain our adaptive thresholding scheme, which will be used on top of described multi-modal density estimator to form our two-step anomaly detection algorithm.

V Anomaly Detection Using Adaptive Thresholding

We construct an algorithm, which thresholds the estimated density function \hat{p}_{t}(\mbox{\mathbf{x}}_{t}) to label the sample vectors. To this end, we label the sample \mbox{\mathbf{x}}_{t} by comparing \hat{p}_{t}(\mbox{\mathbf{x}}_{t}) with a threshold $\tau_{t}$ as

[TABLE]

Suppose at some specific rounds $t\in T_{f}$, after we declared our decision $\hat{d}_{t}$ the true label $d_{t}$ is revealed. We seek to use this information to minimize the total regret defined in (8). However, since we observe the incurred loss only at rounds $t\in T_{f}$, we restrict ourselves to these rounds. Moreover, since the loss function used in (8) is based on the indicator function that is not differentiable, we substitute the loss function defined in (6) with the well known logistic loss function defined as

[TABLE]

Our aim is to achieve the performance of the best constant $\tau$ in a convex feasible set $\mathbf{G}$. To this end, we define our regret as

[TABLE]

We use the Online Gradient Descent algorithm [39] to produce our threshold level $\tau_{t}$. To this end, we choose ${\tau}_{1}\in\mathbf{G}$ arbitrarily. At each round $t$, after declaring our decision $\hat{d}_{t}$, we construct

[TABLE]

where $\alpha_{t}$ is the step size at time $t$ and $\mathbb{P}_{\mathbf{G}}(\cdot)$ is a projection function defined as

[TABLE]

The complete algorithm is provided in Alg. 3.

For the sake of notational simplicity, from now on, we assume that $d_{t}$ is revealed at all time steps. We emphasize that since the rounds with no feedback do not affect neither the threshold in (25), nor the regret in (24), we can simply ignore them in our analysis. The following theorem shows that using Alg. 3, we achieve a regret upper bound of $O(\log T)$, against the best fixed threshold level selected in hindsight.

Theorem 2**.**

Using Alg. 3 with step size

[TABLE]

our anomaly detection regret in (24) is upper bounded as

[TABLE]

where $\mathcal{D}_{\mathbf{G}}=\max\limits_{a,b\in\mathbf{G}}\|a-b\|$ is the diameter of the feasible set $\mathbf{G}$ including $\tau_{t}$ and \hat{p}(\mbox{\mathbf{x}}_{t}). $C_{\max}$ and $C_{\min}$ are the maximum and minimum of $\{C_{1},C_{-1}\}$, respectively.

Proof.

Considering the loss function in (23), we take the first derivatives of $\tilde{l}$ as

[TABLE]

and its second derivative as

[TABLE]

The first derivative can be bounded as

[TABLE]

Similarly, the second derivative is bounded as

[TABLE]

Using Online Gradient Descent [39], with step size given in (27) we achieve the regret upper bound in (28). ∎

VI Experiments

In this section, we demonstrate the performance of our algorithm in different scenarios involving both real and synthetic data. In the first experiment, we have created a synthetic scenario to illustrate how our algorithm works. In this scenario, we sequentially observe samples drawn from a $4$-component distribution, where the probability density function is a convex combination of $4$ multivariate Gaussian distributions. The samples generated from one of the components are considered anomalous. The objective is to detect these anomalous samples. In the second experiment, we have shown the superior performance of our algorithm with respect to the state-of-the-art methods on a synthetic dataset, where the underlying PDF cannot be modeled as a multi-modal Gaussian distribution. The third experiment shows the performance of the algorithms on a real multi-class dataset. In this experiment, the objective is to detect the samples belonging to one specific class, which are considered anomalous.

We compare the density estimation performance of our algorithm ITAN, against a set of state-of-the-art opponents composed of wGMM [40], wKDE [40], and ML algorithms. The wGMM [41] is an algorithm which uses a sliding window of the last $\log{t}$ normal samples to train a GMM using the well known Expectation-Maximization (EM) [41] method. The length of sliding window is set to $\log{t}$ in order to have a fair comparison against our algorithm in the sense of computational complexity. In favor of the wGMM algorithm, we provide to it the number of components that provides the best performance for that algorithm. The wKDE is the well-known KDE [40] algorithm that uses a sliding window of the last $\sqrt{t}$ normal samples to produce its estimate on the density function. The length of sliding window is $\sqrt{t}$ in favor of the wKDE algorithm to produce competitive results. The kernel bandwidth parameters are chosen based on Silverman’s rule [40]. Finally, ML algorithm is the basic Maximum Likelihood algorithm which fits the best single-component density function to the normal samples. We use our algorithm ITAN with the parameters $\beta=2$ and $\xi=0.8$ in the all three experiments. We emphasize that no optimization has been performed to tune these parameters to the datasets.

In order to compare the anomaly detection performance of the algorithms, we use the same thresholding scheme described in Alg. 3 for all algorithms. We use the ROC curve as our performance metric. Given a pair of false negative and false positive costs, denoted by $C_{1}$ and $C_{-1}$, respectively, each algorithm achieves a pair of True Positive Rate (TPR) and False Positive Rate (FPR), which determines a single point on its corresponding ROC curve. In order to plot the ROC curves, we have repeated the experiments $100$ times, where $C_{1}=1$ and $C_{2}$ is selected from the set of $\{{\frac{i}{100}}|i=0,1,...,99\}$. The ROC curves are plotted using the resulting $100$ samples. The Area Under Curve (AUC) of the algorithms are also calculated using these samples as another performance metric.

VI-A Synthetic Multi-modal Distribution

In the first experiment, we have created $10$ datasets of length $1000$ and compared the performance of the algorithms in both density estimation and anomaly detection tasks. Each sample is labeled as “normal” or “anomalous” with probabilities of $0.9$ and $0.1$, respectively. The normal samples are randomly generated using the density function

[TABLE]

while the anomalous samples are generated using

[TABLE]

Fig. 3 shows the samples in one of the datasets used in this experiment to provide a clear visualization.

In order to show how our algorithm learns, we illustrate how the tree splits the observation space, how the density estimations train their single component PDFs and how the combination of single component PDFs models the normal data in the experiment over one of the $10$ datasets. Fig. 4 shows five growth steps of the tree. In each subfigure, the observed samples are shown using black cross signs. The centroids of the 2-means algorithm running over the node that is going to split are shown using two blue and red points. The thicker green line is the new splitting line, while the thiner green lines show previous splitting lines. The splittings shown in this figure result in a tree structure that is shown in Fig. 2. Fig. 6 shows how the single component PDFs defined over nodes are combined to construct our multi-modal density function. In Fig. 6(a) the contour plot of the normal data distribution function is shown. Fig. 6(c) shows the structure of the tree at the end of the experiment, the contour plots of the single component PDFs learned over the nodes, and their coefficient in the convex combination which yields the final multi-modal density function. The contour plot of this final multi-modal PDF is shown in Fig. 6(b). As shown in these figures, the three components of the underlying PDF are almost captured by the three nodes, which are generated in the second level of our tree.

In order to compare the density estimation performance of the algorithms, their averaged loss per round defined by \text{Loss}(t)=\sum_{\tau=1}^{t}l_{P}(\hat{p}_{\tau}(\mbox{\mathbf{x}}_{\tau}))/t, are shown in Fig. 8(a). The loss of all algorithms on the rounds with anomalous observations are considered as [math] in these plots. The anomaly performance of the algorithms are compared in Fig. 8(d). This figure shows the ROC curves of the algorithms averaged over 10 datasets. The averaged log-loss performance, AUC results and running time of the algorithms are provided in Table. I. All results are obtained using a Intel(R) Core(TM) i5-4570 CPU with 3.20 GHz clock rate and 8 GB RAM.

As shown in Figs. 8(a), our algorithm achieves a significantly superior performance for the density estimation task. This superior performance was expected because in the dataset used for this experiment the components are far from each other. Hence, our tree can generate proper subspaces, which contain only the samples from one of the components of the underlying PDF, as shown in Fig. 6. For the anomaly detection task, as shown in Fig. 8(d), our algorithm and wGMM provide close performance, where ITAN performs better in low FPRs and wGMM provides superior performance in high FPRs. However, as shown in Table I, we achieve this performance with a significantly lower computational complexity. Comparing Fig. 8(a) and Fig. 8(d) shows that while satisfactory log-loss performance is required for successful anomaly detection, it is not sufficient in general. For instance, while ML algorithm performs as well as wKDE and wGMM in the log-loss sense, its anomaly detection performance is much worse than the others. In fact, labeling the samples exactly opposite of the suggestions of the ML algorithm provides way better anomaly detection performance. This is because of the weakness of the model assumed by the ML algorithm. This weak performance of the ML algorithm was expected due to the underlying PDF of the normal and anomalous data. It can be also seen from Fig. 3. If we fit a single component Gaussian PDF to the normal samples shown in blue, roughly speaking, the anomalous samples shown in red will get the highest normality score when evaluated using our PDF.

In the next experiment, we compare the algorithms in a scenario, where the data cannot be modeled as a convex combination of Gaussian density functions.

VI-B Synthetic Arbitrary Distribution

In this experiment, we have created $10$ datasets of length $1000$. In order to generate each sample, first its label is randomly determined to be “normal” or “anomalous” with probabilities of $0.9$ and $0.1$, respectively. The normal samples are $2$-dimensional vectors \mbox{\mathbf{x}}_{t}=[x_{t,1},x_{t,2}]^{T} generated using the following distribution:

[TABLE]

where $\mathcal{U}(a,b)$ is the uniform distribution between $a$ and $b$. The anomalous samples are generated using the following distribution:

[TABLE]

Fig. 7 shows the samples in one of the datasets used in this experiment to provide a clear visualization.

Fig. 8(b) shows the averaged accumulated loss of the algorithms averaged over 10 data sets. As shown in the figure, our algorithm outperforms the competitors for the density estimation task. This superior performance is due to the growing in time modeling power of out algorithm. The ROC curves of the algorithms for the anomaly detection task are shown in Fig. 8(e). This figure shows that our algorithm provides superior anomaly detection performance as well. This superior performance is due to the better approximation of the underlying PDF made by our algorithm. The averaged log-loss performance, AUC results and running time of the algorithms in this experiment are summarized in Table. II.

VI-C Real multi-class dataset

In this experiment, we use Vehicle Silhouettes [42] dataset. This dataset contains $846$ samples. Each sample includes a $18$-dimensional feature vector extracted from an image of a vehicle. The labels are the vehicle class among $4$ possible classes of “Opel”, “Saab”, “Bus” and “Van”. Our objective in this experiment is to detect the vehicles with “Van” labels as our anomalies. Fig. 8(c), shows the density estimation loss of the opponents, based on the rounds in which they have observed “normal” samples. Fig. 8(f) shows the ROC curves of the algorithms. As shown in the figures, our algorithm achieves a significantly superior performance in both density estimation and anomaly detection tasks over this dataset. The AUC results and running time of the algorithms are summarized in Table. 2.

Fig. 8(c) shows that the performance of wGMM highly depends on the stationarity of normal samples stream. The intrinsic abrupt change of the underlying model at around round $250$ has caused a heavy degradation in its density estimation performance. However, our algorithm shows a robust log-loss performance even in the case of non-stationarity. Fig. 8(f) shows that our algorithm achieves the best anomaly detection performance among the competitors. Note that ML algorithm outperforms both wGMM and wKDE algorithms in both density estimation and anomaly detection tasks. This is because wGMM and wKDE suffer from overtraining due to the high dimensionality of the sample vectors and short time horizon of the experiment. However, due to the growing tree structure used in our algorithm, we significantly outperform the ML algorithm and provide highly superior and more robust performance compared to the all other three algorithms.

VII Concluding Remarks

We studied the sequential outlier detection problem and introduced a highly efficient algorithm to detect outliers or anomalous samples in a series of observations. We use a two-stage method, where we learn a PDF that best describes the normal samples, and decide on the label of the new observations based on their conformity to our model of normal samples. Our algorithm uses an incremental decision tree to split the observation space into subspaces whose number grow in time. A single component PDF is trained using the samples inside each subspace. These PDFs are adaptively combined to form our multi-modal density function. Using the aforementioned incremental decision tree, while avoiding overtraining issues, our modeling power increases as we observe more samples. We threshold our density function to decide on the label of new observations using an adaptive thresholding scheme. We prove performance upper bounds for both density estimation and thresholding stages of our algorithm. Due to our competitive algorithm framework, we refrain from any statistical assumptions on the underlying normal data and our performance bounds are guaranteed to hold in an individual sequence manner. Through extensive set of experiments involving synthetic and real datasets, we demonstrate the significant performance gains of our algorithm compared to the state-of-the-art methods.