Machine learning for transient discovery in Pan-STARRS1 difference imaging.

Prior to ingesting detections from IPP difference imaging into a MySQL database at Queen’s University Belfast (QUB) we perform pre-ingest cuts based on the detection of saturated, masked or suspected defective pixels within the PSF area. Taking as a typical night 3rd September 2013 (56548 MJD (Modified Julian Date)), the 7 nightly stacks produced 366267 detections ($\sim$52000 detections per stack), the pre-ingest cuts rejected 94.88% of these detections.
The $\sim$18750 detections passing the pre-ingest cuts are associated with transient candidates if there are two or more quality detections within the last seven observations of the field, including detections in more than one filter, and an rms scatter in the positions of $\le$0.5 arcsec. Each quality detection must be of more than 3$\sigma$ significance and have a Gaussian morphology (XYmoments $<$1.2). These post-ingest cuts also include checks for convolution issues, proximity to bright objects and ‘NaN’ values close to the centre of bright PSFs. 63% of the detections that passed the pre-ingest cuts were rejected during the post-ingest cuts. The remaining detections were promoted for human screening, where 37% of the detections were deemed to be real. These real transient candidates are cross-matched with catalogues of astronomical sources in the MDS fields. We use our own MDS catalogue and also extensive external catalogues (e.g. SDSS, GSC, 2MASS, NED, Milliquas²²2http://quasars.org/milliquas.htm, Veron AGN, X-ray catalogues) to make a contextual classification of supernova, variable star, active galactic nuclei or nuclear transient. We also cross-match with the Minor Planet Centre to reject asteroids, though most are removed during the construction of the nightly stacks.

In general supervised learning entails learning a model from a training set of data for which we provide the desired output for each training example. For the purposes of designating a detection as a real transient or a processing artefact, the desired output for each image is discrete. In such cases the problem is a supervised classification task for which there are a vast array of machine learning algorithms. In Section 4 we discuss the algorithms we try; however all such algorithms are trying to learn a model from the training data that will allow them to map the input parameterisation of each training example (see Section 3.2) to the desired output or label, while at the same time ensuring the model performs well on data not seen during the training phase. For building a real-bogus classifier this is an obvious avenue to pursue as we have a large sample of historical data for which we have labels provided by our current cuts and also through human eyeballing.
In Fig. 2 we show a sample of both real and bogus examples drawn at random from the training data. Often bogus detections show a combination of the factors we describe in Section 2.1 and typically this affects the centroiding during the source detection stage.

Over the past 3 years $\sim$1 million potential transients have been catalogued in the MDS by the TSS. Approximately 8000 of these objects have been selected by humans as real transients and promoted as potential targets for spectroscopic follow-up. As of the end of the survey in May 2014, 515 transients had spectroscopic classifications.
The aggregate catalogue information for all objects extracted by IPP and which pass the pre-ingest cuts described in Section 2.2 are stored in a database at QUB. Individual detections are associated with an object if they are spatially coincident within 0.5 arcsec. This information is presented to humans in the form of webpages³³3Similar webpages are made public for the PS1 3Pi Survey at: http://star.pst.qub.ac.uk/ps1threepi/psdb/public/. The webpages show all the photometric points produced by IPP in a multi-colour lightcurve. The number of photometric detections typically ranges from a few to a few dozen depending on the magnitude and timescale of the transient objects (see Rest et al. (2014), McCrum et al. (2014), Chomiuk et al. (2011) for examples of lightcurves). These webpages also present a subset of the image postage-stamps of the detections associated with an object (target image, reference image and difference image). This subset contains the first detections of the object of which there are always at least 2 (see Section 2.2) and up to 5 subsequent detections. Each object is then eyeballed by a human, those that appear to be real transients are promoted as potential targets for scientific follow-up, while artefacts are discarded.
Our training examples are drawn from the subset of detections we choose to present on the human digestible webpages for each object, as detailed above (typically 2-3 but less than 7). The majority of real examples were taken from detections of promoted objects with no spectroscopic classification. There is no guarantee that all detections of a promoted object are necessarily a result of good image subtractions. This prohibits simply assigning a label of real to all individual detections associated with a promoted target. In order to to ensure that we have a secure, reliable and clean set of real detections for training, we inspected and individually labelled 4352 detections (from 1919 different transients) as real, discarding any artefacts from the training set. We augmented this sample of real detections with data from 53 spectroscopically confirmed supernovae (from Dec. 2012 to Jan. 2014) for which we used the complete set of detections ($\sim$31 detections per object on average). These were again manually checked to remove bogus detections. We held out the first detections of all 53 SNe, which we use for testing in Section 5.7 and all detections of PS1-13avb, which we use in Section 5.6. This leaves an additional 1603 real training examples bringing the total to 5955 real detections.
Over the course of the survey approximately 800000 objects have been discarded as artefacts providing on the order of ${10}^6$ examples of bogus detections. We randomly sample from the available bogus examples and aim for 4 times more bogus examples as real, this is similar to the proportions used by Brink et al. (2013). Initial tests with classifiers showed that a significant proportion of the false positives appeared to be clean subtractions. We improved the purity of the bogus sample by examining the randomly selected bogus detections and added any detections that looked like real transient subtractions to the list of real examples (the effect of label contamination is further discussed in Section 5.3). This produced an extra 464 examples for the set of real detections resulting in a final total of 6419. We then selected 4 times as many bogus images from the remainder of the bogus examples we inspected, producing a sample of 25676 bogus detections.
The final training set contains 32095 training examples. We divide the training examples into 2 sets, distributed as follows; 75% for training and cross validation, and 25% for testing. The training examples are randomly shuffled prior to splitting with the caveat that all detections on the same night of a given object are included in the same set. This is to avoid detections with almost identical statistics being in multiple sets and giving a false impression of a classifiers performance. The construction of the data set is summarised in Table 1. The label for each training example is a 1 or 0, with 1 representing a label of real and 0 bogus.
The training set we have constructed is representative, containing examples of detections from different chips, seeing conditions, and filters, with various levels of S/N and examples of all types of processing artefact.

Machine learning algorithms require a 1-dimensional (1-D) vector representation of each training example, where each element of the vector corresponds to some numeric data or feature that may be useful to the algorithm for discerning examples belonging to each class. Previous work in the area of real-bogus classification, has focused on using parameters contained in catalogues generated by the processing pipeline and more complex features derived from that information to represent the detections, see Table 1 from Brink et al. (2013) and Table 1 from Romano, Aragon & Ding (2006).
The catalogue features available to individual surveys depend on the implementation of their image processing pipeline. When applying machine learning for real-bogus classification to a new survey it may not be possible to calculate these features based on the information available in the catalogues. There is also potential to spend a lot of time deriving and testing ways to combine the catalogue information that is available into features that we hope capture the differences between real and bogus detections. Bogus detections are the result of many factors and establishing a set of features that can encapsulate them all is difficult. In contrast simply representing the detections by their pixel intensity values requires no time spent developing or tuning feature extractors. Previous work that relies solely on the pixel data has proven effective for simple visual classification tasks, such as hand written digits (LeCun et al., 1998). For more complex tasks or to boost performance much of this work has been performed by learning a hierarchy of unsupervised features from the pixel data (Coates, Lee, & Ng (2011); LeCun et al. (1998)). Establishing a firm benchmark on the pixel intensity representation allows us to assess the potential gains from applying these more complex methods and is the main focus of this paper. Using this representation we expect the learning algorithm to identify salient relationships between pixels for the classification task. In the next section we discuss our choice of features and continue in the following section by describing the preprocessing steps we apply before training.

Artificial Neural Networks (ANN) comprise a number of interconnected nodes arranged into a series of layers. In this study we limit ourselves to a 3-layer ANN (consisting of an input layer, a hidden layer and an output layer) as those with more than one hidden layer need more careful training and require more computational power (Hinton, Osindero, & Teh, 2006). For our purposes we train feed-forward ANNs with back-propagation and randomly initialised weights, where the activation of each node is calculated with the logistic (sigmoid) function.
By limiting many of the choices for the structure of the ANNs we remove the need to select these hyperparameters during the cross validation phase in Section 4.4.1 which significantly reduces the complexity of the space we have to search. This economy of computation comes at the cost of not testing regions of the parameter space (e.g. other activation functions) and restricting the representational power of the ANNs by requiring a single hidden layer. We are however left with only 2 hyperparamters to choose namely, the number of nodes that make up the hidden layer, $s_2$ and the regularisation parameter, $\lambda$ through which we attempt to prevent overfitting. There is some suggestion (Murtagh, 1991; Geva & Sitte, 1992) that the optimal number of nodes in the hidden layer ($s_2$) is $2n+1$, where n is the number of input features. In our case $n$ is fixed at 400 input features, suggesting we should train ANNs with $s_2=801$ nodes, however training such large networks is beyond the scope of this work and we instead choose to test values of $s_2$ in the range 25-200.
We use our own vectorised implementation of ANNs written in python⁵⁵5https://www.python.org. The code relies on numpy⁶⁶6http://www.numpy.org for efficient array manipulations and scipy⁷⁷7http://docs.scipy.org/doc/scipy/reference/index.html for optimisation of the objective function.

Random Forests (RFs) aim to classify examples by building many decision trees from bootstrapped (sampled with replacement) versions of the training data (Breiman, 2001). Classifications are then assigned based on the average of the ensemble of decision trees. Each individual tree is grown by randomly sampling $k$ features from the $n$ input features and selecting the feature that best separates real examples from bogus as informed by the gini function. We use scikit-learn’s⁸⁸8http://scikit-learn.org/stable/index.html implementation of RFs where we select hyperparameters by assigning values to variables n_estimators, max_features and min_samples_leaf; the total number of trees in the ensemble, the number of features considered at each split and the minimum number of examples that define a leaf, below which no further splitting is allowed. RFs provide the ability to estimate the importance of each feature which we use in Section 5.2.

Support Vector Machines (SVMs) (Cortes & Vapnik, 1995) aim to find the hyperplane in the input feature space that optimally classifies training examples for linearly separable patterns, while simultaneously maximising the margin, the distance between the training examples which lie closest to the hyperplane, known as the support vectors. SVMs can be extended to non-linear patterns with the inclusion of a kernel, where the kernel transforms the original input data into a new parameter space. We again use scikit-learn’s implementation of SVMs where we choose the free parameters namely, the penalty parameter, C (similar to $\lambda$ for ANNs) and the kernel parameter gamma, which controls the local influence that support vectors have on the decision boundary. We only try SVMs with a Radial Basis Function (RBF) kernel, this being the most common choice and again reduces the parameter space that must be searched.

For each algorithm discussed above we need a method to choose the optimal combination of hyperparameters that will achieve the best performance for the classification task. In order to compare the relative performance of the different models we need some Figure of Merit (FoM). We use the FoM of Brink et al. (2013) which captures the essence of the problem we are trying to solve. The FoM is defined as the minimum Missed Detection Rate (MDR) (False Negative Rate) that gives a False Positive Rate (FPR) of 1%. That is, assuming we are willing to accept that 1% of the images deemed real by the classifier and promoted to human scanners will turn out to be bogus, what fraction of the real images would be discarded? With this we can select the model that would discard the least real images while 1% of images classified as real can be expected to be bogus.

As a last step toward boosting performance we investigated a selection of methods to combine the RF, SVM and ANN from Table 2. The predictions of the 3 methods are correlated; a candidate highly ranked by the RF is likely to also be highly ranked by the other 2 classifiers, but there are still detections of real transients that are discarded by only one of the classifiers. From Fig. 9 there are 24 detections labelled as real that only the RF wrongly rejects, it is these examples that we hope to recover by combining classifiers.
We tried only a few of the simplest combination strategies. First we simply classified a detection based on the majority vote of the 3 classifiers. Second we assigned each detection a hypothesis that was the mean of the hypothesis values output by each classifier. This produced a new distribution of mean hypotheses, where we again selected the decision boundary to produce the FoM. Finally we trained a SVM using the 3 hypotheses for each detection as the features representing that detection. In the end none of these methods outperformed the RF classifier, though the performance was comparable (see Table 3).
This result is unsurprising given that the classifiers are highly correlated and there is no guarantee that these methods will outperform the best individual classifier (Fumera & Roli, 2005). The RF is in itself an ensemble of classifiers (the individual decision trees) and may already incorporate much of the gain in performance we can expect from these simple methods.

Random Forests provide a built-in method to estimate the relative importance of each feature to the classification (Breiman, 2001). By inspecting the ‘depth’ at which each feature is used as a decision node we can estimate the relative importance of that feature, as those features used closer to the top of the tree will contribute to the prediction of a larger fraction of the training examples. The fraction of samples for which we expect a feature to contribute to the classification can be used to gauge its relative importance.
Fig. 10 shows the relative importance of each pixel determined from the training set. The relative importance metric is normalised such that it sums to 1. The most important features have the highest values and as would be expected are located in the centre of the image. The pixels on the edges of the images are thought to be important for identifying many of the bogus examples, where the object is not centred in the substamp and often lies at the edge. For reference if features were equally important they would each have a relative importance of $1/400=0.0025$.
Fig. 10 may suggest some redundancy in the features bounding the central pixels. It is expected that omitting these features would have little effect on the performance of our classifier as RFs are thought to be unaffected by the inclusion of noise variables in the feature vector (Biau, 2010). In contrast Brink et al. (2013) find that the MDR for their RF classifier improves by $\sim$4% by omitting noisy features using a backward feature selection method. The effect of feature selection is an interesting area for future work and attempts at optimisation.

We took care to eliminate label contamination in Section 3.1, by visually checking and manually labelling each training example. Nonetheless we expect that there remain some examples with incorrect labels. In this section we employ similar methods to those in Brink et al. (2013) to investigate the effect that label contamination has on our ability to train and test the optimal RF model.
First we investigate the effect of adding label contamination to the training set. We add contamination by randomly selecting a subset of the detections from the training set and flipping their labels. Those labelled as real are now labelled as bogus and vice versa. In Fig. 11 we plot the effect of randomly flipping labels in the training set while leaving the original labels in the test set untouched. The measured MDR appears fairly unaffected up to around 6% contamination. The approach of Brink et al. (2013) is robust to around 10% suggesting our method may be more susceptible to incorrectly labelled training data.
Next we flip labels in the test set, while using the original training set labels as they are. Given that the RF has been trained with correctly labelled data, for the most part we expect it to provide the correct labels for the images in the test set. However, the flipped labels affect our ability to accurately measure the FoM. Although the classifier makes sensible predictions, when we compare these predictions to the flipped labels the otherwise correct predictions are now evaluated as False Positives or Missed Detections. Fig. 11 shows how the FoM is affected as we increase the fraction of flipped labels, we see that even at low proportions labelling noise in the test set can have a significant effect.

To investigate the classifier performance as a function of Signal-to-Noise (S/N), we also follow a similar analysis to Brink et al. (2013). We plot the distribution of magnitudes for each example in the test set labelled as real in Fig. 12. We divide the examples into 11 bins, each spanning 1 magnitude in the range 13 to 24 mag. We then use the classifier to make a prediction for the examples in each bin and calculate the fraction of examples classified as bogus which we take as an estimate of the classifier performance for objects at that level of S/N. For objects with magnitudes $\gtrsim$20 there is a $\sim$6% chance of missing real detections. Counterintuitively the detection performance deteriorates for higher S/N objects. The number of examples of these cases are low as typically these objects result in artefacts from saturation and subsequent masking or unclean subtractions. However, this can also be understood as an effect of our feature representation, where we are learning classifications based on the relative intensity of pixels across the substamp. The tendency to misclassify such detections could stem from a combination of the large relative intensity differences between pixels in these substamps that often characterise artefacts and the low numbers of high S/N images of real transients. This explanation is further supported by both the ANN and SVM, which also misclassify these objects, suggesting the issue is with the data and not a consequence of the realisation of the RF. In the next section we try to identify any relationships in the missed detections.

We inspected the 172 missed detections (see Fig. 13) looking for similarities that may explain why they were rejected. We found that these missed detections are associated with 112 individual transients. Although we took care to limit label contamination during the construction of the training set, we identified some examples of obvious bogus detections mislabelled as real that account for a small fraction ($\sim$1%) of the missed detections.
We also find about 29% of the missed detections appear to be a result of faint galaxy convolution problems (see Section 2.1). These artefacts are difficult to identify by eye and as a result have been incorrectly labelled as real detections significantly contributing to the label contamination of the test set.

In order to demonstrate how we might expect the classifier to perform on a live data stream we first use the classifier to make predictions for the supernova PS1-13avb for which we held out all associated detections from both the test and training set. This object has been spectroscopically classified as a Type Ib SN and has a well sampled lightcurve from about -18 days pre-maximum to around 106 days post maximum, including exposures in all 5 filters ranging in magnitude from around 23 to 20 mags (see Fig. 16 top panel). We selected this object for its high quality lightcurve and magnitude range which represents the majority of objects discovered in the PS1 MDS. In the bottom panel of Fig. 16 we show the hypothesis for each epoch of this target. The plot shows that the hypothesis is consistently above the decision boundary of 0.539 (selected in Section 4.4.2) with the exception of the detection from 56480.406 MJD (Modified Julian Date) which shows the transient at a magnitude of g${}_{P1}=23.12\pm 0.21$ approaching the detection limit in this filter. The detection is displayed as an inset in Fig. 16 with its hypothesis of 0.506, showing the low S/N and deviation from a PSF-like morphology.

One of the major aims of recent supernova searches has been to try to detect the transient as soon after explosion as possible in order to trigger rapid follow-up to spectroscopically study regions of the transients evolution that remain relatively unexplored (Gal-Yam et al. (2014); Cao et al. (2013)). To this end we carry out a simple test by using the classifier to make predictions for the first detections of all 53 classified SNe in our database. Again we held these detections out from the training and test sets. In Table 4 we list the 53 SNe and the details of the first detections along with the hypothesis for each detection. The classifier correctly predicts all detections as real and had it been running on a live data stream would have promoted all objects to humans for follow-up.