ENSEMBLE LEARNING METHOD FOR OUTLIER DETECTION AND ITS APPLICATION TO ASTRONOMICAL LIGHT CURVES

We use a modified version of the standard k-NN methods described in the outlier detection and classification literature (Hodge & Austin 2004). The classical approach is an unsupervised algorithm (type 1) which assumes that anomalies are significantly distant from the "normal" data and thus can be found by computing a distance metric.

In other words, for each point, we compute its distance to its k-nearest neighbors, where the distance in ${{\mathbb{R}}}^{n}$ is the Euclidean distance, or L² norm. Note that this distance metric is generalized and can be replaced by other metrics as required by the data set. The sum of the distances are used as the final score.

K-NN methods are geometric in nature, and thus they are density-based outlier detection methods. This leads to drawbacks for this type of method. For example, two classes that are very similar in density may lead to misclassification if careful tuning of parameters is not performed. Also, the curse of dimensionality strikes when higher dimensions are reached. Because of the increasing sparsity of data in higher dimensions as well as the tendency for pairs of points to become close to equidistant from each other in high dimensional space, the k-nearest neighbor algorithm becomes unstable when the number of dimensions becomes too large. Therefore, dimensionality reduction, such as principal components analysis (PCA), is often suggested before attempting outlier detection via k-nearest neighbors for high dimensional data.

As with the other k-NN algorithm, this version is also discussed in standard classification/outlier detection approaches (Hodge & Austin 2004). K-NN2 also belongs to type 2 algorithms. If m of the k-nearest neighbors, for m < k, lie within a specific distance threshold, d, then the exemplar is considered to lie in a sufficiently dense region of the data distribution and is classified as normal. However, if there are less than m neighbors inside the distance threshold, then the exemplar is an outlier.

The naive computation complexity is $O({n}^{2})$, but can be reduced to $O(n\mathrm{log}n)$ with the use of kd-trees (see details in Section 5.2). The process of labeling outliers is similar to the above: assume that scores are normally distributed, then create a standard deviation cutoff. Such a cutoff can be tuned to select a certain proportion of possible outliers.

Both nearest neighbor methods share the same weaknesses when dealing with high dimensional data.

This is a supervised learning method modified from Nun et al. (2014), and it is a type 3 algorithm. Here, a classifier is trained on data with known-class labels. Then, for each point, a membership probability vector is produced that lists the probability of belonging to each category. The joint probability for the particular combination of membership probabilities can then be calculated. Nun et al. (2014) used the random forest classifier to produce the probabilities and a Bayesian Network to construct the joint probabilities. Outliers are then identified as points that have low joint probabilities because their class membership probability vectors are not been seen often enough in the training data.

An illustration of how this method works is presented in Figure 1. In most unsupervised methods, the red points in the middle will not be considered as outliers because they are in a region with point density that is not separable. In the naivest supervised methods, anything that is outside the boundaries is considered as an outlier. For the example of the outlier class in the middle, the product of the probabilities or the sum of the distances to the known classes may not be adequate as an outlier score, and therefore, the joint probability is a better measure for outliers. This case occurs when the conditional probability is lower than the marginal probability¹ as it can be seen from this simple illustration. The conditional probability shown on the left is smaller than the marginal probability shown on the right. This model will consider those objects as outliers.

Zoom In Zoom Out Reset image size

LoCI is a type 1 method that determines whether or not a point is an outlier by looking at the distribution of distances between pairs of objects. The algorithm involves a multi-granularity deviation factor (MDEF), and then selects a point as an outlier if its MDEF value deviates significantly (more than 3σ) from local averages. Intuitively, the LoCI method finds points that deviate significantly from the density of points in its local neighborhood. This is formalized in the MDEF concept. Let the r-neighborhood of an object p_i be the set of objects within distance r of p_i. Then, the MDEF at radius r for a point p_i is the relative deviation of its local neighborhood density from the average local neighborhood density in its r-neighborhood. An object with neighborhood density that matches the average local neighborhood density will have MDEF 0; outliers will have MDEFs far from 0. More formally, we have:

$$\begin{eqnarray}&&{\rm{MDEF}}({p}_{i},r,\alpha )=1-\displaystyle \frac{n({p}_{i},\alpha r)}{\bar{n}({p}_{i},\alpha ,r)}.\end{eqnarray} \tag{ 1 }$$

Here, n(p_i, αr) is the number of αr-neighbors of p_i; that is, the number of points $p\in {\mathbb{P}}$ such that $d({p}_{i},p)\leqslant \alpha r$, including p_i itself, so $n({p}_{i},\alpha r)\gt 0$. Let $\bar{n}({p}_{i},\alpha ,r)$ denote the average of n(p, αr) over the set of r-neighbors of p_i, and define

$$\begin{eqnarray}&&{\sigma }_{{\rm{MDEF}}}({p}_{i},r,\alpha )=\displaystyle \frac{{\sigma }_{\hat{n}}({p}_{i},r,\alpha )}{\bar{n}({p}_{i},r,\alpha )}\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&{\sigma }_{\bar{n}}^{2}{({p}_{i},\alpha ,r)=\displaystyle \frac{1}{n({p}_{i},r)}\displaystyle \sum _{p\in { \mathcal N }({p}_{i},r)}(n(p,\alpha r)-\bar{n}({p}_{i},r,\alpha ))}^{2}.\end{eqnarray} \tag{ 3 }$$

Then, to determine if a point is an outlier, we use the following algorithm: for each ${p}_{i}\in {\mathbb{P}}$, compute MDEF(p_i, r, α) and ${\sigma }_{{\rm{MDEF}}}({p}_{i},r,\alpha )$. If MDEF > 3σ_MDEF, flag p_i as an outlier. If, for any ${r}_{\min }\leqslant r\leqslant {r}_{\max }$, a point p_i is flagged as an outlier via the aforementioned mechanism, then we consider that point to be an outlier. These cutoffs can be determined on a per-problem basis, but in general we use the following: We set ${r}_{\max }\simeq {\alpha }^{-1}{R}_{p}$ and r_min such that we have ${\bar{n}}_{\min }=20$ neighbors. The number of neighbors can be tuned as necessary for the data set being analyzed. A graphical illustration of the aforementioned quantities can be found in Figure 2.

Zoom In Zoom Out Reset image size

Modified from Eskin (2000), LPD is also a type 1 algorithm. The main goal of this method is to model the underlying generating process for more complicated data as being drawn from a mixture of two simpler distributions. Outlier classification is achieved by first assuming that all points belong to the normal group. Conventionally, the points in the non-outlier group are assumed to be normally distributed, and a T-statistic can be computed as a measure of the probability of observing this grouping given that they are generated from a normal distribution. Then, for every order k, we iterate through the $\left(\genfrac{}{}{0em}{}{n}{k}\right)$ space and move k points at a time into the outlier class, which is usually assumed to follow a uniform distribution. Then, the score of the new grouping is computed, which represents the probability of observing the updated grouping. If the change in score is past a given threshold, the k points are marked as outliers in that order of space.

This means that optimizing for the threshold parameter can only be performed to target a certain number of outliers which again requires additional knowledge. For computational efficiency, the scoring is often based on log scores.

There are three main assumptions which form the basis of this method's efficacy. First, it is assumed that the non-outlier group can be fairly modeled by the assumed non -outlier probability distribution (in our case, a Gaussian distribution). Second, we assume that the anomalous elements are sufficiently distinct from the non-outlier elements. Finally, the model assumes that anomalies are few (<5% of the entire data set), otherwise the model will become distorted.

The greatest weakness of this algorithm is its dependence on the underlying simple models that form the mixture. This means that the normal points' need to follow a known distribution that can lead to computed t-statistics (such as joint distributions, etc.). When in doubt, a Gaussian model can be assumed as the underlying non-outlier distribution, but, if possible, the appropriate distribution should always be chosen.

With the mixture of experts approach, we assume that each outlier detection method performs best within a particular domain of the sample space. In this ensemble method, we combine the results from each method in a smart way so that the diversity of experts model can make up for the deficiencies of the individual methods over particular domains. Therefore, the result of each expert approach is weighted by values generated based on the location of the point in the 64-dimensional space and the approach's strength in this particular part of the feature space. Next, we explain what these parameters represent and the methodology we use to obtain them.

The gating probability g^x_i is the weight assigned to each expert i for data point ${\boldsymbol{x}}$. The weights are calculated using a soft-max gating network. The idea of using the soft-max procedure is to follow a "winner takes all" mechanism for choosing the best expert (Arbib 2003). This is useful in some network models for enforcing competition between different possible outputs of the network. In the soft-max procedure, a weight is assigned to each input so that all weights add to one, and the largest input receives the biggest weight. Let $I=\{{I}_{a}:a=1,\ldots ,N\}$ be the input and β be a positive parameter; the weight w_a(I; β) for each input I_a is defined by:

$$\begin{eqnarray}&&{w}_{a}(I;\beta )=\displaystyle \frac{{e}^{\beta {I}_{a}}}{{\displaystyle \sum }_{b}{e}^{\beta {I}_{b}}}\,.\end{eqnarray} \tag{ 4 }$$

Since the exponential function is monotonically increasing, in the limit of $\beta \longrightarrow \infty$, the weight of the largest input tends to1, and all the other weights tend to 0. On the other hand, performing the optimization of a function that takes the maximum of two or multiple numbers might be difficult to solve. For example f(x, y) = argmax(x, y) has a sharp corner along the line x = y and consequently is not differentiable. Therefore, an alternative to handle this is to use the soft-max function since it is differentiable everywhere.

Applying the soft-max function in the mixture of experts context, the gating probability g^x_i of each expert i for data point ${\boldsymbol{x}}$ is given by the following equation:

$$\begin{eqnarray}&&{g}_{i}^{x}=\displaystyle \frac{\exp ({{\boldsymbol{\eta }}}_{i}^{T}{\boldsymbol{x}})}{{\displaystyle \sum }_{j=1}^{k}\exp ({{\boldsymbol{\eta }}}_{j}^{T}{\boldsymbol{x}})},\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{x}}$ is a 64-dimensional data vector (one data point), and ${\eta }_{i}$ is the gating parameter for each expert i, with a total of e experts.

As can be seen from Equation (5), to obtain the gates g_i^x, it is first necessary to calculate the parameter η_i. η_i is a matrix of dimensions—number of experts by number of features. To determine their values, we perform a gradient descent optimization. Before explaining how this optimization is performed, it is first important to inquire into the meaning of the η_i variable.

The matrix η_i weights experts per feature but not per feature value. In other words, the feature value originates from x, and therefore, there is a linear relationship between the experts model and the feature values. This is not exactly the ideal scenario, and we consider this as the main drawback of the mixture of experts model approach. For example, in one-dimensional case, one expert wins, independently, all of the feature value. In the multidimensional case, there is a linear combination of expert methods; along the axes, one expert model wins regardless of the feature values, but in between the axes there is a linear combination of methods. In the ideal scenario, we would have one expert approach per feature domain or range. Unfortunately, this implies that there is an exponential number of parameters to be learned. An efficient method for this approach might be developed in a future work.

To find the optimal η, we define the following cost function I that we aim to minimize with respect to η:

$$\begin{eqnarray}&&I=\displaystyle \frac{1}{N}\displaystyle \sum _{i}\parallel {L}_{i}-{y}_{i}\parallel ,\end{eqnarray} \tag{ 6 }$$

where y_i represents the objects labels: 0 for non-outliers and 1 for outliers and L_i is the abnormality score for each particular data point i, as

$$\begin{eqnarray}&&{L}_{i}=\displaystyle \sum _{e}{g}_{{ie}}{S}_{{ie}}\end{eqnarray} \tag{ 7 }$$

where S_ie is the abnormality score of expert e for object i.

As we previously mentioned, the objective function is minimized by using a stochastic gradient descent (SGD) approach. The reader can find details about the implementation in Section 5.7 and an extended explanation in the Appendix in this document. In our work, we use a two-level hierarchical mixture of experts model, in which multiple expert networks are gated according to the process described above in order to provide a final recommendation. It is also possible, as seen in Figure 3, to have each expert network be recursively expanded into its own gating network and a set of its own sub-experts. In this way, multiple expert models with an expertise in specific domains can be combined using a hierarchy of gating probabilities to form one comprehensive prediction.

Zoom In Zoom Out Reset image size

As we mentioned in Sections 3.1, and 3.2, when using algorithms such as k-nearest neighbors on high dimensional data, it might be useful to first perform dimensionality reduction. Using the Feature Analysis for Time Series (FATS) library developed by Nun et al. (2015; for more details, please see Section 6), the high-dimensional astronomical light curve data was mapped to 64 astronomically descriptive features per light curve. We tried to apply in this new space other methods such as PCA. However, the principal components are chosen to pick the bulk of the data variance. Since we are looking for a few rare objects, PCA does not fit our goal. Furthermore, given that both k-NN1, and k-NN2 have a high accuracy (>90%) when testing them with known outliers, we decided not to dig in further in dimensionality reduction.

Computing all neighbor distances is naively an ${ \mathcal O }({n}^{2})$ problem. To reduce the running time of k-NN1, k-NN2, and LoCI, we use kd-trees (Friedman et al. 1977).² These structures are binary trees where tree nodes correspond to splits of space, allowing for ${ \mathcal O }(n\mathrm{log}n)$ performance for all neighbor distances. Other spatial partitioning methods exist as well, but performance here is sufficient for our application.

RF + JP is a supervised algorithm and consequently needs to be trained. The remaining four outlier detection methods are unsupervised and thus do not need training. Nevertheless, a base "construction" is necessary for all methods. This refers to making a pass through all the data; since these algorithms are either distance based (k-NN1, and k-NN2), density based (LoCI), or distribution based (LPD), they need to "learn" the data distances, the data densities, or the data distributions.

In the ideal scenario, each algorithm would be run on the whole data set, and each object would be compared to the rest. This would be memory and time inefficient in our case. Given that the data are large enough (roughly 250,000 light curves per field, see Section 6) and that we assume that outliers are a minority, we considered that "constructing" the algorithms with a portion of the data does not affect the final results. To do so, we took 5% of the data at random and built the models with it. In the case of k-NN1, k-NN2, and LoCI, we used it to construct the kd-trees. In the case of LPD, we used it to estimate the normal distribution. Consequently, to make our method scalable when analyzing a new light curve, we did not run each of the methods in the whole data but used the already built algorithms.

As we described in Section 3.4, the original algorithm considers every point that follows ${\rm{MDEF}}\gt 3{\sigma }_{{\rm{MDEF}}}$ as an outlier. Since in our case we do not want a binary score (outlier or non-outlier), but a score that quantifies the degree of abnormality of each object, we define the LoCI's score to simply be the value of ${\rm{MDEF}}$. Non-outliers will have ${\rm{MDEF}}$ values closer to 0, and outliers' values will be far from 0.

When using the LPD method, recall that we are modeling the generative process for the data as a mixture of two models; in particular, we use a Gaussian distribution for the normal class and a uniform distribution for the outlier class. Because of the large size of our data set, we assume for simplicity that moving one element out of the normal class and into the outlier class for testing does not change the parameters of the Gaussian distribution of the normal class. This is a reasonable assumption because of the size of our data set; no one point, even if it is an outlier, will change the distribution of the non-outlier class. Therefore, to find the outliers, the problem is simplified: we must only calculate the log-likelihood of each point belonging to the normal class (modeled as a Gaussian distribution with all the data). Then, we take the points with the smallest likelihoods and flag them as outliers.

As Zimek et al. (2014b) describe in their work, there are several aspects to consider when normalizing and combining scores from different models, but the most important aspect is to make the scores comparable. Since many of the outlier detection methods used in this work do not assign probabilities to outputs, their scores are not equivalent. Consequently, a sigmoid function is applied to the outputted scores of each model.

Specifically, we combined two sigmoid functions for each method: one function for the highest scores (outliers) and one for the lowest scores (non-outliers). This is done to avoid the bimodality of the data; outliers and non-outliers have their separable distribution (see top figure in Figure 5).

To accomplish this, we follow two steps for each model. First, we obtain the ROC (receiver operating characteristic) curve (Metz 1978) of the expert by using the training data set. For a better comprehension of our method, a description of ROC curves and its characteristics is presented next. ROC curves are graphical representations that show the performance of a binary classifier under different discrimination thresholds. Specifically, the ROC curve is a plot of the true positive rate (TPR) versus the false positive rate (FPR) for a range of decision cutpoints or thresholds. Once we obtain the ROC curve, we calculate its maximum Youden's index (Youden 1950). Youden's index is defined for each point in the ROC curve and corresponds graphically to the height of the point to the diagonal (see Figure 4). It is maximum for the optimum cutoff value of the performance test.

Zoom In Zoom Out Reset image size

Youden's index can consequently be used as a threshold to optimally outliers separate from non-outliers: anything above S_o is considered as outlier and anything below S_n as non-outlier. This index is important because it is the index at which the F-score is maximized.

As a second step, we calculate the median absolute deviation (MAD) of the non-outliers' scores, S_n, and the median absolute deviation of the outliers' scores, S_o. These values correspond to the δ parameters (sharpness of the sigmoid function, see below) of our two sigmoid functions.

We choose to use MAD as a variability measure because of its robustness to extreme differences in scale. For some models, we found out that a small number of outliers had extremely high values. Such extreme differences in scale among scores for an expert model leads to undesirable bimodality in our new score functions, a problem alleviated by utilizing the more robust mean absolute deviation metric.

By using the value of Youden's index Y we determine the new scores function (sigmoid function) S2 in the following way:

$$\begin{eqnarray}&&S2({i}_{i\in \{S(i)\gt Y\}})=\displaystyle \frac{1}{1+{e}^{\tfrac{-(S(i)-Y)}{{\delta }_{o}}}}\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&S2({i}_{i\in \{S(i)\lt Y\}})=\displaystyle \frac{1}{1+{e}^{\tfrac{-(S(i)-Y)}{{\delta }_{n}}}}\end{eqnarray} \tag{ 9 }$$

where δ_o is the MAD value for S_o and δ_n is the MAD value for S_n.

Figure 5 shows the difference between the old and the new distribution before and after pre-processing. In particular, it corresponds to the results obtained from testing on a balanced test set, which has the same amount of outliers and non-outliers in order to avoid the accuracy paradox (Zhu 2007). As expected, S2 has values from 0 to 1.

Zoom In Zoom Out Reset image size

SGD is an optimization method developed by Bottou (1998) which seeks for the minimization of a loss function $f(\hat{y},y)$ that measures the cost of predicting $\hat{y}$ when the actual answer is y (Bottou 2012). It is an iterative algorithm in which each step moves in the negative direction of the function gradient. To do so, in each iteration, a set of parameters are updated to minimize the loss function. While a standard gradient descent approach requires running through all the samples in the training set to perform a single update for a parameter in a particular iteration, SGD estimates the gradient using only one randomly picked example from the training set. It can be shown using the Robbins–Siegmund theorem that, for an appropriately chosen (decreasing) learning rate η and under other mild assumptions, SGD almost always converges (Bottou 2012).

During the optimization process of the mixture of expert models, we encountered overflow and underflow limitations. The values of $\exp ({{\boldsymbol{\eta }}}_{i}^{T}{\boldsymbol{x}})$ (see Equation (5)) are either too large or too little to be handled by a 64-bit machine (i.e., they cannot be represented in a floating point number). To deal with this, for every overflow or underflow case, we subtract from $\exp ({{\boldsymbol{\eta }}}_{i}^{T}{\boldsymbol{x}})$ its range (the difference between the largest and smallest value). This process does not change the results of the optimization since we subtract it from the nominator and the denominator in Equation (5) and consequently the operation cancels out.

Next, we present a small selection of the outliers we obtained for Field 77. We visually inspected the top 500 outliers among a total of 100,000 objects and cross-matched them with publicly available astronomical catalogs. Some of these catalogs are collections of known types. For example, LMC Long Period Variables (Fraser et al. 2008) is a collection of long period variables from LMC. On the other hand, catalogs like XMM-Newton (Watson et al. 2009) contain X-ray information, which can be useful to further understand the nature of the candidates. We did this in order to check whether there was any information about these objects and, if so, to determine if they belonged to a rare class. Nevertheless, the scope of this paper is not to analyze each of the detected anomalies. In the next section, we detail the second stage of this work where we aim to elucidate the nature of our outliers.

As an example of our results, we show in Figure 8 the light curves of some of the most interesting objects we detected.

Zoom In Zoom Out Reset image size

MACHO_77.7667.1117, as well as a large portion of our top outliers, belong to the "blue variable class" (Keller et al. 2002). This is a generic class without a unified light curve morphology, features, or underlying physics. Thus, further studies need to be conducted on these objects. We found out that MACHO_77.7552.49 is an Asymptotic Giant Branch Star (He-burning)⁴ and MACHO_77.7546.1541 is a Nova.⁵ MACHO_77.7918.68 is a hot emission-line star (Reid & Parker 2012). MACHO_77.7306.17 is a red giant (Olivier & Wood 2005).

Figure 9 shows some of the outliers that did not present any information when cross-matching them with known catalogs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we previously mentioned, the scope of this paper is not to find and identify every outlier in the MACHO catalog, but rather to show the efficiency and efficacy of our method. As we showed above, some of the objects have no information about their class, and it would take significant effort to find this information. Thus, to better understand the nature of the discovered outliers, we created an online catalog with our results http://iacs-courses.seas.harvard.edu/courses/TSC/OutliersWWW/Home.php. A diagram and screenshot of the website are shown in Figure 10. The goal is to clarify in a collaborative way the identity of unusual, rare, or unknown types of astronomical objects or phenomena that have not been yet classified. A user can contribute with information or possible conjectures about the presented outliers. A user can also provide outliers to be included in the catalog by sending them to any of the contact emails and later receive updates with the posted comments. By collecting this feedback, we also expect to select the most "interesting" outliers for a telescope follow-up observation and possibly solve some of the mysteries behind these anomalies.