PHOTOMETRIC SUPERNOVA CLASSIFICATION WITH MACHINE LEARNING

Although feature extraction generally reduces dimensionality from raw data, the number of features often remains large, making it difficult to determine whether or not a particular feature set separates the classes well. t-SNE (Van der Maaten & Hinton 2008) is a sophisticated technique designed to produce a low-dimensional representation of a high-dimensional space, while clustering similar features together. t-SNE works by first computing, for each pair of points, the probability that the two points are similar based on their Euclidean distance. The aim then is to find the corresponding low-dimensional values for these points that preserves the same probability. Stochastic neighbor embedding (SNE) determines these values by minimizing the sum of Kullback–Leibler divergences (a simple measure of divergence between probability distributions) over all data points, using a gradient descent algorithm with a stochastic element to avoid local minima. t-SNE differs from SNE in the choice of exact probability distribution defining the similarity between points (using a Student-t distribution instead of a Gaussian) and in the exact choice of cost function to minimize (see Van der Maaten & Hinton 2008 for details). If the classes are well-separated on the t-SNE plot, one can generally expect accurate classification. However, the converse is not strictly true and a feature set with poor class-separation may still be well-classified with a sophisticated machine learning algorithm. Each of the feature extraction methods discussed below has an accompanying t-SNE plot for visualization purposes.

Historically, classification of photometric supernovae (generally used to identify candidates for spectroscopic follow-up) has focused on the use of supernova templates constructed from existing data sets of supernovae (Sako et al. 2008). Light curve features such as stretch and color parameters have been used to classify photometric supernovae from SDSS (Sako et al. 2011). We choose the SALT2 (Spectral Adaptive Light curve Template 2) model (Guy et al. 2007) of type Ia supernovae, as it is the most commonly used.

In the SALT2 model, the flux in a given band of a type Ia light curve is given by

$$\begin{eqnarray}F(t,\lambda ) & = & {x}_{0}\times [{M}_{0}(t,\lambda )+{x}_{1}{M}_{1}(t,\lambda )+...]\\ & & \times {\rm{\exp }}[\mathrm{cCL}(\lambda )],\end{eqnarray} \tag{ 1 }$$

where t is the phase (time since maximum brightness in the B-band), λ is the rest frame wavelength, ${M}_{0}(t,\lambda )$ is average spectral sequence (using past supernova data as described in Guy et al. (2007)), and ${M}_{k}(t,\lambda )$ for $k\gt 0$ are higher order components to capture variability in the supernovae. In practice, only ${M}_{0}(t,\lambda )$ and ${M}_{1}(t,\lambda )$ are used. ${CL}(\lambda )$ is the average color correction law. For each object, the redshift z is also used as either a given or fitted parameter. Thus this method has five features: z, t₀ (the time of peak brightness in the B-band), x₀, x₁, and c. We use the implementation of SALT2 in sncosmo ³ (Barbary 2014) and obtain the best fitting parameters using MultiNest ⁴ (Feroz & Hobson 2008; Feroz et al. 2009, 2013) with the pyMultiNest ⁵ (Buchner et al. 2014) software. Table 1 lists the priors used on all parameters.

Table 1.  Priors on the SALT2 Model Parameters

Parameter Name	Prior
z	${ \mathcal U }(0.01,1.5)$
t0	${ \mathcal U }(-100,100)$
x0	${ \mathcal U }(-{10}^{-3},{10}^{-3})$
x1	${ \mathcal U }(-3,3)$
c	${ \mathcal U }(-0.5,0.5)$

Download table as:  ASCIITypeset image

Figures 2(a) and (b) show the t-SNE plots for the SALT2 features both without and with redshift information. The t-SNE plots show a lot of structure, partly due to the relative low-dimensionality of the feature space, which corresponds to a highly constrained 2D plot. The classes are fairly well-separated, especially when redshift is included, but it is clear that it is difficult to distinguish between type Ia and Ibc supernovae.

Zoom In Zoom Out Reset image size

The template-fitting approach relies on a significant amount of prior knowledge about supernovae. An alternative is to use a more general parametric description and use the fitted parameters as features in the machine learning pipeline. We use two different parametric models that have been proposed in the literature as good parameterizations of a supernova light curve and were both used in solutions to the SPCC. Model 1 is the parameterization used in Newling et al. (2011), while Model 2 is proposed in Karpenka et al. (2013).

We use these models because they are both proposed solutions to precisely the problem of supernova classification. While similar, the models have some differences, especially in the treatment of the tail of the light curves, the ability to describe double peaks and the small difference in the number of parameters, which makes it interesting to compare them.

We also considered the linear piecewise model proposed in Sanders et al. (2015) but found that the relatively large number of parameters (11 per filter) made it less robust when fitting to data, and consequently classification with this feature set did not typically perform well.

3.3.1. Model 1

The model used in Newling et al. (2011) describes the flux of a supernova (in any band) by

$$\begin{eqnarray}&&F(t)=A{\left(\displaystyle \frac{t-\phi }{\sigma }\right)}^{k}\exp \left(-\displaystyle \frac{t-\phi }{\sigma }\right){k}^{-k}{e}^{k}+{\rm{\Psi }}(t),\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}{\rm{\Psi }}(t)=\left\{\begin{array}{lc}0 & -\infty \lt t\lt \infty \\ \mathrm{cubic}\ \mathrm{spline} & \phi \lt t\lt \tau ,\\ \psi & \tau \lt t\lt \infty \end{array}\right.\end{eqnarray} \tag{ 3 }$$

and $A+\psi$ is the peak flux, ϕ is the starting time of the explosion, k governs the relative rise and decay times of the light curve, and σ is the width or "stretch" of the light curve. The time of the peak flux, τ, is given by $\tau =k\sigma +\phi$. The function ${\rm{\Psi }}(t)$ is a "tail" function that ensures that the flux tends to a finite value as $t\to \infty$ and the cubic spline is determined to have zero derivative at $t=\phi$ and $t=\tau$.

Each filter band is fitted individually with these 5 parameters (A, ϕ, σ, k, ψ), giving a total of 20 features for each object (or 21 if redshift is included). The parameters are all fitted with MultiNest, as in Section 3.2, and the best fit is taken to be the maximum posterior value for each parameter. The priors on the parameters (which are the same for all four filter bands) are given in Table 2.

Table 2.  Priors on the Model 1 Parameters

Parameter Name	Prior
$\mathrm{log}(A)$	${ \mathcal U }(0,10)$
ϕ	${ \mathcal U }(-60,100)$
$\mathrm{log}(\sigma )$	${ \mathcal U }(-3,4)$
$\mathrm{log}(k)$	${ \mathcal U }(-4,4)$
$\mathrm{log}(\psi )$	${ \mathcal U }(-6,10)$

Download table as:  ASCIITypeset image

Figures 2(c) and (d) show the t-SNE plots for the Model 1 feature set both without and with redshift information, respectively. In both cases, the features do not appear to be easily separable.

3.3.2. Model 2

The proposed parameterization from Karpenka et al. (2013) is similar to Model 1 but allows for the possibility of double peaks in light curves. The flux is given by

$$\begin{eqnarray}&&F(t)=A[1+B{(t-{t}_{1})}^{2}]\displaystyle \frac{\exp [-(t-{t}_{0})/{T}_{\mathrm{fall}}]}{1+\exp [-(t-{t}_{0})/{T}_{\mathrm{rise}}]},\end{eqnarray} \tag{ 4 }$$

where t is assumed to be the earliest measurement in the r-band light curve. We again fit all parameters using MultiNest and take the maximum posterior parameters. There are 6 parameters per filter band, thus forming a feature set of 24 features (or 25 including redshift) for each object. We use the same priors as in Karpenka et al. (2013), given in Table 3.

Table 3.  Priors on the Model 2 Parameters

Parameter Name	Prior
$\mathrm{log}(A)$	${ \mathcal U }(\mathrm{log}({10}^{-5}),\mathrm{log}(1000))$
$\mathrm{log}(B)$	${ \mathcal U }(\mathrm{log}({10}^{-5}),\mathrm{log}(100))$
t0	${ \mathcal U }(0,100)$
t1	${ \mathcal U }(0,100)$
${T}_{\mathrm{rise}}$	${ \mathcal U }(0,100)$
${T}_{\mathrm{fall}}$	${ \mathcal U }(0,100)$

Download table as:  ASCIITypeset image

The t-SNE plots for the Model 2 features can be seen in Figures 2(e) and (f) without and with redshift, respectively. As in Model 1, the features do not appear to be well-separated.

While the model-dependent approaches to feature extraction considered so far are useful for incorporating prior information about SN light curves, they are sensitive to mis-specification of light curve models. It is quite possible that important light curve characteristics cannot be captured efficiently by simple models, such as Models 1 and 2 described above. A model-independent, non-parameteric approach to feature extraction is therefore required as a complementary approach that is not dependent on prior information.

A wavelet decomposition is a good choice for feature extraction since certain wavelet transforms can be constructed that are approximately invariant under translation and stretch. One could then expect the wavelet coefficients of two similar supernovae to be similar, even if the light curves have different explosion times or redshifts.

Recently, Varughese et al. (2015) showed that a wavelet-based approach to photometric supernova classification can be highly effective. Here we consider a substantially different method in the actual wavelet decomposition of the light curves and subsequent classification, but also find wavelets to be an excellent feature extraction method.

The procedure we follow to extract features using wavelets is to first interpolate the light curves onto the same grid of points (using Gaussian processes), then to perform a redundant wavelet decomposition, and lastly to use PCA to reduce the dimensionality of the feature set. We note that there are many alternative dimensionality reduction techniques, such as independent component analysis (Comon 1994), autoencoders (Ballard 1987; Kramer 1991) and Sammon mapping (Sammon 1969), which may improve our classification results, but conclude that for this work PCA is more than sufficient.

3.4.1. Interpolation Using a Gaussian Process Regression

3.4.2. Wavelet Decomposition

3.4.3. Dimensionality Reduction with PCA

The algorithms used in this work are: naive Bayes (a basic probabilistic classifier), k-nearest neighbors (KNN) (a classifier based on clustering of features), a multi-layer perceptron (MLP; the simplest form of artificial neural network; ANN), a support vector machine (SVM; which works by finding a separating hyperplane between classes) and BDTs (an ensemble method that combines a multitude of decision trees, a weaker classifier). What follows is a brief overview of the five algorithms considered, in each case referring the reader to references for more details.

4.1.1. Naive Bayes (NB)

Given a set of features for an object, ${x}_{1},{x}_{2},\ldots {x}_{n}$, the Naive Bayes (NB) algorithm (e.g., Zhang 2004) states that the probability of that object belonging to class y is

$$\begin{eqnarray}&&P(y| {x}_{1},\ldots ,{x}_{n})\propto P(y)\displaystyle \prod _{i}P({x}_{i}| y),\end{eqnarray} \tag{ 5 }$$

assuming independence between features. P(y) can be estimated by determining the frequency of y in the training set. In the implementation of NB that we use, the likelihood, $P({x}_{i}| y)$, is assumed to be Gaussian, the parameters of which are determined from the training data by maximizing the likelihood. The predicted class is simply determined to be the class that maximizes Equation (5).

The advantage of NB is that it is fast and scales very well to high dimensions. The disadvantage is that it assumes independence between features and that the features are Gaussian-distributed. One or both of these assumptions is frequently broken. NB has recently been used in astronomy to classify asteroids to determine targets for spectroscopic follow-up (Oszkiewicz et al. 2014).

4.1.2. k-nearest Neighbors

4.1.3. Artificial Neural Networks

4.1.4. Support Vector Machines

4.1.5. Boosted Decision Trees

There are a variety of metrics for measuring the performance of a machine learning algorithm. Every machine learning algorithm we consider returns a probability for each classification. It is common to produce a set of predicted classifications by simply selecting, for each object, the class associated with the highest probability. However, this is not the full picture, as the probability can be used as a threshold at which an object is considered to belong to a particular class. For instance, a very pure sample of type Ia supernovae can be obtained by demanding a probability of, say, 95% before considering the object a type Ia. This purity will always come at the cost of completeness, however. Many metrics common in machine learning literature (such as purity or completeness) will depend on this subjective threshold and can be optimized for a specific class or goal (for example, obtaining a pure sample of type Ia supernovae for cosmology). Here we use a more general metric, receiver operating characteristic (ROC) curves, to directly compare algorithms.

A confusion matrix is useful for understanding most commonly used machine learning metrics. For a binary classification of two classes, positive and negative, the confusion matrix is shown in Table 4. An ideal classifier would produce only true positives without any contamination from false positives or losses from false negatives. In reality, any classification problem has a trade-off between purity and completeness.

An excellent method of comparing machine learning algorithms and visualizing this trade-off is to use ROC curves (Green & Swets 1966; Hanley & McNeil 1982; Spackman 1989) (see Fawcett (2004) for an introductory tutorial). These allow one to get an overview of the algorithms without having to select an arbitrary probability threshold at which to consider an object as belonging to a particular class. ROC curves can only compare one class (assumed to be the desired class) against all others, but can easily be constructed for each class in turn. In this paper, all ROC curves are plotted assuming the binary classification of type Ia versus non-Ia.

ROC curves are constructed by comparing the true positive rate (TPR) against the false positive rate (FPR), as the probability threshold is varied. The TPR, also called recall, sensitivity, or completeness, is the ratio between the number of correctly classified positive objects and the total number of positive objects in the data set,

$$\begin{eqnarray}&&{\rm{TPR}}=\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}.\end{eqnarray} \tag{ 6 }$$

The FPR, also referred to as the false alarm rate or 1-specificity, is the ratio between the number of objects incorrectly classified as positive and the total number of negative objects in the data set,

$$\begin{eqnarray}&&{\rm{FPR}}=\displaystyle \frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}.\end{eqnarray} \tag{ 7 }$$

It is straightforward to see that a good classification algorithm will maximize the TPR (thus resulting in a more complete data set) while simultaneously minimizing the FPR (thus reducing contamination). The area-under-the-curve (AUC) of a ROC curve provides a straightforward way to directly compare classifiers. An AUC of 1 represents perfect classification, visually represented by a curve that reaches the top left corner. By contrast, a diagonal line would correspond to an AUC of 0.5, meaning the classifier does no better than random. In the machine learning literature, an AUC higher than 0.9 typically indicates an excellent classifier.

In Section 6 we use ROC curves to examine the differences between machine learning algorithms and our three different approaches to feature extraction.

Commonly used metrics for classification include purity and completeness. These are only defined for classified objects, meaning some choice of probability threshold must be made. Usually, an object's class is selected as that with the highest probability. A more pure subset can be obtained by applying a high threshold probability. For any subset chosen, it is useful to determine the purity and completeness to evaluate it.

We can use the confusion matrix in Table 4 to define purity (also known as precision) as

$$\begin{eqnarray}&&{\rm{purity}}=\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}\end{eqnarray} \tag{ 8 }$$

and the completeness, which is equivalent to the TPR, as

$$\begin{eqnarray}&&{\rm{completeness}}=\displaystyle \frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}.\end{eqnarray} \tag{ 9 }$$

In supervised learning, the data are split into a training set, where each object has a known class label, and a test set, with unknown objects on which the classifier is run. Choice of training set is crucial, as a training set that is not representative of the entire data set can dramatically decrease a classifier's performance (see Section 6). Additionally, almost all machine learning algorithms have hyperparameters that need to be optimized. In this work, we use five-fold cross-validation (Kohavi et al. 1995) to determine the values of the hyperparameters (see Table 5 in the Appendix for some example parameter values).

In general, k-fold cross-validation divides the training set into k equally sized "folds." The classifier is trained on $k-1$ folds, leaving the remaining fold as a validation set on which the chosen scoring metric (here we use AUC) is evaluated. The score is averaged over all k validation sets to produce a final score for the full training set. We use a grid search algorithm with cross-validation to find the hyperparameters for each algorithm that maximize AUC.

The effect of redshift is summarized succinctly in Figure 5, comparing all feature extraction methods and machine learning algorithms. Figure 5 shows that providing redshift information plays an important role for SALT2 features, mildly improves the results of the parametric features, and is fairly unimportant to the wavelet features. The notable exception to this is when considering the best machine learning algorithm, BDTs. With this algorithm, including redshift information is not essential, as the algorithm is capable of differentiating between classes without redshift information.

Zoom In Zoom Out Reset image size

This raises the interesting point that with an algorithm with poor performance, such as KNN, redshift information is crucial to classify well with the SALT2 model. However, a better machine learning algorithm eliminates the need for redshift information (under the assumption of a representative training set). This is important since obtaining reliable redshifts for every object in future surveys will be challenging. We have shown that it will be possible to obtain a relatively pure subsample of supernovae without redshift information, on which to focus follow-up observations.

It is also interesting to note that the wavelet method, which uses no prior knowledge of supernovae, is able to perform as well as the current state-of-the-art classification technique, based on supernova templates. This is promising, since the simulated data set used here was constructed using largely the same known templates, so it is expected that the SALT2 model should perform well. However, the wavelet approach requires no previous knowledge and thus should perform well with new data sets. This also suggests that wavelet approaches to classification are likely to be useful for broader transient classification, as also noted in Varughese et al. (2015).

We also investigate the effect of using a representative or non-representative training set. It is well known that most machine learning algorithms underperform when the training set is in any way not representative of the distribution of features. The training set from the SPCC was known to be biased in the way spectroscopic follow-up data sets usually are, in that only the brightest objects tend to be followed up due to telescope time constraints and the difficulty in obtaining good spectra. It is not surprising then that, as shown in Figure 6, all feature extraction methods and all algorithms perform significantly worse with the non-representative training set. It should be noted that redshift information was not provided here, which only mitigates the effect of non-representativeness in the case of the SALT2 model. This is because the SALT2 model is already based on a large training set of supernovae and thus, if redshift information is given, will produce similar sets of features.

Zoom In Zoom Out Reset image size

The clear conclusion is that the more representative the training set, the better any classification method will perform. As the training set was only around 5% of the size of the full data set, we would argue that if spectroscopic follow-up time is limited, it is rather better to obtain fewer, fainter objects of all types, than to get a larger sample of brighter objects (as noted in Richards et al. 2012).

However, the effect of non-representativeness, also called the "domain adaptation problem" or "sample selection bias," can be mitigated using modern techniques in transductive transfer learning (Pan & Yang 2010). These include reweighting techniques to reweight the data in the test set to more closely match that of the training set and techniques to learn new sets of features that minimize the difference between the two data sets. Transfer learning has been successfully used in astronomy in Sasdelli et al. (2015). An in-depth study of transfer learning is beyond the scope of this paper but will likely prove useful for dealing with non-representativeness in the training sets of future supernova surveys.

It is possible to investigate the relative importance of each feature in a feature set using BDTs using the "gini importance" (also known as "mean decrease impurity"), described in Breiman et al. (1984). This is illustrated for each feature set in Figure 7, both with and without redshift. While there are many other methods to determine the relevance of each feature (see Li et al. 2016 for a summary), we only use this simple measure of feature relevance to gain insight into the features and use the full classification pipeline to evaluate the effectiveness of our chosen features, rather than performing initial feature selection based on relevance.

Zoom In Zoom Out Reset image size

For the SALT2 model, the most important features are the shape and color parameters. Note that counterintuitively, the importance of redshift goes down when external redshift information is added. By adding redshift information, the Ia's in the data set are better fit by the SALT2 model as expected, resulting in accurate estimates of the stretch and color parameters. On the other hand, with the redshift fixed the best fits for the non-Ia's result in unphysical stretch and color parameters. Because the distributions of these important features now differ much more between the non-Ia's and Ia's, they increase in importance when classifying. Because all the importances must add up to one, this necessarily dictates a decrease in the importance of other features, including the redshift.

There are a few features that are more important for the parametric models, but none are notable. Interestingly, it is the third and fourth principal components for the wavelet features that are most important, while all others have fairly equivalent importance.