SuperRAENN: A Semisupervised Supernova Photometric Classification Pipeline Trained on Pan-STARRS1 Medium-Deep Survey Supernovae

Approximately 10% of the PS1-MDS transients were spectroscopically observed in real time throughout the survey, without a specific selection function (although brighter objects were more likely to be targeted). For this work, we limit our spectroscopic sample (557 objects) to five potential classes:

• 1.  
  Type I SLSNe (17 objects), which are thought to arise from the birth of highly magnetized neutron stars (Quimby et al. 2007; Chomiuk et al. 2011; Nicholl et al. 2017).
• 2.  
  SNe II (93 objects; including Type IIP and Type IIL SNe²² ), which arise from red supergiant progenitors.
• 3.  
  SNe IIn (24 objects), powered by the interaction the SN ejecta with preexisting circumstellar material (e.g., Smith et al. 2014).
• 4.  
  SNe Ia (404 objects), which are the thermonuclear explosions of white dwarfs.
• 5.  
  SNe Ibc (19 objects), which arise from the core collapse of massive stars that have lost their hydrogen (Ib) and helium (Ic) envelopes. Due to the small sample size, we consider SNe Ib and SNe Ic as a single class.

The SLSN and SN Ia light-curve samples have been previously published in Lunnan et al. (2018) and Jones et al. (2017), respectively. Model fits to the Type II light curves were presented in Sanders et al. (2015). For four objects, the transient spectra yield a reliable redshift but an ambiguous classification. A fifth object, PSc130816, has previously been identified as both an SN IIP/L (Sanders et al. 2015) and an SN IIn (Drout 2016). We do not include these five objects in our spectroscopic sample. An additional 15 objects are spectroscopically identified but do not fall in one of our five classes, including two tidal disruption events (TDEs), a lensed Type Ia, a Type Ibn, a Type Iax, and 10 fast-evolving luminous transients (FELTs). All except the TDEs are included in our photometric sample for training purposes, but not included in our spectroscopic sample. These objects are discussed in more detail in Section 5.1.

Our photometric sample contains 2315 objects with host galaxy spectroscopic redshifts, which are independent of the 557 SNe that are spectroscopically classified, as well as 13 SNe that do not fall in one of the five spectroscopic classes, discussed in Section 5.1. We refer to the union of the photometric and spectroscopic samples (the full set of 2885 events) as the "complete" photometric data set. We summarize the PS1-MDS SN-like objects, their associated hosts, and redshift information in Table 1. We also specify which SNe are used in the supervised/unsupervised portions of our classification algorithm.

Our spectroscopic data set is brighter than our complete photometric data set. As shown in Figure 1, the spectroscopic sample has a median peak r-band magnitude of ∼−21 mag, about 1 mag brighter than the photometric sample. We directly compare the redshift distributions in Figure 2. The spectroscopic sample peaks at a slightly lower redshift compared to the photometric data set (z ≈ 0.27 vs. z ≈ 0.35), with a tail extending to z ≈ 1.0. The lack of confident high-redshift measurements is likely due to the key spectroscopic lines shifting out of the optical and due to the peak absolute magnitudes of most SNe falling below our limiting magnitude. The mismatch between the spectroscopic and photometric samples may translate to biases in our classification pipeline, which we explore in more detail in Section 5. The complete griz light curves of our sample are available on Zenodo (Villar et al. 2020).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We explore the overall data quality of our sample in Figure 3, finding that the majority of events have ∼20 data points across all filters with S/Ns of ≳3. Given a typical SN duration of a month and our typical cadence of a few days, we expect the majority of (but not all) SNe to have fairly complete light curves.

Zoom In Zoom Out Reset image size

We generate and preprocess absolute magnitude light curves before extracting features. We correct each light curve for Milky Way reddening using the extinction map of Schlafly & Finkbeiner (2011). We estimate and normalize the absolute magnitude using the measured host redshift:

$$\begin{eqnarray}\begin{array}{rcl}{M}_{\mathrm{norm}} & = & m-5{\mathrm{log}}_{10}({d}_{{\rm{L}}}/10\mathrm{pc})\\ & & +\ 2.5{\mathrm{log}}_{10}(1+z)-{m}_{\mathrm{lim}}-{A}_{\lambda },\end{array}\end{eqnarray} \tag{ 1 }$$

where m_lim is a chosen limiting magnitude, which we take to be m_lim = 25. This value is dimmer than the 5σ limiting magnitude of PS1-MDS. We choose a dimmer magnitude to ensure that even marginal detections will be included in the light curve. We perform the renormalization so that the GP mean will be zero (i.e., the light curve will be zero when no data are available). Finally, we correct all light curves for time dilation based on the measured redshifts. We do not attempt to make a wavelength-dependent k-correction to the rest-frame data given the complicated, diverse, and time-evolving spectral energy distributions (SEDs) of the various SN types.

We do not correct the SN light curves for host galaxy reddening. The intrinsic reddening of SNe adds an additional scatter in our feature space. Correcting for host galaxy reddening would require estimating both the color excess and dust law, which is not possible given our current data set.

The PS1-MDS light curves are irregularly sampled across the four filters (see Section 2 for the PS1 observing strategy). The architecture of the RAENN does not require uniformly sampled light curves. However, it does require that each observation is made in all four filters. For example, if an observation is made in g band, we need to provide interpolated values for riz bands for that time.

To interpolate the griz light curves, we fit a GP using the open-source Python package George (Ambikasaran et al. 2015). GPs are a nonparametric model that has been previously used to interpolate and classify SN light curves (see, e.g., Lochner et al. 2016; Revsbech et al. 2018; Boone 2019). GPs define a prior over a family of functions, which is then conditioned on the light-curve observations. A key assumption is that the posterior distribution describing the light curve is Gaussian, described by a mean, μ(t), and a covariance matrix, Σ(t), given by Σ_i,j = κ(x_i, x_j) with kernel κ. We use a 2D squared exponential kernel to simultaneously fit all four filtered light curves:

$$\begin{eqnarray}&&\begin{array}{l}\kappa ({{\boldsymbol{t}}}_{i}{{\boldsymbol{t}}}_{j}{{\boldsymbol{f}}}_{i}{{\boldsymbol{f}}}_{j};\sigma ,{l}_{t}{l}_{f})={\sigma }^{2}\exp \ \left[-\displaystyle \frac{{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{l}_{t}^{2}}\right]\\ \qquad \ \times \ \exp \ \left[-\displaystyle \frac{{\left({f}_{i}-{f}_{j}\right)}^{2}}{2{l}_{f}^{2}}\right],\end{array}\end{eqnarray} \tag{ 2 }$$

where f is an integer between 0 and 3 that represents the griz filters, and the parameters l_t and l_f are characteristic correlation length scales in time and filter integer, respectively. This fitting process accounts for the measured data uncertainties, making it robust to low-confidence outliers.

We independently optimize the kernel parameters for each SN using the minimize function implemented in scipy, with initial values of l_t = 100 days and l_f = 1. We find that our choice of initialization values has little effect on the resulting best fit. We find that l_t is typically about 1 week, and l_f is typically 2–3, indicating that the filters are highly correlated. Examples of the GP interpolation for SNe Ia, SNe Ic, and SNe II are shown in Figure 4. The GP is able to produce reasonable interpolated light curves, even in cases with sparse and noisy data, and provide reasonable error estimates.

Zoom In Zoom Out Reset image size

A similar GP method was implemented by Boone (2019) to classify a variety of SN types in the Photometric LSST Astronomical Time-series Classification (PLAsTiCC; The PLAsTiCC team et al. 2018; Kessler et al. 2019) data set. Instead of an integer, Boone (2019) used the rest-frame central wavelength of each filter for each object. We avoid this added layer of complexity because the k-corrections and time-evolving SN SED change the weighted central filter wavelength. However, the simple 2D kernel still allows the four bands to share mutual information.

Our light curves contain several years of data, most of which are nondetections. To limit our input data, we keep data points (of any significance) within 100 days of peak flux (in whichever filter is brightest). For ease of optimization, the light curves need to contain the same number of data points. The data must be a consistent size during the back-propagation step of optimization for the RAENN for each iteration (see next section). Our longest light curve contains 169 data points, so we pad all light curves to match this length. We do so by appending a value dimmer than the estimated absolute limiting flux (we use m_lim = 25) to 100 days after the last detection in the light curve.

We note that using luminosity-based light curves (rather than magnitudes) is an alternative preprocessing choice. Luminosity-based light curves would remove the need to renormalize the light curves to a chosen limiting magnitude. We find that using luminosity-based light curves results in worse performance of the RAENN, likely due to the orders-of-magnitude differences in scale between events.

To extract unique features from the complete (unlabeled and labeled) PS1-MDS photometric sample, we construct an RAENN, inspired by the work of Naul et al. (2018), who uses a similar method to classify variable stars.

Neural networks are a class of machine learning algorithms that use many latent layers to model complex functions. These and other machine learning algorithms are becoming increasingly common in astronomy (see Ntampaka et al. 2019 for an overview). Autoencoders (AEs; Kramer 1991) are a class of neural network architectures that learn a compressed representation of input data. By training an AE to return the original data given a limited set of variables, it learns an "encoded" version of the data.

In astrophysics, AEs have been used for feature learning in galaxy SEDs (Frontera-Pons et al. 2017), image denoising (Lucas et al. 2018; Ma et al. 2019), and event classification (Naul et al. 2018; Pasquet et al. 2019). AEs are also increasingly being used in the astrophysics literature for dimensionality reduction (see, e.g., Ralph et al. 2019 and Portillo et al. 2020 for recent examples).

Here, our model is designed to address several concerns of SN light curves: (1) the temporal irregularity of data, (2) data across multiple filters, and (3) streaming data that update on a given cadence. The last point is not a concern for our PS1-MDS archival data set, but it will become important as LSST comes online and discovers thousands of SNe nightly.

The RAENN uses the GP light curves as input, by codifying the light curves as matrices of size 9 × T₀, where T₀ = 169, as described in the previous section. The nine values are as follows: one time value, relative to maximum (in whichever filter is brightest); four magnitude values (griz) at that time; and four magnitude uncertainties. Recall that the magnitude values are either measured or estimated from the GP. For the uncertainties, we use the 1σ errors for the measured points. For the GP points, we use a large error of 1 mag. We note that the GP produces estimate errors, but we find that, in practice, using this larger error bar leads to better performance. We leave exploration of utilizing error bars to future work. We emphasize that, while T₀ = 169 for training, the RAENN architecture allows a user to input a light curve of any size without needing to pad the light curve.

The RAENN architecture is divided into an encoder and a decoder. Our encoder is a series of fully connected layers that decrease in size until the final encoded layer with size N_E (i.e., the number of neurons used to fully encode the SN light curve). We note that N_E is a free parameter of our model that needs to be optimized. Similarly, the fully connected layer has N_N neurons, where N_N > N_E and is also a tunable parameter. Following the encoded layer, the decoder half of the architecture mirrors the encoder with increasing layer sizes.

A novel feature of our architecture is the inclusion of a repeat layer immediately after our final encoding layer (the layer of size N_E). In this layer, we repeat the encoded version of the light curve T_N times. To each copy, we append the time of each data point, relative to peak brightness in one filter (whichever filter happens to be brightest). One way to view the purpose of this layer is to imagine the autoencoder as two functions. The first function (the encoder) takes in the original data points, including observation times, and outputs a set of N_E values. This is similar to the idea of taking a light curve and fitting it to a model with N_E free parameters. A second function (the decoder) takes in a set of N_E values and T_N times to generate a light curve at the T_N times. This architecture allows us to generate a light curve at different T_N times, e.g., interpolated or extrapolated light curves, which is further explored in Section 5. In this work, we choose T_N = T₀, namely, we repeat the encoded values to match the original light-curve length.

Our autoencoder utilizes gated recurrent neurons (GRUs; Rumelhart et al. 1988; Cho et al. 2014). In addition to the typical hidden weights that are optimized during training, recurrent neurons have additional weights that act as "memory" of previous input. GRUs in particular utilize an update value (called a gate) and a reset gate. The values of these neurons determine how the current and previous inputs affect the value of the output. With each light-curve data point, the gates become updated with new information that informs the next prediction. This class of neurons is useful for our light curves with various numbers of observed data points. Our GRU neurons use the $\tanh$ activation functions with a hard sigmoid for the gate activation function.

Our RAENN is implemented in Keras (Chollet 2015) with a Tensorflow back end (Abadi et al. 2016). A diagram of the architecture is shown in Figure 5 and is outlined as follows:

• 1.  
  Input Layer: Input light curve of size T₀ × 9 with each griz data point labeled with a time (one value) in days relative to light-curve peak (four values) and an uncertainty (four values).
• 2.  
  Encoding Layer: Encoding layer with N_N neurons, where N_N is a hyperparameter.
• 3.  
  Encoded Layer: Encoded light curve with N_E neurons, where N_E is a hyperparameter.
• 4.  
  Repeat Layer: Layer to repeat encoded light curve to match with new time array, with size T₀ × N_E.
• 5.  
  Concatenate Layer: Layer to concatenate new times to encoded light curve, with size T₀ × (N_E + 1).
• 6.  
  Decoding Layer: Decoding layer with N_N neurons.
• 7.  
  Decoded Layer: T₀ × 4 decoded griz light curve.

Zoom In Zoom Out Reset image size

To optimize the free parameters (the weights) of the RAENN model, we must define a loss function. Our loss function is a simple mean square error function:

$$\begin{eqnarray}&&{ \mathcal L }=\displaystyle \sum _{i=0}^{N}\ \displaystyle \frac{{\left[{F}_{i,\mathrm{True}}(t,f)-{F}_{i,\mathrm{Predicted}}(t,f)\right]}^{2}}{N},\end{eqnarray} \tag{ 3 }$$

where F is the SN flux as a function of time t and filter f. Although we feed uncertainties into the network, we find that excluding flux errors in our loss function substantially improved the ability of the RAENN to match the input light curves. We minimize our loss function using the gradient-descent-based optimizer, Adam (Kingma & Ba 2014), finding an optimal learning rate of 10⁻⁴, which is a typical value.

We randomly split our unlabeled data set into training (2/3) and test (1/3) sets. We optimize the number of neurons in both the encoding and decoding layers (fixed to be the same number, N_N) and the number of encoding neurons (N_E) through a grid search, allowing N_N to vary from 20 to 160 in intervals of 20 and N_E to vary between 2 and 24 in intervals of 2. We find that, when optimizing over final classification F1-score (defined below), purity, and completeness, our results are relatively insensitive to N_E and N_N for values of N_E ∼ 10 and N_N ∼ 100. For our final model, we use N_E = 8 and N_N = 120, which is not our optimal model but a representative model. Utilizing our optimal model without creating a test set (in addition to a training and validation set) would likely overestimate performance. Given our limited data set, we are unable to properly optimize our hyperparameters and thus present representative results. Our best-performing model (N_E = 4 and N_N = 160) performs slightly better, with a macro-average completeness and purity of 71% and 69%, respectively. We note that N_N is slightly below the maximum number of data points in our set of light curves (where the longest light curve has 169 observed data points). The number of encoding neurons N_E is similar to the number of free parameters for the analytical model used in Villar et al. (2019) to capture the shape of a single-filter SN light curve.

We contrast our architecture with methods from Naul et al. (2018) and Pasquet et al. (2019), who present similar methodologies. Naul et al. (2018) use a similar GRU-based RAENN to classify variable stars with unevenly sampled light curves in one filter from the All Sky Automated Survey Catalog of Variable stars (Pojmanski 2002). The flux and time since last observation (Δt) are sequentially read into the recurrent layers. The same time array is fed into the decoder for output. In our case, we feed in a time series across four filters and give a time array relative to peak rather than relative to the previous data point. This is more natural in our problem, in which the SNe have a clear beginning and end, versus the periodic signals of variable stars. Additionally, our architecture allows us to give the decoder a different time series to allow for interpolation or extrapolation of the data.

Pasquet et al. (2019) use a semisupervised method to classify simulated SN light curves from the SPCC (Kessler et al. 2010). They use an AE with convolutional layers by transforming the light curves into "light-curve images" (see Pasquet-Itam & Pasquet 2018). Rather than interpolate the light curves, Pasquet et al. (2019) apply a mask to filters that are missing data at a certain time. In contrast, we interpolate our light curves but assign interpolated values a large uncertainty of 1 mag, as explained above. We found that the method of transforming light curves into images and masking across four filters led to unstable training and poorer performance. This is likely due to the large data gaps in the real PS1-MDS light curves, compared to the high-cadence (2 days for each filter) simulated light curves of SPCC. Since the LSST data are expected to more closely resemble the PS1-MDS light curves than the SPCC simulated events, we expect our method to be more robust in a real-life application.

As a final step, we use the encoded light curves as features for a supervised classification method. Following Villar et al. (2019), we train an RF classifier on the PS1-MDS spectroscopically classified SNe, including the RAENN encodings as features.

In addition to the encoding (eight features), we use the following 36 features based on the GP-interpolated light curves:

• 1.  
  The griz rise times in the rest frame, calculated 1, 2, and 3 mag below peak (12 features).
• 2.  
  The griz decline times in the rest frame, calculated 1, 2, and 3 mag below peak (12 features).
• 3.  
  The griz peak absolute magnitudes (four features).
• 4.  
  The median griz slope between 10 and 30 days post-peak in observer frame. This area was chosen by eye to specifically help the model differentiate between SNe II and SNe Ibc (four features).
• 5.  
  The integral of griz light curves (four features).

We measure these values from the GP-interpolated light curves rather than the decoded light curves. The decoded light curves are, at best, approximations of the GP-interpolated light curves. Therefore, using them would only result in noisier features. The decoded light curves are necessary, however, as a means to train the RAENN to extract the N_E encoding neuron values. We note that for some features, e.g., the rise and decline times, the feature values are heavily dependent on the GP extrapolation in cases where there are no measured data. Including GP errors in the supervised portion of our analysis could help capture this intrinsic uncertainty in the underlying light curve, but we leave that exploration to future work.

These features were chosen through trial and error while optimizing classification accuracy. We find that inclusion of all features leads to our optimal classification accuracy, although we do explore how well our classifier performs with the RAENN features alone in the following section.

We pass these features through an RF classifier, utilizing 350 trees in the RF and the Gini information criterion. The number of trees was determined based on trial-and-error optimization. To counteract the imbalance across the five spectroscopic classes, we tested several algorithms to generate synthetic data to augment our training set. Following Villar et al. (2019), we use a Synthetic Minority Over-sampling Technique (SMOTE; Chawla et al. 2002) and a multivariate Gaussian (MVG) fit. We additionally test using a kernel density estimate of the training set, using a Gaussian kernel with bandwidth equal to 0.2 (or 20% of the whitened feature standard deviation). We find that the MVG with a halved covariance matrix performs best. We test our classifier using leave-one-out cross-validation, in which we remove one SN from the sample, oversample the remaining objects by generating new objects using the MVG, and then apply the trained RF to the single, removed event and recording the result. For each object, our RF reports probabilities associated with each class, which are calculated using the fraction of trees that vote for each class. We take the class with the highest probability as the predicted SN type.

We visualize the completeness and purity of the spectroscopic sample using confusion matrices in Figure 6. A confusion matrix compares our RAENN label (horizontal axis, in the case of the completeness matrix) with the spectroscopic label (vertical). Results are shown for leave-one-out cross-validation, in which one event is removed from the sample for training and the trained model is applied to the left-out event. As with Villar et al. (2019), we find that our classifier performs best for SNe Ia (92% completeness), SLSNe (76%), and SNe II (82%) and worst for SNe Ibc (37%). Our class-averaged classification completeness is 69% across the five SN types. This is worse than the performance of Villar et al. (2019), who find a class-averaged completeness of 80%. Our class-averaged purity is 66%, again slightly worse than the average purity of 72% found in Villar et al. (2019). When limiting the sample to only objects in which the classification probability is ≥0.7 (a total of 438 objects), we find that our performance increases, with a class-averaged completeness of 75% and a class-averaged purity of 74% with a loss of 20% of the sample size.

Zoom In Zoom Out Reset image size

Next, we explore the classification confidence reported by our algorithm. The confidence estimates are directly outputted by the RF. With larger data sets, one can calibrate the outputted uncertainties using, e.g., an additional logistic function. Given our small data set, we do not perform any additional calibration. In Figure 7 we plot the cumulative fractions of SNe in our training set, grouped by their spectroscopic and photometric classifications. The majority of high-confidence objects are SNe Ia, with nearly half of the spectroscopic SNe Ia having a confidence (p > 0.98). Similarly, half of the SLSNe have high-confidence identifications (p > 0.8). SNe Ibc and SNe IIn have the lowest confidence on average, with the majority of events having p < 0.5. This is likely reflective of the fact that SNe Ibc and SNe IIn span a wide range of observed properties, including overlap with SNe Ia.

Zoom In Zoom Out Reset image size

Figure 7 also indicates the misclassified objects. Ideally, we want our misclassifications to largely occur in low-confidence objects. This is the case for SLSNe, SNe Ia, and SNe II. For SNe IIn and SNe Ibc, the misclassifications occur even for high-confidence events. This indicates that for SNe Ia, SNe II, and SLSNe, misclassifications are likely tied to data quality. In contrast, misclassifications of SNe IIn and SNe Ibc seem to be due to intrinsic overlap of the classes in feature space with other SNe (mainly SNe Ia).

We additionally attempt to sort events based on the number of data points, rather than classification confidence (see the right column of Figure 7). Our photometric data set has, on average, fewer >5σ data points compared to our spectroscopic data set (∼15 vs. ∼30 data points on average). Because of this mismatch and the lack of a strong correlation between number of points and classification confidence, we do not further explore how cutting sparse light curves affects our final classification accuracy.

We next turn our attention to the performance of our classifier when constrained to only data-driven (RAENN) features. Using the same set of RAENN features without any additional information, we produce the confusion matrices shown in Figure 8. We find a class-averaged completeness of 53% (and purity of 50%), approximately 20% worse than including the additional features. The overall breakdown is similar to our final confusion matrix, with the worst-performing classes being Type Ibc and Type IIn. We find that our RAENN-only classification is more inclined to label events as SNe Ia, likely a bias from the fact that our SN data set used to train the RAENN is highly dominated by SNe Ia. One could correct for this bias by adding some type of prior on the relative observational rates of the SN classes explored here. If we run our classification algorithm without the RAENN features, we find that SuperRAENN performs similarly (slightly worse, with a macro-averaged purity and completeness of 66% and 62%, respectively), implying that the RAENN has not picked up on uniquely helpful features independent from our hand-selected feature set. To be clear, the intent of RAENN is not necessarily to outperform hand-selected features but to create model-independent features in real time. In this work, we determine final classifications with the RAENN and hand-selected features to provide the highest-confidence photometric classifications. Improvements to classifications based solely on RAENN features are left to future work.

Zoom In Zoom Out Reset image size

While not optimized for Type Ia versus non-SN Ia classification, we explore how well our classifier (using the full set of features) performs when we collapse the confusion matrix into just two classes. In Figure 9, we show the completeness and purity confusion matrices for Type Ia versus core-collapse SN (CCSN) classifications, finding ≈90% completeness and >80% purity in both classes.

Zoom In Zoom Out Reset image size

The RF classifier allows us to measure the relative "importance" of the 44 features used to classify the SNe. We define importance as the decrease in the Gini impurity, which accounts for how often a feature is used to split a node and how often a node is reached in the forest (Breiman et al. 1984). We show the importance of each RF feature in Figure 10, along with the measured importance for a normal random variable. The peak magnitudes and decline rates are the most important features for classification. However, the RAENN features also have significant influence on the final classifications, with two RAENN features appearing in the top 10 important features.

Zoom In Zoom Out Reset image size

The feature importance unfortunately loses some quantitative meaning if the features are correlated, which is the case with our features. When two features are highly correlated, one may be arbitrarily measured as more important, so the general trends are more meaningful than precise order. We show the magnitudes of the feature correlations in Figure 11 to better understand the underlying correlations. There are clear correlations between features derived in multiple bands (e.g., the peak magnitude in g band is highly correlated to that in r band). However, we also see correlations between the RAENN features and the more traditional light-curve features. About half of the RAENN features seem strongly correlated with the peak magnitudes, while two others seem well correlated with rise and decline times. A more detailed exploration of the physical interpretation of the RAENN feature space may be worthwhile but is beyond the scope of this work.

Zoom In Zoom Out Reset image size

We apply our trained classification algorithm to the PS1-MDS data set of SN-like transients that pass our quality cuts described in Section 2. We report the probabilities of each class type for each light curve in Table 2. Error bars for each class probability are calculated by running the trained RF classifier 25 times with unique random seeds. We show the class breakdown of the complete photometric set (2885 SNe) in Figure 12. Excluding the spectroscopic sample, we present 2315 new SNe, with 1435 (61.9%) SNe Ia, 459 (19.9%) SNe IIP, 272 (11.7%) SNe Ibc, 112 (4.8%) SNe IIn, and 37 (1.6%) SLSNe. Of these, 1311 are high-confidence (p > 0.7) photometric classifications. A cumulative plot of the confidences grouped by each photometric class is shown in Figure 7; the distribution of these probabilities largely matches the spectroscopic sample.

Zoom In Zoom Out Reset image size

Table 2.  SN Classification

Event Name	pSLSN	pII	pIIn	pIa	pIbc
PSc000001	${0.00}_{0.00}^{0.00}$	${1.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$
PSc000006	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${1.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$
PSc000010	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${1.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$
PSc000011	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${1.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$
PSc000014	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${1.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$
PSc000034	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$	${1.00}_{0.00}^{0.00}$	${0.00}_{0.00}^{0.00}$

Note. Spectroscopically identified SNe have probabilities of 1.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

A sample of SNe from each photometric class is shown in Figure 13, including high- and low-probability examples. For the low-probability examples, it seems that even well-sampled light curves can have low confidence scores, likely because the features of their light curves reside on a region of feature space in which various SN classes reside. This is highlighted in Figure 14 in luminosity–decline space. Several misclassified events (shown as squares) have a high number of data points (represented in misclassified events as marker size).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The redshift distributions of the new, photometrically classified events are roughly consistent with those seen in the spectroscopic sample (Figure 2), with Type Ibc peaking at z ∼ 0.19, Type II peaking at z ∼ 0.21, Type Ia and Type IIn peaking at z ∼ 0.42, and SLSNe peaking at z ∼ 0.58.

We compare the overall photometric breakdown of SN types to that of the ZTF Bright Transient Survey (Fremling et al. 2020), which spectroscopically classified 761 SNe with peak g- or r-band magnitude of <18.5. Fremling et al. (2020) find that their magnitude-limited survey breaks down into 72% SNe Ia, 16% "normal" SNe II (Type IIP/L), 3% SNe IIn (including their Type IIn and SLSN-II category), 5% SNe Ibc, and 1.6% Type I SLSNe. This is a similar breakdown to that found in our spectroscopic sample. Comparing to our photometric set, we find a slightly higher fraction of SNe II and SNe IIn and a lower fraction of SNe Ia (all within ∼20% of the ZTF BTS values), as shown in Figure 12. For our high-confidence (p > 0.7) sample (also shown in Figure 12), our class breakdowns are closer to those of our spectroscopic and the ZTF BTS sample, with a slight overabundance of SNe Ia (≈78%). Based on our understanding of how our classifier performs on the training set, we can understand the biases present (e.g., that some spectroscopic SNe Ibc are classified photometrically as SNe Ia). We can use these known biases, encoded within the confusion matrices, to correct our class breakdown. Mathematically, this is calculated as the dot product of the purity matrix and our original class breakdown. Applying this correction to the photometric data set, the class breakdown is well aligned with the breakdown of our spectroscopic sample, as shown in Figure 12. This study should not be used to rigorously study the observational breakdown of SN classes; however, the fact that our p > 0.7 sample is in relatively good agreement with the ZTF BTS provides some evidence that our photometric sample is correctly labeled.

Our algorithm assumes that every SN belongs in one of five classes: SLSNe, SNe II, SNe IIn, SNe Ia, and SNe Ibc. Yet what does our algorithm do for transients that do not fall in these five classes? Here we address this question for a number of spectroscopically classified extragalactic transients. We summarize the photometric classification for these rare transients in Table 3.

Drout et al. (2014) presented a sample of 10 extragalactic transients discovered with PS1-MDS with redshift measurements that rise too rapidly to be powered solely with 56Ni.²⁴ Following Rest et al. (2018), we refer to these as FELTs. FELTs have a broad range of peak magnitude (−16 ≳ M ≳ −20), which is reflected in the distribution of photometric classifications. Of these 10 objects, six have "high-confidence" (p > 0.7) classifications in one of our five categories, four of which are SNe Ia and two of which are SNe II. The other four objects are classified as low-confidence Type Ia (one object), Type II (two objects), and Type Ibc (one object). As expected, the higher-luminosity objects are those classified as Type Ia, while the lower-luminosity objects are classified as Type II. The majority of objects have Type Ibc as their second-highest classification. Based on this analysis, FELTs are likely a (small) contaminant of both SNe II and SNe Ia in our sample, and our algorithm would need to be retrained to specifically classify FELTs.

Two known TDEs were discovered in PS1-MDS: PS1-10jh (PSc040777; Gezari et al. 2012) and PS1-11af (PSc120170; Chornock et al. 2014). Both objects are classified as SNe IIn with p ∼ 0.8 and p ∼ 0.6, respectively. This makes intuitive sense, as the light curves tend to be long-lived and bright like some SNe IIn. Both objects have Type Ia and Type II as their next most likely classifications. Based on these, it may be possible to search for TDEs in our sample within the photometric Type IIn sample.

We highlight four other SNe that do not fit in our five categories. PS1-10afx (PSc080333) is a lensed SN Ia (Chornock et al. 2013; Quimby et al. 2014), which peaks at −22 mag. We classify PS1-10afx as a high probability (p ∼ 0.9) SLSN. PS1-12sk (PSc370290) is an SN Ibn (Sanders et al. 2013) that peaks at M ∼ −19. We classify PS1-12sk as a low-probability Type Ia (p ∼ 0.6) or Type IIn (p ∼ 0.4). We classify PS1-12sz (PSc370330) as a likely SN IIb using SNID; PS1-12sz peaks at M ∼ −18.5. We photometrically classify this object as a low-probability Type Ibc (p ∼ 0.6). Finally, SN 2009ku (PS0910012) is a spectroscopically identified Type Iax (Narayan et al. 2011) that peaks at M ∼ −18.5. We classify this object as a low-probability Type Ia (p ∼ 0.5) or Type Ibc (p ∼ 0.3).

As discussed in Section 2, our spectroscopic sample is somewhat brighter and at a lower redshift than our test set. This difference may introduce biases in our final classifications, although this effect should be minimal considering the small (∼1 mag) difference between the two sets. De-redshifting the SNe removes some of this bias, by removing knowledge of the underlying redshift as a feature. The relative fractions of SN subtypes may evolve with redshift as host properties change (see, e.g., Graur et al. 2017 for an exploration of the correlations between host properties and SN type). Our spectroscopic and photometric sets differ most greatly at z ≳ 0.5 (see Figure 2). In this redshift range, average host metallicity is not expected to drastically shift (Lilly et al. 2003), implying a small potential bias. A separate bias may arise from the fact that our photometric sample relies on a measured spectroscopic redshift. At higher redshift, our galaxy redshift measurements become increasingly uncertain as dominant emission lines shift out of the optical band and intrinsically dim hosts fall below our observational limits. In contrast, rest-frame UV features of SNe (especially SLSNe) remain in the optical band, making it easier to confidently measure a distance from SN spectra. In the future, this problem can be mitigated with photometrically derived host galaxy redshifts.

As expected, the relative observed fraction of SN subtypes evolves with redshift owing to the magnitude limit of the survey. We trace this evolution in Figure 15. We show the cumulative fraction (integrating from z = 0) of each subclass as a function of redshift. Each subclass peaks in order of luminosity function. The dimmest subclass, SNe II, dominates the sample for z < 0.3, peaking near z ∼ 0.

Zoom In Zoom Out Reset image size

Using the high-redshift (z > 0.75) sample, we can test whether redshift information is playing an unwanted role in our training. The spectroscopic sample at z ≳ 0.75 is solely made up of SLSNe; however, we do not expect all high-z objects to be SLSNe. Given a typical limiting magnitude of m_r,lim ∼ 23.3, the corresponding absolute magnitude is ∼−20 at z = 0.75. At this sensitivity, we expect to find SLSNe, SNe IIn, and potentially bright SNe Ia (if the limiting magnitude is slightly deeper). For z > 0.75, we find that our photometric sample (a total of 28 SNe) is 68% SLSNe, 18% SNe IIn, and 14% SNe Ia (with all SNe Ia occurring at z < 0.85), implying that our classifier has not learned to simply classify all high-z events as SLSNe. The high-z SNe Ia, in particular, have noisy light curves that peak at M ∼ −20.

We first compare our results to H20, which extends the work of Villar et al. (2019) to classify the PS1-MDS photometric sample using features extracted from analytical fits to the light curves. Overall, H20 (and Villar et al. 2019) achieve better performance at the cost of a more computationally expensive feature extraction method. While our algorithm takes roughly seconds to run on new events after training, Superphot takes roughly several minutes to fit an analytical model to the data set. We agree with 74% of the photometric classifications of H20. If we compare the top two labels, the algorithms agree on 95% of classifications. Indeed, often the top two classification choices are flipped for either algorithm, occurring most often with SNe II/Ibc and SNe IIn/Ia. We find stronger agreement if we exclude objects with low classification confidence, namely, using only p > 0.7 in both algorithms, our classifications agree 84% of the time (with 1597 objects remaining after the cut, i.e., a loss of ∼50% of the sample). The agreement increases further for even higher probability cuts of p > 0.8 (>0.9), with 88% (92%) agreement with 1249 (888) objects remaining. Most classification disagreements lie in the Type Ibc/IIn categories, which have low-confidence classifications. We find that our algorithm is more likely to classify SNe as Type Ia, likely a bias built into the unsupervised step of training on the complete data set (which is dominated by SNe Ia). A more detailed comparison of these two results is provided in H20.

Villar et al. (2019) discuss the difficulty in comparing our results to the broader literature. In short, previous works have largely focused on Type Ia versus CCSN classification (e.g., Ishida & de Souza 2013; Jones et al. 2017; Brunel et al. 2019) or have been trained and tested on simulated data (e.g., Kessler et al. 2010; Muthukrishna et al. 2019a; Möller & de Boissière 2020). In the case of Type Ia versus CCSN classification, we achieve an accuracy of ≈92%, similar to (but somewhat worse than) specialized classifiers (Jones et al. 2017; Brunel et al. 2019). When comparing to works based on simulated data, we caution that not all simulated data sets are suitable for multiclass SN classification. In particular, the SNPCC training set (Kessler et al. 2010) lacks the SN diversity necessary to accurately train classifiers and will lead to artificially promising results. PLAsTiCC (The PLAsTiCC team et al. 2018; Kessler et al. 2019) is better suited for this task, and we encourage future work to be built on this data set or the PS1-MDS data set presented here.

We next compare our results to Jones et al. (2017), who presented a PS1-MDS sample of 1169 likely SNe Ia, focusing on Type Ia versus non-Ia classification. Jones et al. (2017) used four classification algorithms: the template-matching algorithm PSNID (Sako et al. 2011), a nearest-neighbor approach using the PSNID templates, an algorithm based on fitting light curves to SALT2 templates (Guy et al. 2007), and a method, GALSNID (Foley & Mandel 2013), that only utilizes host galaxy properties. Jones et al. (2017) similarly removed objects with unreliable host redshifts and potential active galactic nucleus hosts, but unlike our analysis, they removed objects at z > 0.75. Of their 1169 identified SNe Ia, only 1046 SNe pass our quality cuts to be classified in this work. For these, we find 95% agreement. Of the remaining 48 SNe, we identified Type Ia as the second-highest choice in 24 cases. Of the remaining 24 cases, 15 have low Type Ia probabilities (p < 0.8 from Jones et al. 2017) or classification probabilities based entirely on host galaxy. It is worth noting that our classifier, similar to Jones et al. (2017), achieves 96% purity in SNe Ia, making it likely usable for cosmological studies (Jones et al. 2018).

We compare our results to those trained on PLAsTiCC—in particular, Boone (2019), Muthukrishna et al. (2019a), and Gabruseva et al. (2020). These classifiers present average completenesses of ≈0.88 for SLSNe (higher than our score), ≈0.5 for SNe II/IIn (lower than our averaged Type II/IIn score), ≈0.92 for SNe Ia (similar to our score), and ≈0.46 for SNe Ibc (similar to our score given low-number statistics). These results are based on simulated data that lack the complexity of real data, so it is encouraging that our algorithm performs similarly or outperforms these works. It would be interesting and useful to the community to know how these algorithms perform on the PS1-MDS data set, but we leave this for future work.

We now turn to the architecture of the RAENN itself and its use in future surveys. The recurrent neurons allow our neural network to generate light-curve features that can be updated in real time, in addition to extrapolating and interpolating light curves. We highlight the accuracy of the RAENN light-curve model as a function of light-curve completeness in Figure 16, parameterized by the rms error. In particular, we train the RAENN model on the available observations up to a certain epoch and then compare the complete model estimate to the complete light curve. We find that SuperRAENN performance drastically improves post-peak, but that it can provide accurate light-curve estimates somewhat before peak. We leave the optimization of SuperRAENN to classify SNe in real-time data streams to future work. To explore why SuperRAENN improves near peak, we track how the RAENN features change as the light curves evolve. In Figure 17, we plot the values of representative encoding values of an SN Ia; we find similar behavior for other types of SNe. The encodings vary smoothly until settling on the final values ∼10 days post-peak. For CCSNe (Type II and Type Ibc), the encoding values also vary smoothly but trend more slowly to the final values.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The ability of the RAENN to extrapolate light curves without built-in physical assumptions allows it to search for anomalous events in real time for the purpose of spectroscopic and multiwavelength follow-up. Given the millions of events expected from LSST, it is essential to search for unexpected or previously unknown physical effects. One concern is that our algorithm is potentially not robust to noisy live-streaming data; in other words, our algorithm must be able to distinguish between anomalous data and noisy data. We check the stability of our encoded values as a function of scaled white noise by adding white noise to a light curve. We then use our RAENN to encode the noisy light curve and record the scatter of the encoded values. We report the results of this test in Figure 18, in which we show the scaled scatter of the encoded values as a function of the magnitude of the injected noise. The scatter grows linearly with noise; however, even with 1 mag of scatter added to the light curve, the overall scatter of the encoded values only increases to 30% of the overall spread of class's features. This implies that the RAENN is largely robust to noise.

Zoom In Zoom Out Reset image size

Several steps need to be taken to allow our architecture to work on streaming data. First, we use phases relative to maximum light, which will be unavailable during the rise of the SN. A shift to a time measurable early in the light curve, like time of first detection, will allow the RAENN to otherwise perform as designed. Similarly, the features utilized during the supervised portion of our classifier rely on the full light curve being available. All features can be estimated from extrapolated RAENN light curves, or a new set of features may be used on streaming data. Finally, although not necessary, our RAENN could output uncertainties on the SN light curves by converting the network into a variational AE, which is designed to simultaneously find an encoding space and uncertainties on the encoded data. This more complex architecture would likely require a larger training set to be reliable. Finally, we note that an algorithm like RAENN could be used in conjunction with an active oracle (a software that recommends new observations to improve classification) such as REFITT (Sravan et al. 2020), in order to actively optimize classification accuracy.

Finally, we note that we build the architecture and SuperRAENN pipeline assuming that measured redshift will be available. In practice, LSST will build an archive of redshift estimates as template images become more sensitive and photo-z algorithms become more accurate. We test how reliant our method is on redshift estimates by training a model (with the same hyerparameters as our final pipeline) in which all estimates are assumed to be equal (equivalent to working in flux space rather than luminosity space). We find that our algorithm performs, unsurprisingly, worse without redshift information, obtaining a macro-averaged completeness of 0.50 and purity of 0.49 across the five classes. The highest-performing classes are Type Ia and Type Ibc (which are similar to our redshift-dependent pipeline). The pipeline performs significantly worse on SLSNe and SNe II without redshift information.