PELICAN: deeP architecturE for the LIght Curve ANalysis

TL;DR

PELICAN is a deep learning architecture designed for light curve analysis, capable of classifying supernovae with high accuracy despite data sparsity and non-representative training sets, applicable to large astronomical surveys.

Contribution

This paper introduces PELICAN, a novel deep learning model that effectively classifies light curves without additional features, handling irregular sampling and non-representative training data.

Findings

Achieved 81.1% accuracy on Supernova Photometric Classification Challenge.

Detected 87.4% of supernovae Ia with over 98% precision in simulated LSST data.

Reached 86.8% accuracy on real SDSS data using simulated training data.

Abstract

We developed a deeP architecturE for the LIght Curve ANalysis (PELICAN) for the characterization and the classification of light curves. It takes light curves as input, without any additional features. PELICAN can deal with the sparsity and the irregular sampling of light curves. It is designed to remove the problem of non-representativeness between the training and test databases coming from the limitations of the spectroscopic follow-up. We applied our methodology on different supernovae light curve databases. First, we evaluated PELICAN on the Supernova Photometric Classification Challenge for which we obtained the best performance ever achieved with a non-representative training database, by reaching an accuracy of 0.811. Then we tested PELICAN on simulated light curves of the LSST Deep Fields for which PELICAN is able to detect 87.4% of supernovae Ia with a precision higher than 98%, by considering a non-representative training database of 2k light curves. PELICAN can be trained on light curves of LSST Deep Fields to classify light curves of LSST main survey, that have a lower sampling rate and are more noisy. In this scenario, it reaches an accuracy of 96.5% with a training database of 2k light curves of the Deep Fields. It constitutes a pivotal result as type Ia supernovae candidates from the main survey might then be used to increase the statistics without additional spectroscopic follow-up. Finally we evaluated PELICAN on real data from the Sloan Digital Sky Survey. PELICAN reaches an accuracy of 86.8% with a training database composed of simulated data and a fraction of 10% of real data. The ability of PELICAN to deal with the different causes of non-representativeness between the training and test databases, and its robustness against survey properties and observational conditions, put it on the forefront of the light curves classification tools for the LSST era.

Figures

Tables

Table 1. Table 1 : Statistics obtained for SPCC challenge by PELICAN and BDTs results in parenthesis. The first part reports results for a non-representative training database and the second part for a representative training database. We consider both cases by adding or not the redshift values inside the training.

SPCC	Redshift	Training database	Accuracy	AUC
Non-representative training database (SALT2+BDTs)	no	1103   spec	0.856 (0.705)	0.934 (0.818)
	yes	1103     spec	0.863 (0.713)	0.939 (0.855)
Representative training database (SALT2+BDTs)	no	1103
mix of spec and phot	0.911 (0.843)	0.970 (0.905)
	yes	1103
mix of spec and phot	0.917 (0.878)	0.971 (0.948)

Table 2. Table 2 : Statistics for various training configurations on the DDF survey (first part) and the WFD survey (second part), with BDTs results in parenthesis. The different metrics are defined in Section 7.1 .

	Total (LC)	Training database (Spec only)	test database (Phot only)	Accuracy	${\text{Recall}}_{Ia}$ ${\text{Precision}}_{Ia}>0.95$	${\text{Recall}}_{Ia}$ ${\text{Precision}}_{Ia}>0.98$	AUC
LSST DDF	2,000	$500$	$1,500$	0.895 (0.795)	0.769 (0.382)	0.602 (0.183)	0.966 (0.885)
		$1,000$	$1,000$	0.917 (0.888)	0.827 (0.505)	0.687 (0.353)	0.973 (0.898)
	4,000	$1,000$	$3,000$	0.928 (0.850)	0.890 (0.490)	0.764 (0.246)	0.979 (0.895)
		$2,000$	$2,000$	0.938 (0.813)	0.927 (0.594)	0.806 (0.256)	0.984 (0.905)
	10,000	$2,000$	$8,000$	0.946 (0.796)	0.944 (0.496)	0.886 (0.284)	0.989 (0.891)
		$3,000$	$7,000$	0.950 (0.809)	0.950 (0.548)	0.903 (0.285)	0.990 (0.905)
		$5,000$	$5,000$	0.971 (0.818)	0.981 (0.510)	0.959 (0.315)	0.996 (0.910)
	24,000	$2,000$	$22,000$	0.942 (0.792)	0.940 (0.477)	0.874 (0.209)	0.986 (0.890)
		$3,000$	$21,000$	0.945 (0.797)	0.937 (0.474)	0.891 (0.254)	0.986 (0.892)
		$5,000$	$19,000$	0.968 (0.805)	0.978 (0.485)	0.957 (0.228)	0.996 (0.898)
		$10,000$	$14,000$	0.971 (0.790)	0.983 (0.465)	0.965 (0.260)	0.997 (0.888)
LSST WFD	17,000	DDF Spec : $2,000$	WFD : $15,000$	0.965 (0.620)	0.982 (0.041)	0.905 (0.008)	0.992 (0.703)
	43,000	DDF Spec : $3,000$	WFD : $40,000$	0.976 (0.623)	0.995 (0.018)	0.964 (0.000)	0.996 (0.711)
	90,000	DDF Spec : $10,000$	WFD : $80,000$	0.978 (0.620)	0.995 (0.046)	0.973 (0.000)	0.997 (0.709)

Table 3. Table 3 : Statistics obtained on real SDSS data. The first line reports results obtained by training PELICAN only on simulated data. By including only 80 SDSS light curves in the training results are significantly improved (second line).

Training database	test database	Accuracy	AUC
SDSS simulations : 219,362	SDSS-II SN confirmed : 582	0.462	0.722
SDSS simulations : 219,362 SDSS-II SN confirmed : 80	SDSS-II SN confirmed : 582	0.868	0.850

Equations

$$L_{auto}=||X-g(f(X))||,$$

$$L_{denoising}=||X-g(f(\tilde{X}))||,$$

$$KL(\rho ||\widehat{{\rho }_j})=\rho \text{log}\frac{\rho }{\widehat{{\rho }_j}}+(1-\rho )\text{log}\left(\frac{1-\rho }{1-\widehat{{\rho }_j}}\right).$$

$$L_{sparse}=||X-g(f(X))||+\alpha {\Omega }_L,$$

$${\Omega }_L=\sum _{j}KL(\rho ||\widehat{{\rho }_j})=\sum _j\rho \text{log}\frac{\rho }{\widehat{{\rho }_j}}+(1-\rho )\text{log}\frac{1-\rho }{1-\widehat{{\rho }_j}},$$

$$L(D,(Y,\text{@vec}{Z}_1,\text{@vec}{Z}_2))=(1-Y)L_S(D)+Y{L}_D(D)$$

$$L_{contrastive}=(1-Y)||\text{@vec}{Z}_1-\text{@vec}{Z}_2||+Y\text{max}{\left(0,m-\sqrt{||\text{@vec}{Z}_1-\text{@vec}{Z}_2||}\right)}^2$$

$$L_{auto}=||X-g(f(X))\odot M(X)||$$

$${\rho }^{\prime }(l)={\rho }_a{n}_l+{\rho }_b$$

$${\Omega }_L^{\prime }(l)=\sum _j{\rho }^{\prime }(l)\text{log}\frac{{\rho }^{\prime }(l)}{\widehat{{\rho }_j}}+(1-{\rho }^{\prime }(l))\text{log}\frac{1-{\rho }^{\prime }(l)}{1-\widehat{{\rho }_j}}$$

$$L=$$

$$\frac{1}{2}\text{max}{\left(0,m-\sqrt{||LC2-LC3||}\right)}^2+\frac{1}{2}\text{max}{\left(0,m-\sqrt{||LC2-LC4||}\right)}^2+$$

$$||LC1-LC2||+||LC3-LC4||$$

$$Acc=\frac{TP+TN}{(TP+FP+TN+FN)}$$

$$R(\text{or }TPR)=\frac{TP}{(TP+FN)}$$

$$P=\frac{TP}{(TP+FP)}$$

$$FPR=\frac{FP}{(FP+TN)}$$

Full text

1 Introduction

A major challenge in cosmology is to understand the observed acceleration of the expansion of the universe. A direct and very powerful method to measure this acceleration is to use a class of objects, called standard candles due to their constant intrinsic brightness, which are used to measure luminosity distances. Type Ia supernovae (SNe Ia), a violent endpoint of stellar evolution, is a very good example of such a class of objects as they are considered as standardizable candles. The acceleration of the expansion of the universe was derived from observations of several tens of such supernovae at low and high redshift (Perlmutter et al., 1999; Riess et al., 1998). Then, several dedicated SN Ia surveys have together measured light curves for over a thousand SNe Ia, confirming the evidence for acceleration expansion (e.g. Betoule et al., 2014; Scolnic et al., 2018).

The future Large Survey Synoptic Telescope (LSST, LSST Science Collaboration et al., 2009) will improve on past surveys by observing a much higher number of supernovae. By increasing statistics by at least an order of magnitude and controlling systematic errors, it will be possible to pave the way for advances in precision cosmology with supernovae.

A key element for such analysis is the identification of type Ia supernova. But the spectroscopic follow-up will be limited and LSST will discover more supernovae than can be spectroscopically confirmed. Therefore an effective automatic classification tool, based on photometric information, has to be developed to distinguish between the different types of supernovae with a minimum contamination rate to avoid bias in the cosmology study. This issue was raised before and has led to the launch of the Supernova Photometric Classification Challenge in 2010 (SPCC, Kessler et al.) to the astrophysical community. Several classification algorithms were proposed with different techniques resulting in similar performance without resolving the problem of non-representativeness between the training and test databases. Nonetheless, the method developed by Sako et al. (2008, 2018) based on template fitting, shows the highest average figure of merit on a representative training database, with an efficiency of 0.96 and an SN Ia purity of 0.79.

Since then, several machine learning methods were applied to classify supernovae light curves (e.g. Richards et al., 2012; Ishida & de Souza, 2013; Karpenka et al., 2013; Varughese et al., 2015; Möller et al., 2016; Lochner et al., 2016; Dai et al., 2018). They showed interesting results when they are applied on a representative training dataset but the performance dramatically decreases when the learning stage is made on a non-representative training subset, which represents however the real scenario.

We propose to explore in this paper a new branch of machine learning, called deep learning, proved to be very efficient for image and time series classification (e.g. Szegedy et al., 2015; He et al., 2016; Schmidhuber et al., 2005). One of the main difference with the classical machine learning methods is that the raw data are directly transmitted to the algorithm that extracts by itself the best feature representation for a given problem. In the field of astrophysics, deep learning methods have shown better results than the state of the art applied to images for the classification of galaxy morphologies (Domínguez Sánchez et al., 2018), the classification of transients (du Buisson et al., 2015; Gieseke et al., 2017) and the estimation of photometric redshifts (Pasquet et al., 2018) to name a few. This method have also showed impressive performance for the classification of light curves (Mahabal et al., 2017; Pasquet-Itam & Pasquet, 2018) and especially the classification of supernovae (Charnock & Moss, 2017; Brunel et al., 2019).

In this work we develop a complex Deep Learning architecture to classify light curves. We apply our study to the classification of light curves of supernovae. Unlike the other studies, our method overcome the problem of non-representativeness between the training and the test databases, while considering a small training database. We apply our method on the SPCC challenge, then on LSST simulated data including a biased and small training database. We also validate our method on real data coming from the Sloan Digital Sky Survey (SDSS) data. The paper is organized as follows. In Section 2, we explain the different issues for the classification of light curves. In section 3, we introduce deep learning concepts that we used and developed in this work. In Section 4, we present our architecture named PELICAN (deeP architecturE for the LIght Curve ANalysis). In section 5, we describe the different datasets used in this study. In Section 6 we present the experimental protocol that we adapted to make PELICAN robust against the differences of sampling and noise. In Section 7 we present our results for different databases. In Section 8 we analyze the behaviour of PELICAN with respect to a number of light curve properties and observational conditions. Finally we conclude and expose perspectives in Section 9.

2 Light curve classification issues

Light curves are fundamental signals to measure the variability of astrophysical objects. They represent the flux of an object along time in different photometric bands (e.g. ugriz system). Due to the observational strategy and conditions, light curves have an irregular sampling, often sparse. Therefore a sampling with two different observational cadences present several differences. Figure 1 shows, as an example, two simulated light curves with two different cadence models (see Section 5.3). Compared to an image, such light curves have incomplete and non-continuous information, thus imposing dedicated training algorithms.

The other issue is the non-representativeness between the training and the test databases. As the spectroscopic follow-up, used to label the lightcurves, is limited, the coverage of the training database in brightness, redshift and number is different from the test database as shown on Fig. 2. Moreover the absolute magnitude of SN Ia is correlated with two quantities. First, brighter SN Ia have wider, slower declining light curves. This variability can be described as a timescale stretch of the light curve (Phillips, 1993). In addition brighter SN Ia are bluer and a color correction has to be applied to standardize them (van den Bergh, 1995; Tripp, 1998). So due to these correlations, the lower-stretch and redder supernovae are fainter and tend to have small recovery efficiency (Malmquist bias) and so are under-represented in the training database which is limited in brightness. The non-representativeness of the databases, which is a problem of mismatch, is critical for machine learning process.

In general, machine learning methods require a sufficient large number of training data in order to correctly classify, so the small size of the training database involves another difficulty.

To provide a solution for each of these issues, we have designed a specific architecture. First, light curves from the test database are trained with a non-supervised model without using the knowledge of labels. This allows to reduce the mismatch between the training and the test databases and provides a solution to the small training dataset, by extracting features from the larger test database.

To reduce again the problem of non-representativeness we performed a second training step to minimize the distances in the feature representation space between bright and faint supernovae of same labels and maximize distances of supernovae with different labels.

Finally we integrated a regularization term into the training to adapt the model to the sparsity of data.

The resulting deep architecture, dedicated to the characterization and classification of light curves, is presented in Section 4.

3 Deep Learning Model

In this section, we present the main deep learning concepts which are used to build our network architecture. Namely the convolution layer network which describes the basics, the autoencoder which is a non-supervised module where we detail the notion of loss function and finally the contrastive approach where an adapted loss function is defined to improve the clustering of same label entries.

3.1 Convolutional Neural Network

Convolutional Neural Network (CNN) is a special type of multilayered neural network that is made up of neurons that have learnable weights and biases. The architecture of a CNN is designed to take advantage of the 2D structure of an input image. It takes as input a $h\times w\times c$ image where $h$ is the height and $w$ is the width of the image and $c$ is the number of channels. As a light curve is a 1D signal, we transformed it into a 2D ”light curve image” (LCI) as we did in Pasquet-Itam & Pasquet (2018). The width of the LCI is the temporal axis (expressed in days) and the height is the number of photometric bands. For example, if we consider a light curve measured in ugriz bands in 200 days, the corresponding LCI is a 2D array of dimension ($5\times 200$) pixels. As a light curve is not a continuous signal, the corresponding array is composed of many blank cells that we filled with zero values. A harmful consequence is the overfitting of the position of missing data, which could dramatically degrade the performance. In this case, the model learns the exact position of missing data of the training light curves and is not able to generalize to the unseen test data. To prevent this overfitting and make the model invariant against the position of missing data we have integrated several techniques that are explained in Section 6.

3.1.1 Convolution layers

In a convolution layer, each neuron applies a convolution operation to the input data using a 2D map of weights used as a kernel. Then resulting convolved images are summed, a bias is added and a non-linearity function is applied to form a new image called a feature map. In the first convolution layer, the convolution operation is realized between the input LCI and the set of convolution kernels to form feature maps that are then convolved with convolution kernels in the next convolution layer. For the non-linearity function, we mainly use the most commonly used activation function : the ReLU (Rectified Linear Unit, Nair & Hinton, 2010) defined by $f(x)=max(x,0)$. The weights of the kernels are updated during the training by back-propagation process.

3.1.2 Pooling layers

The network can be composed of pooling layers which quantify the information while reducing the data volume. The two most used methods consist in selecting only the maximum or the average value of the data in a local region.

3.1.3 Fully connected layers

Finally, fully connected layers are composed of neurons that are connected to every neuron of the previous layer and perform the classification process with features from previous convolution layers.

More details on CNN models can be found in Pasquet-Itam & Pasquet (2018); Pasquet et al. (2018).

3.2 Autoencoder

To benefit from the information of the light curves of the test database and so reduce the mismatch between the training and test databases, we have adapted a non-supervised autoencoder. An autoencoder is an unsupervised learning algorithm that tries to learn an approximation to the identity function such as output should mimic the input. As a consequence the internal layers exhibit a good representation of the input data. The input $X\in\mathbb{R}^{D}$, with $D=h\times w$ is transformed into an embedding $Z\in\mathbb{R}^{K}$, often such as $K<<D$. The mapping from $X$ to $Z$ is made by the encoder, noted $f$, that could perform a dimension reduction to finally get a good representation of the data in a compressed format. The reconstruction of the original signal, $X^{\prime}\in\mathbb{R}^{D}$ , is obtained by the decoder, noted $g$, that uses the compressed embedding representation $Z$ (see Fig. 3). The objective of the autoencoder is to minimize the distance function (for example L2 distance), named loss function, between each input $X$ and each output $X^{\prime}$. The learning process of the autoencoder consists in iteratively refining its internal parameters such that the evaluation of the loss function, on all the learning set, is reduced. The loss function associated to the autoencoder, noted $L_{auto}$ is defined as follow:

[TABLE]

with $X$ represents the input signal and ——.—— symbolizes the L2 distance.

The minimization of the loss function that maps from $X$ to $Z$ in order to obtain $X^{\prime}$ does not guarantee the extraction of useful features. Indeed the network can achieve a perfect reconstruction by simply ”copying” the input data and thus obtain a minimal mapping error. Without any other constraints, the network can miss a good representation of the input data. A strategy to avoid this problem is to constrain the reconstruction criterion by cleaning or denoising partially corrupted input data with a denoising autoencoder (Vincent et al., 2008).

The methodology consists to first apply noise, for example an additive Gaussian noise, on input data to corrupt the initial input $X$ into $\widetilde{X}$. Then the autoencoder maps $\widetilde{X}$ to $Z$ via the encoder $f$ and attempt to reconstruct $X$ via the decoder $g$ (see Fig. 4). Although $Z$ is now obtained by applying the encoder on corrupted data $\widetilde{X}$, the autoencoder is still minimizing the reconstruction loss between a clean $X$ and its reconstruction from $\widetilde{X}$ with the following loss function:

[TABLE]

where $\widetilde{X}=X+\epsilon$, with $\epsilon$ an additive uniform noise with the same dimension of $X$ (see Fig. 4).

As the number of parameters of the autoencoder is high, light curves are sparse and the size of the training database that we will have is small, the overfitting can seriously affect the performance of the network in our case. A solution is to introduce sparsity to the learned weights to avoid learning ”noisy” patterns. The sparsity concept does the assumption that only a few neurons are required in order to reconstruct a complete signal. Let us consider that a neuron is active if its output value is close to 1 and inactive if its output value is close to 0. We want a very sparse representation such that the neurons of a given layer should be inactive most of the time and so obtain an average activation of neurons close to 0. Therefore the output signal is reconstructed with a very limited number of activated neurons.

We note $\hat{\rho_{j}}$ the average activation of hidden unit $j$ over the training set. To force neurons to be inactive, we enforce the constraint: $\hat{\rho_{j}}=\rho$, with $\rho$ a sparsity parameter that is initialized close to zero. Thus the average activation of each hidden neuron has to be close to a small value and so most of the neurons have to be inactive to satisfy this constrain.

In practice, a penalty term is added to the loss function to penalize $\hat{\rho_{j}}=\rho$ deviating significantly from $\rho$. One of the main frequently used penalty term is the Kullback-Leibler (KL) divergence (Hinton, 2002) defined as:

[TABLE]

KL divergence measure how two distributions are different from one another. It has the property that $KL(\rho||\hat{\rho_{j}})=0$ if $\hat{\rho_{j}}=\rho$ and otherwise it increases monotonically as $\hat{\rho_{j}}$ diverges from $\rho$.

The loss function is now integrating this penalty term as:

[TABLE]

with $\alpha\in\mathbb{R}$ a scalar that weights the sparse regularization term $\Omega_{L}$ which is the sum of the KL divergence term for neurons of a chosen layer numbered $L$ defined as:

[TABLE]

$j\in\{1,...,N_{j}\}$ with $N_{j}$ the number of neurons of layer $L$.

3.3 The contrastive loss function

The contrastive loss function was introduced to perform a dimensionality reduction by ensuring that semantically similar examples are embedded close together (Hadsell et al., 2006). It was shown that this method provides invariance to certain transformations on images. The contrastive loss function is computed over pairs of samples unlike traditional loss functions which is a sum over all the training database. Let ($\@vec{Z_{1}}$, $\@vec{Z_{2}}$) a pair of input data and $Y$ a binary label assigned to this pair. If $\@vec{Z_{1}}$ and $\@vec{Z_{2}}$ have the same label, $Y=0$, otherwise $Y=1$. The distance function between $\@vec{Z_{1}}$ and $\@vec{Z_{2}}$ is learned as the euclidean distance: $D=||\@vec{Z_{1}}-\@vec{Z_{2}}||$. Thus the loss function tends to maximize the distance $D$ if they have dissimilar labels and minimize $D$ if they have similar labels. So we can write the loss function as:

[TABLE]

with $(Y,\@vec{Z_{1}},\@vec{Z_{2}})$ a labeled sample pair and $L_{S}$ the partial loss function for a pair of similar labels, $L_{D}$ the partial loss function for a pair of dissimilar labels. To get low values of $D$ for similar pairs and high values of D for dissimilar pairs, $L_{S}$ and $L_{D}$ must be designed to minimize $L$. We introduce the margin $m>0$ which defines a minimum distance between $(\@vec{Z_{1}},\@vec{Z_{2}})$. Dissimilar pairs contribute to the loss function only if their distance is below this minimum distance so that pairs who share a same label will be bring closer, and those who does not share a same label will be drive away if their distance is less than $m$. The final contrastive loss function is defined as:

[TABLE]

4 The proposed architecture

We developed the method named deeP architecturE for the LIght Curve ANalysis (PELICAN) to obtain the best feature-space representation from light curves and perform a classification task. In this work, we apply PELICAN for the classification of supernovae light curves but it can be extended to the classification of other variable or transient astrophysical objects.

PELICAN is composed of three successive modules (see Fig. 7). Each of them has a specific purpose with a loss function associated. The first module learns a deep representation of light curves from the test database under an unsupervised autoencoder method. The second module optimizes a contrastive loss function to learn invariance features between the bright and fainter supernovae from the training database. Finally, the third module performs the classification task

In this section we explain in more details the different mechanisms and objectives of the operations related to each module.

4.1 The autoencoder branch

To deal with the low number of examples in the training database that leads to overfitting and mismatch between the spectroscopic and photometric distributions (see Fig. 2), we propose to train an unsupervised sparse autoencoder method on the test database. In this way we can benefit from information of light curves in the test database without knowing the label associated to each object.

The autoencoder takes as input a batch of LCIs of size $h\times w$ from the test database, that are encoded and decoded through a convolutional neural network architecture (CNN). To extract useful features, we applied an uniform noise, which affects differently each magnitude on the light curve by adding a random value $\in[-0.5,0.5]\,\textrm{mag}$, before passing through the encoder (see Fig. 4).

In the first part of the CNN, the encoder which is composed of 9 convolution layers (conv 1 to conv 9 in Fig. 7) and 4 pooling layers (Pool 1,4, 6 and 8), converts the input noisy LCIs into an embedding representation. Then, the decoder reconstructs the original LCIs from the embedding representation through two fully connected layers (FC10 and FC11) of 5000 neurons. So the output of the autoencoder is a reconstructed LCI with the same size than input, $h\times w$. The loss function associated to the autoencoder, called Autoencoder Loss on Fig. 7, minimizes the difference between the original LCI and the reconstruction one. However, we have to pay attention at the overfitting of missing data on the light curve. An illustration of the overfitting problem and the solution we propose is given on Fig. 5. By construction, the input LCI is composed of many zero values (see Section 3.1), and the autoencoder model can ”learn” the position of missing data and just replace them at the same position on the reconstructed LCI (see case 1 on Fig. 5). To avoid this overfitting, we introduce a regularization term. In this case, the autoencoder learns to interpolate zero values (case 2 on Fig. 5). But the computation of the autoencoder loss can not take into account these interpolated values as they are compared to zero values on the initial LCI, and so lead to a divergence of the loss function. Therefore we propose to define a mask with the same size as the considered original light curve, filling with 1 if there is an observation on the light curve, and 0 otherwise. The reconstructed LCI is then multiplied by the mask before the minimization of the loss function (case 3 on Fig. 5). Equation 1 becomes:

[TABLE]

with $M(X)$ the mask related to the input light curve $X$.

Finally, we compute the penalty term as defined in equation 5, in the second fully connected layer, FC11 and called it Sparsity loss. It depends on two hyperparameters: the sparsity parameter $\rho$ and the weight of the sparse regularization $\alpha$. To determine the best values of $\rho$ and $\alpha$, we searched the best combination using a 2D grid search among values into the following finite sets: $\{10^{-5},5\times 10^{-4},5\times 10^{-3},10^{-3},5\times 10^{-2},10^{-2},5\times 10^{-1},10^{-1}\}$ and $\{10^{-3},5\times 10^{-2},10^{-2},5\times 10^{-1},10^{-1}\}$ respectively.

However the regularization term does not take into account the number of observations on each light curve which varies significantly. It may cause overfitting as the number of active neurons is then always the same whatever the number of data points on each light curve. So the number of active neurons has to be adapted depending on the number of observations in all filters. Thus, we propose to express the sparsity parameter, $\rho$, as a linear function depending on the number of observation for each light curve. This contribution allows to increase the number of active (inactive) neurons when the light curve is densely (poorly) populated with observations. We define a new sparsity parameter $\rho^{\prime}(l)$ for the specific light curve noted $l$ as follow:

[TABLE]

with $n_{l}$ the number of observations on the light curve $l$, $\rho_{a}$ and $\rho_{b}$ are two hyper-parameters. They are determined as the same time as $\alpha$ using a 3D grid search among the same values as $\rho$.

In this case, the sparse regularization term (see equation 5) of our autoencoder module take the form :

[TABLE]

4.2 The contrastive branch

Once the autoencoder training has converged on the test database the weights of its convolution and fully connected layers are fixed. Then, the output of a chosen layer of the encoder part is given as input to the contrastive branch. This second module is designed to reduce the mismatch between the training (higher magnitudes) and the test (lower magnitudes) databases. This requires a specific contrastive loss function that is minimized through a CNN architecture. So we propose a loss function that minimizes the variations of intra-class light curves and maximizes the variations of inter-class light curves. In this way, we split the training database in four subsets following a cut magnitude $m^{c}$ in the i-band magnitude.

If we note $m^{Ia}(l)$ the i-band median magnitude of a type Ia light curve and $m^{non-Ia}(l)$ the i-band median magnitude of a non-Ia type light curve, a given light curve can belong to one of the four following subsets:

• •

LC1 : type Ia light curves with $m^{Ia}(l)<m^{c}$,

• •

LC2 : type Ia light curves with $m^{Ia}(l)>m^{c}$,

• •

LC3 : non-Ia type light curves with $m^{non-Ia}(l)<m^{c}$,

• •

LC4 : non-Ia type light curves with $m^{non-Ia}(l)>m^{c}$

Therefore the goal is to define a loss function that minimizes the variation between intra-class light curves, i.e between the LC1-LC2 and LC3-LC4 sets; and maximizes the variation between inter-class light curves, i.e between LC1-LC3, LC1-LC4, LC2-LC3, and LC2-LC4 sets.

Equation 7 becomes:

[TABLE]

We introduce $\frac{1}{2}$ terms into the formula to weight the inter-class distances so that the inter-class and the intra-class distances have the same weight in the computation of the loss function.

In practice this means that the encoder is fed with sets of four light curves from the training database, with one light curve from each subset. The learning of the contrastive branch (see Fig. 7) is done without updating the training weights of the autoencoder, that have been adjusted during the non-supervised step on the test database. This step allows also to solve the problem of asymmetry that exists between the classes as this module takes as input both light curves of type Ia and non-Ia supernovae at the same time. Features from the seventh convolution (conv 7 on Fig. 7) are given as input to the contrastive branch where the training weights are updated. Therefore the minimization of the contrastive loss is made only on the training database. The choice of the seventh convolution layer as input to the contrastive branch was made for several reasons. First of all, as the encoder part of the first module is dedicated to extract relevant features from the test light curves to characterize them precisely, while the decoder part is designed to reconstruct the original light curve, we decided to extract features from the first part of the autoencoder to reduce the mismatch between the training and the test databases. Figure 6, which represents the t-SNE111The t-distributed stochastic neighbor embedding (t-SNE, van der Maaten & Hinton, 2008) is a nonlinear dimensionality reduction technique well-suited for embedding high-dimensional data for visualization in a low-dimensional space of two or three dimensions. projections of features, offers means of better understanding. If the projection of features from the first fully connected layer (FC 10) of the autoencoder part shows a better separation of type Ia and non-Ia supernovae, than from the seventh convolution layer, the extraction of these features for the contrastive branch degrades the performance. This means that it is preferable to consider a feature representation space of light curves of high abstraction level rather than a representation apparently more suited for classification in the autoencoder layers, as it allows a significant reduction of the mismatch between the training and the test databases.

The last layers of the contrastive module (conv 9c and FC 11c) mark a clear separation between type Ia and non-Ia supernovae (bottom panel of Fig. 6).

4.2.1 The classification branch

The last module of PELICAN is composed of three fully connected layers (see Fig. 7). It takes as input features from the two first modules to perform the classification step. Indeed to make the final classification, this part needs information of the first module that fully characterizes light curves of the test database and so give a large variety of features that allows to reduce the mismatch between the training and test databases. However, this time we extract features from the decoder part as it was shown that it is able to make a separation of the classes that is relevant for this final step (see Fig. 6). Then the classification branch must benefit from features of the second contrastive branch, and particularly the fully connected layer (FC11c) that reduce again the mismatch while marking a separation between classes.

Finally to fight against the overfitting of missing data, the third module takes also as input, features from the ninth and tenth convolution layers of the contrastive branch (conv 9c and conv 10c). We apply a specific operation, called a global pooling, which allows to transform a 2D output feature vector of a convolution layer into a 1D feature vector given as input to a fully connected layer. We choose to apply a global max pooling that will select only the maximum value on the 2D output feature maps from the convolution layers, excluding zero values and so missing data.

We also make use of dropout technique (Srivastava et al., 2014) on the two fully connected layers FC13 and FC14 to fight against overfitting.

5 Light curve data

We test and adapt our method on three different databases. First we evaluate the techniques on simulated data from the Supernovae Photometric Classification Challenge (SPCC, Kessler et al., 2010a, b) then on simulated LSST light curves for the main survey and the deep fields. Finally we explore the possibility to make the learning step on simulated light curves and then to test on real data. We apply this last work on SDSS supernovae light curves (Frieman et al., 2008; Sako et al., 2008).

5.1 SNANA simulator

Light curves have been simulated using the SNANA simulator (Kessler et al., 2009). It is an analysis package for supernovae light curves that contains a simulation, a light curve fitter and a cosmology fitter. It takes into account actual survey conditions and so generate realistic light curves by using the measured observing conditions at each survey epoch and sky location. First the supernovae properties are generated by choosing a shape-luminosity and color parameters, that are used in addition to other internal model parameters to determine the rest-frame magnitude at each epoch. Then K-corrections are applied to transform rest-frame model magnitudes to observed magnitudes in the telescope system. Finally the ideal above atmosphere magnitudes are translated into observed fluxes and uncertainties. Type Ia supernovae light curves are simulated from SALT2 (Guy et al., 2007) or MLCS models (Jha et al., 2007; Kessler et al., 2009). But there are no such models for non-Ia types. So the simulations uses a library of spectral templates that give the supernovae flux as a function of epoch and wavelength. Only well-sampled photometric light curves are used because spectral templates are interpolated to cover all wavelengths and epochs. The current library contains composite and individual templates for types Ib, Ibc, IIn and IIP.

5.2 SPCC data

The SPCC dataset is composed of simulated light curves of supernovae in griz filters of the Dark Energy Survey (DES). The dataset were subdivided in a spectroscopically confirmed subset of 1,103 light curves, which should constitute the training dataset, and a test dataset of 20,216 light curves. However the training dataset is small and highly biased as it is not representative in brightness and in redshift compared to the test set.

5.3 Simulated LSST data

As LSST will observe a large amount of supernovae, the photometric classification of supernovae types from multi-band light curves is necessary. There will be two main kind of cadences. The first one dedicated to the main survey is called the Wide-Fast-Deep (WFD). It will scan a very large area of the sky. The second one, called Deep Drilling Fields (DDF), will focus on small part of the sky with a higher cadence and deeper images. Thus this will correspond to well measured light curves (see Fig. 1) but for a smaller sample.

To validate our method in the context of the future LSST data we simulated light curves of supernovae as observed with the WFD and DDF observational strategies, with the minion 1016 baseline model (Biswas et al., 2017). The simulation was realized in the ugriz filters of the LSST. We assume a $\Lambda$CDM cosmology with $\Omega_{M}=0.3156$, $\Omega_{\Lambda}=0.6844$ and $w_{0}=-1$. Simulations are made in a finite redshift range, $z\in[0.05,1.20]$. We consider an efficiency for the image subtraction pipelines reaching 50% around Signal-to-Noise Ratio (SNR) $\sim 5$. Each object must have 2 epochs in any band. For the simulation of type Ia light curves the color and the light curve shapes parameters vary in the following intervals: $c\in[-0.3,0.5]$, $x_{1}\in[-3,2]$. The simulation of non-Ia types is based on a library of spectral templates for types Ib, Ibc, IIn and IIP.

Our simulation includes a spectroscopically confirmed sample from the DDF survey. It is based on observations from a 8 m class telescope with a limiting i-band magnitude of 23.5. In this work we assume different allocating time for the spectroscopic follow-up. A reasonable scenario allows a spectroscopic follow-up of 10% of observed light curves in DDF, i.e 2k spectroscopically confirmed light curves of supernovae. But we also consider a most restrictive case by assuming that only 500 then 1,000 light curves are spectroscopically confirmed. Moreover we explore two ideal cases for which 5k then 10k supernovae have been followed up. Finally we also consider different number of photometric observations of light curves as it is interesting to classify light curves before 10-years observation of LSST. All the configuration are summarized on Table 7.3.1.

5.4 SDSS data

As simulated data do not reproduce perfectly the real data, it is interesting to test our method on real data. The ideal strategy is to simulate light curves that corresponds to SDDS survey to train the model and then test on real data. This is a challenging methodology as there is a severe mismatch between the training and the test databases. However making a model able to remove this kind of mismatch is crucial for the future surveys where the spectroscopic follow-up is limited. Therefore we simulated light curves of supernovae that corresponds to SDSS data. Then, we extracted light curves in ugriz filters from the SDSS-II Supernova Survey Data (Frieman et al., 2008; Sako et al., 2008). The SDSS-II SN data were obtained during three month campaigns in the Fall of 2005, 2006 and 2007 as part of the extension of the original SDSS. The Stripe 82 region was observed with a rolling cadence. Some spectroscopic measurements were performed for promising candidates depending on the availability and capabilities of telescopes (Sako et al., 2008). A total of 500 SN Ia and 82 core collapse SN were spectroscopically confirmed.

6 Experimental protocol

In this section, we explained the protocol and the different techniques used for the training process.

6.1 Data augmentation

In this classification context, data augmentation is a crucial step. Indeed, in order to make PELICAN robust against the differences between the training and the test databases (i.e sampling, mean magnitude, noise…), it is essential to use different data augmentation techniques. Moreover when light curves that composed the training and test databases are measured with different observational strategies the difference in sampling is increased and the data augmentation has to be reinforced. It is the case in the context of LSST if we compare light curves from the WFD survey on the one hand, and light curves from the DDF survey on the other hand. To make PELICAN able to learn on DDF light curves and generalize on WFD light curves, the data augmentation has to be adapted.

Finally as supernovae from the test database are often fainter, errors on their fluxes are often bigger. Therefore the data augmentation needs also to be applied to the errors.

Thus, in addition to the uniform noise applied differently on each magnitude of light curves given as input to the denoising autoencoder, we add two other kind of noise on magnitudes of light curves:

• •

an uniform constant noise $\in[-1.0,1.0]\,\textrm{mag}$ which is added to all the magnitudes of the light curve,

• •

an uniform noise $\in[0.93,1.07]$ which is multiplied by all the magnitudes of the light curve,

Then we randomly remove one or several magnitudes or/and all magnitudes for a given band. This process is particularly effective for the classification of light curves of supernovae observed with a WFD strategy based on a training on supernovae light curves from the DDF survey.

Finally, to prevent PELICAN model from learning the missing value positions on each light curve, we perform random time-translations keeping all the data points but varying their positions in time. So the learning becomes invariant to the position of points.

6.2 Setting learning parameters

We used the Adam optimizer (Kingma & Ba, 2014) for all the training steps in different modules with a learning rate decreasing by a 10 factor after 25 000 iterations. The batch size during the learning is fixed to 96 light curves.

For the autoencoder learning, we optimized the values of the sparsity parameters over one validation base and use them for all the different configurations, as it is not sensitive to the database. We set $\rho_{a}$ and $\rho_{b}$ equal to $5\times 10^{-4}$ and $0.0$ respectively and $\alpha$ to $0.01$.

The cut parameter $m^{c}$ in the i-band magnitude from the second module, depends on the database. We choose its value in order to have enough examples on both sides of the cut in magnitude. We set $m^{c}$ to 23.5 mag, 24.5 mag and 22.5 mag for the SPCC, LSST and SDSS databases respectively. The values of these parameters are not sensitive and a small variation of them did not change the results.

6.3 Ensemble of classifiers

To increase the performance, we trained an ensemble of classifiers as it was shown to be more accurate than individual classifiers (e.g. Polikar, 2006). Moreover the generalization ability of an ensemble is usually stronger than that of base learners. This step involves training N times one model with the same training database but a different initialization of the weights. We chose N=7 and the individual decisions were then averaged out to obtain the final values of probabilities. This step allows to increase the accuracy of $2\%$ on average.

7 Results

In this section we present the results that we obtained for each dataset.

7.1 Metrics

To evaluate the performance of PELICAN in different contexts we use several commonly used statistic metrics that are the Accuracy (Acc), the recall (R) or true positive rate (TPR), the precision (P) and the false positive rate (FPR). They are defined from the following confusion matrix:

[TABLE]

[TABLE]

As a graphical performance measurement we use the ROC (Receiver Operating Characteristic) curve which is plotted with TPR on y-axis against FPR on x-axis. It gives an estimation of the performance of a classifier at different thresholds settings. The best possible prediction method would yield a point in the upper left corner or coordinate (0,1) of the ROC space, representing the lack of false negatives and false positives. A random guess would give a point along a diagonal line from the left bottom to the top right corners.

From the ROC graphic, we can extract the value of the AUC (Area Under the Curve) which captures the extent to which the curve is up in the upper left corner. The score has to be higher than 0.5 which is no better than random guessing.

7.2 SPCC

The first evaluation of PELICAN is made on the SPCC dataset. We trained the model with two different training datasets: a representative training database and a non-representative training dataset. The representative training database is a simplified theoretical scenario, in which there is no limitation in brightness and redshift of the spectroscopic follow-up. It is built by selecting randomly 1,103 light curves from the whole dataset. The non-representative training database, which represents the real scenario, is the spectroscopically confirmed subset of 1,103 light curves that was proposed for the challenge. This last one is non-representative of the test dataset in brightness and redshift.

As shown in Lochner et al. (2016, noted L16 hereafter) the best average AUC is obtained by extracting SALT2 features and using boosted decision trees (noted BDTs hereafter) as classifier. Therefore we compared the performance of PELICAN with BDTs algorithm that take as input SALT2 features. We test both methods with and without the information of the redshift inside the training.

The ROC curves for both methods are represented on Fig. 8 (on left panels) and values of statistics are reported in Table 1.

By considering a non-representative training database, without including the redshift during the training, PELICAN obtains an accuracy of 0.856 and an AUC of 0.934 which outperforms BDTs method which reaches 0.705 and 0.818. If we train both methods on a representative training database, as expected, the performance increases. The accuracy and the AUC become 0.911 and 0.970 with PELICAN, against 0.843 and 0.905 with BDTs algorithm. It is interesting to note that the gain in statistics obtained with PELICAN, is lower than of BDTs values, which means that PELICAN is able to better deal with the problem of mismatch. This ability will be confirmed by the promising results obtained with a non-representative training database composed of LSST light curves simulations (see Section 7.3).

The performance of PELICAN does not change by adding the redshift information during the training, which is not the case for BDTs algorithm, for which the accuracy and the AUC are slightly increased. This might means that PELICAN is able to extract by itself the redshift information during its training. Figure 8 shows the accuracy as a function of redshift, with the corresponding redshift distributions and BDTs results for comparison. If PELICAN is trained on a representative training database, the accuracy tends to decrease at low redshifts and at redshift above 1.0, as the corresponding redshift bins are poorly populated at these extreme values. A further trend is observed for BDTs, except at redshift above 1.0, and only if redshift values are included as an additional feature for the training. By considering a non-representative training database, the accuracy significantly decreases at high redshift for both methods. As the addition of redshift inside the training does not change the tendency obtained by PELICAN, this trend in function of redshift is likely due to the too small number of examples at high redshifts.

7.3 LSST simulated light curves

The next step is to evaluate PELICAN on simulated LSST light curves under realistic conditions. In this way, we consider for all the tests a non-representative spectrocopically confirmed database from the DDF survey as represented on Fig. 2. We consider different configurations of the training and test databases. We constrain the number of spectroscopically confirmed light curves to vary between 500 light curves to 10k. Even if the upper bound corresponds to an ideal scenario in which roughly more than 40% of light curves in the DDF have been following up, it is interesting to compare the performance of PELICAN with large training sample.

We simulated light curves with the minion 1016 cadence model. This model includes a WFD survey and five DDF (see Fig. 9)
It is not certain that a spectroscopic follow-up will be performed on supernovae light curves in WFD fields. So we use a different approach which consists to train PELICAN on DDF light curves and then adapt the pre-trained model to classify supernovae light curves observed in WFD survey. This strategy allows to consider the possibility to benefit from SN Ia candidates from WFD fields to constrain the cosmological parameters, without any spectroscopic follow-up of the main survey.

7.3.1 Classification of DDF light curves

The results of the different configurations are reported in Table 7.3.1 and the ROC curves for some of these configurations on Fig. 10. In addition to the values of the accuracy and the AUC, we compare the recall of SN Ia by constraining the precision to be higher than $95\%$ and $98\%$. Such a level of contamination becomes competitive with spectroscopy contamination. Again we compared the performance of PELICAN with the best method highlighted in L16, that is BDTs and SALT2 features. Even if this method was not designed for such training configurations, it allows to compare a feature-based machine learning method to PELICAN.

If we consider the most constraining configuration composed of only 500 spectrocopically confirmed light curves for the training and 1,500 light curves for the test database, PELICAN reaches an accuracy of 0.895, and an AUC of 0.966. Moreover PELICAN is then able to detect 76.9% of SN Ia with a precision higher than 95% and 60.2% with a precision higher than 98%. These results are quickly improved by considering more examples on both training and test databases. The number of light curves inside the test database is important, especially if the number of examples in the training database is small, as the autoencoder is trained on the test database. Indeed there is about a 8% improvement factor of the recall by going from 1k to 3k light curves in the test database with a fixed number of examples in the training database of 1k light curves. But this factor becomes negligible if the number of spectroscopic confirmed light curves is sufficient, i.e from 5k examples, with an improvement of around $0.2\%$.

The configuration that seems reasonable after 10 years of observation includes a spectroscopic follow-up of 10% of the observed light curves, i.e 2k light curves of supernovae, and a test database of 22K light curves. For this realistic scenario, PELICAN reaches an accuracy of 0.942 and is able to correctly classify 87.4% of SN Ia with a precision higher than 98%, which constitutes a major result of our study meaning that a large fraction of SN Ia are well-classified by PELICAN, with a precision comparable to a spectroscopy measurement. By considering 10k light curves in the training database, the number of detected SN Ia is then increased by 9%. All results obtained by PELICAN outperform those obtained by BDTs (the BDTs values are listed in parenthesis in Table 7.3.1).

The right panel of Fig. 10 shows the accuracy as a function of redshift, with the corresponding redshift distributions on both training and test databases, and the BDTs results for comparison. The accuracy of PELICAN does not depend on redshift until 1.0 where it slightly decreases. This tendency is likely due to the small number of training examples at redshift higher than 1.0. BDTs method shows the same behaviour at high redshifts.