A Comparative Study of Machine-learning Methods for X-Ray Binary Classification

All the data used in this paper were obtained with the Gas Slit Camera (GSC; Mihara et al. 2002; Tomida et al. 2011) on board MAXI (Matsuoka et al. 2009). In operation since 2009, MAXI is the first astronomical mission to be operated on the International Space Station (ISS; Matsuoka et al. 2009; Sugizaki et al. 2011). MAXI has higher sensitivity and higher energy resolution than any other all-sky X-ray-monitors flown to date (Mihara et al. 2002; Tomida et al. 2011). GSC has a 1 day sensitivity of 9 mCrab (3σ) compared to 15 mCrab with RXTE/ASM (Levine et al. 1996) and 16 mCrab with Swift/BAT (Krimm et al. 2013). The GSC on MAXI covers the energy band from 2 to 30 keV and contains six units of Xe-gas proportional counters, which are assembled to cover wide fields of view of 15 × 160° (Sugizaki et al. 2011). In addition, MAXI provides on-demand processing to extract lightcurves in user-specified energy bands. ⁸ We specify the energy bands 2–3 (Low), 3–5 (Medium), and 5–12 keV (High), which are close to those used by VB13 and GP15 and demonstrated by Islam et al. (2021) to be effective in separating classes of XRBs.

In this paper, the soft colors (SC) are defined as

where Low and Medium refer to the counts in the 2–3 keV and 3–5 keV energy bands, respectively. The hard colors (HC) are defined as

where High refers to the counts in the 5–12 keV energy band. Park et al. (2006) found that using this fractional difference ratio works well both in high-count and low-count regimes, whereas simple ratios or logarithms of ratios tend to fail in low-count regimes. After we compute SC and HC, we rescale them to match the scale of the relative intensity, which ranges from 0 to 1.

The relative intensity is computed by summing the counts in each of the energy bands (2–3 keV, 3–5 keV, and 5–12 keV) and then normalizing each source by the average of the top 0.01% of the counts to ensure the relative intensity is scaled from 0 to 1. The relative intensity (RelInt) can thus be written as

We only used detections with at least 3σ significance for each XRB observed by MAXI. For a data point to be statistically valid, we require that it has at least 3σ significance in the sum of the counts in the three energy bands. We achieve this by requiring that the counts in the individual energy bands be detected at the accuracy of such that the total is summed in quadrature:

From each source, we also remove data points that deviate 10σ or greater from the mean as these outliers may be unphysical in origin (e.g., these outliers can appear when the reflection from solar panels comes into view). We demonstrate the effect of this 10σ cutoff in Figure 1. For some sources in the MAXI data set, we found that very few to no data points remain once checked for 3σ significance. This can happen because a source may vary in brightness and become fainter over time. For proper statistics, we considered only sources with at least 100 data points with 3σ significance (Figure 2).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Of the 58 sources from MAXI that we checked, only 12 BHs, 20 NPNSs, and 12 pulsars matched our criteria for statistical significance. These 44 sources are listed in Table 1. The CCI diagram of all 44 sources is plotted in Figure 3, and their corresponding 2D projections are plotted in Figure 4. Because most of the BHs and some NPNS are transients, the sources that went into outburst during 16 yr of RXTE/ASM may or may not go into outburst during the 10 yr of MAXI operation. However, we have used new transients (which had not been active during the 16 yr of RXTE/ASM) that have been discovered by MAXI.

Zoom In Zoom Out Reset image size

In our selection of data points, observations that are less bright are inherently less likely to be included. Because XRB systems of the same type can exhibit considerable variability in terms of which regions they occupy in CCI diagrams (Figure 3), this bias in data selection can be a limitation. We address this by subsampling where the probability that a particular observation is included in the training set is inversely proportional to the total number of observations of its system in the entire training set. An additional motivation for subsampling the data is for a computational reason: Gaussian process models can be particularly computationally expensive due to the inversion of a large matrix. The time required to perform this matrix inversion scales with O(N³), where N is the number of data points in the data set.

We run the subsampling algorithm 10 times to create 10 independent subsets of the data on which we can run our algorithms. The histograms in Figure 5 show the number of observations between various systems before subsampling and five examples of the distribution after subsampling. For each independent run, we sample 20% of the training data without replacement, resulting in a data set of 10,314 observations.

Zoom In Zoom Out Reset image size

BGP methods have been extensively used in ML for the classification and quantification of prediction uncertainty (Rasmussen & Williams 2005). In our work, we are using the BGP described in Williams & Barber (1998) as implemented in the R library Kernlab (Karatzoglou et al. 2004).

To use a BGP to solve our classification problem, we must choose a covariance kernel. In ML, capturing the similarity between data points is essential; on a basic level, we assume that points with inputs x that are close are likely to have similar target y values. The covariance function determines how this proximity or similarity is defined, and in this paper, we consider squared-exponential (i.e., Gaussian), Laplacian, and ANOVA radial basis function kernels as implemented in Karatzoglou et al. (2004). All three kernels are variants of radial basis functions, meaning that they are both stationary and isotropic kernels. Stationary kernels are invariant to translations in the input space and are a function of . If a kernel is also isotropic, it is more specifically a function of , making it invariant to rigid motions (Rasmussen & Williams 2005). We explicitly define the Gaussian and Laplacian Kernels below:

where refers to the Euclidean distance (i.e., the two-norm) between two inputs x and , σ is the characteristic length scale. To see a full derivation of the ANOVA kernel, please refer to Hofmann et al. (2008). In Section 4.2.2, we discuss how we compare and select from these three kernels for our specific classification problem.

We determine the median probability of belonging to a certain class and estimate the errors in our BGP model by undertaking five independent runs of the best BGP model (Table A1) as described further in Section 4.3. Although we did this for 10 runs for each of our other models, computational costs limited us to 5 runs for the BGP. Table A1 lists the median and standard deviation (SD) from the five runs.

A simpler ML algorithm called KNN (Altman 1992) is also used to predict XRB classifications. To implement our KNN algorithm, we customized code from the R caret package (Kuhn 2008) and the R class package (Venables & Ripley 2002). For any given input CCI data point, the KNN model first calculates the Euclidean distance between that data point and every other point in the training set. Certain variables might have an unusually large impact on this distance because their magnitude is larger without normalization. For example, the relative intensity spans a larger range of values than SCs and HCs, meaning that any distances along the SCs and HCs could be overshadowed by distances along the relative intensity. To adjust for this, we normalized the data set by the average of the top 0.01% of the counts before we began. Normalizing the training data has been shown to improve the KNN algorithm’s predictive accuracy dramatically (Hastie et al. 2001; Piryonesi & El-Diraby 2020).

After calculating the distances, we find the k data points closest to the input point p, which are the k-nearest neighbors. Then, we classify the point p based on the classifications of those neighbors. For example, if the majority of the k closest data points are labeled BH, then we classify the point as a BH. In Figure 6, KNN for a simple 2D case is illustrated where the two classes are BHs and NPNSs. Here, k = 5 and for the unknown point p, 60% of the k-nearest neighbors are NPNSs while 40% are BHs. Thus, the source has a 60% probability of being an NPNS and a 40% probability of being a BH.

Zoom In Zoom Out Reset image size

We choose k with a cross-validation procedure described in Section 4.2. After choosing k, we determine the median probability of belonging to a certain class and estimate the errors in our KNN model by undertaking 10 independent runs. Independence is achieved by generating a different subsample for each run (see Figure 5). Table A2 lists the median and SD from the 10 runs.

Finally, we compare the use of the SVM (Cortes & Vapnik 1995) algorithm to make predictions on the compact object type of an XRB. To implement our SVM algorithm, we customized code from the LIBSVM library (Chang & Lin 2011). In an SVM, given the input and output data for each source, the goal is to find the hyperplane that creates the largest separation between the three classes of data. This is most clearly explained by examining a simple two-dimensional example.

Consider Figure 7(a), where data points are graphed with respect to variables x and y and colored purple or orange depending on whether they belong to the BH or NPNS class. Given these inputs, a linear SVM model would calculate the two-dimensional line that provides the largest possible separation between the two classes. In this case, the points border the two dotted lines and are a distance D apart. The separator between the two classes is the solid line between them. New points are classified based on their location with respect to this line. Note that to create this model, only the points near the boundary determine the location of the separator. Thus, points far away from the boundary could be removed from the training set and the resulting model would be identical.

Zoom In Zoom Out Reset image size

However, this method has the obvious problem that the points may not always be linearly separable. To resolve this, a nonlinear SVM model uses a method known as the “kernel trick” where the points are projected into a space with even more dimensions, allowing the separator line to become a nonlinear curve. An example of this can be seen in Figure 7(b), where a nonlinear curve separates the BHs from NPNSs.

In higher dimensions, the logic is similar: the points are graphed in an n-dimensional space, and the optimal (n – 1)-dimensional dividing surface is then calculated. In the case of a linear SVM, this is a linear hyperplane; in the case of a nonlinear SVM, which we will use here, this is a nonlinear surface.

We computed both the probability of belonging to a class and error estimates. SVMs can produce class probabilities as outputs by fitting a sigmoid function

to the decision values f of the binary SVM classifiers, where A and B are estimated by minimizing the negative log-likelihood function. For this multiclass case, all three binary classifier probability outputs can be combined as described in Wu et al. (2009). This method is equivalent to fitting a logistic regression model to the estimated decision values (Karatzoglou et al.2006).

We determine the median probability of belonging to a certain class and estimate the errors in our SVM model by undertaking 10 independent runs. Independence is achieved by generating a different subsample for each run (see Figure 5). Table A3 lists the median and SD from the 10 runs.

For each example in our training data set, the input vector consists of the coordinates in the CCI diagram and the output vector is the class (1 = BH, 2 = NPNS, or 3 = pulsar). The model is trained with these examples and asked to predict the probability that the input vector belongs to each of the classes. The output vector thus takes the form of Y_pred = (0.12, 0.73, 0.15) for an example that has a 12%, 73%, and 15% probability of belonging to the BH, NPNS, and pulsar classes, respectively.

For the MAXI data set, these examples can be described as two-tuple (X_train, Y_train), where X_train is the N_train by 3 matrix and Y_train is the classification corresponding to the compact object type 1, 2, or 3. Our inputs can thus be represented as

where SC₁ stands for the soft color, HC₁ stands for the hard color, and RelInt₁ stands for the relative intensity of the first data point. The output Y_train is a vector of length N, where N refers to the number of statistically significant data points in the training set. For a prediction for a single data point of an unclassified system, the ML algorithms use three features (SC, HC, and RelInt). Thus, for each row of SC, HC, and RelInt values in Equation (8), the ML models output a probability distribution across the BH, NPNS, pulsar classes. All predictions belonging to the observations for a particular source are then aggregated to generate a mean probability distribution per source as described in Section 4.3. Thus, for an XRB source with M observations, the number of features used to compute the overall probability distribution is 3 × M.

To ensure each model generalizes to XRB examples it has never seen before, we use a cross-validation method described in Section 4.2. This method checks whether the model can be applied to predict out-of-sample data rather than only being capable of memorizing and reproducing the input training set.

A fundamental challenge in ML is that algorithms must perform well on novel, previously unseen inputs—not just examples on which the model was trained (Goodfellow et al. 2016). Being able to perform well on novel, previously unseen inputs is often referred to as generalization.

To estimate the ability of each ML algorithm to generalize to out-of-sample data, we use a cross-validation procedure. At each iteration of the procedure, we remove the CCI points corresponding to a particular source and use the remaining sources as training data, which comprises between 97% and 98% of the entire data set depending on the number of observations belonging to the removed source. Then, we predict the classification of the left-out source. We iterate this procedure for all of the sources in the data set, and the total number of correct classifications out of the 44 sources is indicative of the predictive accuracy on out-of-sample data. This cross-validation procedure is similar to k-fold cross-validation except that the data set is chunked based on the sources, as opposed to k equally sized chunks typically used in cross-validation. Chunking the data by source is more realistic because the ML algorithms will ultimately be used to predict the classification of a new source that the algorithms were not trained on.

4.2.1. BGP Parameter Tuning

4.2.2. KNN Hyperparameter Tuning

4.2.3. SVM Hyperparameter Tuning

For the SVM cross-validation, we chose to sample 100 values across the parameter space of the hyperparameters gamma and cost (C). Because these parameters span a large range of possible ideal values (2⁻¹⁵–2³ for gamma, Chang & Lin 2011; 2⁻¹⁵–2⁵ for C, Chang & Lin 2011), they required an efficient sampling method. We implemented Latin hypercube sampling (Carnell 2020), which is a statistical technique that achieves near-random sampling across parameters (Stein 1987) and is commonly used for sampling large parameter spaces.

Gamma can be thought of as defining the radius of influence of the support vectors. When gamma is very small, the model is too constrained and cannot capture the features or complexity in the data. However, making gamma too large corresponds to the radius of the support vectors only including the support vectors themselves, resulting in overfitting on the training data and large generalization errors.

The other parameter we simultaneously optimized was C, which is a regularization term in the Lagrange formulation of SVMs. C presents a trade-off between the correct classification of training examples and against maximizing the margin of the decision function. Increasing C means that a smaller margin will be deemed acceptable and result in better classification of training points. Lower values of C will create a larger margin and a simpler decision function, which may prevent overfitting but also decreases training accuracy. In this way, C can prevent the SVM model from overfitting, decrease the generalization error, and serve as a regularization parameter.

The best-performing hyperparameter combinations are summarized in Table 4, and the performance across all 100 samples is illustrated in Figure 8. Ultimately, the best-performing models with ∼86.4% accuracy had values of C between ∼0.06 and ∼0.1 and for gamma between ∼0.57 and ∼2.3. For these values of gamma and C, we found that C = 0.655 and gamma = 0.585 provide predictions with the highest probabilities of belonging to the correct class across all three source types. Thus, we chose those hyperparameter values for our final model. We ran our final model 10 times and computed the median and SD across those 10 runs (Table A3).

Zoom In Zoom Out Reset image size

Table 4. SVM Cross-validation Results

Model Parameters			Percentage Classified as X per Source Type Xa			Overall Accuracy
C	Gamma	Time (minutes)	BH	NPNS	Pulsar	Total
0.040	0.484	∼20	50.00%	95.00%	100.00%	84.09%
0.116	0.318	∼20	50.00%	95.00%	100.00%	84.09%
0.184	0.406	∼20	50.00%	95.00%	100.00%	84.09%
0.231	6.437	∼20	50.00%	95.00%	100.00%	84.09%
0.371	0.270	∼20	50.00%	95.00%	100.00%	84.09%
0.397	1.898	∼20	50.00%	95.00%	100.00%	84.09%
0.870	2.859	∼20	50.00%	95.00%	100.00%	84.09%
1.969	0.051	∼20	50.00%	95.00%	100.00%	84.09%
2.162	0.126	∼20	50.00%	95.00%	100.00%	84.09%
0.020	5.141	∼20	58.33%	90.00%	100.00%	84.09%
0.021	7.729	∼20	58.33%	90.00%	100.00%	84.09%
0.062	2.314	∼20	58.33%	95.00%	100.00%	86.36%
0.093	1.572	∼20	58.33%	95.00%	100.00%	86.36%
0.100	0.579	∼20	58.33%	95.00%	100.00%	86.36%
0.655	0.585	∼20	58.33%	95.00%	100.00%	86.36%

Notes. To compute how the model accuracy changed based on values of the hyperparameters C and gamma, we implemented Latin hypercube sampling across the parameter space (C = 2⁻¹⁵–2⁵, gamma = 2⁻¹⁵–2³). The 15 best-performing combinations of C and gamma are summarized here.

^a Here X can be a BH, NPNS, or pulsar. For example, for the BH column, the percentage listed is the percentage of BH sources that are classified as BHs. For the NPNS column, the percentage listed is the percentage of NPNS sources that are classified as NPNSs, etc.

Download table as:  ASCIITypeset image

Because the ML algorithms treat each observation independently and compute their predicted probability distribution, we aggregate the predictions per source to compute an overall probability distribution. Specifically, for a given XRB system with M observations, we compute the overall probability distribution and corresponding SDs after performing the cross-validation by taking the following two steps:

• 1.  
  First, we compute the mean probability distribution across the M observations belonging to that system.
• 2.  
  We then estimate the errors in our models by undertaking 10 independent runs (5 for the BGP due to computational constraints). Independence is achieved by generating a different subsample for each run (see Figure 5). Across those 10 independent runs, we repeat step 1, aggregate those results, and compute the median and sample SD across the 10 probability distributions.

Resampling data to estimate errors has a long history in statistics through methods such as bootstrapping (Efron & Tibshirani 1993) and has also been shown to reliably produce error estimates when applied to ML algorithms (Stracuzzi et al.2017).

In this way, we aggregate the predictions to get a median and SD for each XRB source in Table 1.

For some of the XRB sources, the median probability of belonging to two or more classes is equal (within 1σ). In this case, we chose to assign the sources an indeterminate classification because the algorithms cannot distinguish them within the error.

To examine the role of the quiescent observations in the misclassification of BHs as NPNSs, we ran two independent analyses where we included only the burster observations when sources were in an “on” or “burst” state and all observations for other source types. In these two analyses, we changed how we determined which points were considered in quiescence versus in a “burst” state.

In the first analysis, we applied a uniform cutoff where all points with a relative intensity below 0.4 were removed from burster sources in the training set. The number 0.4 was chosen because most “on” or “burst” observations fall above this threshold. We found that this uniform cutoff significantly increased the BH classification accuracy, leading to twice as many BHs being classified correctly for all three ML methods. In addition, pulsars and NPNSs are correctly classified with rates of 100% and 90%, respectively.

In the second analysis, we changed the exact cutoff based on each burster source’s unique variation in relative intensity. In particular, we computed the mean and SD of the relative intensity for each burster source. Then, we define the cutoff as

where RelInt is the relative intensity. We assume all points below this cutoff to be in quiescence and remove them from the training set. Similar to the 0.4 cutoff analysis, we find that this cutoff_2σ method also leads to twice as many BHs being classified correctly for all three ML methods. Pulsars and NPNSs are correctly classified with rates of 100% and 85%, respectively. Overall, the results across these two quiescence analyses support our hypothesis that significant parameter overlap during the quiescence of BHs and NPNS could play a role in the misclassifications we observe for BHs.

The accuracy of each of the ML algorithms can also be affected by the choice of a specific probability threshold. In the results reported throughout this paper, the class with the largest probability is assigned as the predicted class, regardless of the exact value of the largest probability. However, in some ML applications, a probability threshold is chosen instead to assign the predicted class. We investigated how choosing a probability threshold rather than choosing the class with the largest predicted probability would affect the accuracy of our models by using Receiver Operator Characteristic (ROC) curves.

ROC curves plot the True Positive Rate (TPR) and False Positive Rate (FPR) for various probability thresholds. The TPR and FPR are defined in the following way:

ROC curves can allow us to choose a threshold where the TPR → 1 while the FPR → 0. For a perfect classifier (TPR = 1, FPR = 0), the “elbow” of the ROC curve tightly fits in the upper-left corner. This perfect classifier is plotted in each of the subplots of Figure A1 along with the ROC curves of the ML methods. We can also compute the area under the curve (AUC) for each of these ROC curves. AUC is a threshold-sensitive measure of model accuracy. A perfect classifier has an AUC = 1, whereas a completely random classifier would have an AUC = 0.5.

We listed the AUC values for each XRB class and each ML model in Figure A1. The AUC values for column 3 are the highest, ranging from 0.97 to 0.98, meaning that all three algorithms are nearly perfect classifiers for pulsars. For NPNS, the SVM model outperforms the others with AUC = 0.9 compared to AUC = 0.88 for the KNN and BGP models. Across BHs, the differences between model accuracy for the methods are marginally smaller (BGP AUC = 0.86, SVM AUC = 0.86, KNN AUC = 0.85).

Users of our software may use these ROC curves to choose an appropriate threshold for their own specific purposes. If a user’s primary concern is reducing the FPR to zero, a higher-probability threshold could be chosen based on the ROC curves. For instance, say a user primarily wants to minimize BH false positives with the KNN method, Figure A1(g) indicates that a probability threshold of ∼0.40 would keep the FPR at zero while providing them with a TPR of ∼0.6. Depending on a user’s exact purposes, they may instead consider a few false positives acceptable and would prefer to increase the TPR. In that type of scenario, a lower-probability threshold may be more appropriate. Overall, these ROC curves allow users to choose a threshold that balances the FPR and TPR in a manner that is suitable for their own specific purposes.

Animated three-dimensional versions of Figures 3, 13, 14, and 15 are available. Each of these figures rotates from 0° to 360° and illustrates the spatial patterns in the data.