On the cosmological performance of photometrically classified supernovae with machine learning

TL;DR

This paper evaluates machine learning methods for photometric supernova classification to improve cosmological measurements, demonstrating effective feature selection and thresholding techniques on simulated data.

Contribution

It introduces a robust threshold selection method based on Bias-Variance tradeoff for supernova classification, enhancing cosmological constraint accuracy.

Findings

Achieved strong cosmological constraints with ~75% information retention using SALT2 features.

Demonstrated the effectiveness of decision tree ensembles and automated search pipelines.

Proposed a robust thresholding approach to optimize classification performance.

Abstract

The efficient classification of different types of supernova is one of the most important problems for observational cosmology. However, spectroscopic confirmation of most objects in upcoming photometric surveys, such as the The Rubin Observatory Legacy Survey of Space and Time (LSST), will be unfeasible. The development of automated classification processes based on photometry has thus become crucial. In this paper we investigate the performance of machine learning (ML) classification on the final cosmological constraints using simulated lightcurves from The Supernova Photometric Classification Challenge, released in 2010. We study the use of different feature sets for the lightcurves and many different ML pipelines based on either decision tree ensembles or automated search processes. To construct the final catalogs we propose a threshold selection method, by employing a \emph{Bias-Variance tradeoff}. This is a very robust and efficient way to minimize the Mean Squared Error. With this method we were able to get very strong cosmological constraints, which allowed us to keep $\sim 75\%$ of the total information in the type Ia SNe when using the SALT2 feature set and $\sim 33\%$ for the other cases (based on either the Newling model or on standard wavelet decomposition).

Tables

Table 1. Table 1: Confusion matrix for a boolean classification problem.

		True Class
		P	N
Predicted class	P	True positive (TP)	False positive (FP)
	N	False negative (FN)	True negative (TN)

Table 2. Table 2: Uniform prior ranges on the SALT2 model parameters.

Parameter	Range
$z$	$(0.01,1.5)$
$t_0$	$(-60,100)$
$x_0$	$(-{10}^{-3},{10}^{-3})$
$x_1$	$(-3,3)$
$c$	$(-0.5,0.5)$

Table 3. Table 3: Uniform prior ranges on the Newling model parameters.

Parameter	Range
$\log \left(A\right)$	$(0,10)$
$\varphi$	$(-60,100)$
$\log \left(\sigma \right)$	$(-3,4)$
$\log \left(k\right)$	$(-4,4)$
$\log \left(\psi \right)$	$(-6,10)$

Table 4. Table 4: Uniform prior ranges on the Karpenka model parameters.

Parameter	Range
$\log \left(A\right)$	$(-3,3)$
$\log \left(B\right)$	$(-3,2)$
$t_0$	$(0,100)$
$t_1$	$(0,100)$
$T_{\text{rise}}$	$(0,100)$
$T_{\text{fall}}$	$(0,100)$

Table 5. Table 5: Comparison between the best AUC scores in this work and the literature. (*): this information is absent in the original paper. The comparison is made between Accuracy values (defined as ( TP + TN ) / N TP TN N \mathrm{(TP+TN)/N} )

	AUC	AUC ($+z_{\mathrm{host}}$)
This work	0.969	0.985
Lochner et al. (2016)	0.969	0.979
Charnock & Moss (2017)	0.901	0.940
Möller & Deboissì Ere (2018)	0.980	(*)

Table 6. Table 6: The effective completeness (in percentage) of each method, for both fitted models. This quantity illustrates how much information is preserved after the photometric classification. The quoted error bars are the standard deviation in this quantity among our 10000 catalogs.

	$\Lambda CDM$				$wCDM$
	Binned Purity		Constant Purity		Binned Purity		Constant Purity
	AP	AUC	AP	AUC	AP	AUC	AP	AUC
Karpenka	$19\pm 5$	$19\pm 6$	$20\pm 4$	$16\pm 5$	$21\pm 6$	$19\pm 6$	$23\pm 4$	$20\pm 5$
Newling	$24\pm 7$	$32\pm 9$	$23\pm 6$	$30\pm 4$	$22\pm 7$	$30\pm 10$	$24\pm 5$	$27\pm 5$
Wavelets	$36\pm 13$	$16\pm 8$	$22\pm 4$	$22\pm 4$	$30\pm 13$	$10\pm 8$	$12\pm 3$	$13\pm 3$
SALT2	$66\pm 15$	$75\pm 15$	$48\pm 14$	$54\pm 17$	$62\pm 16$	$75\pm 15$	$37\pm 12$	$44\pm 16$

Table 7. Table 7: Functions used in pipelines built by the tpot genetic optimization.

Library	Module	Function	Website
copy		copy	docs.python.org/3/library/copy.html
xgboost		XGBClassifier	xgboost.readthedocs.io
tpot	builtins	StackingEstimator	epistasislab.github.io/tpot
		ZeroCount
sklearn	pipeline	make_union	scikit-learn.org
	decomposition	FastICA
	linear_model	LogisticRegression
	svm	LinearSVC
	neighbors	KNeighborsClassifier
	naive_bayes	GaussianNB
	kernel_approximation	Nystroem
		RBFSampler
	ensemble	ExtraTreesClassifier
		GradientBoostingClassifier
		RandomForestClassifier
	feature_selection	RFE
		SelectFwe
		f_classif
		VarianceThreshold
	preprocessing	FunctionTransformer
		MaxAbsScaler
		MinMaxScaler
		PolynomialFeatures
		Normalizer
		RobustScaler

Equations

$$\mathrm{TP}+\mathrm{FP}+\mathrm{TN}+\mathrm{FN}=\mathrm{N},$$

$$\mathcal{C}=\mathrm{TPR}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}.$$

$$\mathrm{FPR}=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}.$$

$$\mathcal{P}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FP}}.$$

$${\mathcal{P}}_{\min }=\frac{\mathrm{TP}+\mathrm{FN}}{\mathrm{N}},$$

$$\mathrm{AP}\equiv \sum _n({\mathcal{C}}_n-{\mathcal{C}}_{n-1}){\mathcal{P}}_n,$$

$$\mathrm{b}=⟨\Delta \mu ⟩,\mathrm{MSE}=⟨\Delta {\mu }^2⟩.$$

$$\mathrm{Var}[\Delta \mu ]=⟨\Delta {\mu }^2⟩-⟨\Delta \mu {⟩}^2=\mathrm{MSE}-{\mathrm{b}}^2,$$

$$\mathrm{MSE}=\mathrm{Var}+{\mathrm{b}}^2.$$

$$\varphi (p,\lambda ;x_0,x_1,c)=x_0[M_0(p,\lambda )+x_1{M}_1(p,\lambda )]\exp [cCL(\lambda )],$$

$$F^Y(p(1+z))=(1+z)\int \varphi (p,{\lambda }^{\prime })T^Y({\lambda }^{\prime }(1+z))d{\lambda }^{\prime },$$

$$F(t)=A{\left(\frac{t-\varphi }{\sigma }\right)}^k\exp \left(-\frac{t-\varphi }{\sigma }\right)k^{-k}{e}^k+\Psi (t),$$

$$\Psi (t)=\{\begin{array}{ll}0& -\infty <t<\varphi \\ \text{cubic spline}& \varphi <t<\tau \\ \psi & \tau <t<\infty \end{array}.$$

$$F(t)=A\big{[}1+B(t-t_{1})^{2}\big{]}\,\frac{\exp[-(t-t_{0})/T_{\text{fall}}]}{1+\exp[-(t-t_{0})/T_{\text{rise}}]},$$

$${\chi }^2$$

$$C_{ij}=\left(\frac{{\mathrm{Var}}_i}{N_{SN}^i}\right){\delta }_{ij}+{\mathrm{b}}_i{\mathrm{b}}_j.$$

$${\mathcal{C}}^{\mathrm{eff}}\equiv \frac{N_{\mathrm{Ia}}^{\mathrm{eff}}}{N_{\mathrm{Ia}}}={\left(\frac{{\sigma }_{Ia}}{\sigma }\right)}^2.$$

$$\mu =m_B-M_B=m_B^{\ast }+\alpha x_1-\beta c-M_B,$$

$$m_B^{\ast }=-2.5{\log }_{10}\left[x_0\int M_0(t=t_0,{\lambda }^{\prime })T^B({\lambda }^{\prime })d{\lambda }^{\prime }\right].$$

Full text

††pubyear: 2019††pagerange: On the cosmological performance of photometrically classified supernovae with machine learning–B

1 Introduction

Current and next generation of optical surveys which cover large areas with good cadence, like the Dark Energy Survey (DES, Abbott et al., 2016), The Rubin Observatory Legacy Survey of Space and Time (LSST, LSST Science Collaboration et al., 2009), and the Zwicky Transient Facility (ZTF, Bellm, 2014) will not provide spectra of the majority of observed objects. Alternatively, the Javalambre Physics of the Accelerating Universe (J-PAS) survey will provide effectively very coarse spectra with the use of 56 narrow band filters (Benitez et al., 2014), and other smaller surveys follow similar strategies, albeit with less filters (Cenarro et al., 2019; Mendes de Oliveira et al., 2019). This has pushed the community to develop tools which optimize the amount of science that can be gained from the purely photometric data. Most effort has been focused on how to determine the redshift (see, for instance, Abdalla et al. (2011); Dahlen et al. (2013); Benitez et al. (2014); Sadeh et al. (2016); Gomes et al. (2018) and references therein). For the cosmological use of supernovae (SNe), determining the redshift is not enough, as it is crucial to correctly classify them in order to isolate the type Ia SNe (henceforth SNeIa, which are standardizable candles) from all other (core-collapse) types of SNe, for which distance estimations are not yet competitive when compared to the ones with SNeIa.

For the flat $w$CDM cosmological model, in which dark energy has a constant equation of state, taking into account all systematic uncertainties, SNeIa alone are able to constrain $\Omega_{m0}$ and $w$ to $23\%$ and $20\%$, respectively (Scolnic et al., 2018). And since SNe data imply different parameter degeneracies than CMB (Cosmic Microwave Background) and BAO (Baryon Acoustic Oscillations), the final uncertainties get much smaller when combining these datasets. A proper classification and redshift determination of future SNeIa would allow for unprecedent precision tests on the $\Lambda$CDM model as the number of events would grow from the current number of $\sim 10^{3}$ (Betoule et al., 2014; Scolnic et al., 2018) to over $\sim 10^{6}$ (LSST Science Collaboration et al., 2009).

Such an increase in statistics opens up many new avenues for SN data exploration beyond the conventional Hubble-Lemaître diagram. For instance one can directly measure: anisotropies in the Hubble-Lemaître diagram (Koivisto et al., 2011; Colin et al., 2011; Soltis et al., 2019); the Hubble constant for SNeIa in nearby galaxies (Riess et al., 2016; Burns et al., 2018); non-Gaussianities in the scatter of the Hubble-Lemaître diagram due to lensing (Quartin et al., 2014; Castro & Quartin, 2014; Macaulay et al., 2017); correlations in the same scatter induced by SNeIa peculiar velocities (Gordon et al., 2007; Castro et al., 2016; Howlett et al., 2017; Garcia et al., 2020); the Hubble function combining SNeIa velocities and number densities (Amendola & Quartin, 2019); and time-delays between multiple images of strong lensed SNeIa (Oguri & Marshall, 2010; Bonvin et al., 2017; Birrer et al., 2019; Huber et al., 2019).

The SN classification scheme does not need to be boolean. Techniques have been proposed in Kunz et al. (2007); Newling et al. (2012) to make use of SN catalogs for which each event has a set of probabilities of belonging to different types. Nevertheless, the higher the classification probability of each object, the better these kind of techniques can work. Current SN catalogs have employed a SN classification scheme based on an empirical template fitting (Sako et al., 2018; Jones et al., 2018). But in recent years alternative classifications scheme have been proposed (Ishida & de Souza, 2013) and a couple of large classification challenges have been held (Kessler et al., 2010a; Malz et al., 2019). In particular, the use of machine learning (ML) techniques in SN classification is starting to be studied in more detail (Lochner et al., 2016; Charnock & Moss, 2017; Moss, 2018). Ishida (2019) indeed claims that the use of ML methods will be inevitable for the future of supernovae cosmology and automated transient classification.

The historical standard technique to perform photometric classification is template matching. It consists in constructing the photometric light curves as a combination of spectra from a template set, and comparing it with the real one. The coefficients are then determined performing a numerical optimization (say, by a $\chi^{2}$ minimization); see for instance Sako et al. (2018). This is still the most common technique in Supernovae Surveys, and form the basis of SNe compilations such as the SDSS-II (Kessler et al., 2009), JLA (Betoule et al., 2014) and Pantheon (Scolnic et al., 2018) catalogs. In this decade some groups started to use Machine Learning techniques to improve this classification. In SDSS-II the authors applied the Nearest Neighbors (NN) algorithm to see if it could improve their Photometric SN IDentification (PSNID) template fitting software (Sako et al., 2018). The Pan-STARRS team (Jones et al., 2018) compared the NN with 3 different template fitting codes: PSNID, GalSNID (Foley & Mandel, 2013) and Fitprob (the fit probability from the SALT2, Spectral Adaptive Lightcurve Template 2, light curve fit multiplied by a redshift-dependent SN type prior). Recently Villar et al. (2019) employed machine learning methods to photometrically classify Pan-STARRS SNe. The first public release of the scikit-learn111scikit-learn.org (Pedregosa et al., 2011) library in 2010 increased the popularity of the ML methods.222scikit-learn provides sophisticated machine learning models in a simple implementation, making it easy to write efficient codes with few lines. In Lochner et al. (2016) (henceforth L16) the authors performed the first thorough comparison between many scikit-learn models (including the NN) and found that an ensemble of decision-trees , the AdaBoost classifier, reached the best AUC score, which is the area under the Receiver Operating Characteristic (ROC) curve (see below).

In this work we use the simulations performed as part of the Supernova Photometric Classification Challenge (SNPCC), discussed below in Section 3, in order to perform our ML training and tests. We start by performing a comparison similar to L16 but expanding to other ensemble Decision Trees provided by scikit-learn: Adaptive Boosting (ADA), Random Forest (RF), Extra-Trees (EXT), Gradient Boosting (GB), and the Extreme Gradient Boosting (XGB) provided by the XGBoost library333xgboost.ai and testing also an alternative metric to characterize the performance of the classifier: the Average Purity (AP), while L16 used only the AUC. We also performed an automated machine learning training by using TPOT444epistasislab.github.io/tpot (Tree-Based Pipeline Optimization Tool), which allowed us to find the best pipelines for each feature set, while L16 used a simple pipeline for all cases, a StandardScalar transform followed by the machine learning model. We then explore how a ML-classified SN catalog performs in cosmological parameter estimation, quantifying the degradation due to the imperfect completeness and purity of the resulting catalogs. These results should be directly applicable to the SNe in DES due to the similar numbers of both SNe with only photometry and of SNe with spectra that can be used to train the algorithm.

This paper is structured as follows: in section 2 we summarize the SNPCC and its specifications; in section 3 we revise the supervised learning approach paying special attention to the Bias-Variance tradeoff method 3.1, with more details on decision trees 3.2 and automated machine learning 3.3; in section 4 discuss the different features sets. Finally we present our results and discuss them in section 5 and the conclusions in 6.

2 Large photometric surveys and the Supernova Photometric Classification Challenge

One of the greatest challenges in photometry is to infer values of observable quantities which are traditionally computed from spectra. The most recurrent example is the photometric determination of galaxy redshifts. In this work we focus solely on a separate photometric problem: supernova classification. While the different SN types are defined by the presence of specific spectral features, photometric classification must rely on their different lightcurve shapes and colors.

The Supernova Photometric Classification Challenge (SNPCC) was proposed by Kessler et al. (2010b) in order to stimulate the development of tools to address this problem. They have publicly released a blinded mix of simulated SNe, with the different types selected in proportion to their expected rate and as if they had been measured with the griz filters of the DES with realistic observing conditions (sky noise, point spread function and atmospheric transparency) based on years of recorded conditions at the DES site. Simulations of non-Ia SNe are based on spectroscopically confirmed light curves that include non-Ia samples donated by the Carnegie Supernova Project (CSP), the Supernova Legacy Survey (SNLS), and the Sloan Digital Sky Survey-II (SDSS–II).

The SNPCC catalog contains 21319 supernova light curves corresponding to five seasons of observations, simulated by the SNANA software. This set is composed by 5086 SNeIa and 16231 core collapse SNe, distributed among 7 different types: II, IIn, IIp, IIL, Ib, Ib/c and Ic. The type Ia supernovae were simulated using an equal mix of the models MLCS2k2 (Jha et al., 2007) and SALT2 (Guy et al., 2007) with an additional random color variation. An extinction correction, MLCS-U2 was used in order to make the models agree in the ultraviolet (Kessler et al., 2009). The non-Ia supernovae light curves were based on spectroscopically confirmed observed ones and constructed by warping a standard spectral template to match the photometry. The SN rates were based on the results from Dilday et al. (2008) and Bazin et al. (2009). Ten groups addressed the challenge using different algorithms including template matching and ML techniques. The results of the challenge were published in Kessler et al. (2010a): the template matching code PSNID from Sako et al. (2008) was found to be the best one based on the figure-of-merit employed (which was a function of purity and completeness).

After the end of the SNPCC the catalog was unblinded, and it now serves as a very useful tool to train and test ML algorithms relying on supervised training, in which a subset of the data has known classification a priori. In our case in particular this corresponds to assuming we will have a subset of the SNe observed spectroscopically.

We discuss some technicalities involving the SALT2 lightcurve fitting analysis pertaining the simulated SNPCC supernovae in Appendix A.

3 The supervised machine learning approach

The objective of a supervised learning algorithm is to construct rules, according to a method, which relate a set of feature (observable) variables to a target one. If the rules are too simple it leads to underfitting, which returns poor predictions and leads to poor performance. If the rules are instead too complex they may result in overfitting, which results in a model that reproduces peculiar and unimportant behavior present in the training set data but not in the test data or on the real data. A balance must be found and the hyperparameters (parameters of the model) chosen in order to optimize a score function. For a linear model, for instance, the hyperparameters are simply the intercept and the slope, usually set to minimize the mean squared error.

In classification problems, like the Supernovae classification, the target variables are classes, represented by integer numbers. In this work we are interested in the cosmological use of type Ia supernovae so we consider only two classes: Type Ia (Ia), the positive class (1); and Core Collapse (CC) supernovae, the negative class (0). The output of a classification is the probability that an input corresponds to a positive class, Type Ia in our case. In order to perform an univocal classification we must choose a probability threshold, i.e. a probability value above which we assume to correspond to the positive class. The methods used to determine this value are referred to as the selection criteria.

The quality of a classifier is completely summarized by four quantities: (i) The number of true positives (TP), i.e. the positive values which were correctly classified; (ii) The number of true negatives (TN), i.e. the negative values which were correctly classified; (iii) The number of false positives (FP), i.e. the positive values which were not correctly classified; (iv) The number of false negatives (FN), i.e. the negative values which were not correctly classified. These are the components of the confusion matrix, shown in Table 1.

All these quantities vary with the probability threshold, but always respect the constraint

[TABLE]

where $N$ is the length of the data set. So we must choose a value that maximizes the true predictions (TP and TN) and minimizes the false predictions (FP and FN). However we can not be sure that the best threshold found in the training process is the best for predictions. To circumvent this problem during the optimization process we choose a score function which depends of all threshold values (between 0 and 1). The score functions are derived from the confusion matrix, and in this work we use two different functions: the area under the receiver operating characteristic (ROC) curve score, which is the same used in L16, and the Average Precision. To define these scores we must first define three other quantities.

Completeness ($\mathcal{C}$) or True Positive Rate (TPR): the fraction of real positives, which are correctly classified, used in final set. Also referred sometimes in the ML literature just as “recall” or “sensitivity”, it is given by

[TABLE]

False Positive Rate (FPR): the fraction of real negatives which are not correctly classified, which contaminate the final set. Also referred as false alarm rate, it is given by

[TABLE]

Purity ($\mathcal{P}$): The fraction of real positives in the final set. Sometimes referred to in the ML literature as just “precision”, it is given by

[TABLE]

With the above definitions we may define the aforementioned scores.

Area Under the ROC Curve (AUC): both TPR and FPR are, by definition, quantities between 0 and 1. They are used to construct the ROC curve, which is just a plot of TPR versus FPR, and is used as a metric for evaluating the performance of classification algorithms. Imperfect classifiers typically have monotonically increasing ROC curves, which often tend to zero in the limit in which FPR goes to zero and the area under this curve is naturally less than one, as shown in Figure 1. A perfect classifier has FN$=$FP$=0$ and thus ${\cal C}={\cal P}=1$ and TPR $=1$. Its ROC is technically ill-defined, but can be assumed to correspond to the limit in which the ROC tends to a step-function with TPR $=1$ and unity area. The area under this curve is our first score function.

Average Purity (AP): similarly we may define another curve as the plot of Purity versus Completeness. This curve is know as Completeness-Purity (Markel & Bayless, 2019; Kgoadi et al., 2019; Villar et al., 2019). In contrast with the ROC, this curve is decreasing. For the minimum threshold, where all objects are positively classified, the purity equates to the total fraction of the positives in the raw, unclassified data (Saito & Rehmsmeier, 2015). For higher thresholds the purity is expected to increase, so usually the minimum purity is simply

[TABLE]

as shown in Figure 1. For a perfect classifier, the Purity is equal to one and the AP curve corresponds to the unity square with ${\cal C}=1$, similar to the case of the ROC curve. However the area under this curve is not a reliable score thanks to the numerical noise. To avoid this problem we use the Average Purity (dubbed Average Precision on scikit-learn and many ML references) as our second score function:

[TABLE]

where the summation is over the thresholds.

Typically a ML problem consists in constructing an optimized pipeline of numerical preprocessing and machine learning methods that can give the best scores, where the score function is chosen according to the problem. As mentioned above, in this work we test the overall classification performance when optimizing one of two different score functions. The first, which follows the SNPCC approach, is the AUC – see also Swets et al. (2000); Fawcett (2004) for accessible discussions on the subject. Since this metric is one of many which can be derived from the confusion matrix, we explore an alternative metric as well, which relies on the area under the curve of completeness versus purity. This is motivated by the fact that in Astronomy one is usually more interested in high-purity catalogs. And since purity is not a linear combination of TPR and FPR, a code that maximizes this AP will not necessarily maximize the AUC. In addition, with the Completeness-Purity curve we can easily infer the incompleteness of the final sample, which is important to manage possible biases in SN cosmology.

Our algorithm is depicted in Figure 2. It was constructed as follows: from the 21319 supernovae, we select randomly 1100 as our training set and leave the remaining 20219 as a testing set. We further subdivide the training set into 5 subsamples in order to perform a jack-knife approach: we use 4 out of 5 subsamples in the training process and cross-validate on the fifth subsample. This procedure is repeated for all 5 combinations internally by the scikit-learn function RandomizedSearchCV, and the result with the best score in the validation set is kept. This preliminary procedure is called 5-fold cross validation split in the ML literature, and is very important to avoid overfitting. It was also used by L16. Each training step consists in drawing a set of 30 points in the ML model hyperparameter space and finding which one maximizes the score in the testing set.

This whole procedure is repeated as many times as needed in order to get an accurate evaluation of the performance of the algorithm. In our case, we repeated this procedure at first 200 times for each ML method, each time using a different training set of 1100 supernovae, which is enough to get accurate results for the mean and variance of the AUC and AP. We later extended this number to 10000 repetitions for the best performing method in order to also get accurate results for the Bias-Variance tradeoff performance (discussed below).

We also store a catalog of all SNeIa from each test set, which we identify as the result of a perfect classifier, i.e. the one which always returns $\mathcal{C}=1$ and $\mathcal{P}=1$. This corresponds to approximately 4800 SNeIa, for any probability threshold. This will be useful when we compare the cosmological performance of the photometric classified SNeIa with the ideal case in Section 5.2.

3.1 The Bias-Variance tradeoff

We define the bias as the average displacement of the estimated distance modulus from the fiducial value, $\Delta\mu=\bar{\mu}-\bar{\mu}_{\rm fid}(z_{\rm sim})$, and the Mean Squared Error (MSE) as the average squared displacement:

[TABLE]

We will also make use of the Root Mean Squared Error (RMSE) throughout the text. Note that there is a difference between two averages: the one denoted by a bar over the quantity (such as $\bar{\mu}$) represents the mean in a given bin of a given catalog, while the average denoted by $\langle...\rangle$ represents the ensemble average over the 10,000 catalogs. Using the definition of variance

[TABLE]

we find the relation between these quantities

[TABLE]

In general an increase in purity leads to a decrease in the bias at the cost of an increased variance. The Bias-Variance tradeoff consists in determining which hyperparameters values minimize the full MSE, instead of solely the variance. The typical Bias-Variance tradeoff, as presented in most textbooks and tutorials, is enunciated for regression problems, where the balance between the variance and bias depends on the model complexity. However, since in this work we investigate how a classification problem affects the cosmological distance estimator, the tradeoff depends on the probability threshold, which can be mapped into the target purity of the catalog, as shown in Figure 3. In this plot, the green line depicts the absolute value of the Bias and the dashed segment represents negative values. This sign flip can be very important for the optimization process, and in fact for our purposes it is, as we show in Section 5.2.

3.2 Decision trees and ensembles

The first approach that we have investigated is the ensemble of decision trees. The decision tree method consists in constructing a decision flow structure, known as tree, which splits the feature space in a rectangular grid. Each grid cell is then associated with a class (although a given class may end up mapped into various cells). Figure 4 shows an example of a small decision tree for the SALT2 feature set. The main parameters of a typical tree are the maximal depth and the minimum sample per leaf. The leaf is the final node, which returns a class.

An ensemble of trees is a combination of many different trees, randomly trained. These methods are very useful for analysis of light curves, and some of then were already investigated in recent papers (Lochner et al., 2016; Markel & Bayless, 2019; Kgoadi et al., 2019; Villar et al., 2019). In this work we perform a broad comparison between five different ensemble methods. There are two different ways to build an ensemble of decision trees:

Averaging Methods: the driving principle is to build several independent weak estimators and then use the average of their predictions as the final prediction. On average, the ensemble estimator is usually better than any single estimator because its variance is reduced. In this work we use two such methods.

• i

Random Forest (RF), which constructs a distribution for the classes of the individual trees and take the mode of that as the final class.

• ii

Extremely Randomized Trees or Extra-Trees (EXT), which includes a special kind of Decision Tree, known as Extra-Tree, that differs from classic decision trees in the way they are built. The thresholds are sorted randomly for each candidate feature and the ones with best results are used as the splitting rule.

• 2.

Boosting Methods: weak estimators, typically slightly better than a completely random choice, are built sequentially and one tries to reduce the bias of the combined estimator. In this work we used three different boost methods.

• i

Adaptive Boost, or AdaBoost (ADA), which is a weighted average of the individual predictions, from weak learners, where the weights are computed iteratively.

• ii

Gradient Boost (GB), which uses the gradient descendent method to minimize the loss function.

• iii

Extreme Gradient Boost(XGB), which is an optimized implementation of the GB available in xgboost.ai.

3.3 Automated Machine Learning

The automated Machine Learning approach aims to determine what is the best pipeline to train the machine learning model by an automated process, i.e. a code that can test several pipelines and assess which one is the best for the problem. The TPOT(Tree-based Pipeline Optimization Tool) uses the evolutionary computation framework provided by the DEAP (Distributed Evolutionary Algorithms in Python)555deap.readthedocs.io library to choose the best machine learning pipeline.

A typical TPOT pipeline is composed by three main steps: (i) a combination of one or more transformed copies of the features; (ii) a recursive feature elimination; and (iii) a classifier (or a regressor), which could be a single method or a stacking of many.

One of the most simple pipelines in our work was found when we optimized a classifier for Karpenka features, with AUC as metric:

⬇

GradientBoostingClassifier(

learning_rate=0.1, max_depth=8,

max_features=0.85, min_samples_leaf=4,

min_samples_split=11, n_estimators=100,

subsample=0.9

)

This means that steps (i) and (ii) are reduced to the identity, and (iii) becomes a single GradientBoostingClassifier. But when we include $z_{\rm host}$ as a new feature we find a completely different pipeline:

⬇

make_union(

FunctionTransformer(copy),

make_pipeline(

    VarianceThreshold(threshold=0.2),

    FastICA(tol=0.0)

)

),

XGBClassifier(

learning_rate=0.1, max_depth=9,

min_child_weight=2, n_estimators=100,

nthread=1, subsample=0.8

)

Step (i) is now to take two copies of the input data and transform one applying a VarianceThreshold and FastICA in this order, step (ii) is a union of these two copies and (iii) is a XGBClassifier. All the TPOTpipelines are shown in Appendix B.

4 Features Extraction

4.1 Template fitting

The template fitting666Note that Template Fitting refers to methods which use templates of Type Ia light curves to fit other light curves, while Template Matching refers to ones which use templates of different types to classify any light curve. approach consists in reconstruct the light curve of a type Ia supernova, based on templates constructed from previous data sets of SNeIa confirmed by spectroscopic surveys. We fit the light curves computing the flux of the SALT2 model by using the implementation provided by the SNCosmo777sncosmo.readthedocs.io library (sncosmo.readthedocs.io), and sample 1000 points using the multinest method provided by the PyMultinest888github.com/JohannesBuchner/PyMultiNest python library.

In SALT2, the rest-frame specific flux at wavelength $\lambda$ and phase (time) $p$ is modeled by ($p=0$ at $B$-band maximum)

[TABLE]

where $M_{0}(t,\lambda)$ is the average spectral sequence (using past supernova data, Guy et al. (2007)), and $M_{k}(t,\lambda)$, $k>0$, are higher order components related to the supernova variability (terms with $k>1$ are negligible). Finally, $CL(\lambda)$ is the average color correction.

The observer-frame flux in passband $Y$ is calculated as

[TABLE]

where $T^{Y}(\lambda)$ is the transmission curve of the observer-frame filter $Y$, and $z$ is the redshift. The free parameters to be fitted are $x_{0}$, $x_{1}$ and $c$ which are the SED sequence normalization, the stretch and the color parameters respectively, and also $t_{0}$ (the time of B band maximum), and the supernova redshift $z$. So in total we have 5 free parameters, for which we use the priors given in Table 2.

4.2 Parametric models

This approach consists in fitting a functional shape for the light curves for each filter. The advantage over the previous method is that we make less prior assumptions about supernovae. In this work we use the fitted parameters as features for our machine learning algorithms. To fit the light curves we use two models proposed for the SNPCC (Kessler et al., 2010a, b): (a) The Newling model (Newling et al., 2011); (b) the Karpenka model (Karpenka et al., 2012). The parameters for both parametric models were determined by sampling 1000 points with PyMultiNest.

4.2.1 Newling model

The function form proposed by the Newling model is given by

[TABLE]

where

[TABLE]

The peak flux is equal to $A+\psi$, $\phi$ is the explosion starting time, $k$ is related to the relative rise and decay times and $\sigma$ is the light curve width, which gives 5 free parameters. The $\Psi(t)$ is a tail function which ensures that the flux tends to a finite value for $t\gg\tau$, where $\tau=\phi+k\sigma$ is the peak time. The cubic spline in equation (13) is determined to null the derivative in $t=\phi$ and $t=\tau$. The model is fitted separately for each filter, giving a total of 20 features, for the SDSS $griz$ filter system used in the SNPCC. The parameters ranges (which are the same for all bands) are given in Table 3.

4.2.2 Karpenka model

The Karpenka model is similar to Newling’s, generating a light curve with a peak and a tail which tends to a finite value, but can also describe a situation where a second flux peak can happen. The function is given by

[TABLE]

where $t=0$ is taken as the earliest measurement in the r-band light curve. There are 6 free parameters $\{A,B,t_{0},t_{1},T_{\text{fall}},T_{\text{rise}}\}$ for each filter, totalizing 24 free parameters. The ranges of the priors are shown in Table 4, which are the same for all bands.

4.3 Wavelet decomposition

The wavelet decomposition is a very useful tool for signal processing and time series analysis, which is the case of supernovae light curves. Following L16 we perform the à trous wavelet transform, which achieves dyadic scale-invariance, and use the symlet family of wavelets, which are a more symmetric version of the widely used Daubechies family of wavelets, by the implementation sym2 family provided by the PyWavelets library.999pywavelets.readthedocs.io Before using the wavelet decomposition we construct an interpolated version of the light curve, for each band, using the Gaussian Random Process provided by the George library.101010We choose the exponential radial kernel optimized following this procedure: george.readthedocs.io/en/latest/tutorials/hyper

We choose 100 points, evenly spaced in the time range of the light curve, for the Gaussian process curve and a two-level wavelet transform, returning 400 coefficients per filter, or 1600 coefficients per object. To avoid redundancy we perform a dimensionality reduction with Principal Components Analysis (PCA) restricting to only the first 20 components, which we use as our feature set and refer only by Wavelets instead Wavelet decomposition.

4.4 Host galaxy photometric redshift ($z_{\rm host}$) as a new feature

Following the SNPCC we also investigated the effects of introducing the host galaxy photometric redshift, $z_{\rm host}$, as a new feature. So for each feature set we performed the classification twice, both with and without $z_{\rm host}$ in the feature set. For Wavelets, for instance, we first train using the features $({\rm PCA}_{1},{\rm PCA}_{2},...,{\rm PCA}_{20})$, and then a second time using $({\rm PCA}_{1},{\rm PCA}_{2},...,{\rm PCA}_{20},z_{\rm host})$. The only exception is the SALT2 case, where the features are the best fit parameters of two different light curve fits and the redshift can simply be a free parameter in the fit. First, we let the redshift be free and determined only by the light curve, thus using the feature set $(z,t_{0},x_{0},x_{1},c)$. To include $z_{\rm host}$ we refit the light curve fixing $z=z_{\rm host}$ and fitting $(t_{0},x_{0},x_{1},c)$ again, in which case the feature set becomes $(z_{\rm host},t_{0},x_{0},x_{1},c)$.

The value of $z_{\rm host}$ is the photometric one, both for training and predicting. We have tested using the true spectroscopy redshift to train both the models and the photometric prediction. But this resulted in smaller scores than the ones shown below. Even though $z_{\rm host}$ is not exactly the spectroscopic one, the latter will not be measurable. This downgrade in performance when using the exact redshift is likely due to the fact than in ML methods, the closer your training set is to the real data, the better the performance.

5 Results and discussion

5.1 Comparison of ML methods with AUC and AP

Let us first analyze the results of the decision trees ensembles. In Figure 5 we compare the ROC curves and the AUC value (given in per mille in the plots) for each ML model and feature set. The curves and shaded regions around them are the mean and the standard deviation in the catalog set. The AUC values and uncertainties are likewise obtained from this set. We notice that the boosting methods (ADA, GB, and XGB) reach higher scores than the averaging methods (EXT and RF) for any set of features. It is also clear that the different boosting methods exhibit very similar performance for all features. TPOT was not able to provide a significant improvement, except for two cases, wavelets and SALT2 in the case in which we include the $z_{\rm host}$. The Completeness-Purity curves are in turn depicted in Figure 6. We note in general the same difference between boosting and averaging methods, although the latter exhibit broader distributions. An interesting feature which is made clear in this score is that even methods with high AP may result in considerably low completeness values for high purities. This happens in general for Wavelets, Karpenka and Newling, but in some cases using TPOT this issue was alleviated. This means that we may lose a lot of candidates if we need to impose a restrictive purity threshold, which in turn increases the statistical error. For Wavelets (ADA, GB, XGB), for example, if we impose a 90% threshold for the purity, we get a catalog with less than $10\%$ completeness. In contrast, for SALT2, we get around 70% of the real positives even when we impose a 90% purity.

Figure 7 compiles the scores for all cases, where we can confirm the ranking between the methods described above, and also compare the features more clearly. The ranking, from low to high score, is: Karpenka, Newling, Wavelets and SALT2. We can also confirm that the introduction of $z_{\rm host}$ as feature increases the score as expected and already reported in L16, except for the case of Wavelets in which it only made any difference when we used TPOT.

Using AP as score, we found that the TPOT method was able to increase significantly the scores, as can be seen clearly in Figure 6. It was particularly successful for the Newling model and, when $z_{\rm host}$ is included, and also for the Wavelets method. Summarizing the scores for all cases, we find the following.

• •

Comparing the feature set for different methods we find from best to worst: 1. SALT2; 2. Newling and Wavelets (tied); 3. Karpenka (although when using the averaging ML methods Newling performed better). We stress that since the SALT2 model was used to create the lightcurves the results may be biased toward it. Nevertheless as we discuss below the performance of Newling and Wavelets are still very competitive.

• •

Comparing the ML methods for given feature sets, from best to worst: 1. TPOT; 2. XGB; 3. GB; 4. ADA; 5. RF; 6. EXT.

• •

Including $z_{\rm host}$ as a new feature: the scores increased in general, as expected, except for Wavelets where it only increased when using TPOT. However, one should note that the increase in performance obtained by including the host redshift was comparable to the performance difference among ML methods and feature sets. Therefore using TPOT without redshift information still yielded very good results, which is a promising result given that future photometric surveys will only have limited spectroscopic follow-ups.

When we look at the Completeness-Purity curves themselves in Figure 6 and not just the AP score, we note some important details. When we compare the curves for Newling and Wavelets, both without $z_{\rm host}$ as a feature and looking only at the TPOT results we see that if we demand a purity of $95\%$, for instance, Newling reaches a completeness of around $30\%$, significantly larger than that of Wavelets, which lies around $5\%$. Including $z_{\rm host}$ the difference is much smaller but still significant. The final score of Newling is however very similar to the one of Wavelets. This difference in completeness can affect the final cosmological analysis, as we show below. From now on, we focus the analysis only on the catalogs generated by the TPOT pipelines with $z_{\rm host}$ as a feature.

In Table 5 we show a comparison between the AUC scores, in our best case, and other values quoted in the literature. We focused the comparison only between cases where the authors performed the classification on the SNPCC data set and trained their models with 1100 light curves as in this work. Of the 3 quoted studies, only the one conducted by L16 uses the same features as ours. All studies use different classifiers. Only Möller & Deboissì Ere (2018) (using Recurrent Neural Networks) reaches higher scores. One advantage of our method is that Decision Tree Ensembles are easier to train than a Neural Network: it takes much less time and can be done by CPUs without the need for GPUs.

5.2 Analysis of the Bias-Variance tradeoff

Figure 8 shows the evolution of bias forall feature sets in comparison to a completely random classifier and a perfect classifier (which returns all SNeIa with $\mathcal{P}=\mathcal{C}=1$), considering two purity cuts (80% and 90%). A common feature for all cases is that the bias is low for intermediate redshifts ($0.4<z<0.9$), and high for both low and high redshifts ($z<0.4$, and $z>0.9$). One may note that the bias values decrease as the purity increases (from the top panel to the bottom one), while the error bars increase due to the loss in the number of classified objects. This general trend does not always hold because the purity in each bin is not necessarily the same, since we fixed only the value of this quantity for the whole sample. The balance is set by the Bias-Variance tradeoff, described in Section 3.1. Interestingly though, for high purities the bias becomes negative, so the tradeoff actually penalises very high purities not only due to the low completeness, but also due to the higher absolute values of the bias. One would naively expect the bias to be a monotonic function of the purity, but this a priori unexpected behavior of the bias is explained below.

In Figure 9 we depict how each SN type contributes to the magnitudes by showing the bin counts of fitted individual absolute magnitudes ($M_{B,i}\equiv m_{B,i}-\mu_{\rm fid}(z_{i,{\rm sim}})$) for four different purity values (70%, 80%, 90% and 95%). The average of all $M_{B,i}$ is the mean absolute magnitude $M_{B}$ of the SNeIa (see Appendix A). The overall trend is that the type II SN are the first to be excluded as we increase the purity threshold, followed by the type Ibc. This means that the algorithms are more efficient at distinguishing between types I and II than between Ia and CC. Since types II and Ibc are on opposing $M_{B}$ extremes, and type II has a larger magnitude range, as the purity threshold increases the magnitude bias due to CC SNe changes sign, instead of monotonically approaching zero. Interestingly, in the case of SALT2 and Wavelets, for large purity thresholds our classifier ends up mis-classifying mostly the CC SNe which have similar magnitudes to the ones of type Ia. This means that for these feature sets the CC contaminants for high purity thresholds may end up mimicking the type Ia standard candles. This peculiar behavior of the bias is fundamentally due to the ML trying to reduce the 8 classes of SN to a simple binary classification: Ia versus CC.

The dependence of the RMSE with the purity is shown in Figure 10. The green line is the contribution due to the variance, the blue one is the contribution due to the bias and the red thick line is the total RMSE. The dotted gray lines are the RMSE computed in each redshift bin. Note that, close to the minimum, the variance is almost constant while the bias changes sign. I.e. the point at which the bias changes sign corresponds approximately to the purity at which RSME reaches a minimum.

We use this dependence to choose what is the best probability threshold, once it is univocally determined by the purity. The simplest choice is the probability which minimizes the total MSE, a constant purity cut, however we may get very high values of RMSE in some redshift bins. We thus understand that the best way of selecting the sample is to choose the purity values in each redshift bin that minimize the RMSE, which we dub binned purity. Figure 11 shows explicitly that the bias decreases substantially and becomes smaller than $0.03$ for all cases in a large range of redshifts ($0.4<z<1.1$), when we choose this selection criteria.

We therefore construct a new catalog in which each bin is analyzed enforcing the purity threshold that minimizes the MSE of that bin. To find the best method we compute how much each choice affects the final cosmological constraints. In this context the aim is finding the model that fits the data by minimizing

[TABLE]

where $\boldsymbol{\mu}=(\bar{\mu}_{1},\bar{\mu}_{2},...,\bar{\mu}_{11})$, $\boldsymbol{z}=(\bar{z}_{1},\bar{z}_{2},...,\bar{z}_{11})$, and $\boldsymbol{C}$ is the covariance matrix. Defining the error as $\sigma_{i}^{2}={\rm MSE}_{i}=\langle\Delta\bar{\mu}_{i}^{2}\rangle$, where as before $\Delta\bar{\mu}\equiv\bar{\mu}-\bar{\mu}_{fid}(z_{\rm sim})$, we generalize it for correlated bins as $C_{ij}=\langle\Delta\bar{\mu}_{i}\Delta\bar{\mu}_{j}\rangle$ so the covariance for a single catalog is $C_{ij}=\Delta\bar{\mu}_{i}\Delta\bar{\mu}_{j}$, therefore

[TABLE]

In order to assess how the different ML techniques behaved in a cosmological analysis, we computed fits for two models: (i) a flat $\Lambda$CDM model with the current matter density parameter $\Omega_{m0}$ as a free parameter; and (ii) a flat $w$CDM model with the dark energy equation of state parameter $w$ running free while fixing $\Omega_{m0}=0.3$. In both cases we marginalize over the overall distance scale which is a combination of the unknown SN absolute magnitude and the Hubble parameter $H_{0}$. The constraints on $\Omega_{m0}$ and $w$ are shown in Figure 12, where we compare the selection criteria (constant or binned purity thresholds) and the different ML feature sets. We also indicate the result for a perfect classifier (i.e., composed of all the 4800 SNeIa) in red and the fiducial values used in SNPCC ($\Omega_{m0}=0.3$ and $w=-1$) in black. All catalogs classified by ML show small biases when compared to the perfect classifier or the fiducial values (less than half $\sigma$). We also note that using the binned purity yielded in general slightly improved constraints, but the difference was small. Finally, the SALT2 feature set gives the best results in all cases, approaching the ideal classifier in the binned purity case. On the other hand on average the Karpenka model lead to the worst results, with error bars over twice as large as the ideal case.

We remark that since the small bias on the perfect classifier is smaller than the error bars, it is consistent with random noise in the simulation data. As it is well-known in the literature (Kessler & Scolnic, 2017; Kessler et al., 2019), the combination of lightcurve fitting, selection effects and the intrinsic scatter of SNeIa introduces bias in the cosmological analysis that must be corrected for. Bias removal is thus a very important aspect of supernova cosmology, especially in magnitude-limited surveys. Here however we remove the light-curve fitting biases in a simple way, which we explain in more detail in Appendix A.

If we assume that all Type Ia supernovae contribute with the same amount of information, the effective statistical error on a cosmological parameter should scale as $\sigma_{\mathrm{Ia}}=\sigma_{\mathrm{I}}/\sqrt{N_{\mathrm{Ia}}}$, where $\sigma_{\mathrm{I}}$ is the standard deviation of the Ia supernovae. We extrapolate on this idea to define $N_{\mathrm{Ia}}^{\rm eff}$, the effective number of SNeIa in a final photometric classified catalog. The effective statistical error for this catalog is thus $\sigma=\sigma_{\mathrm{I}}/\sqrt{N_{\mathrm{Ia}}^{\rm eff}}$. We likewise define the effective completeness $\mathcal{C}^{\rm eff}$ as

[TABLE]

This quantity is computed in Table 6 for different methods. SALT2 preserves up to $75\%$ of the information in the best scenarios, Newling and Wavelets around $30-35\%$ and Karpenka only $20\%$. This means that even if one drops SALT2 since our SNe may be biased toward it, one still retains around a third of all the information in the SNeIa data without spectroscopy.

To check the validity of our fit we also computed the value of the so-called $R^{2}$ score (or coefficient of determination111111scikit-learn.org/stable/modules/generated/sklearn.metrics.r2_score.html). In all cases $R^{2}$ was found to be very close to unity ($|1-R^{2}|\sim 10^{-5}$).

6 Conclusions

In this work we used simulated data from the Supernova Photometric Classification Challenge (SNPCC) to assess the performance of a series of machine learning methods applied to the problem of supernova classification. We considered some flavors of ensemble Decision Trees: Adaptive Boosting (ADA), RandomForest (RF), ExtraTrees (EXT), Gradient Boosting (GB), and the Extreme Gradient Boosting (XGB) provided by the XGBoost library. We also tested an alternative metric to characterize the performance of the classifier, the Average Purity (AP), and compared the results with the AUC used in L16. In addition we performed an automated machine learning training by using TPOT, in order to find the best pipelines for each feature set. Regarding the features, we tested four different models: Karpenka, Newling, wavelets and SALT2. We then explored how a ML-classified SN catalog performs in cosmological parameter estimation, quantifying the degradation due to the imperfect completeness and purity of the resulting catalogs.

Among the feature sets, the SALT2 features gave the best overall results. Considering that this model was used to generate part of the simulated data, it is not clear to what extent the data itself was possibly biased toward this model, and this question certainly deserves to be investigated in more detail in the future. If one thus neglects the SALT2 feature set, both Wavelets and Newling models performed similarly. The Karpenka model performed consistently the worst.

Regarding the scores, we found that the AP was more effective to discriminate between the ML methods as opposed to the AUC, which gave essentially the same results for all boosting methods and TPOT. The boosting methods were considerably better than the averaging ones for both metrics. We showed that, considering the AP score, TPOT was capable to provide a pipeline more efficient than the boosting methods tested in this work.

We tested two ways for selecting the purity probability threshold: either imposing a constant purity for the whole sample or choosing different purity levels in each redshift bin in order to optimize the MSE. Both choices yielded good results for a flat $\Lambda$CDM model. Employing the Bias-Variance tradeoff we were able to consistently produce photometric classified catalogs that were within half-$\sigma$ of the fiducial model and of a perfect classifier.

One of the most important open questions for future photometric SN surveys is how much cosmological information is preserved after classification. We thus tested this using each feature set with the best ML method, TPOT. We found out that SALT2 is able to preserve up to $75\%$ of the information, while the other second-best alternatives, based on either the Newling model or Wavelet decomposition, around $30-35\%$. This is a very important quantity (it can be interpreted as an effective completeness), as future surveys like LSST will have a huge number of observed SN, but the majority without any spectra. Table 6 summarizes how much information each fit preserves, compared to a perfect classifier.

It is important to stress that an analysis based on the Bias-Variance Trade-off assumes that the bias cannot be modeled. In fact, as CC SNe lightcurves become better understood it may become possible to predict their impact on the distance estimation and thus mitigate the bias they introduce. This would allow even better performance of the ML photometric classification techniques.

Automated machine learning training techniques are an interesting and very efficient way of attacking the problem of supernova photometric classification for cosmological uses, and we recommend their application in future ML studies in astronomy.

Acknowledgements

This work made use of the CHE cluster, managed and funded by COSMO/CBPF/MCTI, with financial support from FINEP and FAPERJ, and of the Milliways cluster at UFRJ, funded by FAPERJ and CNPq. We would like to thank Michelle Lochner for sharing her code, Martin Makler for access to the CHE cluster and Luca Amendola and Luciano Casarini for useful suggestions. MVS is supported by the Brazilian research agency FAPERJ. MQ is also supported by FAPERJ and by the Brazilian research agency CNPq.

Data availability

No new data were generated or analysed in support of this research.

Appendix A Distance estimation by SALT2

In order to make the residual plots and perform the cosmological constraints we use the SALT2 (Guy et al., 2007) as distance estimator, in which the distance modulus for a single supernovae is given by

[TABLE]

where $m_{B}^{*}$, $x_{1}$ and $c$ are light curve fitting parameters. In particular, $m_{B}^{*}$ and $x_{0}$ are related by

[TABLE]

The parameters $\alpha$, $\beta$ must be fitted along with the cosmological ones. $M_{B}$ is often left as a free parameter to be marginalized together with $H_{0}$. In our analysis we fixed the parameters to the fiducial values used in the SNPCC simulations, to wit: $\alpha=0.11$, $\beta=3.2$. We also assumed throughout the simple fiducial model ($\mu_{fid}$) used in SNPCC: flat $\Lambda$CDM with $\Omega_{m0}=0.3$ and $H_{0}=70\,\mathrm{km/(s\,Mpc)}$.

Based on the assumption that SNeIa are standard candles, $M_{B}$ in principle should also be a constant and assigned the value used by SNPCC ($M_{B}=-19.365$). However, as discussed in detail by Kessler & Scolnic (2017); Kessler et al. (2019), in any magnitude-limited survey the light-curve fitting itself introduces distance biases, which must be corrected for. These are much larger and independent of usual biases such as the Malmquist bias, and are also independent of selection criteria. The SNPCC data also has additional sources of scatter introduced in the simulations. Together, these effects introduce important biases in the Hubble diagram.

In this work, as we are not particularly concerned with the usual light-curve fitting biases or the particularities of the SNPCC simulation, we chose not to implement any sophisticated bias removal methods such as the BBC (BEAMS121212BEAMS : “Bayesian Estimation Applied to Multiple Species” (Kunz et al., 2007) with Bias Corrections) method of Kessler & Scolnic (2017) or a detailed modelling of extra scatter sources. Instead we simply remove by hand a priori the bias in $\mu$ by using different values of $M_{B}$ in each of the $\Delta z=0.1$ redshift bins. This is achieved by making use of the pure SNeIa catalog and the $\mu(z)$ values of the fiducial cosmological model to measure $M_{B}(z)$ in each redshift bin using (18). The difference between the measured $M_{B}(z)$ value and the fiducial value of $M_{B}$ used in SNPCC is thus the bias in each bin, which we subsequently remove. Figure 13 illustrates the fitted values of $M_{B}$, together with their uncertainties. This procedure guarantees that any cosmological bias in the classified catalogs would be a result of the ML classification code, or of problems in the Bias-Variance tradeoff implementation.

Appendix B TPOT pipelines

Here we provide the best pipelines obtained with TPOT for each feature set and optimization criterion. They are ready to use as argument of the function make_pipeline from sklearn.pipeline module. The used functions are from the libraries listed in table 7.

• •

Karpenka, optimized by AP:

⬇

make_union(

make_pipeline(

    PolynomialFeatures(

        degree=2, include_bias=False,

        interaction_only=False

    ),

    SelectFwe(

        score_func=f_classif,

        alpha=0.033

    )

),

FunctionTransformer(copy)

),

GradientBoostingClassifier(

learning_rate=0.1, max_depth=9,

max_features=0.25, min_samples_leaf=8,

min_samples_split=12, n_estimators=100,

subsample=0.9

)

• •

Karpenka, optimized by AUC:

⬇

GradientBoostingClassifier(

learning_rate=0.1, max_depth=8,

max_features=0.85, min_samples_leaf=4,

min_samples_split=11, n_estimators=100,

subsample=0.9

)

• •

Karpenka $+z_{\rm host}$, optimized by AP:

⬇

ZeroCount(),

GradientBoostingClassifier(

learning_rate=0.1, max_depth=10,

max_features=1.0, min_samples_leaf=8,

n_estimators=100, subsample=0.75

)

• •

Karpenka $+z_{\rm host}$, optimezed by AUC:

⬇

make_union(

FunctionTransformer(copy),

make_pipeline(

    VarianceThreshold(threshold=0.2),

    FastICA(tol=0.0)

)

),

XGBClassifier(

learning_rate=0.1, max_depth=9,

min_child_weight=2, n_estimators=100,

nthread=1, subsample=0.8

)

• •

Newling, optimized by AP:

⬇

make_union(

make_union(

    FunctionTransformer(copy),

    FunctionTransformer(copy)

),

make_pipeline(

    Normalizer(norm="l2"),

    RFE(estimator=ExtraTreesClassifier(

        criterion="entropy",

        max_features=0.75,

        n_estimators=100

        ), step=0.15

    ),

    RFE(

        estimator=ExtraTreesClassifier(

            criterion="gini",

            max_features=0.05,

            n_estimators=100

        ),

        step=0.95

    ),

    FastICA(tol=0.55),

    RBFSampler(gamma=0.8),

    RFE(

        estimator=ExtraTreesClassifier(

            criterion="entropy",

            max_features=0.15,

            n_estimators=100

        ),

        step=0.4

    )

)

),

ExtraTreesClassifier(

bootstrap=False, criterion="entropy",

max_features=0.95, min_samples_leaf=2,

min_samples_split=5, n_estimators=100

)

• •

Newling, optimized by AUC:

⬇

Normalizer(norm="l2"),

StackingEstimator(

estimator=LinearSVC(

    C=5.0, dual=False,

    loss="squared_hinge",

    penalty="l1", tol=0.01

)

),

GradientBoostingClassifier(

learning_rate=0.1, max_depth=9,

max_features=0.85, min_samples_leaf=14,

min_samples_split=9, n_estimators=100,

subsample=0.85

)

• •

Newling $+z_{\rm host}$, optimized by AP:

⬇

make_union(

FunctionTransformer(copy),

FunctionTransformer(copy)

),

MinMaxScaler(),

RFE(

estimator=ExtraTreesClassifier(

    criterion="entropy",

    max_features=0.3,

    n_estimators=100

),

step=1.0

),

PolynomialFeatures(

degree=2, include_bias=False,

interaction_only=False

),

SelectFwe(

score_func=f_classif, alpha=0.004

),

ExtraTreesClassifier(

bootstrap=False, criterion="entropy",

max_features=0.25, min_samples_leaf=2,

min_samples_split=3, n_estimators=100

)

• •

Newling $+z_{\rm host}$, optimized by AUC:

⬇

RFE(

estimator=ExtraTreesClassifier(

    criterion="entropy",

    max_features=0.75,

    n_estimators=100

),

step=0.2

),

ExtraTreesClassifier(

bootstrap=False,

criterion="gini",

max_features=0.95,

min_samples_leaf=1,

min_samples_split=2,

n_estimators=100

)

• •

SALT2, optimized by AP:

⬇

StackingEstimator(

estimator=LinearSVC(

    C=20.0, dual=False,

    loss="squared_hinge",

    penalty="l2", tol=0.0001

)

),

StackingEstimator(

estimator=LinearSVC(

    C=10.0, dual=True,

    loss="squared_hinge",

    penalty="l2", tol=0.1

)

),

ZeroCount(),

RandomForestClassifier(

bootstrap=False, criterion="entropy",

max_features=0.35, min_samples_leaf=3,

min_samples_split=9, n_estimators=100

)

• •

SALT2, optimized by AUC:

⬇

make_union(

FunctionTransformer(copy),

make_union(

    make_union(

        make_pipeline(

            RBFSampler(gamma=0.35),

            Normalizer(norm="l1")

        ),

        FunctionTransformer(copy)

    ),

    FunctionTransformer(copy)

)

),

RandomForestClassifier(

bootstrap=False, criterion="entropy",

max_features=0.2, min_samples_leaf=2,

min_samples_split=7, n_estimators=100

)

• •

SALT2 $+z_{\rm host}$, optimized by AP:

⬇

make_union(

FunctionTransformer(copy),

make_union(

    make_pipeline(

        FastICA(tol=0.65),

        FastICA(tol=0.0),

        ZeroCount()

    ),

    FunctionTransformer(copy)

)

),

RandomForestClassifier(

bootstrap=False, criterion="entropy",

max_features=0.4, min_samples_leaf=1,

min_samples_split=4, n_estimators=100

)

• •

SALT2 $+z_{\rm host}$, optimized by AUC:

⬇

make_union(

FunctionTransformer(copy),

FastICA(tol=0.05)

),

StackingEstimator(

estimator=LogisticRegression(

    C=0.5, dual=False,

    penalty="l1"

)

),

RandomForestClassifier(

bootstrap=True, criterion="entropy",

max_features=0.2, min_samples_leaf=1,

min_samples_split=2, n_estimators=100

)

• •

Wavelets, optimized by AP:

⬇

PolynomialFeatures(

degree=2, include_bias=False,

interaction_only=False

),

StackingEstimator(

estimator=ExtraTreesClassifier(

    bootstrap=True,

    criterion="gini",

    max_features=0.25,

    min_samples_leaf=7,

    min_samples_split=17,

    n_estimators=100

)

),

StackingEstimator(

estimator=LogisticRegression(

    C=1.0, dual=False, penalty="l1"

)

),

MaxAbsScaler(),

RFE(

estimator=ExtraTreesClassifier(

    criterion="entropy",

    max_features=0.75,

    n_estimators=100

),

step=0.1

),

StackingEstimator(estimator=GaussianNB()),

MinMaxScaler(),

XGBClassifier(

learning_rate=0.1, max_depth=4,

min_child_weight=1, n_estimators=100,

nthread=1, subsample=0.6

)

• •

Wavelets, optimized by AUC:

⬇

make_union(

make_pipeline(

    Normalizer(norm="l2"),

    FastICA(tol=0.7)

),

FunctionTransformer(copy)

),

StackingEstimator(

estimator=ExtraTreesClassifier(

    bootstrap=False,

    criterion="entropy",

    max_features=0.9,

    min_samples_leaf=3,

    min_samples_split=6,

    n_estimators=100

)

),

RobustScaler(),

ZeroCount(),

LogisticRegression(

C=0.5, dual=True,

penalty="l2"

)

• •

Wavelets $+z_{\rm host}$, optimized by AP:

⬇

make_union(

Nystroem(

    gamma=0.15, kernel="linear",

    n_components=10

),

StackingEstimator(

    estimator=make_pipeline(

        StackingEstimator(

            estimator=LogisticRegression(

                C=25.0,

                dual=False,

                penalty="l1"

            )

        ),

        PolynomialFeatures(

            degree=2,

            include_bias=False,

            interaction_only=False

        ),

        LogisticRegression(

            C=0.001, dual=False,

            penalty="l2"

        )

    )

)

),

GradientBoostingClassifier(

learning_rate=0.1, max_depth=8,

max_features=0.8, min_samples_leaf=4,

min_samples_split=10, n_estimators=100,

subsample=0.9

)

• •

Wavelets $+z_{\rm host}$, optimized by AUC:

⬇

make_union(

StackingEstimator(

    estimator=KNeighborsClassifier(

        n_neighbors=21, p=1,

        weights="uniform"

    )

),

StackingEstimator(

    estimator=make_pipeline(

        RBFSampler(gamma=0.1),

        ExtraTreesClassifier(

            bootstrap=False,

            criterion="entropy",

            max_features=0.25,

            min_samples_leaf=18,

            min_samples_split=7,

            n_estimators=100

        )

    )

)

),

MaxAbsScaler(),

GradientBoostingClassifier(

learning_rate=0.1, max_depth=7,

max_features=0.9, min_samples_leaf=12,

min_samples_split=9, n_estimators=100,

subsample=0.85

)