Classifying Kepler light curves for 12,000 A and F stars using supervised feature-based machine learning

Selecting an adequate training set to train a feature-based classifier for all 200,000 stars in the Kepler field is the most challenging and time-consuming aspect of developing a general classification algorithm. Many classes of variable stars are rare, while others remain poorly understood, and still others have not yet been identified. Indeed, new classes are occasionally proposed (Debosscher et al. 2011; Bass & Borne 2016; Pietrukowicz et al. 2017). For this reason, a good training set must be representative of a wide range of variable stars to construct a suitably general feature-space representation of Kepler light curves. Rather than attempting to compile a collection of all known classes of variable and non-variable stars, as attempted with limited success by Debosscher et al. (2011), we focused our research on a subset of seven well-studied classes, within the temperature range $6500\,\text{K}\leq T_{\mathrm{eff}}\leq 10,000\,\text{K}$, as an initial demonstration. The seven chosen classes are $\delta$ Scuti stars, $\gamma$ Doradus stars, RR Lyrae stars, rotational variables, contact eclipsing binaries, detached eclipsing binaries, and non-variable stars. Typical light curves and power spectra for each class are included in Fig. 1. These classes are commonly represented in the wider Kepler data (with the exception of the RR Lyrae class), are interesting stars that we wish to study in further detail, and are not intermediate or hybrid classes. We excluded hybrid stars to avoid having light curves in the training data that belong to multiple classes.

Two of our seven classes are subclasses of eclipsing binary (EB) systems, which were catalogued in Kepler data by Kirk et al. (2016). EBs can be classed as one of: detached binaries, where the two stars are far from each other to give highly separated eclipses; contact binaries, where there is no space between the two stellar envelopes, producing almost sinusoidal eclipse patterns; and semi-detached binaries, the intermediate class of the two extremes. In accordance with our decision to exclude hybrid classes, only the detached and contact subclasses have been included in our training set.

Rotational variables are most commonly found among cooler stars, whose star spots present darker patches on the surface that rotate in and out of view, but rotational variability is also seen by Kepler across the effective temperature range of our sample (Balona, 2013; Sikora et al., 2020). Typically, the spots have lifetimes not much longer than the rotation period, and they may occur at different latitudes, so the variability is only quasi-periodic (Nielsen et al., 2013; McQuillan et al., 2014). The $\alpha^{2}$ CVn stars are hotter stars whose strong dipolar magnetic fields concentrate certain elements into spots. These also rotate with the star, but occur near the magnetic poles and are much longer lived, leading to light curves that do not change rapidly in period, amplitude, or shape (Wolff, 1983). In our chosen temperature range, examples of both are found.

Three classes of pulsating variable star were included (for a recent review of pulsating stars, see Kurtz 2022). RR Lyr variables are bright, evolved stars burning helium in their cores. As they traverse the horizontal branch, they cross the instability strip and pulsate periodically with a characteristic phase curve. Their use as standard candles has allowed measurements of the distance to the Galactic centre and to globular clusters (Oort & Plaut, 1975; Walker, 1992). The two other pulsating star classes, $\gamma$ Doradus and $\delta$ Scuti variables, both comprise A or F-type stars on or near the main sequence, and embody two distinct types of oscillation: g modes, or buoyancy-driven modes sensitive to the near-core region of a star; and p modes, pressure-driven modes most sensitive to the envelope. $\gamma$ Doradus stars are multiperiodic g-mode pulsators with periods between approximately 0.3 d and 3 d (Kaye et al., 1999). Despite having periods similar to the RR Lyr variables, the multiperiodic $\gamma$ Doradus stars do not have simple phased curves. There are several hundred in the Kepler field (Li et al., 2020), and they have seen substantial recent attention because of their ability to probe internal rotation (Van Reeth et al., 2018; Ouazzani et al., 2019), diffusive mixing (Bouabid et al., 2013), and core overshooting (Mombarg et al., 2019).

Finally, $\delta$ Scuti stars are the most common class of pulsating star at A and F spectral types, with approximately 2000 known in Kepler data alone (Murphy et al., 2019; Guzik, 2021). These stars are p-mode oscillators, and with periods between 18 min and 8 hr they are the highest-frequency variables in our sample. For unknown reasons, even in the middle of the $\delta$ Scuti instability strip only half of the stars pulsate as $\delta$ Scuti stars (Murphy et al., 2019), hence we include a non-variable class in this work. We note that some $\delta$ Scuti stars are known to lie outside of the instability strip (e.g., Bowman & Kurtz, 2018), but our classifications are based only on the Kepler light curves and not on parameters such as effective temperature.

We created a training set across all seven classes by hand-picking 1319 stars from candidate lists according to specific criteria. Stars were restricted to the temperature range $6500\,\text{K}\leq T_{\mathrm{eff}}\leq 10,000\,\text{K}$ using effective temperatures from Mathur et al. (2017). We examined one quarter of long-cadence Kepler photometry for each star to prepare the training data. Quarter 9 (Q9) was chosen because it has no prolonged gaps in observation, such as those arising from telescope safe mode events, and no anomalies in data quality. We used light curves made with simple aperture photometry (SAP), downloaded from the Kepler Asteroseismic Science Operations Center (KASOC) website (Data Release 25).¹¹1http://kasoc.phys.au.dk/ The choice to examine a single quarter was made to reduce computation time, but this also precludes the analysis of variability on timescales longer than a typical 90-d quarter. While we certainly recommend the investigation of four-year data in future research focused on a wider range of Kepler variables, this will not have a great effect on the stars chosen for our investigation. Of the seven classes, only a handful of detached binaries are known to have a period greater than 90 days (Kirk et al., 2016), and we did not include these in the training set.

The training stars were chosen from lists of possible candidates for each of the seven classes by visually inspecting their light curves and Fourier transforms. This laborious process embodies the motivation for automated variable star classification, and was a necessary task to ensure that the training data were accurate and would not mislead the automated feature-selection process. In the following paragraphs, we describe the selection of stars in training sets for each class. Table 1 summarises the class-specific number of stars in the training set. The full list of stars is provided as supplementary material, with a sample shown in Table 2.

Eclipsing binary systems were selected from the Kepler Eclipsing Binary Catalog (Kirk et al., 2016), restricted to periods $<$90 d and a morphology index of $0\leq c\leq 0.5$ (detached binaries) or $0.75\leq c\leq 1.0$ (contact binaries), as recommended by Matijevič et al. (2012).

Our selection of $\delta$ Scuti stars began with 2405 stars manually identified as variable at frequencies above 7 d^-1 from a preliminary version of the Murphy et al. (2019) catalogue. We randomly selected 1000 of these, and further refined this list to remove any stars that were also $\gamma$ Doradus stars (i.e. $\gamma$ Doradus/$\delta$ Scuti hybrids) by manual inspection. From the same source, we also chose 500 stars that were not $\delta$ Scuti pulsators, and removed stars with low-frequency variability to arrive at the 201-star non-variable class.

We selected the $\gamma$ Doradus sample from the Debosscher et al. (2011) catalogue by choosing stars with a label confidence of $>$95% and an effective temperature in the appropriate range. While the Debosscher et al. catalogue is known to have errors, this approach was taken due to a lack of an available list of $\gamma$ Doradus stars exhibiting a broad range of oscillatory behaviours characteristic of the class — that is, a sample not restricted to neat and well-studied $\gamma$ Doradus stars from which scientific inference has been made (references in Sec. 2.1). The addition of rigorous manual inspection ensured that the $\gamma$ Doradus stars included in the final sample were significantly more likely to be correctly classified than in the Debosscher et al. catalogue, and that hybrid pulsators were removed.

Unlike the other classes, RR Lyr variables are not common in the Kepler data set. Of the 47 Kepler RR Lyr stars we found in the literature (Molnár et al., 2018; Nemec et al., 2013; Murphy et al., 2018), only 25 were observed in Q9. We admitted all 25 of these in the hope that we might discover additional RR Lyrae variables when classifying the remainder of the Kepler field (we did not).

The rotational variables were selected after trialling a preliminary version of our classifier, trained on the other six classes, on a test sample of Kepler stars. When visually inspecting the classification results across the six classes, we found that rotational variables constituted a considerable fraction of stars (approximately $15-20\%$). From these, we added a list of 166 rotational variables to the training set after a second manual verification.

Starting with SAP fluxes from Q9 light curves, we processed the data to remove instrumental variability by eliminating long-period trends in the light curve of each star. Such variability can arise from physical drift of the telescope, causing changes in the flux levels falling in the aperture mask, as well as other instrumental effects distinct from stellar variability. Our processing involved division by a smoothed version of each light curve (smoothed using a Savitzky-Golay filter), removal of single-point outliers more than $3\sigma$ from the mean of the smoothed light curve, and converting units to magnitudes.

Any gaps in the data of more than an hour (corresponding to two 29.45-min integrations) were padded with either the mean of the time series for long gaps of four or more integrations, or the mean of the points on either side of smaller gaps. Even in high-quality quarters, long gaps arise from standard telescope operations such as the data downlinks that happen for approximately 24 hrs twice every quarter, while small gaps may be caused by e.g. cosmic ray events. Most machine-learning tools operate as functions of array index rather than explicit functions of time, hence it is imperative that these gaps are filled.

The task of selecting relevant properties of a time series for a given application, like light-curve classification, is commonly a manual one performed by a given researcher (e.g., Pashchenko et al. 2017). An alternative approach, termed ‘highly comparative time-series analysis’ (Fulcher et al., 2013; Fulcher & Jones, 2014), is to include a large and comprehensive candidate set of possible time-series features, and take a data-driven approach to selecting those that are most relevant to the task at hand. To extract features from a light curve, we used a comprehensive candidate set of over 7000 time-series features from the hctsa software package (v0.96) (Fulcher & Jones, 2017).²²2https://github.com/benfulcher/hctsa The hctsa feature set encompasses a wide range of time-series analysis methods, from properties of: the distribution of time-series values, linear and nonlinear autocorrelation, entropy and complexity measures, stationarity, time-series model fits, wavelet and Fourier basis-function decompositions, and others (Fulcher et al., 2013). This approach allowed us to represent a set of $L$ light curves as an $L\times F$ matrix, where $F$ is the number of features; applying hctsa to our training data set yielded a $1319\times 7873$ light curve $\times$ feature matrix, where each row is labelled according to one of the seven classes listed in Table 1. After performing feature extraction, we excluded features that contained special values (NaN, Inf), returned an error, or produced near-constant outputs (within $10\times$ machine precision) across all 1319 time series, resulting in 6492 features after filtering. As a preprocessing prior to classification, feature values were normalized to the unit interval using a scaled, outlier-robust sigmoidal transformation (Fulcher et al., 2013).

In modern applications of machine learning, complexity is often introduced at the level of the classifier. In this work we instead focused on selecting from a large candidate set of complex features, but using simple classifiers. This has the advantage of yielding features that can provide clear scientific interpretation, and follows the approach of Timmer et al. (1993): “The crucial problem is not the classificator function (linear or nonlinear), but the selection of well-discriminating features. In addition, the features should contribute to an understanding”. For classification, we used a Gaussian Mixture Model (GMM) (McLachlan & Peel, 2000) on a labelled time series $\times$ feature matrix (described above). We fitted a single Gaussian component to each of the seven training classes in feature space, combining them with equal prior probabilities to form a seven-component Probability Density Function (PDF). While all classes are not equally common, equal priors are the simplest choice without knowing the true distribution of variable stars in the Kepler field. Classification was performed by evaluating the (posterior) probability of a star belonging to each class using the trained PDF, and selecting the class with highest probability. This GMM approach was substantially faster (by factors of approximately 10–100) than alternative algorithms such as nearest-neighbour clustering or support vector machines, but achieved similar classification performance on our training set.

We evaluated classification performance as the average balanced accuracy computed using 10-fold stratified cross-validation Hastie et al. (2009). Balanced accuracy, $C_{\mathrm{bal}}$, accounts for class imbalance (the unequal number of observations in each class) in our data set and is defined as:

where $m$ is the number of classes, $t_{i}$ is the number of successfully identified time series in the $i$th class, and $c_{i}$ is the total number of time series in this class.

To extract a small number of hctsa features that are most informative of the class labels, we used greedy forward feature selection (Hastie et al., 2009; Fulcher & Jones, 2014). This simple algorithm iteratively builds a feature set, one feature at a time, with the objective of maximising the balanced classification accuracy, $C_{\mathrm{bal}}$, at each iteration. That is, at iteration $k$, the algorithm searches across all individual features for the feature that maximizes $C_{\mathrm{bal}}$ when used in combination with the features selected in the $k-1$ previous iterations.

hctsa was developed to encompass a comprehensive sample of the interdisciplinary time-series analysis literature, and thus contains groups of features with highly correlated behavior (Fulcher et al., 2013; Henderson & Fulcher, 2021). When multiple features exhibit similar classification performance, we implemented a simple heuristic constraint to favour the selection of faster-to-compute features: at each iteration, of the features with an accuracy within a margin of 1% of the best-performing feature, the feature with the fastest computation time was selected. The iterative procedure was terminated when the improvement in training-set $C_{\mathrm{bal}}$ from adding another feature dropped below 1%. Note that applying this algorithm to the full training data set has the potential to overfit, since the selection step at each iteration (despite using cross-validation for each feature) uses the training set itself to select the best-performing feature. Accordingly, we evaluate the performance of our reduced feature set on an independent test set in Section 4.

We first investigated the structure of the seven labelled classes of 1319 Kepler stars in the 6492-dimensional hctsa feature space. We found that the hctsa feature space is able to capture characteristic properties of the seven labelled classes of stars, obtaining a high mean 10-fold cross-validated balanced accuracy of 95.9% (using a linear support vector machine (Hastie et al., 2009), compared with a chance rate for seven classes of 14.3%). This indicates that each type of star displays distinctive dynamics in ways that can be detected by the features included in hctsa. To better understand the structure of light curves in the high-dimensional hctsa feature space, we inspected a two-dimensional $t$-SNE visualization ($t$-distributed stochastic neighbour embedding; Van Der Maaten & Hinton 2008). The result is shown in Fig. 2, where each point is a light curve, and light curves with similar features tend to be positioned closely in the space. While $t$-SNE is an unsupervised technique (Fig. 2 was constructed without knowledge of the class labels), stars are meaningfully organized according to their labelled class, with most stars clustering with other stars of the same type. Consistent with the high classification results reported above, this indicates that the hctsa feature space captures distinctive dynamical properties of the light curves corresponding to the seven different types of stars. The plot also reveals scientifically meaningful structure between classes, such as the continuum from non-variable (light orange) stars to rotational-variable (light blue) stars. We also see a small overlap between RR Lyr stars and contact EBs, which reflects the similar morphologies of their light curves (see Fig. 1).

The results above demonstrate that time-series properties in hctsa can capture differences in light-curve dynamics between different types of stars. But which types of individual time-series features are most informative of these differences? To address this question, we aimed to construct a reduced set of hctsa features that display strong classification performance using greedy forward selection (see Sec. 3.3 for details). The cross-validated balanced misclassification rate on the training set is shown as a function of the number of selected features in Fig. 3. This plot reveals that strong in-sample classification performance can be obtained with a relatively small set of well-chosen time-series features, e.g., a balanced accuracy of 95.2% with just three features. According to our termination criterion—when an additional feature provides $<1$% marginal improvement in balanced accuracy—we obtained an informative five-dimensional feature space in which to represent Kepler light curves.

To visualize how stars are organized in the reduced feature space, we plotted the training set in the space corresponding to three of the selected features in Fig. 4 (left). Despite a dramatic dimensionality reduction of each time series—from the 4767 data points in a typical Q9 time series to just three extracted summary statistics—the space meaningfully organizes all seven training classes in this low-dimensional feature representation, with each occupying a characteristic region of the space. Much like the $t$-SNE construction in Fig. 2, the relative positions of each class are consistent with what we would intuitively expect from their light curves and power spectra in Fig. 1. For example: detached binaries are highly separated from the other classes, as their light curves are the most distinct; non-variable stars blend with rotational variables when the rotations are weak and difficult to distinguish by eye, such that the light curves are almost non-varying; the $\gamma$ Doradus and $\delta$ Scuti stars lie on opposite sides of the space, reflecting their contrasting low- and high-frequency pulsations; and the contact binaries are close to the $\gamma$ Doradus and RR Lyrae clusters, which all characteristically exhibit regular low-frequency variability.

In this section, we use the low-dimensional feature space learned from the training set and validated on the test set to classify variable stars across the entire Kepler field. Our full classification catalogue is provided as supplementary material (Table 4).

When our paper was in the final stages of preparation, a new classification of Kepler light curves was published by Audenaert et al. (2021, hereafter Aud21). Their work was done as part of efforts to design an automated classification algorithm for the TESS mission. Given the complementary nature of Aud21 and our own study, especially given that both were based on Q9 data, it is worthwhile to carry out a brief comparison. We should keep in mind that the emphasis in Aud21 was on providing a high-performance classification pipeline from existing methods, whereas ours involved designing an interpretable classifier from a rich library of time-series features.

The classification by Aud21 included about 167,000 Kepler Q9 light curves, regardless of effective temperature, whereas our work is restricted to about 12,000 stars with $6500\,\text{K}\leq T_{\mathrm{eff}}\leq 10,000\,\text{K}$. We have compared our classifications with Aud21 in Fig. 10. There is an obvious mapping between most of our classes and those used by Aud21, with the following differences:

• •

  Aud21 combined contact eclipsing binaries and rotational (spotted) variables into a single class.

• •

  Aud21 included $\delta$ Scuti stars in a class with $\beta$ Cephei stars. These have similar light curves but the $\beta$ Cep pulsators have higher effective temperatures that lie outside the range of our sample. Similarly, Aud21 combined $\gamma$ Doradus stars with SPBs (Slowly Pulsating B stars), which are also hotter than our sample.

• •

  Aud21 included a class for solar-like oscillators, which should not appear in our sample because they occur in stars whose effective temperatures fall below our range.

• •

  Aud21 also included a class for aperiodic variables.

We see from the confusion matrix in Fig. 10 that there is generally excellent agreement between our results and those of Audenaert et al. (2021). We briefly discuss the areas with the greatest disagreement:

1. 1.

   3094 stars that our classifier labelled as non-variable were classified by Aud21 as rotational/contact EB. We inspected 200 of these light curves (and their Fourier amplitude spectra) and found that most are non-variable, with some showing a weak rotation signal.

2. 2.

   637 stars that we labelled as contact EBs or rotational variables were classified by Aud21 as $\gamma$ Doradus pulsators. Inspection of 200 light curves shows that most are indeed $\gamma$ Doradus stars. This may be a shortcoming of our specific feature space and classifier, particularly when considering Fig. 7, where the same disagreement occurs between our GMM classifications and our independent test set.

3. 3.

   153 stars were labelled by us as $\gamma$ Doradus stars and by Aud21 as contact EBs or rotational variables. Inspection of these shows that many are indeed $\gamma$ Doradus stars, although it is sometimes difficult to be sure.

4. 4.

   139 stars in our sample were labelled by Aud21 as having solar-like oscillations, which is not a class that we considered because these oscillations occur in stars below our temperature range. Our classifications for these light curves were mainly as rotational variables, contact EBs or non-variable. We inspected all 139 light curves and found that our classifications were mostly correct.

5. 5.

   106 stars were labelled by us as contact EBs or rotational variables, and by Aud21 as $\delta$ Scuti stars. We inspected all light curves and found that most have $\delta$ Scuti pulsations, but many also have low-frequency variability.

Finally, we note KIC 10024862, which is one of two stars listed by Aud21 as non-variable and by our algorithm as a detached binary. In fact, Kawahara & Masuda (2019) identified this as a Jupiter-sized exoplanet in a long-period orbit that has only one transit during the 4-year Kepler mission, which happened to be in Q9. This suggests that it might be worthwhile to look in more detail at groups for which classification methods are in disagreement for a small number of stars.

Much like our approach, Aud21 assigned labels to each star according to the class with the highest posterior probability from their classifier. Figure 10 therefore contains samples where either classifier may be confused — for example, a given light curve may have probabilities of 0.34, 0.35, 0.01 split between three classes and the maximum probability (0.35) is relatively low. Not surprisingly, we found that by restricting to stars with a high maximum probability in both samples, the agreement increased between our classification labels and those from Aud21. A detailed comparison of the two catalogues goes beyond the scope of this paper and would require a measure of label confidence from the Aud21 classifier similar to the probability density heuristic from Section 4.3.2.

In general, we conclude that the two approaches produce results that generally agree well. The difference in point (i) reflects the subjectivity in drawing the line between variables and non-variables (and perhaps also different amounts of filtering applied to the light curves). Points (ii) and (iii) reflect the difficulty—especially with short data sets—in deciding whether low-frequency variability is due to pulsation or rotation (e.g., Briquet et al., 2007; Lee, 2021; Kurtz, 2022).