Physics
Computer Science-2021.gz:143140
Computer Science
abe89c5be897b20f8ecd3e39f8b9a857bebfa72b
SCONE: Supernova Classification with a Convolutional Neural Network

We present a novel method of classifying Type Ia supernovae using convolutional neural networks, a neural network framework typically used for image recognition. Our model is trained on photometric information only, eliminating the need for accurate redshift data. Photometric data is pre-processed via 2D Gaussian process regression into two-dimensional images created from flux values at each location in wavelength-time space. These"flux heatmaps"of each supernova detection, along with"uncertainty heatmaps"of the Gaussian process uncertainty, constitute the dataset for our model. This preprocessing step not only smooths over irregular sampling rates between filters but also allows SCONE to be independent of the filter set on which it was trained. Our model has achieved impressive performance without redshift on the in-distribution SNIa classification problem: $99.73 \pm 0.26$% test accuracy with no over/underfitting on a subset of supernovae from PLAsTiCC's unblinded test dataset. We have also achieved $98.18 \pm 0.3$% test accuracy performing 6-way classification of supernovae by type. The out-of-distribution performance does not fully match the in-distribution results, suggesting that the detailed characteristics of the training sample in comparison to the test sample have a big impact on the performance. We discuss the implication and directions for future work. All of the data processing and model code developed for this paper can be found in the SCONE software package located at github.com/helenqu/scone.


INTRODUCTION
The discovery of the accelerating expansion of the universe (Perlmutter et al. 1999;Riess et al. 1998) has led to an era of sky surveys designed to probe the nature of dark energy. Type Ia supernovae (SNe Ia) have been instrumental to this effort due to their standard brightness and light curve profiles. Building a robust dataset of SNe Ia across a wide range of redshifts will allow for the construction of an accurate Hubble diagram that will enrich our understanding of the expansion history of the universe as well as place constraints on the dark energy equation of state.
Modern-scale sky surveys, including SDSS, Pan-STARRS, and the Dark Energy Survey, have identified thousands of supernovae throughout their operational lifetimes (Frieman et al. 2008;Chambers et al. 2016;Smith et al. 2020). However, it has been logistically challenging to follow up on most of these detections spectroscopically. The result is a low number of spectroscopically confirmed SNe Ia and a large photometric dataset of SNe Ia candidates. The upcoming Rubin Observatory Legacy Survey of Space and Time (LSST) is projected to discover 10 7 supernovae (LSST Science Collaboration et al. 2009), with millions of transient alerts each observing night. As spectroscopic resources are not expected to scale with the size of these surveys, the ratio of spectroscopically confirmed SNe to total detections will continue to shrink. With only photometric data, distinguishing between SNe Ia and other types can be difficult. A reliable photometric SNe Ia classification algorithm will allow us to tap into the vast potential of the photometric dataset and pave the way for confident classification and analysis of the ever-growing library of transients from current and future sky surveys.
Significant progress has been made in the past decade in the development of such an algorithm. Most approaches involve lightcurve template matching (Sako et al. 2011), or feature extraction paired with either sequential cuts (Bazin et al. 2011;Campbell et al. 2013) or machine learning algorithms (Möller et al. 2016;Lochner et al. 2016;Dai et al. 2018;Boone 2019). Most recently, the spotlight has been on deep learning techniques since it has been shown that classification based on handcrafted features is not only more time-intensive for the researcher but is outperformed by neural networks trained on raw data (Charnock & Moss 2017;Moss 2018;Kimura et al. 2017). Since then, many neural network architectures have been explored for SN photometric classification, such as PELICAN's CNN architecture (Pasquet et al. 2019) and SuperNNova's deep recurrent network (Möller & de Boissière 2020).
Several photometric classification competitions have been hosted, including the Supernova Photometric Classification Challenge (SPCC) (Kessler et al. 2010) and the Photometric LSST Astronomical Time Series Classification Challenge (PLAsTiCC) (The PLAsTiCC team et al. 2018). These have not only resulted in the development of new techniques, such as PSNID (Sako et al. 2011) and Avocado (Boone 2019), but have also provided representative datasets available to researchers during and after the competition, such as the PLAs-TiCC unblinded dataset used in this paper.
In this paper we present SCONE, a novel application of deep learning to the photometric classification problem. SCONEis a convolutional neural network (CNN), an architecture prized in the deep learning community for its state-of-the-art image recognition capabilities (LeCun et al. 1989(LeCun et al. , 1998Krizhevsky et al. 2012;Zeiler & Fergus 2014;Simonyan & Zisserman 2014). Our model requires raw photometric data only, precluding the necessity for accurate redshift approximations. The dataset is preprocessed using a lightcurve modeling technique via Gaussian processes described in Boone (2019), which alleviates the issue of irregular sampling between filters but also allows the CNN to learn from information in all filters simultaneously. The model also has relatively low computational and dataset size requirements without compromising on performance -400 epochs of training on our ∼ 10 4 dataset requires around 15 minutes on a GPU.
We will introduce the datasets used to train and evaluate SCONE its computational requirements, as well as the algorithm itself in Section 2. Section 3 will focus on the performance of SCONE on both binary and categorical classification, and Section 4 presents an analysis of misclassified lightcurves and heatmaps for both modes of classification.

Datasets
The PLAsTiCC training and test datasets were originally created for the 2018 Photometric LSST Astronomical Time Series Classification Challenge.
The PLAsTiCC training set includes ∼ 8000 simulated observations of low-redshift, bright astronomical sources, representing objects that are good candidates for spectroscopic follow-up observation. This dataset will be referred to in future sections as the "spectroscopic dataset". We use this dataset to evaluate the outof-distribution performance of SCONE in section 3.1.2.
The PLAsTiCC test set consists of 453 million simulated observations of 3.5 million transient and variable sources, representing 3 years of expected LSST output (Kessler et al. 2019). The objects in this dataset are generally fainter, higher redshift, and do not have associated spectroscopy. Note that most of the results presented in this paper are produced from this dataset alone and will be referred to as "the dataset", "the main dataset", or "the PLAsTiCC dataset" in future sections.
All observations in both datasets were made in LSST's ugrizY bands and realistic observing conditions were simulated using the LSST Operations Simulator (Delgado & Schumacher 2014). While PLAsTiCC includes data from many other transient sources, we are using only the supernovae in the datasets. We selected all of the type II, Iax, Ibc, Ia-91bg, Ia, and SLSN-1 sources (corresponding to true target values of 42, 52, 62, 67, 90, and 95) from this dataset and chose only well-sampled lightcurves by restricting ourselves to observations simulating LSST's deep drilling fields (ddf=1). peak mjd values, the modified Julian date of peak flux for each object, were calculated for both the main dataset and the spectroscopic dataset by taking the signal-to-noise weighted average of all observation dates. peak mjd is referred to as t peak in future sections. The total source count for the spectroscopic set is 4,556, and the total source count for the main dataset is 32,087. A detailed breakdown by type is provided in Table 1.
The categorical dataset was created using SNANA (Kessler et al. 2009) with the same models and redshift distribution as the main dataset in order to perform categorical classification with balanced classes. 2,000 examples of each type were randomly selected to constitute a class-balanced dataset of 12,000 examples.

Quality Cuts
In order to ensure that the model is learning only from high-quality information, we have instituted some additional quality-based cuts on all datasets. These cuts are based on lightcurve quality, so all metrics are defined for a single source. The metrics evaluated for these cuts are as follows: • number of detection datapoints (n detected ): number of observations where the source was detected. We chose a detection threshold of S/N > 5, based on Fig. 8 of Kessler et al. (2015) • cumulative signal-to-noise ratio (CSN R): cumulative S/N for all points in the lightcurve where f represents the flux measurements from all observations of a given source, σ f represents the corresponding uncertainties, t represents the timestamps of all observations of a given source, t first is the time of initial detection, and t last is the time of final detection. Our established quality thresholds require that: • n detected ≥ 5 • CSN R > 10 • t active ≥ 30 days for lightcurve points in the range t peak −50 ≤ t ≤ t peak + 130. 1,150 out of 4,556 sources passed these cuts in the spectroscopic dataset and 12,611 out of 32,087 sources passed these cuts in the main dataset. The makeup of these datasets is detailed in Table 2. The categorical dataset was created from sources that already passed the cuts, so the makeup is unchanged.

Class Balancing
Maintaining an equal number of examples of each class, or a balanced class distribution, is important for machine learning datasets. Balanced datasets allow for an intuitive interpretation of the accuracy metric as well as provide ample examples of each class for the machine learning model to learn from. As shown in Table 2, the natural distribution of the spectroscopic dataset is more abundant in Ia sources than non-Ia sources. Thus, all non-Ia sources were retained in the class balancing process for binary classification and an equivalent number of Ia sources were randomly chosen. SNIax and SNIa-91bg sources were labeled as non-Ia sources for binary classification. The class-balanced spectroscopic dataset has 496 sources of each class for a total of 992 sources.
In contrast, the natural distribution of the main dataset is more abundant in non-Ia sources than Ia sources. Thus, all Ia sources were retained in the class balancing process for binary classification and an equivalent number of non-Ia sources were randomly chosen. The random selection process does not necessarily preserve the original distribution of non-Ia types. The classbalanced dataset has 6,128 sources of each class for a total of 12,256 sources.
The categorical dataset of 2,000 sources for each of the 6 types was created explicitly for the purpose of retaining balanced classes in categorical classification, as mentioned in Section 2.1.
All datasets were split by class into 80% training, 10% validation, and 10% testing. Splitting by class ensures balanced classes in each of the training, validation, and test sets.

Heatmap Creation
Prior to training, we preprocess our lightcurve data into heatmap form. First, all observations are labeled with the central wavelength of the observing filter according to Table 3, which was calculated from the filter functions used by The PLAsTiCC team et al. (2018). We then use the approach described by Boone (2019) to apply 2-dimensional Gaussian process regression to the raw lightcurve data to model the event in the wave-  length (λ) and time (t) dimensions. We use the Matern 32 kernel with a fixed 6000Å characteristic length scale in λ and fit for the length scale in t. Once the Gaussian process model has been trained, we obtain its predictions on a λ, t grid and call this our "heatmap". Our choice for the λ, t grid was t peak − 50 ≤ t ≤ t peak + 130 with a 1-day interval and 3000 < λ < 10, 100Å with a 221.875Å interval. The significance of this choice is explored further in Section 3. The result of this is a n λ ×n t image-like matrix, where n λ and n t are the lengths of the wavelengths array and times array, respectively, given to the Gaussian process. We also take into account the uncertainties on the Gaussian process predictions at each time and wavelength, producing a second image-like matrix. We stack these two matrices depthwise and divide by the maximum flux value to constrain all entries to [0,1]. This matrix is our input to the convolutional neural network. Figure 1 shows the raw lightcurve data in the t peak − 50 ≤ t ≤ t peak + 130 range for each filter in blue, the Gaussian process model in gray, and the resulting flux and uncertainty heatmaps at the bottom. Figure 2 shows a representative example of a heatmap for each supernova type.

Convolutional Neural Networks
Convolutional neural networks (CNNs) are a type of artificial neural network that makes use of the convolution operation to learn local, small-scale structures in an image. This is paired with an averaging or sub-sampling layer, often called a pooling layer, that reduces the resolution of the image and allows the subsequent convolutional layers to learn hierarchically more complex and less localized structures.
In a convolutional layer, each unit receives input from only a small neighborhood of the input image. This use of a restricted receptive field, or kernel size, resembles the neural architecture of the animal visual cortex and allows for extraction of local, elementary features such as edges, endpoints, or corners. All units, each corresponding to a different small neighborhood of the image, share the same set of learned weights. This allows them to detect the presence of the same feature in each neighborhood of the image. Each convolutional layer often has several layers of these units, each of which is called a filter and extracts a different feature. The output of a convolutional layer is called a feature map.
Pooling layers reduce the local precision of a detected feature by sub-sampling the each unit's receptive field, often 2 × 2 pixels, based on some rule. Average-pooling, for example, extracts the average of the 4 pixels, and max-pooling extracts the maximum value. Assuming no overlap in the receptive fields, the spatial dimensions of the resulting feature maps will be reduced by half.
Dropout ) is a commonly used regularization technique in fully connected layers. A dropout layer chooses a random user-defined percentage of the input weights to set to zero, improving the robustness of the learning process.
A convolutional neural network typically consists of alternating convolutional layers and pooling layers, followed by a series of fully-connected layers that learn a mapping between the result of the convolutions and the desired output.

SCONE Architecture
The relatively simple architecture of SCONE, shown in without compromising on performance. It has a total of 22,606 trainable parameters for categorical classification and 22,441 trainable parameters for binary classification when trained on heatmaps of size 32×180×2 (h×w×d).
As mentioned in Section 2.4, each heatmap is divided by its maximum flux value for normalization. After receiving the normalized heatmap as input, the network pads the heatmap with a column of zeros on both sides, bringing the heatmap size to 32 × 182 × 2. Then, a convolutional layer is applied with h filters and a kernel size of h × 3, which in this case is 32 filters and a 32 × 3 kernel, resulting in a feature map of size 1 × 180 × 32. We reshape this feature map to be 32 × 180 × 1, apply batch normalization, and repeat the above process one more time. We have now processed our heatmap through two "convolutional blocks" with an output feature map of size 32 × 180 × 1.
We apply 2 × 2 max pooling to our output, reducing its dimensions to 16 × 90 × 1, and pass it through two more convolutional blocks, but this time h = 16.
We pass our output through a final 2 × 2 max pooling layer, resulting in a 8 × 45 × 1 feature map. This is subsequently flattened into a 360-element array and passed through a 50% dropout layer. A 32-unit fully connected layer followed by a 30% dropout layer feeds into the final layer. For binary classification, this is a node with a sigmoid activation that returns the model's predicted Ia probability. For categorical classification, the final layer contains 6 nodes with softmax activations that return the respective probabilities of each of the 6 SN types. Both of these versions of SCONE are shown in Figure 3.
The model is trained with the binary crossentropy loss function for binary classification and the sparse categorical crossentropy loss function for categorical classification. Both classification modes use the Adam optimizer (Kingma & Ba 2017) at a constant 1e-3 learning rate for 400 epochs.

Evaluation Metrics
Accuracy is defined as the number of correct predictions divided by the number of total predictions. We also evaluated the model on a number of other performance metrics: purity, efficiency, and AUC.
Purity and efficiency are defined as: purity= TP TP+FP ; efficiency= TP TP+FN where TP is true positive, FP is false positive, and FN is false negative.
AUC, or area under the receiver operating characteristic (ROC) curve, is a common metric used to evaluate binary classifiers. The ROC curve is created by plotting the false positive rate against the true positive rate at various discrimination thresholds, showing the sensitivity of the classifier to the chosen threshold. A perfect classifier would score a 1.0 AUC value while a random classifier would score a 0.5.

Computational Requirements
Due to the minimal number of trainable parameters and dataset size, the time and hardware requirements for training and evaluating with SCONE are relatively  low. The first training epoch on one NVIDIA V100 Volta GPU takes approximately 2 seconds, and subsequent training epochs take approximately 1 second each with TensorFlow's dataset caching. Training epochs on one Haswell node (with Intel Xeon Processor E5-2698 v3), which has 32 cores, take approximately 26 seconds each.

Binary Classification
Our model achieved 99.93 ± 0.06% training accuracy, 99.71 ± 0.2% validation accuracy, and 99.73 ± 0.26% test accuracy on the Ia vs. non-Ia binary classification problem performed on the class-balanced dataset of 12,256 sources. Figure 4 shows the confusion matrices for binary classification, and Table 4 shows the model's performance on all the evaluation metrics described in Section 2.7. Both Figure 4 and Table 4 were created with data from five independent training, validation, and test runs of the classifier. Unless otherwise noted, the default threshold for binary classification is 0.5, where classifier confidence equal to or exceeding 0.5 counts as an SNIa classification, and vice versa. Altering this threshold produces Figure 5, the ROC curve of one of these test runs.

Heatmap Dimensions
We explored the binary classification performance of different heatmap dimensions in both the time (width)  and wavelength (height) axes. For our 7100Å wavelength range (3000Å< λ < 10, 100Å), we chose intervals of 221.875Å, 443.8Å, 887.5Å, and 1775Å, resulting in heatmaps with 32, 16, 8, and 4 wavelength "bins", respectively. For our 180 day range (t peak − 50 ≤ t ≤ t peak + 130 ), we chose intervals of 1, 2, 6, and 18 days, resulting in heatmaps with 180, 90, 30, and 10 time bins. Figure 6 shows the training, validation, and test accuracies for each choice of wavelength and time dimensions. Our classifier seems relatively robust to these changes, showing minimal performance impacts for wavelength bins ≥ 16 and impressive performance for the smaller sizes as well. Test accuracy drops a maximum of 2.67% between the best-and worst-performing variants, 32 × 180 and 4 × 10. This is noteworthy as the 32 × 180 heatmaps contain 144 times the number of pixels of the 4 × 10 heatmaps.
The performance seems to unilaterally improve as expected as the number of wavelength bins increase, but increasing the number of time bins seems to yield varying, though likely not statistically significant, results. We have reported all of our SCONE results using one of the best performing variants, the 32 × 180 heatmaps.

Out-of-Distribution Results
Preliminary exploration into the out-of-distribution task of training on the spectroscopic dataset and testing on the main dataset yielded 80.6% test accuracy. The full results of training on 85% of the spectroscopic dataset, validating on the remaining 15%, and testing on the full main dataset are shown in Table 5. Since the redshift distribution of the spectroscopic dataset is skewed toward lower redshifts, testing on class-balanced low-redshift subsets of the main dataset yielded 83% test accuracy for z < 0.4 and 87% test accuracy for z < 0.3. Boone (2019) introduced redshift augmentation to mitigate this effect. The mismatch between the test and training sets, however, comes in other forms. For example, if the training data is generated using models rather than real data, differences in the characteristics of the spectral surfaces can have an impact on classification. The unknown relative rates of of each type of event also affect the overall performance. Since out-of-distribution robustness is an integral part of the challenge of photometric SNe classification, improving the performance of SCONE on these metrics will be the topic of a future paper. Class-Balanced Test Set

Categorical Classification
In addition to binary classification, SCONE is able to perform categorical classification and discriminate between different types of SNe. We performed 6way categorical classification with the same PLAsTiCC dataset used for binary classification as well as the class-balanced dataset described in Section 2.1. Our model differentiated between SN types Ia, II, Ibc, Iax, SN-91bg, and SLSN-1. On the PLAsTiCC dataset (not class-balanced), it achieved 99.26 ± 0.18% training accuracy, 99.13 ± 0.34% validation accuracy, and 99.18 ± 0.18% test accuracy. The confusion matrices in Figure 8 show the average by-type breakdown for five independent runs.
On the balanced dataset, it achieved 97.8 ± 0.32% training accuracy, 98.52±0.28% validation accuracy, and 98.18 ± 0.3% test accuracy. The confusion matrices in Figure 7 show the average and standard deviations of the by-type breakdown for 5 independent runs. It is worth noting that we trained and tested with the PLAsTiCC dataset even though it is not class-balanced for this task to try to evaluate the model's performance on a dataset emulating the relative frequencies of these events in nature.

Binary Classification
According to the data presented in Figure 4, the model seems to mispredict about the same number of Ia's as non-Ia's. An average of 3.33 Ia's are mispredicted in the training set compared with 3.5 non-Ia's. For the validation set, 2 Ia's are mispredicted compared to 1.5, and 1.33 Ia's and 2 non-Ia's are mislabeled for the test set.
The misclassified set summarized in Table 6 is the result of one of the five runs represented by the data in Figure 4. In this example, the model missed 5 examples total during testing -1 SNIa, 2 SNII, and 1 SNIax. It was > 90% confident about all of these misclassifications, which is certainly not the case for categorical classification. This could be due to the fact that there are more examples of each type for binary classification than categorical (∼ 6, 000 and 1,200 per type for training, respectively).

Categorical Classification
The data in Figures 7 and 8 show some level of symmetry between misclassifications. SNIax and SLSN-1 seem to be easily distinguishable across the board, for example, with 0's in all relevant cells in both figures except one. In Figure 7, SNIbc and SLSN-1 are seemingly very similar, as they are misclassified as one another at similarly high rates.
There are notable differences between Figures 7 and 8, however. In Figure 7, representing the classifier's performance on class-balanced categorical classification, the model mispredicts SNIa's as other types at a similar rate as non-Ia's mispredicted as Ia's. An average of 3.2 Ia's are mispredicted, whereas an average of 4 non-Ia's are misclassified as Ia. In the confusion matrix shown in Figure 8, significantly more non-Ia's were mispredicted as Ia. An average of 1 Ia was mispredicted compared to an average of 4.2 non-Ia's misclassified as Ia. The rate of SNIbc's mispredicted as SLSN-1 is also significantly lower for the PLAsTiCC dataset than for the balanced dataset. These observations further reinforce the impact of imbalanced classes in classification tasks.
The misclassified set summarized in Table 7 is the result of one of the five class-balanced categorical classification runs represented by the data in Figure 7.
One point of interest is the lack of symmetry between misclassifications, in contrast with the analysis of Figures 7 and 8. This is clear in the significantly larger number of SLSN-1 misclassified as SNIbc (7) compared with the number of SNIbc misclassified as SLSN-1 (1). SNIax is also more often misclassified as other types (3 as SNIa, 3 as SNIbc, and 1 as SNIa-91bg) than non-Iax misclassified as Iax (1 SNIa and 1 SNIbc). The more symmetric Figure 7 suggests that the asymmetry of this table is due to randomness and would be corrected with data from other runs.
The distribution of misclassified examples across the confidence spectrum is non-uniform. In this table, confidence refers to the probability assigned to the predicted type by the classifier. Confidence near 100% for a misclassified example is potentially more insightful than one near 50%. 9 out of the 21 misclassified heatmaps were misclassified at > 90% confidence, 6 at 90-70% confidence, and 6 at 70-50%. Surprisingly, the classifier is confidently wrong almost half the time. One particularly interesting example is a SNIbc "misclassified" as SLSN-1, but the classification probabilities for both SLSN-1 and SNIbc were 50%.

Limitations and Future Work
As stated in section 2.1, it is important to note that the metrics reported in this paper are in-distribution results since the training, validation, and test sets are mutually exclusive segments of the main dataset. The out-of-distribution performance of SCONE, as evaluated in section 3.1.2, is noticeably diminished from the > 99% in-distribution test accuracy. The high in-distribution test accuracy shows that SCONE is robust to previously unseen data, but the lower out-of-distribution test accuracy demonstrates SCONE's sensitivity to variations in the parameters of the dataset, such as the redshift distribution, relative rates of different types of SNe, small variations in the SN Ia model, as well as telescope characteristics. Generalizing SCONE to become robust to these variations will be the subject of a future paper.

CONCLUSIONS
In this paper we have presented SCONE, a novel application of deep learning to the photometric supernova classification problem. We have shown that SCONE has achieved unprecedented performance on the indistribution Ia vs. non-Ia classification problem and impressive performance on classifying SNe by type without the need for accurate redshift approximations or handcrafted features.
Using the wavelength-time flux and error heatmaps from the Gaussian process for image recognition also allows the convolutional neural network to learn about the development of the supernova over time in all filter bands simultaneously. This provides the network with far more information than a photograph taken at one moment in time. Our choice of an h × 3 convolutional kernel, where h is the number of wavelength bins, supplements these benefits by allowing the network to learn from data on the full spectrum of wavelengths in a sliding window of 3 days.
As future large-scale sky surveys continue to add to our ever-expanding transients library, we will need an accurate and computationally inexpensive photometric classification algorithm. Such a model can inform the best choice for allocation of our limited spectroscopic resources as well as allow researchers to further cosmological science using minimally contaminated SNIa datasets. SCONE can be trained on tens of thousands of lightcurves in minutes and confidently classify thousands of lightcurves every second at > 99% accuracy.
Although SCONE was formulated with supernovae in mind, it can easily be applied to classification problems with other transient sources.
The documented source code has been released on Github (github.com/helenqu/scone) to ensure reproducibility and encourage the discovery of new applications.