A MACHINE LEARNS TO PREDICT THE STABILITY OF TIGHTLY PACKED PLANETARY SYSTEMS

In order to train our algorithms, we generated a data set of 5000 N-body integrations of three-planet systems over 10⁷ orbits of the innermost body. We focus on three-planet systems since there exists an analytic criterion for the case of two planets (Gladman 1993), and systems with more planets exhibit qualitatively similar behavior (Chambers et al. 1996). The number of simulations and length of integration were chosen to generate a data set at limited computational cost (∼1000 CPU hours) and assess the value of investing significant computing time to train classifiers on more astrophysically relevant timescales (∼10⁹ orbits). Because we expect that instability timescales of 10⁷ and 10⁹ orbits are both physically driven by Chirikov diffusion due to the overlap of mean-motion resonances (see Figure 2 in Obertas et al. 2016), we expect the performance of models trained on this data set to be comparable to that of similar algorithms trained on longer data sets.

With a view toward applying these models to Kepler discoveries, we adopted a solar-mass star, 5 ${M}_{\oplus }$ (Earth-mass) planets, and drew the innermost planet's semimajor axis randomly between 0.04 and 0.06 au. We note, however, that our results are strictly scale-free and can be applied to comparable systems with masses, orbital periods, and semimajor axes expressed in terms of the star's mass, innermost planet's orbital period, and semimajor axis, respectively. The second planet's semimajor axis was separated from the first by a number of mutual Hill radii drawn from a uniform distribution in the range [5, 9]. The third planet's separation was then independently drawn from the same distribution, yielding unevenly spaced planets. The particular range of Hill-radius separations was chosen to capture the regime of interest on our adopted timescale of 10⁷ orbits and roughly generate a balanced data set of stable and unstable systems (this yielded 1479 stable systems out of 5000). Eccentricities and inclinations were drawn independently for each planet from uniform distributions over [0, 0.02] and [0, 1°], and the remaining angles were drawn randomly over [0, 2π].

Because one only needs to integrate until the first Hill sphere crossing (once this happens, strong scatterings happen quickly; e.g., Gladman 1993), we use the efficient WHFAST integrator (Rein & Tamayo 2015) in the open-source REBOUND N-body package (Rein & Liu 2012). REBOUND is written in C99 and comes with an optional python interface. We adopted a conservative time step of 1% of the innermost planet's orbital period and classified systems as unstable if any pair of planets came within 1 Hill radius of each other during the simulation.

Binary classifiers are often evaluated on their precision (here the fraction of systems that are actually stable when the model predicts stability) and recall (the fraction of systems the model predicts are stable out of the truly stable cases). For typical algorithms that predict probabilities of class membership, one can trade off between precision and recall by varying the threshold probability for classification. For example, a conservative model that only classifies systems as stable if it predicts a probability of stability greater than 0.99 will be right most of the times that it predicts stability (high precision), but will miss all the stable systems that were assigned slightly lower probabilities (low recall). The appropriate threshold depends on the application (e.g., if predicting DNA matches for crime cases, one might set a high threshold as above to have confidence in predicted matches).

When comparing two models, one can plot pairs of precision and recall scores for all possible thresholds to generate a precision–recall curve. A common scalar metric for comparing classifiers is the area under this curve (AUC), which would be unity for a perfect model.

After experimenting with several machine-learning algorithms (random forest and support-vector machine implementations in the Python scikit-learn library), we found that gradient-boosted decision trees (GBDT⁶ ) XGBoost v0.6 (Chen & Guestrin 2016) yielded significantly higher AUC values for our data set.

A recurring theme in machine learning is that of "overfitting," an algorithm's tendency to latch onto irrelevant idiosyncrasies in the training set that cause it to predict poorly on unseen examples. Different algorithms therefore try to penalize overly complicated models in an effort to retain only the broad features that are likely to generalize well. In practice, the user navigates this balance between simplicity and complexity empirically, by tuning an algorithm's "hyperparameters" that mediate this tradeoff, training it, and checking performance (Section 2.2) on unseen data; this process, together with trying different features to maximize performance, is known as cross-validation. In order to rule out the possibility of (sometimes subtle) mistakes in cross-validation yielding overly optimistic performance metrics, it can be good practice to assign a subset of the data to a holdout (test) set that is never seen by the algorithm during cross-validation. Evaluation of the final model on the holdout set therefore provides robust metrics of the trained algorithm's expected performance on unseen examples; consistency between the cross-validation and test scores also suggests a robust cross-validation methodology. In our case, we randomly assigned 1500 systems to a holdout test set and used only the remaining 3500 for cross-validation.

A typical technique to reduce statistical fluctuations when comparing the performance of different sets of hyperparameters is k-fold cross-validation. One begins by splitting the training examples into k evenly sized groups; then, for each group, one trains the model on the remaining k − 1 chunks and uses the remaining group to evaluate performance. The scores from the k folds are then averaged, reducing the variance in the estimate. Finally, it is generally good practice to use stratified cross-validation, whereby one ensures that each of the k folds is assigned an approximately equal number of samples from each class (stable and unstable).

XGBoost has several hyperparameters, so we sequentially performed grid searches through two-dimensional cuts through the parameter space, evaluating performance through the precision–recall AUC (Section 2.2) using stratified, five-fold cross-validation on the training set of 3500 examples. Values of the final adopted hyperparameters for the algorithm are discussed in Sections 3.1 and 3.2 and can be found in Table 1.

We begin by considering as features the initial orbital elements and periods of each planet, as well as the interplanetary separations between adjacent planets in units of mutual Hill radii (23 features). We then trained an XGBoost classifier on these features (Section 2.3), allowing us to predict the stability of each system in the test set in the form of an estimated probability.

As discussed in Section 2.2, a threshold probability is required for classifying a system as stable/unstable and is a subjective choice that depends on the desired qualities of the classifier. For our purposes we argue it is logical to adopt a conservative threshold, in the sense that if the model predicts stability, there is a strong likelihood that the system is actually stable (high precision). This follows from the fact that it is computationally much faster to verify that a system is unstable (on short timescales) than it is to check that it is stable (on long timescales). We choose to require a precision of 90% on our test data set, which corresponded to the model only classifying systems as stable if their predicted stability probability is larger than a 0.785 threshold.⁷

Since previous works have identified the Hill separations between adjacent planets as important features (Chambers et al. 1996; Marzari 2014), we plot the performance of the model projected onto this 2D plane (Figure 1).

Zoom In Zoom Out Reset image size

Looking at the dashed, black line in Figure 1, one can see that to first order, the model's prediction boundary roughly obeys the relation Δ₁ + Δ₂ > 16.1. This is in fact the form of the simple criterion suggested by Lissauer et al. (2011), who quote Δ₁ + Δ₂ > 18 for stability on timescales of 10⁹ orbits. Because we consider stability on shorter timescales, the threshold number of Hill radii should be adjusted for a fair comparison. We term models of the form ${{\rm{\Delta }}}_{1}+{{\rm{\Delta }}}_{2}\gt x$ "Lissauer-family" models, and find 16.1 is the threshold for Lissauer-family models that yields a precision of 90% on this data set.

As stated above, a conservative probability threshold has the disadvantage that the model will misclassify many stable systems as unstable (low recall). This is easily seen by considering different Lissauer-family models (dashed black line in Figure 1), i.e., imposing different threshold values than 16.1 and generating lines parallel to the one plotted. Larger threshold values ensure that a larger fraction of the systems to the right of the boundary (i.e., those predicted stable) are in fact stable (blue), leading to higher precision. However, this lowers the recall, since now fewer of the stable systems (blue) are predicted to be stable by the model (i.e., lie to the right of the line). For a fixed precision of 90%, the machine-learning model has a significantly higher recall (52%) than the Lissauer-family model (30%). This is because the machine-learning model can use information in the features not visible in this 2D projection to make better predictions. The uncertainty in the recall should be comparable to the rms variations in the AUC calculated across the k folds during cross-validation, or ≈1% (Table 1).

An important factor determining the performance of a machine-learning algorithm is the quality of the features it is provided for each of the training examples. To this end, we improved upon the previous model (Section 3.1) by generating new features from short N-body integrations. To create the new features, we performed simulations over 10⁴ orbits (0.1% of the simulation timescale probed) for each of the 5000 systems in the data set and recorded each planet's orbital elements and Lyapunov time every 5 orbits.

Because we suspect that the instability is driven by overlapping mean-motion resonances, we first generated features that capture semimajor-axis variations, which would approximately vanish if the dynamics were purely secular (Murray & Dermott 1999). In particular, we generated features for the standard deviation and maximum value of each planet's semimajor axis over the 10⁴ orbits, normalized to the mean value over the same period (std_ai and max_ai, where i denotes the planet number). We also generated features for these quantities over only the first 50 orbital periods to capture variations on orbital timescales (std_window_ai and max_window_ai). Furthermore, we capture any drifts by generating slope features from linear fits to each of the three planets' semimajor axes, normalized to the mean semimajor axis divided by the integration time (slope_ai). For the eccentricities, we took the mean, standard deviation, maximum, and minimum values over the full 10⁴ orbits, normalized to the eccentricity the planet would require to reach its nearest neighbor's semimajor axis (avg_ei, std_ei, max_ei, min_ei). For the Lyapunov time we generated a single feature corresponding to the value measured at the end of the integration, normalized to the innermost orbital period (LyapTime). Finally, we added features for the two pairs of initial Hill-radius separations (daOverRH1, daOverRH2), and for the minimum and maximum initial Hill-radius separations (mindaOverRH, maxdaOverRH). We experimented with features describing orbital inclinations, but they did not significantly improve the models.

A summary of the features can be found in Table 2, ordered by their importance. We quantify the importances through the Gain value recorded by XGBoost, which corresponds to the gain in accuracy that a given feature provides when it is introduced into the underlying decision trees used by the algorithm. The units are normalized so that the gains sum to 100. We find that the variations in the middle planet's semimajor axis are the most informative in this sense. This suggests that the instabilities in these closely packed systems are driven by the overlap of mean-motion resonances (which change the semimajor axes), rather than secular effects (which would keep the semimajor axes constant).

Finally, we compare the performance of this "short-integration" model to the previous "initial-condition" model. Figure 2 shows that while both models often assign unstable systems (blue bins) a low predicted probability, the "initial-condition" model assigns stable systems (green bins) a wide range of predicted probabilities, translating to a lower recall. In contrast, the "short-integration" model more confidently assigns high predicted probabilities to stable systems, better separating the two classes. Again, setting the predicted probability threshold so as to obtain 90% precision, the recall improves to 68%.

Zoom In Zoom Out Reset image size

In summary, generating improved features from short integrations and adopting XGBoost as our algorithm provided our largest AUC gains (Table 1).