Searching for Exoplanets Using Artificial Intelligence

TL;DR

This paper introduces a convolutional neural network approach for detecting exoplanets in large datasets, outperforming traditional methods and enabling more accurate, automated analysis of stellar light curves.

Contribution

The paper presents a novel deep learning method that learns to recognize exoplanet features directly from data, improving detection accuracy over existing algorithms.

Findings

Deep nets outperform least-squares methods in exoplanet detection.

The neural network generalizes well across different time series data.

Periodic transits are detected accurately without explicit model fitting.

Abstract

In the last decade, over a million stars were monitored to detect transiting planets. Manual interpretation of potential exoplanet candidates is labor intensive and subject to human error, the results of which are difficult to quantify. Here we present a new method of detecting exoplanet candidates in large planetary search projects which, unlike current methods uses a neural network. Neural networks, also called "deep learning" or "deep nets" are designed to give a computer perception into a specific problem by training it to recognize patterns. Unlike past transit detection algorithms deep nets learn to recognize planet features instead of relying on hand-coded metrics that humans perceive as the most representative. Our convolutional neural network is capable of detecting Earth-like exoplanets in noisy time-series data with a greater accuracy than a least-squares method. Deep nets are highly generalizable allowing data to be evaluated from different time series after interpolation without compromising performance. As validated by our deep net analysis of Kepler light curves, we detect periodic transits consistent with the true period without any model fitting. Our study indicates that machine learning will facilitate the characterization of exoplanets in future analysis of large astronomy data sets.

Tables

Table 1. Table 1: Summary of Parameters

Training Parameters		Values	Unit
Transit Depth	R${}_{}{}^2{}_p{}^{}$/R${}_{}{}^2{}_s{}^{}$	200,500,1000,2500,5000,10000	ppm
Orbital Period	P	2,2.5,3,3.5,4	days
Inclination	i	86,87,90	degrees
Noise Parameter	${\sigma }_{tol}$	1.5,2.5,3.5,4.5
Phase Offset	$\varphi$	0, $\pi$/3, 2$\pi$/3, $\pi$
Wave Amplitude	A	250,500,1000,2000	ppm
Wave Period	$\omega$	6./24, 12./24, 24./24	days
Amplitude Variability Period	$P_A$	-1, 1, 100	days
Wave Variability Period	$P_{\omega }$	-3, 1, 100	days
Sensitivity Test		Values	Unit
Noise Parameter	${\sigma }_{tol}$	0.25,0.5,0.75,1,1.25,1.5,1.75,2.25,2.5,2.75,3
Fixed Transit Parameters
Scaled Semi-major Axis	a/Rs	12
Eccentricity	$e$	0
Arg. of Periastron	$\Omega$	0
Linear Limb Darkening	$u_1$	0.5
Mid Transit Time	$t_{mid}$	1
Window Size		360	minutes
Cadence	dt	2	minutes
Minimum Time	$t_{min}$	0.875	days
Maximum Time	$t_{max}$	1.125	days
Training Samples	311040
Test Samples	933120

Table 2. Table 2: Summary of Classifiers

	BLS	SVM	MLP	CNN 1D	Wavelet MLP
Input features	180	180	180	180	182
Trainable Params	3	181	13,937	17,293	14,169
Layers		1	4	5	4
Total Neurons		1	105	109	105
Neural Connections		1	2492	2544	2494
Training Accuracy (%)	73.51	91.08	99.72	99.60	99.77
Training False Pos.(%)	22.34	3.05	0.08	0.21	0.08
Training False Neg.(%)	4.10	5.85	0.20	0.19	0.15
Sensitivity Test (%)	63.14	83.10	88.73	88.45	91.50
Test False Pos. (%)	31.58	2.92	0.29	0.25	0.29
Test False Neg. (%)	5.37	13.98	10.97	11.29	11.22

Equations

$$t^{\prime }=t-t_{min}$$

$$A(t^{\prime })=A+A\sin \left(\frac{2\pi t^{\prime }}{P_A}\right)$$

$$\omega (t^{\prime })=\omega +\omega \sin \left(\frac{2\pi t^{\prime }}{P_{\omega }}\right)$$

$$F_{transit}(t)\ast \mathcal{N}\left(\frac{R_p^2}{R_s^2}/{\sigma }_{tol}\right)\ast \left(1+A(t^{\prime })\sin \left(\frac{2\pi t^{\prime }}{\omega (t^{\prime })}+\varphi \right)\right)$$

$$\mathcal{N}\left(\frac{R_p^2}{R_s^2}/{\sigma }_{tol}\right)\ast \left(1+A(t^{\prime })\sin \left(\frac{2\pi t^{\prime }}{\omega (t^{\prime })}+\varphi \right)\right)$$

$$\mathrm{Output}=0\mathrm{if}\sum _j{w}_j{x}_j+b\le 0$$

$$\mathrm{Output}=\sum _j{w}_j{x}_j+b\mathrm{if}\sum _j{w}_j{x}_j+b>0$$

$$P(\hat{y}|X)={\hat{y}}^y(1-\hat{y}{)}^{(1-y)}$$

$$-\log P(\hat{y}|X)=-y\log \hat{y}-(1-y)\log (1-\hat{y})$$

$$Loss=-\frac{1}{M}\sum _i^M{y}_i\log {\hat{y}}_i+(1-y_i)\log (1-{\hat{y}}_i)$$

$$\begin{array}{rl}\text{ Output Activation: }& f(\text{w}\cdot \text{X}+b)=\frac{1}{1+e^{-(\text{w}\cdot \text{X}+b)}}\end{array}$$

$$W^{i+1}=W^i-ϵ(\nabla Los{s}_W^i+\eta \nabla Los{s}_W^{i+1})$$

$$s(t)=\sum _{a}x(t-a)w(a)$$

Full text

††pubyear: 2017††pagerange: Searching for Exoplanets using Artificial Intelligence–2

1 Introduction

Transiting exoplanets provide a remarkable opportunity to detect planetary atmospheres through spectroscopic features. During primary transit, when a planet passes in front of its host star, the light that transmits through the planet’s atmosphere reveals absorption features from atomic and molecular species. Currently 3,513 exoplanets have been discovered from space missions (Kepler (Borucki et al., 2010), K2 (Howell et al., 2014) and CoRoT (Auvergne et al., 2009)) and from the ground (HAT/HATnet (Bakos et al., 2004), SuperWASP (Pollacco et al., 2006), KELT (Pepper et al., 2007) ). Future planet hunting surveys like TESS, PLATO and LSST plan to increase the thresholds that limit current photometric surveys by sampling brighter stars at faster cadences and over larger field of views (LSST Science Collaboration et al. 2009; Ricker et al. 2014; Rauer et al. 2014). Kepler’s initial four-year survey revealed $\sim$15$\%$ of solar type stars have a 1–2 Earth-radius planet with an orbital period between 5–50 days (Fressin et al. 2013; Petigura et al. 2013). The detection of such small Earth-sized planets are difficult because the transit depth, $\sim$100 ppm for a solar type star, reaches the limit of current photometric surveys and is below the average stellar variability (see Figure 1). Stellar variability is present in over 25$\%$ of the 133,030 main sequence Kepler stars and ranges between $\sim$950 ppm (5th percentile) and $\sim$22,700 ppm (95th percentile) with periodicity between 0.2 and 70 days (McQuillan et al., 2014). The analysis of data in the future needs to be both sensitive to Earth-like planets and robust to stellar variability.

Common techniques to find planets maximize the correlation between data and a simple transit model via a least-squares optimization, grid-search, or matched filter approach (Kovács et al. 2002; Jenkins et al. 2002; Carpano et al. 2003; Petigura et al. 2013). A least-squares optimization aims to minimize the mean-squared error (MSE) between data and a model. Since the transit parameters are unknown a priori, a simplified transit model is constructed with a box function. Least-square optimizers are susceptible to finding local minima (see Figure 2) when trying to minimize the MSE and, thus, can result in inaccurate transit detections unless the global solution can be found. When individual transit depths are below the scatter, as is the case for Earth-like planets currently, constructively binning the data can increase the signal-to-noise (SNR). Grid-searches utilize binning by performing a brute-force evaluation over different periods, epochs and durations to search for transits either with a Least-squares approach (Kovács et al., 2002); or matched-filter (Petigura et al., 2013). A matched filter approach tries to optimize the signal of a transit by convolving the data with a hand-designed kernel/filter to accentuate the transit features.

The SNR of a transit detection can be maximized when convolving data with the optimal filter. However, solving for the optimal filter in the case of varying transit shapes cannot be done analytically, so kernels are hand-made to approximate what the human user thinks is best. Deep learning with convolutional neural networks (CNN) have previously been used to solve similar kernel optimization problems (Krizhevsky et al., 2012). Recently, these CNNs have been able to beat human perception at object recognition (He et al. 2015; Ioffe & Szegedy 2015). Other heuristics at automating planet finding have been developed using Random Forests (McCauliff et al. 2015; Mislis et al. 2016), Self-Organizing Maps (Armstrong et al., 2016), and k-nearest neighbors (Thompson et al., 2015). These algorithms are used to remove a substantial fraction of false positive signals prior to detecting the planet’s signal. The fraction of false positives can be reduced by requiring at least three self-consistent transit detections instead of a single one (Petigura et al. 2013; Burke et al. 2014; Rowe et al. 2015; Coughlin et al. 2016).

Complications can arise if the input data is not quiescent where the instrumental systematics or stellar variability are marginal compared to the signal of the transit and change with time. De-trending the data from these effects requires a prior information, where an assumption is made about the behavior of the effects perturbing the data. Techniques to de-trend light curves from instrumental or stellar systematics have been successful in the past using Gaussian Processes (e.g. Gibson et al. 2012; Aigrain et al. 2015; Crossfield et al. 2015), PCA (e.g. Zellem et al. 2014), ICA (e.g. Waldmann et al. 2013; Morello et al. 2015), wavelet based approaches (e.g. Carter & Winn 2009), or more classical aperture photometry de-trending (e.g. Armstrong et al. 2014). Typically these methods involve modeling the systematic trend and subtracting it from the data before characterizing the exoplanet signal. This disjointed approach may allow the transit signal to be partially absorbed into the best-fit stellar variability or instrument model. In this way, the methods may be over-fitting the data and making each transit event appear shallower. Foreman-Mackey et al. (2015) propose an alternate technique that models the transit signal simultaneously with the systematics and find a large improvement on planet detection but at the cost of computational efficiency.

The ideal algorithm for detecting planets should be fast, robust to noise and capable of learning and abstracting highly non-linear systems. A neural network trained to recognize planets with simulated data provides the ideal platform. Deep learning with a neural network is a computational approach at modeling the biological way a brain solves problems using collections of neural units linked together (Rosenblatt 1958; Newell 1969;). Deep nets are composed of layers of “neurons”, each of which are associated with different weights to indicate the importance of one input parameter compared to another. Our neural network is designed to make decisions, such as whether or not an observation detects a planet, based on a set of input parameters that treat, e.g. the shape and depths of a light curve, the noise and systematic error, such as star spots. The discriminative nature of our deep net can only make a qualitative assessment of the candidate signal by indicating the likelihood of finding a transit within a subset of the time series. Within a probabilistic framework, this is done by modeling the conditional probability distribution $P(y|x)$ which can be used for predicting $y$ (planet or not) from $x$ (photometric data). The advantage of a deep net is that it can be trained to identify very subtle features in large data sets. This learning capability is accomplished by algorithms that optimize the weights in such a way as to minimize the difference between the output of the deep net and the expected value from the training data. Deep nets have the ability to model complex non-linear relationships that may not be derivable analytically. The network does not rely on hand designed metrics to search for planets, instead it will learn the optimal features necessary to detect a transit signal from our training data.

Neural networks have been used in a few planetary science applications including multi-planet prediction and atmospheric classification. Kipping & Lam (2017) train a neural network to predict multi-planet systems. While this neural network does not focus on directly detecting transits, it models the correlation between orbital period, planetary radius and mass to predict the presence of additional bodies in a planetary system. Data from TESS is limited to short period planets (i.e. $P$ < 13.7 days) but predictions about longer period planets can be modeled from correlations in Kepler data. Neural networks provide an effective tool to model such non-linear systems when no analytic formula is known or even possibly derivable. Another application of deep learning in exoplanet science predicts atmospheric compositions based on emission spectra (Waldmann, 2016). Additionally, George & Huerta 2017 introduces a new method for time-series classification and regression of highly noisy gravitational waves using raw time-series inputs in a 1D convolutional neural network.

In this paper we design various deep learning algorithms to recognize planetary transit features from a training data set. In Section 3, we explain the architecture of our deep learning algorithms. In Section 4, we explore the sensitivities of each algorithm at detecting planets in noisy data. Section 5 covers the time series evaluation of multiple transits and validates the detection algorithm against known planets in Kepler data. In section 6 and 7, we discuss our results and summarize our findings in the conclusion.

2 Implementation and Training

Simulated training data is used to teach our deep nets how to predict single planetary transits in noisy photometric data. The simulated data is similar to what we would expect from a real planetary search survey. After the deep nets are trained, we use the network to assess the likelihood of potential planetary signals in data it has not seen before.

2.1 Training Data Set

We generate a total of 311040 transits and non-transit data to train our deep nets. The training data are computed from a discretely sampled 9 dimensional hypergrid in our parameter space (see Table 1). The parameters limit transit duration to 30 minutes and no longer than 4 hours, which spans 1./12 to 2./3 of the time domain. We chose our parameters to mimic data from real search surveys by encompassing many possible systematic shapes and transit sizes. The equation used to generate our noisy data with a quasi-periodic systematic trend is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $F_{transit}(t)$ is the transit function given by Mandel & Agol 2002, $A$ is the amplitude of our simulated stellar variability, $\omega$ is the period of oscillation, $\phi$ is the phase shift, and $\mathcal{N}$ is a Gaussian distribution used to generate random numbers with a mean of 1 and standard deviation of ($R_{p}^{2}$/R${}_{s}^{2}$)/ $\sigma_{tol}$ and $R_{p}^{2}$/R${}_{s}^{2}$ is the normalized radius ratio between the planet and star. Non-transit noisy data is generated in a similar manner, but without the transit signal, $F_{transit}(t)$.

The simulated data has a quasiperiodic systematic trend, like signals found in Kepler data (Aigrain et al., 2015). The simulated variability represents a sinusoid with varying frequency and amplitude given our choice of parameters, $P_{A}$ and $P_{\omega}$. We generate a standard sinusoidal model with no variation in amplitude or frequency when $P_{A}$ and $P_{\omega}$ are equal to 1000 days. The amplitude of the sinusoidal function will double over the span of 6 hours when $P_{A}$ is equal to 1 day. The amplitude of the sinusoidal function will diminish to zero when $P_{A}$ is equal to -1 day. It is physically unrealistic to have a negative period, however we use it as a mathematical tool for achieving the shape we want. The oscillation period will double when $P_{\omega}$=1 day and will diminish by half when $P_{\omega}$=-3. A shape analysis for each varying parameter in our systematic trend is shown in Figure 14.

[TABLE]

Each transit light curve has a non-transit sample using the same systematic parameters but newly generated noise of the same distribution shape and size. This allows our deep net to differentiate between transit and non-transit. The synthetic data are normalized to unit variance and have the mean subtracted off prior to input in the deep nets. Various light curves and systematic trends are shown in Figure 3. The variable light curve shapes and systematic trends are shown in appendix Figure 4.

2.2 Neural Network Architecture

An artificial neural network is a computational approach at modeling the biological way a brain solves problems. It uses a collection of neural units connected to many others and has links/synapses which can be enforced or inhibited through the activation state. Networks are composed of layers of “neurons” (see Figure 5), generally referred to as Restricted Boltzmann Machines (RBM), each of which has a set of input parameters, each of which are associated with a weight to indicate the importance of one input parameter compared to another. The circles in the graph represent one unit or neuron in the network and $h$ represents a transformation the neuron performs on the input data to abstract the information for the next layer. A neuron also has a bias, which acts like a threshold number for the final decision of that neuron, whether a yes/no decision or a probability. A fully connected neural network is one where each neuron in a higher layer uses the output from every neuron in the previous layer as the input. We implement a fully connected multi-layer perceptron for each of our deep nets.

For example, a perceptron is a type of neuron that takes several inputs, $x_{1}$, $x_{2}$,…, (e.g. photometric data) and produces a single output. Weights, $w_{1}$, $w_{2}$,…, are introduced to express the importance of the respective inputs to the output. The neuron’s output is determined by whether the weighted sum $\sum_{j}w_{j}x_{j}$ is less than or greater than some threshold value, which is quantified by the bias, $b$. Just like the weights, the threshold is a real number, which is a parameter of the neuron. To put it in more precise algebraic terms:

[TABLE]

[TABLE]

Using simplified notation, $\sum_{j}w_{j}x_{j}$ = $\bf{w}\cdot\bf{X}$. The neurons transform the input data with a weighted sum and then use a non-linear function as the “activation” that defines the neuron output. The activation function we use is the rectified linear unit (reLU) and it zeros the output of a neuron if it has a negative contribution to the next layer (Nair & Hinton, 2010). The ReLU function helps solve the “exploding/vanishing gradient” problem because it is a non-saturated activation function. Saturated activation functions (e.g. sigmoid, tanh) have an undesirable property of producing a derivative close to zero when the activation saturates either tail of 0 or 1. This small gradient will reduce information from being propagated back to the previous layers’ weights and prevent the network from learning efficiently (Krizhevsky et al., 2012). Additionally, we initialize the weights following a method in He et al. (2015) found to help networks (e.g. 30 convolutional/fully connected layers) converge and prevent saturation. The weight initialization draws random numbers from a Gaussian distribution centered on 0 with the standard deviation equal to $\sqrt{2/N}$ where $N$ is the number of input features.

The advantage of a neural network is that it can be trained to identify subtle features inherent in a large data set. This learning capability is accomplished by allowing the weights and biases to vary in such a way as to minimize the difference (i.e. the cross-entropy) between the output of the neural network and the expected or desired value from the training data. The cross-entropy is used in our classifier over the mean squared error because it better represents the error in a binary classifier (transit vs. non-detection). We are interested in separating the input data by a line (or plane) which allow us to make predictions based on whether the input data falls above or below the dissecting surface. Suppose we have a model that predicts 2 classes with a ground truth (correct) label as $y$ and predicted probability $\hat{y}$ which is from the output of the neural network. The probability of predicting the correct value from our data for a given sample based on our two classes would be

[TABLE]

Taking the logarithm and changing sign yields

[TABLE]

Summing Equation 4 over all our samples (i.e. training samples) and then dividing by the number of samples, $M$, yields the cross-entropy for our system,

[TABLE]

In effect we are minimizing the Loss Function, which tells us how well we are achieving our goal. The cross entropy is related to the expectation of the logarithmic difference between predictions and true labels where the expectation is taken using the true probabilities. As the $Loss$ goes to 0, the closer we are to predicting the true labels of the data, $y$.

The last layer in our network is made of only one neuron that is slightly different than the previous layer neurons because it classifies the input data into a transit or non-transit detection. The output of the last layer is modified with a sigmoid function (see equation 6 and Figure 6). The sigmoid function is an “activation function”, $f$, for the last neuron only and it will output a probability that the input falls into one of the two classes (transit or not). The input into the sigmoid function is the weighted sum of the input parameters (i.e. the output of the previous layer’s neurons).

[TABLE]

Optimization of the weights, $w_{j}$, for our deep net is done by minimizing the loss function of our system, the cross-entropy. The weights in each layer are optimized using a backward propagation scheme (Werbos, 1974). The weights are randomly initialized and after the loss function is computed, a backward pass propagates it from the output layer to the previous layers providing each weight parameter with an update value meant to decrease the loss. The new weight value is computed using the method of stochastic gradient descent such that

[TABLE]

where $i$ is the iteration step, $\eta$ is the momentum, $\epsilon$ is the learning rate with a value larger than 0 and \nabla$$Loss^{i}_{W} is the gradient of the loss function at the current iteration with respect to the current weights. We employ the use of Nesterov momentum to modify the weight update by predicting the gradient at a new position and correcting the gradient at the current position (Nesterov, 1983). We use the common technique of stochastic gradient descent (SGD), (Bottou, 1991) whereby we determine the gradient of the loss function using subsets of the training data (here 128 samples). Training on batches of data is more efficient for weight optimization if the total number of training samples is much larger than the batch size. Additionally, when we train, we cycle through all of the training data 30 times but at each epoch the samples in the batches are randomized. We employ the use of dropout with 25$\%$ of our neurons on the first layer of each network. While training, dropout helps prevent the network from over-fitting by randomly zeroing a neuron’s output (Srivastava et al., 2014). The use of dropout is to first-order equivalent to an L2 regularizer because it adaptively distorts the neuron data to control over fitting (Wager et al., 2013).

The neural network relies on a handful of parameters that define the architecture (e.g. the hidden layer size and learning rate) which affect the performance. Tuning of these parameters was accomplished from a grid search where we trained over 1000 different neural works and chose the configuration that yielded the best performance. We use a hidden layer size of 64, 32, 8, 1, where each is the number of neurons in a respectively deeper layer of the network. The optimal parameters are a regularization weight of 0, learning rate of 0.1, momentum of 0.25 and a decay rate of 0.0001. The decay rate corresponds to the learning rate decreasing by $\epsilon^{i+1}=\epsilon^{i}/(1.+decay*i)$, where i is the iteration step. We tried different regularization terms between 10 – 10e-6 at factors of 10 and found they only decreased our network’s performance. A smaller learning rate usually gave us a smaller accuracy and higher loss over the same number of training epochs. Our optimal deep net uses 13,937 trainable weights across 105 neurons with 2492 neural connections or synapses. We refer to this algorithm as MLP in Table 2.

2.3 Wavelet

Wavelets allow us to represent a signal as a series of components where we can discard the "least significant" pieces of that representation and keep the original signal largely intact. Training a neural network on all the wavelet components will allow it to learn the most significant pieces while ignoring lesser ones that do not define the signal (e.g. the noise). We design an another fully connected (fc) deep net using the same structure as above except the input data undergoes a wavelet transform. We compute a discrete wavelet transform using the second order Daubechies wavelet (i.e. four moments) on each whitened light curve (Daubechies, 1992). The approximate and first detailed coefficients are appended together and then used as the input for our learning algorithm, Wavelet MLP.

2.4 Convolutional Neural Network

The photometric measurements from a light curve are correlated to one another through time. We can make use of convolutions (conv) to compute local features from time-ordered input data. Convolutional neural networks (CNN 1D) can be thought of as generating new input data that has been convolved with a specific filter. The weights within each filter are optimized in a similar manner as a fully connected layer. Using fully connected layers where the every input feature has a dedicated weight per neuron can quickly grow the number of trainable parameters a model has. CNNs utilize convolutions and down sampling to compute local properties of the data when they are correlated to one another. The input data are discretely convolved with a filter or kernel as such

[TABLE]

where the $t$ is the time index, $x$ is an input data array, $a$ loops over each element in the filter and $w$ are the weights in the filter. After the data has been convolved with a filter we down sample by averaging every three data points together to reduce the number of features for the next layer. We use an average pooling layer to help reduce the scatter from sources of noise. The average pooling layer mimics binning observational measurements in time. We also tested taking the maximum value within bins of 3 (max pooling layer) but it performed with less accuracy than an averaging layer. Planet finding techniques in the past have used convolutions via a matched filter approach however the filters are hand designed and only one is used. Our CNN 1D uses 4 filters each containing 6 weights that are optimized using the training data. After the input data are convolved and down sampled, we concatenate each light curve and use it as the input for a fully connected network with a layer size of 64,32,8,1. We use a ReLU activation function after convolving the data and before the average pooling layer. Additionally, 30 training epochs are used to teach the network. The training time per epoch for CNN 1D (our most expensive model) is 18 seconds on an Intel i7-7500U using TensorFlow (Abadi et al., 2015). The code for each network is provided online111https://github.com/pearsonkyle/Exoplanet-Artificial-Intelligence.

3 Cross-Validation and Algorithm Comparison

Each models uses 311040 training samples in batches of 128 over 30 epochs to learn the transit features (See Figure 7). The validation/test data set is not used to train any model. The validation set consists of 933120 samples that span a larger range of noise than the training data (see 1).

3.1 Multilayer Perceptron

The large amount of weights in any deep net can create a non-convex loss function where there exists more than one local minimum. Therefore different random weight initializations can lead to different validation accuracies. We find variations in the accuracies of our deep learning algorithms on the order of $\sim$0.1$\%$ suggesting our use of a SGD solver is robust.

3.2 Support Vector Machine

A SVM is a non-probabilistic linear classifier which constructs a hyperplane that optimally separates data points of two clusters given the labels (Vapnik & Lerner, 1963). A SVM is inherently designed to separate linear data but can be extended to separate non-linear data by transforming the input with a kernel (e.g. RBF). We use the method of SGD to optimize the weights from batches of 128 samples in our SVM. The method of SGD has been successfully applied to large-scale problems with more than 105 training samples (Bottou, 2010). A downside, similar to our MLP, is that the algorithm is sensitive to feature scaling so the data is whitened before training (mean=0 with unit variance). Based on a grid search of parameters we determined the optimal classifier has a cross-entropy loss function, linear kernel and an L2 regularization coefficient of 0.01. Regularization helps the weights from growing too large by introducing an additional term on the loss function that takes into account the size of each weight.

3.3 Least-Squares Box Fitting

We compare our machine learning techniques to a more traditional model for detecting single transit events by fitting a box function to data (see Figure 2). However, using a step function at the boundary of the transit does not allow for sub-pixel precision when finding the duration or mid transit because a step function is not continuous. This presents a problem because the optimization function (i.e. $\chi^{2}$) has a Jacobian with zero at multiple locations suggesting the current parameters require no updating and thus convergence has been reached. Instead, we use a 1 pixel slope to define the boundary of our transit such that the depth at the 1 pixel slope is half that of the full transit depth. Adding a slope allows us to evaluate our transit function with sub-pixel precision because it is a continuous function. We optimize our parameter estimation using a simplex algorithm, Nelder-Mead, and minimize the mean squared difference between the data and model (Nelder & Mead, 1965). Random parameter initializations often lead least-squares optimizations away from the global solution within noisy data. The local solution will depend on the initial conditions and optimization algorithm. We were generous by initializing the BLS routine with a mid transit value 20 minutes different (10 data points) than the actual mid transit value and all other parameters fixed to the true value. The BLS algorithm is not a probabilistic classifier like the MLP or CNN so we create a score function such that it can be compared in the receiver operating characteristic (ROC) plot. The score of the transit fit is calculated as a combination of accuracies for mid transit and transit depth. The accuracy score makes sure the mid transit is within our time series and close to the center of the data and that the transit depth is greater than 0 and in most cases, greater than the noise.

3.4 Algorithm Comparison

We use a receiver operating characteristic (ROC) plot to compare the results of each transit finding algorithm (See Figure 8). The ROC plot illustrates the performance of a classifier as its discrimination threshold is varied (Hanley & McNeil 1982). A classifier outputs a probability for a sample pertaining to a specific class and it is often just rounded up to signify a transit or non-transit detection (0 or 1). The ROC plot changes the classification of each model by varying where the probability data gets rounded. The true positive rate (TPR) is known as the probability of detection and taking 1-TPR yields the false negative rate. The ROC plot shows the cumulative distribution of the TPR as the discrimination threshold is varied against the cumulative distribution of false positive rates (FPR). The true negative rate can be calculated as 1-FPR. Ideally, the area under the curve should be close to unity which indicates a perfect classifier. The ROC plot shows the accuracy of each algorithm applied to test data the algorithms have not seen before.

4 Planet Detection Sensitivity

We explore the sensitivity of our algorithms in order to understand the detection limits and robustness of finding new planets. We use our test data set to explore the accuracy under varying amounts of noise (see sensitivity test in Table 1). Figure 9 shows the accuracy of our deep learning algorithms on data with varying amplitudes of noise. All of the initial parameters were kept the same as in Table 1 except for the noise parameter, this leads to 77760 samples in each noise bracket of Figure 9. By whitening the data we remove the scale of the transit depth from the observations such that the limiting detection factor is the noise parameter. As the transit depth approaches the noise the accuracy of detection goes down. Figure 9 provides us with a means to estimate how much binning is required to reduce the scatter to be able to detect small planets.

4.1 Time Series Variation

The nature of our deep net allows us to evaluate data from a different time series than the original training data using interpolation. The flexibility in the evaluation of different data stems from the networks ability to analyze 180 input features. Creating a new deep net for more instrument specific cadences is not necessary and we demonstrate this in Figure 10. Waldmann (2016) argues that one can interpolate a lower cadence signal onto a higher cadence grid. The interpolation would incur an extra noise penalty but still provide enough information for his deep net to make a significant classification. We tested that hypothesis by creating a new data set of varying resolutions compared to the original data and tested the performance of each algorithm. Since the new data have a different sampling rate for the same sized window in time, we linearly interpolate the data back onto the original grid size (i.e. 180 input features). We find that the convolutional neural network has the smallest performance drop when evaluating down sampled data. The accuracy of detection remains the same when evaluating data from a higher resolution grid. Arguably, data from unevenly sampled grids could also be interpolated onto a more uniform grid prior to input.

Our deep nets were originally trained on data within a 6 hour time window at a 2 minute cadence with the transit duration being between at least 30 minutes to 4 hours. The deep net only knows about 180 input features and nothing regarding the time domain, except that the input is time ordered (an important property for the CNN 1D algorithm). This allows us to stretch the boundaries of our data interpretation to different time domains. We can use the same deep net to detect planets within the Kepler data even though it has a time cadence of 30 minutes. The longer cadence limits the time domain to a 90 hour window with detectable transits being between 7.5 hours to 60 hours. These conditions might be beneficial for some planets but it will not find those with shorter transit durations. Given our findings above, data with a lower cadence could be super-sampled and then evaluated with a minimal decrease in performance ($\sim$2$\%$). However, for the purposes of this pilot test we use the native resolution of the Kepler data during our analysis.

4.2 Feature Loss

Observations are often subject to less than desirable conditions that suffer from instrumental malfunctions or weather induced systematics which can yield an incomplete set of measurements. We explore the capability of our deep learning algorithms to evaluate data that are missing features. Figure 11 shows the evaluation of our algorithms with missing data points. We randomly remove chunks of data from each light curve up to a certain extent. For each value, $i$, between 0 to 60$\%$ we select a random integer of chunks between 1 and $i$/10 to take out of the data. The positions of the chunks are randomized and chosen such that there are no overlaps which conserve the amount of data removed, $i$. The convolutional neural network has the best performance with less data because it pools local information together. Test data with a noise parameter above 1 was used in this analysis.

5 Time Series Evaluation

Real mission data is seldom as ideal as the training set. To test the ability of our deep nets we apply them to longer time series data from Kepler. In order to evaluate data from a time series thats longer than the training data (i.e 180 input features) requires separating data into bite sized chunks for the neural network. We do this by creating overlapping patches from the time series that are offset by a small value. Algorithm 1 in the appendix highlights the specifics of performing our evaluation on longer time series data. After the input light curve is broken up, we compute a detection probability for each chunk and superimpose the probabilities in time. The cumulative probability time series or just probability times series, $PTS$, is normalized between 0 and 1 for the whole time series (See bottom left subplot in Figure 11). An estimate on the periodicity of the planetary transit can be made by finding the average difference in maxima within the $PTS$. Smoothing the $PTS$ plots using a Gaussian filter helps to find the maxima.

The search for ever smaller planets depends on beating down the noise enough to detect a signal. When an individual planet signal cannot be found in the data we use a brute force evaluation method for estimating periodic transits (Algorithm 2 in appendix). For a range of periods we compute the phase folded light curve and bin it to a specific cadence. In our test case, we simulate Kepler data and therefore bin the phase folded data to a 30 minute cadence. Next we follow the same steps as the time series evaluation above; break the data up into overlapping light curves and evaluate the cumulative probability phase series, $PPS$. However, we do this four different times for the same period but change the reference position/time for computing the phase. The second for loop in Algorithm 2 shows the step where we compute four different phase positions for the period. The phase curve is shifted to account for the possibility of the transit being at the edges of the data. Statistical fluctuations which might cause false positive detections are also mitigated with this cycling approach. Each phase curve is ran through Algorithm 1 and produces a $PPS$. The variance is computed at each phase value across all four $PPS$ arrays producing a probability variance phase series, $PVPS$. Each phase value in the $PVPS$ has a variance associated with it. The mean value of the $PVPS$ is then plotted with the corresponding period (e.g. see Period Analysis in the subplot of Figure 11).

This method only works well for transits below the noise level. Larger transit signals will deconstructively add into the phase folded data but to an extent that can be greater than the scatter. Smaller transits do not have this problem and can stay hidden beneath the noise even when phase folded at the incorrect period. Only will the true signal reveal itself when the correct period is found and the signal is constructively added together to produce something above the noise. This method merely provides us with an orbital period estimate which we can then use as a prior for transit modeling routines.

5.1 Kepler Data

The Kepler exoplanet survey acquired over 3 years worth of data on transiting exoplanets that ranged in size from Earth to Jupiter-like. We only use one quarter of the Kepler data without any pre-processing to validate our neural network. Due to constraints from the time span the transit duration is limited to $\sim$7–15 hours with an orbital period greater than 90 hours. Detecting individual transit events with our algorithms requires a certain amount of features in-transit to yield an accurate prediction. We chose random targets with a noise parameter ($R_{p}^{2}$/R${}_{s}^{2}$)/$\sigma_{tol}$) larger than 1.5 and a transit period greater than 90 hours so that we can acquire multiple transit events in a single quarter of Kepler data. Figure 13 shows the results from a small analysis of Kepler targets. The probability plots are computed using the same method highlighted above for the time series evaluation. From the probability-time plot we can estimate the period of the planet by finding the average difference between peaks. The estimated orbital periods can deviate from the truth if discrepancies arise from multiple planet systems or strong systematics.

6 Discussion

Neural networks are an extremely versatile tool when it comes to pattern recognition. As demonstrated in Waldmann (2016) deep nets can be trained to characterize planetary emission spectra as a means of narrowing the initial parameter space for atmospheric retrieval codes. Automated detection and characterization will pave the way for future planet finding surveys by eliminating human interaction which can bottleneck analyses and introduce error.

Observations of exoplanets from different platforms contain separate observational limitations that can produce an incomplete set of measurements or sample the data differently. Accounting for each is best done using a convolutional neural network. Our CNN 1D algorithm achieved the highest performance on our interpolation and feature loss tests. The CCN 1D network computes local properties in the data using a filter convolution prior to input in a fully connected network. Additionally, the CNN 1D algorithms uses a down sampling technique by averaging which can help reduce some scatter in the data. Our findings indicate that data interpolated from a different resolution grid only suffered a 2$\%$ decrease in accuracy at most. The flexible nature of interpolation and network performance on different sized grids opens up the window of time for evaluating short and long period planets.

The link between statistical properties of exoplanet populations and planetary formation models must account for the observational detection biases of exoplanet discovery. The transiting exoplanet survey satellite, TESS, will be flying next year and estimates of the planet discovery rate for our deep net are estimated in Figure 14. We use the average photometric precision of TESS shown in Figure 8 of Ricker et al. (2014) and our CNN 1D detection accuracy for various transit depths and noise levels. Figure 14 shows the detection accuracy for potential planets assuming a solar type star with the same radius as the Sun. While TESS will target mostly M stars, which enhance the planet to star contrast, a solar type star places a lower limit on the transit depth. An Earth-sized exoplanet with a transit depth of 100 ppm around a 8th magnitude G star is 87$\%$ detectable with a noise parameter of 1.11.

The BLS method for detecting single transit events without a priori knowledge of the transit location and stellar variability is not adequate at finding transits. The BLS algorithm on its own has trouble finding the correct solution because the optimization function has multiple minima. Different parameter initializations can guide the optimization algorithm into a local minimum. Foreman-Mackey et al. 2015 find modeling the instrumental systematics along with the transit can greatly improve transit detection but at the cost of computational efficiency. Foreman-Mackey et al. 2015 generate a synthetic data set similar to our own, with transit depths ranging between 400–10,000 ppm and inject noise on the order of 25–500 ppm based on Kepler Magnitude to photometric precision relations. This corresponds to a noise parameter, $R_{p}^{2}$/R${}_{s}^{2}$)/$\sigma_{tol}$, between 0.8 and 16 for the smallest transit depth tested in their sample. In all of our test cases the worst accuracy achieved with a noise parameter of 0.8 was $\sim$60%. The worst detection accuracy in Foreman-Mackey et al. 2015 was on the order of 1% because shallower transits at longer periods were harder to detect for fainter stars. We predict our algorithm can detect small planets with large periods based on our findings from interpolating data from a super sampled grid. The average detection accuracy for a span of Kepler magnitudes in Foreman-Mackey et al. 2015 was on the order of $\sim$60% where we find our accuracies to be at least 60% or above for noise parameters greater than 0.8. Figures 10, 11 and 12 from Foreman-Mackey et al. 2015 were used to make this comparison.

Performing a Monte Carlo simulation like a Markov chain or nested sampler can help find global solutions because they sample large portions of the parameter space. However, these algorithms require a long computational time and may not be a suitable means of analysis for big data sets. Current detection methods rely heavily on the uncertainty estimations of transit fitting algorithms when determining a significant signal. The uncertainty estimation is non-trivial but significant transits can be verified with three individual measurements to help decrease false positive signals.

7 Conclusion

In the era of “big data” manual interpretation of potential exoplanet candidates is a labor intensive effort and difficult to do with small transit signals (e.g. Earth-sized planets). Exoplanet transits have different shapes, as a result of, e.g. the stellar activity. Thus a simple template does not suffice to capture the subtle details, especially if the signal is below the noise or strong systematics are present. We use an artificial neural network to learn the photometric features of a transiting exoplanet. Deep machine learning is capable of processing millions of light curves in a matter of seconds. The discriminative nature of neural networks can only make a qualitative assessment of the candidate signal by indicating the likelihood of finding a transit within a subset of the time series. For planet signals smaller than the noise we devise a method for finding periodic transits using a phase folding technique that yields a constraint when fitting for the orbital period. Neural networks are highly generalizable allowing data to be evaluated with different sampling rates by interpolating the data to a standard grid. We validate our deep nets on light curves from the Kepler mission and detect periodic transits similar to the true period without any model fitting. Additionally, we test various methods to improve our planet detection rates including 1D convolutional networks and feature transformations such as wavelets and find significant improvement with CNNs. Machine learning techniques provide an artificially intelligent platform that can learn subtle features from large data sets in a more efficient manner than a human.

The next generation of automation in data processing will be become more adaptive, self-learning and capable of optimizing itself with little to no auxiliary user input. The program will understand what it is looking at, make a qualitative pre-selection followed by a quantitative characterization of the exoplanet signal. In the future we would like to explore the use of more deep learning techniques (e.g. long short-term memory or PReLU) to increase the detection robustness to noise. Additionally, active research is currently being done in machine learning to optimize the network architecture and have it adapt to specific problems. Adding a pre-processing step has the potential to greatly improve the transit detection performance by removing systematics from the time series (e.g. Aigrain et al. 2017).

A special thanks goes out to the anonymous reviewer for providing helpful comments and references, this manuscript would not be same without you.

Appendix A Algorithms for Time Series Evaluation