A Novel Geometric Approach for Outlier Recognition in High Dimension

Besides the aforementioned existing results, many other methods for outlier recognition/anomaly detection were developed previously and the readers can refer to several excellent surveys [22, 8, 18].

In computational geometry, a core-set [1] is a small set of points that approximate the shape of a much larger point set, and thus can be used to significantly reduce the time complexities for many optimization problems (please refer to a recent survey [28]). In particular, a core-set can be applied to efficiently compute an approximate MEB for a set of points in high-dimensional space [5, 24]. Moreover, [4] showed that it is possible to find a core-set of size $\lceil 2/\epsilon\rceil$ that yields a $(1+\epsilon)$-approximate MEB, with an important advantage that the size is independent of the original size and dimensionality of the dataset. In fact, the algorithm for computing the core-set of MEB is a Frank-Wolfe style algorithm [15], which has been systematically studied by Clarkson [9].

The problem of MEB with outliers also falls under the umbrella of the topic robust shape fitting [19, 2], but most of the approaches cannot be applied to high-dimensional data.  [35] studied MEB with outliers in high dimension, however, the resulting approximation is a constant $2$ that is not fit enough for the applications proposed in this paper.

Actually, our idea is inspired by a recent work about removing outliers for SVM [12], where they proposed a novel combinatorial approach called Random Gradient Descent (RGD) Tree. It is known that SVM is equivalent to finding the polytope distance from the origin to the Minkowski Difference of the given two labeled point sets. Gilbert algorithm [17, 16] is an efficient Frank-Wolfe algorithm for computing polytope distance, but a significant drawback is that the performance is too sensitive to outliers. To remedy this issue, RGD Tree accommodates the idea of randomization to Gilbert algorithm. Namely, it selects a small random sample in each step by a carefully designed strategy to overcome the adverse effect from outliers.

As mentioned before, we model outlier recognition as a problem of MEB with outliers in high dimension. Here we first introduce several definitions that are used throughout the paper.

From Definition 2 we can see that the major challenge is to determine the subset of $P$ which makes the problem a challenging combinatorial optimization. Therefore we relax our goal to its approximation as follows. For the sake of convenience, we always use $P_{\text{opt}}$ to denote the optimal subset of $P$, that is, $P_{\text{opt}}=\arg_{P^{\prime}}\min\{$ the radius of $MEB(P^{\prime})\mid P^{\prime}\subset P,\left|P^{\prime}\right|\geq(1-\gamma)n\}$, and $r_{\text{opt}}$ to denote the radius of $MEB(P_{\text{opt}})$.

When both $\epsilon$ and $\delta$ are small, the bi-criteria approximation is very close to the optimal solution with only a slight violation on the number of covering points and radius.

The rest of the paper is organized as follows. We introduce our main algorithm and the theoretical analyses in Section 2. The experimental results are shown in Section 3. Finally, we extend our method to handle multi-class inliers in Section 4.

We present our $(\epsilon,\delta)$-approximation algorithm for MEB with outliers in this section. Although the outlier recognition problem belongs to unsupervised learning, we can estimate the fraction of outliers in the given data before executing our algorithm. In practice, we can randomly collect a small set of samples from the given data, and manually identify the outliers and estimate the outlier ratio $\gamma$. Therefore, in this paper we assume that the outlier ratio is known.

To better understand our algorithm, we first illustrate the high-level idea. If taking a more careful analysis on the previously mentioned core-set construction algorithm [4], we can find that it is not necessary to select the farthest point to the center of $MEB(S)$ in each step. Instead, as long as the selected point has a distance larger than $(1+\epsilon)r_{\text{opt}}$, the minimal extent of improvement would always be guaranteed [10]. As a consequence, we investigate the following approach.

We denote the ball centered at point $c$ with radius $r>0$ as $Ball(c,r)$. Recall that $P_{\text{opt}}$ is the subset of $P$ yielding the optimal MEB with outliers, and $r_{\text{opt}}$ is the radius of $MEB(P_{\text{opt}})$ (see Section 1.2). In the $i$-th step, we add an arbitrary point from $P_{opt}\setminus Ball(c_{i},(1+\epsilon)r_{\text{opt}})$ to $S$ where $c_{i}$ is the current center of $S$. Based on the above observation, we know that a $(1+\epsilon)$-approximation is obtained after at most $\lceil{2/\epsilon}\rceil$ steps, that is, $\left|P\cap Ball(c_{i},(1+\epsilon)r_{\text{opt}})\right|\geq(1-\gamma)n$ when $i\geq\lceil{2/\epsilon}\rceil$.

However, in order to carry out the above approach we need to solve two key issues: how to determine the value of $r_{\text{opt}}$ and how to select a point belonging to $P_{opt}\setminus Ball(c_{i},(1+\epsilon)r_{\text{opt}})$. Actually, we can implicitly avoid the first issue via replacing the radius $(1+\epsilon)r_{\text{opt}}$ by the $k$-th largest distance from the points of $P$ to $c_{i}$, where $k$ is some appropriate number that will be determined in our following analysis. For the second issue, we have to take a small random sample instead of a single point from $P_{opt}\setminus Ball(c_{i},(1+\epsilon)r_{opt})$ and try each of the sampled points; this operation will result in a tree structure that is similar to the RGD Tree introduced by [12] for SVM. We present the algorithm in Algorithm 1 and place the detailed parameter settings, proof of correctness, and complexity analyses in Sections 2.2 & 2.3.

We illustrate Step 2(2)(a-c) of Algorithm 1 in Fig. 1.

We denote the tree constructed by Algorithm 1 as $\mathbb{H}$. The following theorem shows the success probability of Algorithm 1.

Before proving Theorem 2.1, we need to introduce several important lemmas.

In the remaining analyses, we always assume that such a root-to-leaf path $\mathcal{P}^{r}_{u}$ described in Lemma 3 exists and only focus on the nodes along this path. We denote $\hat{R}=(1+\epsilon)r_{\text{opt}}$ where $r_{\text{opt}}$ is the optimal radius of the MEB with outliers. Let $Ball(c_{v},r_{v})$ be the MEB covering $P\setminus P_{v}$ centered at $c_{v}$, and the radii of $MEB(\mathcal{P}^{r}_{v})$ and $MEB(\mathcal{P}^{r}_{v^{\prime}})$ be $\tilde{r}_{v}$ and $\tilde{r}_{v^{\prime}}$ respectively. Readers can refer to Fig. 2. The following lemma is a key observation for MEB.

The following lemma is a key for proving the quality guarantee of Algorithm 1. As mentioned in Section 2.1, the main idea follows the previous works [4, 10]. For the sake of completeness, we present the detailed proof here.

Now we prove Theorem 2.1 by the idea from [4]. Suppose no node in $\mathcal{P}^{r}_{u}$ makes the first event of Lemma 5 occur. As a consequence, we obtain a series of inequalities for each pair of radii $\tilde{r}_{v^{\prime}}$ and $\tilde{r}_{v}$ (see (5)). For ease of analysis, we denote $\tilde{r}_{v}=\lambda_{i}\hat{R}$ if the height of $v$ is $i$ in $\mathbb{H}$. By Inequality (5), we have

Combining the initial case $\lambda_{1}=0$ and (10), we obtain

by induction [4]. Note that the equality in (11) holds only when $i=1$, therefore,

Then, $\tilde{r}_{u}=\lambda_{h}\hat{R}>r_{\text{opt}}$ (recall that $u$ is the leaf node on the path $\mathcal{P}^{r}_{u}$), which is a contradiction to our assumption $\mathcal{P}^{r}_{u}\subset P_{\text{opt}}$. The success probability directly comes from Lemma 3. Overall, we obtain Theorem 2.1.

We analyze the time and space complexities of Algorithm 1 in this section.

Time Complexity. For each node $v$, we need to compute the corresponding approximate $MEB(\mathcal{P}^{r}_{v})$. To avoid computing the exact MEB costly, we apply the approximation algorithm proposed by [4]. See Algorithm 2 for details.

For Algorithm 2, we have the following theorem.

From Theorem 2.2, we know that a $(1+\varepsilon)$-approximation for MEB can be obtained when $N=1/\varepsilon^{2}$ with the time complexity $O\left(\frac{\left|{Q}\right|d}{\varepsilon^{2}}\right)$. Suppose the height of node $v$ is $i$, then the complexity for computing the corresponding approximate $MEB(\mathcal{P}^{r}_{v})$ is $O\left(\frac{id}{\varepsilon^{2}}\right)$. Further, in order to obtain the point set $P_{v}$, we need to find the pivot point that has the $(n-k)$-th smallest distance to $c_{v}$. Here we apply the PICK algorithm [6] which can find the $l$-th smallest from a set of $n$ ($l\leq n$) numbers in linear time. Consequently, the complexity for each node $v$ at the $i$-th layer is $O\left({\left({n+\frac{i}{\varepsilon^{2}}}\right)d}\right)$. Recall that there are $\left|{S_{v}}\right|^{i-1}$ nodes at the $i$-th layer of $\mathbb{H}$. In total, the time complexity of our algorithm is

If we assume $1/\varepsilon$ is a constant, the complexity $T=O(Cnd)$ is linear in $n$ and $d$, where the hidden constant $C=\left({\left({1+\frac{1}{\delta}}\right)\ln\frac{h}{\mu}}\right)^{h-1}$. In our experiment, we can carefully choose the parameters $\delta,\epsilon,\mu$ so as to keep the value of $C$ not too large.

Space Complexity. In our implementation, we use a queue $\mathcal{Q}$ to store the nodes in the tree. When the head of $\mathcal{Q}$ is popped, its $\left|S_{v}\right|$ child nodes are pushed into $\mathcal{Q}$. In other words, we simulate breadth first search on the tree $\mathbb{H}$. Therefore, $\mathcal{Q}$ always keeps its size at most $C=\left({\left({1+\frac{1}{\delta}}\right)\ln\frac{h}{\mu}}\right)^{h-1}$. Note that each node $v$ needs to store $\mathcal{P}^{r}_{v}$ to compute its corresponding MEB, but actually we only need to record the pointers to link the points in $\mathcal{P}^{r}_{v}$. Therefore, the space complexity of $\mathcal{Q}$ is $O(Ch)$. Together with the space complexity of the input data, the total space complexity of our algorithm is $O(Ch+nd)$.

By Theorem 2.1, we know that with probability at least $(1-\mu)(1-\gamma)$ there exists an $(\epsilon,\delta)$-approximation in the resulting tree. However, when outlier ratio is high, say $\gamma=0.5$, the success probability $(1-\gamma)(1-\mu)$ will become small. To further improve the performance of our algorithm, we introduce the following two boosting methods.

1. 1.

   Constructing a forest. Instead of building a single tree, we randomly initialize several root nodes and grow each root node to be a tree. Suppose the number of root nodes is $\kappa$. The probability that there exists an $(\epsilon,\delta)$-approximation in the forest is at least $1-(1-(1-\gamma)(1-\mu))^{\kappa}$ which is much larger than $(1-\gamma)(1-\mu)$.

2. 2.

   Sequentialization. First, initialize one root node and build a tree. Then select the best node in the tree and set it to be the root node for the next tree. After iteratively performing the procedure for several rounds, we can obtain a much more robust solution.

In our experiment, we test the algorithms on two random datasets and two benchmark image datasets. In terms of the random datasets, we generate the data points based on normal and uniform distributions under the assumption that the inliers usually locate in dense regions while the outliers are scattered in the space. The benchmark image datasets include the popular MNIST [25] and Caltech [14].

To make our experiment more convincing, we compare our algorithm with three well known methods for outlier recognition: angle-based outlier detection (ABOD) [23], one-class SVM (OCSVM) [31], and discriminative reconstructions in an autoencoder (DRAE) [34]. Specifically, ABOD distinguishes the inliers and outliers by assessing the distribution of the angles determined by each 3-tuple data points; OCSVM models the problem of outlier recognition as a soft-margin one-class SVM; DRAE applies autoencoder to separate the inliers and outliers based on their reconstruction errors.

The performances of the algorithms are measured by the commonly used $F1$ score $=\frac{2*\text{Precision}*\text{Recall}}{\text{Precision}+\text{Recall}}$, where precision is the proportion of the correctly identified positives relative to the total number of identified positives, and recall is the proportion of the correctly identified positives relative to the total number of positives in the dataset.

We validate our algorithm on the following two random datasets.

A toy example in $2$D. To better illustrate the intuition of our algorithm, we first run it on a random dataset in 2D. We generate an instance of $10,000$ points with the outlier ratio $\gamma=0.4$. The inliers are generated by a normal distribution; the outliers consist of four groups where the first three are generated by normal distributions and the last is generated by a uniform distribution. The four groups of outliers contains $800$, $1200$, $800$, and $1200$ points, respectively. See Fig. 3. The red circle obtained by our algorithm is the boundary to distinguish the inliers and outliers where the resulting $F1$ score is $0.944$. From this case, we can see that our algorithm can efficiently recognize the densest region even if the outlier ratio is high and the outliers also form some dense regions in the space.

High-Dimensional Points. We further test our algorithm and the other three methods on high-dimensional dataset. Similar to the previous 2D case, we generate $20,000$ points with four groups of outliers in $\mathbb{R}^{100}$; the outlier ratio $\gamma$ varies from $0.1$ to $0.5$. The $F1$ scores are displayed in Table 1, from which we can see that our algorithm significantly outperforms the other three methods for all the levels of outlier ratio.

In this section, we evaluate all the four methods on two benchmark image datasets.

From Section 3.2 and Section 3.3 we know that our method achieves the robust and competitive performances in terms of accuracy. In this section, we compare the time complexities of all the four algorithms.

ABOD has the time complexity $O(n^{3}d)$. In the experiment, we use its speed-up edition FastABOD which has the reduced time complexity $O((n^{2}+nk^{2})d)$ where $k$ is some specified parameter.

OCSVM is formulated as a quadratic programming with the time complexity $O(n^{3})$.

DRAE alternatively executes the following two steps: discriminative labeling and reconstruction learning. Suppose it runs in $N_{1}$ rounds; actually the two inner steps are also iterative procedures which both run $N_{2}$ iterations. Thus, the total time complexity of DRAE is $O(N_{1}N_{2}hdn)$, where $h$ is the number of the hidden layer nodes that can be generally expressed as $d/m$ ($m$ is a constant); then the total time complexity becomes $O(\tilde{C}nd^{2})$ where $\tilde{C}$ is a large constant depending on $N_{1}$, $N_{2}$ and $m$.

When the number of points $n$ is large, FastABOD, OCSVM, and DRAE will be very time-consuming. On the contrary, our algorithm takes only linear running time (see Section 2.3) and usually runs much faster in practice. For example, our algorithm often takes less than $1/2$ of the time consumed by the other three methods in our experiment.