# Sparse logistic regression for RR Lyrae vs binaries classification

Piero Trevisan Department of Physics, Università di Roma Sapienza, Piazzale
Aldo Moro 2, 00185 Rome, Italy  Ciela, Institute for Computation and
Astrophysical Data Analysis , Montreal, Quebec, Canada  INAF-Osservatorio
Astronomico di Roma, via Frascati 33, 00040 Monte Porzio Catone, Italy  Mario
Pasquato Physics and Astronomy Department Galileo Galilei, University of
Padova, Vicolo dell’Osservatorio 3, I–35122, Padova
Département de Physique, Université de Montréal, Montreal, Quebec H3T 1J4,
Canada
Mila - Quebec Artificial Intelligence Institute, Montreal, Quebec, Canada
Ciela, Institute for Computation and Astrophysical Data Analysis , Montreal,
Quebec, Canada  Gaia Carenini Department of Computer Science, ENS, PSL
Research University, Paris, France
Scuola Universitaria Superiore IUSS Pavia, Palazzo del Broletto, piazza della
Vittoria 15, I-27100 Pavia, Italy
Nicolas Mekhaël Département de Physique, Université de Montréal, Montreal,
Quebec H3T 1J4, Canada  Vittorio F. Braga INAF-Osservatorio Astronomico di
Roma, via Frascati 33, 00040 Monte Porzio Catone, Italy  Space Science Data
Center, via del Politecnico snc, 00133 Roma, Italy  Giuseppe Bono Department
of Physics, Università di Roma Tor Vergata, via della Ricerca Scientifica 1,
00133 Roma, Italy  INAF-Osservatorio Astronomico di Roma, via Frascati 33,
00040 Monte Porzio Catone, Italy  Mohamad Abbas Center for Astro, Particle and
Planetary Physics (CAP3), New York University Abu Dhabi
New York University Abu Dhabi, Physic Department, PO Box 129188, Abu Dhabi,
United Arab Emirates

(Accepted 24 April 2023)

###### Abstract

RR Lyrae (RRL) are old, low-mass radially pulsating variable stars in their
core helium burning phase. They are popular stellar tracers and primary
distance indicators, since they obey to well defined period-luminosity
relations in the near-infrared regime. Their photometric identification is not
trivial, indeed, RRL samples can be contaminated by eclipsing binaries,
especially in large datasets produced by fully automatic pipelines.
Interpretable machine-learning approaches for separating eclipsing binaries
from RRL are thus needed. Ideally, they should be able to achieve high
precision in identifying RRL while generalizing to new data from different
instruments. In this paper, we train a simple logistic regression classifier
on _Catalina Sky Survey (CSS)_ light curves. It achieves a precision of $87\%$
at $78\%$ recall for the RRL class on unseen CSS light curves. It generalizes
on out-of-sample data (ASAS/ASAS-SN light curves) with a precision of $85\%$
at $96\%$ recall. We also considered a L1-regularized version of our
classifier, which reaches $90\%$ sparsity in the light-curve features with a
limited trade-off in accuracy on our CSS validation set and -remarkably- also
on the ASAS/ASAS-SN light curve test set. Logistic regression is natively
interpretable, and regularization allows us to point out the parts of the
light curves that matter the most in classification. We thus achieved both
good generalization and full interpretability.

Statistical Learning; Automated Classification; Variable Stars;

##  1 Introduction

RR Lyrae (RRL) are old ($t>10$ Gyr), low-mass ($M=0.6-0.8M_{\odot}$) variable
stars. They typically pulsate in two main radial modes: fundamental and first
overtone. A small fraction ($5$-$10\%$) pulsates simultaneously in both modes.
The reader interested in a more detailed discussion concerning the possible
occurrence of non-radial modes in RRL is referred to Smolec (2005). RRL obey
to tight Period-Luminosity (PL) relations for wavelengths longer than the red
end of the visible spectrum (Catelan et al., 2004). Moreover, the slope of the
PL relation becomes steadily steeper and the standard dispersion steadily
decreases when moving from shorter to longer wavelenghts. Near infrared PL
relations provide individual distances with an outstanding accuracy of
$\approx$ 3% (Bono et al., 2003; Braga et al., 2015; Beaton et al., 2018).
Hence, they are a fundamental stepping stone in the calibration of secondary
distance indicators, in turn constraining possible systematics in the cosmic
distance scale (Beaton et al., 2016). Compared to other distance indicators
and stellar tracers, RRL have several advantages: they are ubiquitous, having
been identified both in gas-poor and in gas-rich stellar systems; their
lifetimes of the order of 108 yrs ensure that they outnumber other evolved
radial variables such as type II cepheids, anomalous cepheids, and classical
cepheids (Bono, Braga, Pietrinferni, ARA&A submitted); they cover a very broad
metallicity distribution ranging from [Fe/H]$\sim$-3.5 to super-solar iron
abundance (Crestani et al., 2021a, b; Fabrizio et al., 2021), so they can be
easily identified not only in the MW halo, but also in metal-rich stellar
systems such as the MW Disk.

The main issue with relying on RRL to accurately measure distances is due to
their periods and amplitudes partially overlapping with eclipsing binaries,
leading to spurious contamination in RRL samples. In the period-amplitude
plane, also called Bailey diagram, the location of first overtone RRL overlaps
with $\delta$ Scuti variables and eclipsing contact binaries. Furthermore the
so-called W Ursae Majoris (W UMa) variables have sinusoidal light curves
similar to RRc variables, leading to overlap in diagnostics based on Fourier
parameters, especially in cases where the sampling of the light curve is
either poor or the photometry is noisy so that the typical secondary eclipse
of W UMa light curve can be barely detected. While eclipsing binaries can have
similar periods and similar light curve shapes to RR Lyrae, they obviously do
not follow the same period-luminosity relation, thus degrading the quality of
the measured distances. This is the main reason why RRc variables are only
partially used in distance and in metallicity estimates (Contreras Ramos et
al., 2018).

Consequently, accurate classification of RRL versus eclipsing binary stars has
strong implications for our ability to avoid biases in distance determination.
This has long been a priority in catalog building efforts, as in the work by
Drake et al. (2014a), who used a robust version of the skewness of the
magnitude distribution introduced by Kinemuchi et al. (2006) to separate the
two classes in Catalina Sky Survey (CSS) data. It is thus natural to employ
machine learning methods to address this task. Since the large-scale diffusion
of deep learning, more and more classification tasks in astrophysics have been
entrusted to various neural network architectures (see e.g. Hinners et al.,
2018; Lanusse et al., 2018; Gillet et al., 2019; Bialopetravičius et al.,
2019; Trevisan et al., 2020; Chardin & Bianchini, 2021; Lalande & Trani,
2022). Unfortunately, the excellent performance of such tools is often
accompanied by insufficient theoretical guarantees as well as a lack of human
interpretability. Consequently, making use of deep networks for extracting
high-reliability patterns to be interpreted as physical models is challenging.
While physics-informed deep learning (see e.g. Raissi et al., 2019, 2020;
Haghighat & Juanes, 2021) holds the promise to mitigate this issue, in this
paper we choose to favor a statistical foundation of the classifier that
ensures its reliability. Moreover, in the context of scientific research and
astronomy in particular, the ability to interpret our classifier is crucial.
It has recently been convincingly argued that natively interpretable models
are thus preferable to black boxes even when the latter are explained a
posteriori using relevant techniques because explanations cannot possibly be
faithful everywhere or they could simply be used in place of the underlying
black box (Rudin, 2019).

Several works have already appeared in the literature attempting to address
the problem of automatic classification of variable stars, but in general the
main priorities of these were different from interpretability. In a recent
investigation, Pantoja et al. (2022) computed a wide range of features
(including the Feature Analysis for Time Series (FATS features) introduced by
Nun et al., 2015) on light curve data and apply dimensionality reduction and
clustering. They also leverage the results to perform semi-supervised learning
achieving a precision of $90$% for main variable classes with a support-vector
machine, thus reducing the need for human assigned labels. On the other
hand,Castro et al. (2018) focussed their attention on modelling errors, so
they apply Gaussian process regression to light curves -as we do in the
following- and use it as a basis for bootstrapping and calculate FATS
features, using them to train a random forest. The goal of a variety of works
employing a range of neural network architectures is arguably speed and
accuracy of classification, e.g. Becker et al. (2020), Szklenár et al. (2022),
Zhang & Bloom (2021), Becker et al. (2020), Aguirre et al. (2019), and
Debosscher et al. (2007).

The paper is structured as follows. In Section 2, we discuss the dataset,
point out some of its properties obtained by a preliminary data analysis and
describe the subsample involved in the experiments; afterwards, in Section 3,
we present the statistical methods involved in our classification protocol and
provide justification concerning the implementation choices; subsequently in
Section 3.2, we give an overview of the data analysis by specifying the input
representation, the validation scheme and the performance measures; then in
Section 4, we present the experimental results that are discussed further in
Section 6 where we describe, as well, the possible directions for further
work.

##  2 Data

###  2.1 Observations

The training data for this study comes from the Catalina Sky Survey (CSS;
Drake et al., 2013a). A brief summary of the characteristics of the CSS is
given in Pasquato et al. (2022), which shares part of the light curve sample
with this study. A catalog of variable stars observed by CSS was published in
a series of papers (e.g., Drake et al. 2013a, b, 2014b; Torrealba et al. 2015;
Drake et al. 2017). Variable sources were identified based on the Welch-
Stetson variability index ($I_{WS}$; Welch & Stetson 1993). We initially adopt
the nominal period found by the authors through Analysis of Variance (AoV;
Schwarzenberg-Czerny 1989); however we discuss below how this may implicitly
encode information on the assumed classification of a given variable star, and
the consequent pre-processing decisions we took. The combined CSS variable
catalog contains $\sim$ 110,000 variable sources, each of which was observed
on average at $\sim$ 200 epochs. As pointed out also in Pasquato et al.
(2022), CSS does not cover the sky in the proximity of the Galactic plane,
which in principle may introduce a bias in our results. Our choice of CSS is
however justified by its large number of relatively homogeneous light curves,
which is a crucial prerequisite for training most ML models.
Furthermore, to improve the quality of our variable star sample and remove any
potential false positives, we undertook several steps. One of these steps
involved utilizing the RUWE values obtained from Gaia (Gaia Collaboration et
al., 2016). To do this, we selected a circle with a radius of 3 arcsec,
centered at the positions provided by CSS, and cross-matched it with the Gaia
DR3 database (Gaia Collaboration et al., 2022). We then removed all stars with
RUWE values greater than 1.4. This value was chosen as a threshold, as RUWE
values larger than 1.4 are indicative of potential eclipsing binaries, which
are not single stars and therefore not suitable for our RRc sample. By
implementing this methodology, we were able to eliminate potential false
positives and increase the efficiency of our variable star sample. We point
out that CSS images were collected without filters. The conversion to the
V-band was performed by setting up a sample of comparison stars (G0-G8 dwarfs)
from the 2MASS Point Source Catalog (Skrutskie et al., 2006), converting to
the $JHKs$ bands (first step) and then to the $V$ band (second step) (see
Larson et al., 2003; Skrutskie et al., 2006, for details).

Table 1: Adopted data set summary. The initial number of stars (second column) is shown broken down by variable type (first column). The third column shows the number of stars with a good Gaussian process fit according to visual inspection; only these stars were included in the final adopted data set. Median period and amplitude for each variable type are shown as well. Variable Type | N | $N_{gf}$ | Median Period | Median Amplitude
---|---|---|---|---
RRc | 1429 | 1298 | 0.33 | 0.40
EW | 8803 | 7994 | 0.35 | 0.40
EA | 2509 | 2446 | 1.03 | 0.47

We ran our analysis on the subsample of light curves described in Tab. 1; this
includes both intrinsically variable RR Lyrae and eclipsing binaries of two
subclasses: Algol type stars (EA) and W Uma type (EW). Even though in the
initial Catalina data set the light curves of RRab stars were available, we
did not include them in our analysis (except for showing a sample light curve
in Fig. 2) because they are easier to distinguish from eclipsing binary
interlopers due to their characteristic shape. The Bailey diagram for the full
sample is shown in Fig. 1. It is clear from the figure that achieving a
separation between RRc and eclipsing variables is impossible solely based on
period and amplitude information. Typical light curves for different variable
types are shown as red points with error bars in Fig. 2. The superimposed
shaded areas are obtained by fitting the curves with a Gaussian process, as
described in Sect. 3.
To test the performance of our classifier outside of the sample it was trained
on, we relied on a test set comprised of $614$ stars from the All Sky
Automated Survey data (ASAS; Pojmanski, 1997, 2002) and All-Sky Automated
Survey for Supernovae (ASAS-SN; Shappee et al., 2014; Jayasinghe et al.,
2019). This dataset differs from our CSS dataset because it was taken by a
different instrument and also displays a different class balance: RRL are
$12\%$ of our CSS sample and $35\%$ of our ASAS/ASAS-SN test sample.

Figure 1: Period-amplitude plot of our dataset. Only the region up to 1.5 mag
in amplitude and 1 day in period is shown. The stars labeled as RRc are shown
in red, the eclipsing variables are shown in different shades of grey.
Substantial overlap between the two groups is readily evident.

###  2.2 Inclusion criteria

We selected RRc, EW and EA stars only, because RRab stars are less likely to
be contaminated by misclassified binary stars given the distinctive shape of
their light curves. We visualize our final adopted sample in the space of the
first few principal components and directly inspected the light curves of any
outliers. The limited number of outliers corresponds to light curves where the
Gaussian process interpolation (see below) failed to give periodic results and
this was evident upon visual inspection. We speculate that this may be related
to underestimated error bars in the raw data. The final adopted sample
contains $11738$ stars and is described in Tab. 1.

##  3 Pre-processing and classification algorithm

###  3.1 Kernelized Gaussian Process

Gaussian processes have proved to be a very versatile mathematical tool for
interpolation and are progressively becoming more widespread in astronomy.
Formally, a Gaussian process $\mathcal{GC}$ is a distribution over functions
$f(\textbf{x})$:

$f(\textbf{x})\sim\mathcal{GC}(m(\textbf{x}),k(x_{i},x_{j}))$ (1)

where x is the _function value_ , $(x_{i},x_{j})$ are all possible pair in the
input domain, $m(\textbf{x})$ is the _mean function_ and $k(x_{i},x_{j})$ is a
_positive definite covariance function_ also known as the _kernel_ function.
For any finite subset $\mathcal{X}=\\{x_{1},\dots,x_{n}\\}$ of the domain of
x, the marginal distribution is a multivariate Gaussian distribution:

$f(\mathcal{X})\sim\mathcal{N}(m(\mathcal{X}),k(\mathcal{X},\mathcal{X}))$ (2)

with mean vector $\mu=m(\mathcal{X})$ and covariance matrix
$\Sigma=k(\mathcal{X},\mathcal{X})$. Roughly speaking, each input to this
function is a variable correlated with the other variables in the input
domain, as defined by the covariance function. Since functions can have an
infinite input domain, the Gaussian process can be interpreted as an infinite
dimensional Gaussian random variable.
Sampling functions from a Gaussian process boils down to fixing the mean and
covariance functions. In particular, the covariance function $k(x_{i},x_{j})$
models the joint variability of the Gaussian process random variables and
returns the modelled covariance between each pair in $x_{i}$ and $x_{j}$. Its
specification force properties of the distribution over functions $f(x)$, i.e.
by choosing a specific kernel function, it is possible to set prior
information on this distribution. For the purpose of this work, we use a
combination of four distinct kernels. The first of them belongs to the family
of _matern kernels_ and forces the smoothness of the resulting function
allowing for better interpolation. This kernel is given by:

$k(x_{i},x_{j})=\frac{1}{\Gamma(\nu)2^{\nu-1}}\left(\frac{\sqrt{2v}}{l}d(x_{i},x_{j})\right)^{\nu}K_{\nu}\left(\frac{\sqrt{2v}}{l}d(x_{i},x_{j})\right)$
(3)

where $\nu$ is a parameter controlling the smoothness of the resulting
function, $d(.,.)$ is the Euclidean distance, $K_{\nu}(.)$ is a modified
Bessel function and $\Gamma(.)$ is the gamma function. The second is a _white
kernel_ explaining the noise of the signal as independently and identically
normally-distributed. The third one is the one described by the following
equation:

$k(x_{i},x_{j})=\exp\left(-\frac{2\sin^{2}(\pi
d(x_{i},x_{j})/p)}{\ell^{2}}\right)$ (4)

where $\ell$ is the length scale of the kernel, $p$ its periodicity and
$d(.,.)$ the Euclidean distance. This kernel is commonly used to model
functions which repeat themselves exactly. Last, we also added a _constant
kernel_ to establish the mean of the Gaussian process.

#### Application

To obtain equally spaced data by interpolation, we applied the above specified
Gaussian Process to the phase folded light curves, which was phased according
to the nominal period from CSS for EA and EW variables, and to double the
nominal period for the RRc variables. This allows us to compare the two on an
equal footing, because a miscalculated period may result in mis-
classification. In particular an eclipsing binary folded to half its actual
period will be much more likely to be classified as an intrinsic variable and
vice-versa. We repeated the resulting phased light curve three times to
enforce periodicity. The kernel function selected for this study included the
exponential component with a length scale $l=2$ and a periodicity $p=3$,
combined with the _constant kernel_ and a _white kernel_ components with noise
levels of 5000. We also included the _matern kernel_ with a length scale $l=5$
and a $\nu=3/2$. It is worth noting that the kernel parameters were kept
constant throughout the study, and no optimization was performed. This means
that only one kernel function with fixed parameters was used for all stars. We
used the Python library scikit-learn (Pedregosa et al., 2011) based on
Rasmussen (2003) to fit our Gaussian process model. The resulting interpolated
light curve is sampled at $300$ equally spaced points in phase, normalized
between $0$ and $1$. Finally, we phase-shifted to make the maximum coincide
with the beginning of the period and then, we selected 100 consequent points.
By operating in this way we avoided any possible edge effect that could affect
our final dataset. Not only, we want our classifier to rely only on the light
curve topology, we select one period in order to have all the data in a box
$[0,1]\times[0,1]$. Thus each curve is represented by $100$ features
corresponding to normalized values of the magnitude at given phases. The
results of the Gaussian process fit are shown in Fig. 2.

Figure 2: Example phase-folded light curves of different variable types, shown
as red points with error bars: RRab (top left panel), RRc (top right panel),
EW (bottom left panel) and EA (bottom right panel). The curves are repeated
for three periods and were fit with a Gaussian process following the procedure
described in Sect. 3. The solid blue line is the Gaussian process fit with the
corresponding variance (light blue area).

For purposes of visualization, we calculated the first three principal
components of our data set in feature space. These are shown in Fig. 3, where
the eclipsing binaries (including both EA and EW) and the RR Lyrae are shown
in different colors. Principal components have been calculated on normalized
light curves, so they depend only on the shape of the light curve and not on
amplitude and period. Fig. 3 thus shows that the light curves of RRc and
eclipsing binaries differ in shape enough that the two groups appear as
clusters in principal component space. Even though these clusters are not
easily separated by eye, their degree of overlap is less than that witnessed
in the Bailey diagram in Fig. 1. In principle one could build a classifier by
placing a cut along one of the principal components shown in Fig. 3; this
would result in a separating hyperplane in feature space. However, having
access to the full feature space rather than to just the first three principal
components, a classification algorithm is likely to perform substantially
better than any fiducial cut we could place by hand. Moreover, as we will
explain below, the algorithm we chose actually learns a separating hyperplane
in feature space, which can be intuitively understood as an optimal cutoff on
a linear combination of the features.

Figure 3: Distribution of our data in principal component space. The first
three components are shown. The RRc variables are shown in orange, the
eclipsing binaries in light blue.

###  3.2 (Regularized) Logistic Regression

Logistic Regression is a widely used classification algorithm. Contrary to
linear regression, there is no closed-form solution and one needs to solve it
thanks to iterative convex optimization algorithms. Given a input
$X_{i}\in\mathbb{R}^{d}$, logistic regression classify it, i.e. returns a
boolean label $Y_{i}\in\\{0,1\\}$, thanks to classification rule
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ such that $f(X_{i})\geq 0\Rightarrow
Y_{i}=1$ and $f(X_{i})<0\rightarrow Y_{i}=0$. In the linear case, $f$ is a
linear operator of the form $f_{\beta}:x\rightarrow x^{T}\beta$.
Intuitively, this algorithm seeks to identify the best linear separation for
input data. However, this is achieved by a hyperplane that in general is not
orthogonal with any of the coordinates of feature space, resulting in a model
that has no non-zero coefficients. To obtain a smaller model by forcing some
coefficients to zero we considered regularizing our classifier. This has the
bonus of mitigating overfitting, even though this is not a major concern of
ours at this stage; the main reason for inducing sparsity is to further
promote interpretability. There are several possible choices for the
regularizer; however, when looking for a succinct description of the
separating hyperplane, $L_{1}-$regularization performs better. Roughly
speaking in fact it introduces a penalty upon the number of relevant features
considered for the classification.

#### Application

In our case, regularization allows us to limit the number of features included
in the model, allowing us to understand which parts of a given light curve
really matter for classification. We applied the logistic regression with the
scikit-learn python package.

##  4 Results

We trained our logistic regression model on a training set comprising $70\%$
of the CSS sample and evaluated its performance on a validation set containing
the remaining $30\%$. The assignment to validation was stratified using
_StratifiedShuffleSplit_ function from scikit-learn. We measured the
performance according to a variety of metrics, as follows: accuracy, ratio of
stars correctly classified over the total; precision, the fraction of stars
classified as members of a given class who were actually members of that
class; recall, the fraction of actual members of a given class that were
classified as such; and F-score, the harmonic mean of precision and recall. We
also obtained and plotted the precision-recall curve for our classifiers,
which shows precision as a function of recall when the threshold used to
convert the classifier’s probabilistic prediction to a hard classification is
varied. With the exception of accuracy, the metrics discussed above are class-
dependent, i.e. they are in general different if one considers e.g. the
fraction of actual RRc among the stars classified as RRc (precision for the
RRc class) or the fraction of actual eclipsing binaries among the stars
classified as eclipsing binaries (precision for the eclipsing binary class). A
given decrease in precision for one class does not in general have the same
cost as the same decrease for the other; polluting the sample of RRc on which
distance measurements are based with eclipsing binaries is worse than the
opposite. We thus report these metrics for each class.

The first model we present is a plain logistic regression, fit without a
penalty term. The confusion matrix achieved on our test set is shown in Fig.
4. The performance metrics discussed above are reported in Tab. 2 for a
standard threshold of 0.5. A more detailed view of the classifier performance
at various thresholds is available in Fig. 5 through a precision-recall curve.

Table 2: Performance metrics broken down by class (first column) of our non-penalized (second column) and L1 regularized model (third column) on the validation set. Metrics | $\text{Score}$ | $\text{Score}_{92}$
---|---|---
Accuracy | 0.96 | 0.95
Recall (RRc) | 0.78 | 0.70
Recall (EB) | 0.98 | 0.99
Precision (RRc) | 0.87 | 0.87
Precision (EB) | 0.97 | 0.96
F-score (RRc) | 0.82 | 0.77
F-score (EB) | 0.98 | 0.97

As discussed above, an advantage of logistic regression over other methods is
the interpretability of the model. Compared to e.g. a deep neural network that
can have millions of weights, our model is entirely specified by the values of
the coefficients associated to each feature. Since each feature is a magnitude
at a given phase we can visualize the coefficients as a function of the
associated phase, interpreting them as weights on the light curve to be
classified. This is shown in Fig. 6. A positive coefficient corresponds to an
increased probability of being classified as an RR Lyrae when a given feature
is increased, and vice-versa for a negative coefficient. This is similar to
interpreting a classifier using a saliency map as in e.g. Peruzzi et al.
(2021) but with the enormous advantage that we do not need to compute the
gradient locally for each instance, since our model is linear and the
direction of the gradient is the same for every instance. In the absence of
regularization, most coefficients differ from zero, but the coefficients are
far from independent from each other, leading to contiguous regions of the
light curve that are weighted similarly by the classifier. In the bottom panel
of Fig. 6 such regions appear.

|
---|---
Figure 4: Left: confusion matrix of our non-penalized classifier on the
validation dataset (comprised of CSS light curves not seen in training).
Right: confusion matrix of our L1 regularized classifier on the same
validation dataset.

##  5 Discussion

In Fig. 6 we show the coefficients of our non-penalized logistic regression
model (blue bars, top panel) together with the mean RRc light curve (orange
points) and the mean eclipsing binary light curve (green points). The
qualitative interpretation of the coefficients is as follows: increasing the
value of the corresponding feature (moving one of the points in the light
curve downwards) increases the probability of being classified as a RRc if the
coefficient is positive, decreases it otherwise. Most coefficients are nonzero
in the non-penalized model, making interpretation harder. The penalized model
drastically reduces the number of non-zero coefficients, showing which parts
of the curve are actually important for the classifier. These coincide roughly
with the approach towards the minima of the light curve. This is showing that
the classifier relies on the quicker rise (or fall) of the RRc light curve.

The results obtained in the experimental analysis may seem at a first glance
to be unsatisfactory when compared to results of other machine learning
methods: in fact, as a baseline comparison we trained a random forest (with
the default parameters of sklearn) obtaining a precision on the RRc of 0.91
and a recall of 0.96. However, it should be noted that there are trade-offs
between the accuracy and interpretability of a model (Bell et al., 2022;
Sarkar et al., 2016; van der Veer et al., 2021). As there is no unambiguous
definition for interpretability, there is no canonical quantification for this
trade-off; however, important findings from the field of statistics confirm
the existence of a statistical cost of interpretability Dziugaite et al.
(2020). In our investigation, we strengthen the interpretability side, in fact
we do not simply choose a classical statistical model such as regression, but
increase its sparsity. Augmenting sparsity has an interest per se, for a
survey around sparsity see Hoefler et al. (2021).

---
Figure 5: Precision-recall curve for our non-penalized model (blue) and for
the L1 regularized version of our model (red) for the RRc and EB classes
respectively on the left and right panel. Only the region with recall greater
than $0.5$ is shown. The curve was calculated on the validation set (comprised
of CSS light curves not seen in training). Note the difference on the y-axis
between the two panels.

Adding a sparsity-inducing penalty to the usual loss function optimized by
logistic regression allows us to increase the number of coefficients that
equal zero. This makes the classifier easier to understand and reduces the
risk of overfitting. We chose to use a $L^{1}$ penalty, and we varied the
weight assigned to the penalty term in the resulting loss function, leading to
different regularization strengths. The effect on sparsity and on the
performance of our classifier is shown in Fig. 7; this figure was used to
select a regularization strength around the knee of the performance curve. Our
choice resulted in a model with $92\%$ of the coefficients equal to zero. The
relevant metrics are summarized in Tab. 2; regularization reduced our recall
on the RRc class from $0.78$ to $0.70$ but did not impact precision.

In order to see how well our model predict on an unseen data from a different
catalogue, we tested both classifiers on ASAS/ASAS-SN stars. The relevant
metrics are summarized in Tab. 3 and the confusion matrices are shown in Fig.
8.

---
Figure 6: Top: values of the linear coefficients of the logistic regression
for the non-penalized model (blue bars). A positive coefficient means that an
increased value of the magnitude at the corresponding phase would increase the
confidence of classification as RRc. Meanwhile, a negative coefficient
increase the confidence of classification as EB. Mean light curve of all RRc
in the whole data set (orange), mean light curve of all EW (green) and EA
(red) are shown as dots. Bottom: same for the L1 regularized model: the
regularization strength was chosen so that $92\%$ of the coefficients became
zero. Figure 7: Performance metrics of our classifier (y axis) as a function
of the fraction of coefficients that equal exactly zero (sparsity, x-axis). At
over $90\%$ we reach a knee where the recall for the RRc class drops
precipitously but for lower sparsity the three metrics are approximately
constant.
---
Figure 8: Left: confusion matrix of our non-penalized classifier on the test dataset; right: confusion matrix of our L1 regularized classifier on the test dataset. Table 3: Performance metrics broken down by class (first column) of our non-penalized (second column) and L1 regularized model (third column) on the test set. Metrics | $\text{Score}$ | $\text{Score}_{92}$
---|---|---
Accuracy | 0.92 | 0.92
Recall (RRc) | 0.96 | 0.95
Recall (EB) | 0.91 | 0.97
Precision (RRc) | 0.85 | 0.83
Precision (EB) | 0.98 | 0.90
F-score (RRc) | 0.90 | 0.89
F-score (EB) | 0.94 | 0.93

##  6 Summary and conclusions

Distance determination relying on RR Lyrae stars as standard candles is
negatively affected by contamination due to eclipsing binaries. In this paper
we have shown that a simple logistic regression classifier trained on CSS
light curves can separate RRc variables from eclipsing binaries reaching
$0.87$ precision at $0.78$ recall on a validation sample of unseen CSS light
curves. On ASAS/ASAS-SN light curves recall becomes $0.96$ and precision
becomes $0.85$, showing very good generalization even to data taken by a
different instrument. It has to be noted that the fraction of RRcs in the
ASAS/ASAS-SN sample is different from the CSS sample ($12\%$ VS $35\%$). Thus
good generalization is due to our deliberate choice of a simple model which is
unlikely to overfit the data.
Since our features are normalized magnitudes at different phases, our
classifier is making use only of the shape of the light curve, without
receiving any information about its amplitude, absolute mean magnitude, or
period. We are also able to visualize the logistic regression coefficients
corresponding to the normalized magnitude at each phase point. Through our
sparsification approach, we are able to identify how the classifier is using
the light curve shape to cast its prediction. The features with non-zero
coefficients correspond to the phase approaching the minima, suggesting that
the steepness of the fall is used by the classifier. In terms of physical
explanation, non-contact eclipsing binaries are in fact characterized by much
flatter maxima (both the components are visible for an extended period of
time) broken by steep minima (when one of the two components is occluded).
We provide a precision-recall curve for the full classifier as well as for the
penalized one, which in principle allows for a selection of the classification
threshold based on the relative cost of false positives (eclipsing binaries
classified as RR Lyrae) and false negatives (vice versa).
We conclude with a few caveats and perspective improvements for the future.
Our classifiers rely on turning the observed light curve, which is sampled at
irregular intervals, into an equally spaced time series by means of
interpolation. We carried out this interpolation based on Gaussian processes,
which depend on several parameters to fully specify the covariance structure
in terms of a kernel function. In the work described in this paper we fixed
those parameters once and for all, even though this resulted in a handful of
bad fits which were later discarded. There is room for improvement in this
pre-processing of the data, especially if it is to be adapted to different
data-sets. In this content it is worth mentioning that current Zwicky
Transient Facility and near future Vera C. Rubin Observatory (Masci et al.
2018, Ivezić et al. 2019) long-term variability surveys (will) provide well
sampled, long-baseline light curves in several optical- and near-infrared
bands (see e.g. LSST Science Collaboration et al. 2017). Thus providing the
unique opportunity to use physically rooted diagnostics (amplitude in
different photometric bands) to properly separate intrinsic from geometrical
variables.

Finally, our model should be compared to more expressive ones, such as a 1D
convolutional neural network, which might attain better accuracy at the price
of substantially increased model complexity.

## Acknowledgments

This paper is supported by the Fondazione ICSC , Spoke 3 Astrophysics and
Cosmos Observations, National Recovery and Resilience Plan (Piano Nazionale di
Ripresa e Resilienza, PNRR) Project ID CN$00000013$ Italian Research Center on
High-Performance Computing, Big Data and Quantum Computing funded by MUR
Missione $4$ Componente $2$ Investimento $1.4$: Potenziamento strutture di
ricerca e creazione di campioni nazionali di R$\&$S (M4C2-19 ) - Next
Generation EU (NGEU). This project has received funding from the European
Union’s Horizon research and innovation program under the Marie Skłodowska-
Curie grant agreement No. $896248$. This research was made possible by a
generous donation from Eric and Wendy Schmidt by recommendation of the Schmidt
Futures program

We would like to thank the anonymous referee for their insightful comments and
suggestions, which greatly improved the content and the readability of an
early version of this manuscript.

## References

  * Aguirre et al. (2019) Aguirre, C., Pichara, K., & Becker, I. 2019, MNRAS, 482, 5078, doi: 10.1093/mnras/sty2836
  * Beaton et al. (2016) Beaton, R. L., Freedman, W. L., Madore, B. F., et al. 2016, ApJ, 832, 210, doi: 10.3847/0004-637X/832/2/210
  * Beaton et al. (2018) Beaton, R. L., Bono, G., Braga, V. F., et al. 2018, Space Sci. Rev., 214, 113, doi: 10.1007/s11214-018-0542-1
  * Becker et al. (2020) Becker, I., Pichara, K., Catelan, M., et al. 2020, MNRAS, 493, 2981, doi: 10.1093/mnras/staa350
  * Bell et al. (2022) Bell, A., Solano-Kamaiko, I., Nov, O., & Stoyanovich, J. 2022, in 2022 ACM Conference on Fairness, Accountability, and Transparency, FAccT ’22 (New York, NY, USA: Association for Computing Machinery), 248–266, doi: 10.1145/3531146.3533090
  * Bialopetravičius et al. (2019) Bialopetravičius, J., Narbutis, D., & Vansevičius, V. 2019, A&A, 621, A103, doi: 10.1051/0004-6361/201833833
  * Bono et al. (2003) Bono, G., Caputo, F., Castellani, V., et al. 2003, MNRAS, 344, 1097, doi: 10.1046/j.1365-8711.2003.06878.x
  * Braga et al. (2015) Braga, V. F., Dall’Ora, M., Bono, G., et al. 2015, ApJ, 799, 165, doi: 10.1088/0004-637X/799/2/165
  * Castro et al. (2018) Castro, N., Protopapas, P., & Pichara, K. 2018, AJ, 155, 16, doi: 10.3847/1538-3881/aa9ab8
  * Catelan et al. (2004) Catelan, M., Pritzl, B. J., & Smith, H. A. 2004, ApJS, 154, 633, doi: 10.1086/422916
  * Chardin & Bianchini (2021) Chardin, J., & Bianchini, P. 2021, MNRAS, 504, 5656, doi: 10.1093/mnras/stab737
  * Contreras Ramos et al. (2018) Contreras Ramos, R., Minniti, D., Gran, F., et al. 2018, ApJ, 863, 79, doi: 10.3847/1538-4357/aacf90
  * Crestani et al. (2021a) Crestani, J., Fabrizio, M., Braga, V. F., et al. 2021a, ApJ, 908, 20, doi: 10.3847/1538-4357/abd183
  * Crestani et al. (2021b) Crestani, J., Braga, V. F., Fabrizio, M., et al. 2021b, arXiv e-prints, arXiv:2104.08113.  https://arxiv.org/abs/2104.08113
  * Debosscher et al. (2007) Debosscher, J., Sarro, L. M., Aerts, C., et al. 2007, A&A, 475, 1159, doi: 10.1051/0004-6361:20077638
  * Drake et al. (2013a) Drake, A. J., Catelan, M., Djorgovski, S. G., et al. 2013a, ApJ, 763, 32, doi: 10.1088/0004-637X/763/1/32
  * Drake et al. (2013b) —. 2013b, ApJ, 765, 154, doi: 10.1088/0004-637X/765/2/154
  * Drake et al. (2014a) Drake, A. J., Graham, M. J., Djorgovski, S. G., et al. 2014a, ApJS, 213, 9, doi: 10.1088/0067-0049/213/1/9
  * Drake et al. (2014b) —. 2014b, ApJS, 213, 9, doi: 10.1088/0067-0049/213/1/9
  * Drake et al. (2017) Drake, A. J., Djorgovski, S. G., Catelan, M., et al. 2017, MNRAS, 469, 3688, doi: 10.1093/mnras/stx1085
  * Dziugaite et al. (2020) Dziugaite, G. K., Ben-David, S., & Roy, D. M. 2020, Enforcing Interpretability and its Statistical Impacts: Trade-offs between Accuracy and Interpretability, arXiv, doi: 10.48550/ARXIV.2010.13764
  * Fabrizio et al. (2021) Fabrizio, M., Braga, V. F., Crestani, J., et al. 2021, ApJ, 919, 118, doi: 10.3847/1538-4357/ac1115
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al. 2016, Astronomy & Astrophysics, 595, A2
  * Gaia Collaboration et al. (2022) —. 2022, Astronomy & Astrophysics
  * Gillet et al. (2019) Gillet, N., Mesinger, A., Greig, B., Liu, A., & Ucci, G. 2019, MNRAS, 484, 282, doi: 10.1093/mnras/stz010
  * Haghighat & Juanes (2021) Haghighat, E., & Juanes, R. 2021, Computer Methods in Applied Mechanics and Engineering, 373, 113552, doi: 10.1016/j.cma.2020.113552
  * Hinners et al. (2018) Hinners, T. A., Tat, K., & Thorp, R. 2018, AJ, 156, 7, doi: 10.3847/1538-3881/aac16d
  * Hoefler et al. (2021) Hoefler, T., Alistarh, D., Ben-Nun, T., Dryden, N., & Peste, A. 2021, Sparsity in Deep Learning: Pruning and growth for efficient inference and training in neural networks, arXiv, doi: 10.48550/ARXIV.2102.00554
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, The Astrophysical Journal, 873, 111
  * Jayasinghe et al. (2019) Jayasinghe, T., Stanek, K. Z., Kochanek, C. S., et al. 2019, MNRAS, 485, 961, doi: 10.1093/mnras/stz444
  * Kinemuchi et al. (2006) Kinemuchi, K., Smith, H. A., Woźniak, P. R., McKay, T. A., & ROTSE Collaboration. 2006, AJ, 132, 1202, doi: 10.1086/506198
  * Lalande & Trani (2022) Lalande, F., & Trani, A. A. 2022, arXiv e-prints, arXiv:2206.12402.  https://arxiv.org/abs/2206.12402
  * Lanusse et al. (2018) Lanusse, F., Ma, Q., Li, N., et al. 2018, MNRAS, 473, 3895, doi: 10.1093/mnras/stx1665
  * Larson et al. (2003) Larson, S., Beshore, E., Hill, R., et al. 2003, in AAS/Division for Planetary Sciences Meeting Abstracts, Vol. 35, AAS/Division for Planetary Sciences Meeting Abstracts #35, 36.04
  * LSST Science Collaboration et al. (2017) LSST Science Collaboration et al. 2017, arXiv preprint arXiv:1708.04058
  * Masci et al. (2018) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2018, Publications of the Astronomical Society of the Pacific, 131, 018003
  * Nun et al. (2015) Nun, I., Protopapas, P., Sim, B., et al. 2015, arXiv e-prints, arXiv:1506.00010.  https://arxiv.org/abs/1506.00010
  * Pantoja et al. (2022) Pantoja, R., Catelan, M., Pichara, K., & Protopapas, P. 2022, MNRAS, 517, 3660, doi: 10.1093/mnras/stac2715
  * Pasquato et al. (2022) Pasquato, M., Abbas, M., Trani, A. A., et al. 2022, ApJ, 930, 161, doi: 10.3847/1538-4357/ac5624
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Peruzzi et al. (2021) Peruzzi, T., Pasquato, M., Ciroi, S., et al. 2021, A&A, 652, A19, doi: 10.1051/0004-6361/202038911
  * Pojmanski (1997) Pojmanski, G. 1997, Acta Astron., 47, 467
  * Pojmanski (2002) —. 2002, Acta Astron., 52, 397.  https://arxiv.org/abs/astro-ph/0210283
  * Raissi et al. (2019) Raissi, M., Perdikaris, P., & Karniadakis, G. E. 2019, Journal of Computational Physics, 378, 686, doi: 10.1016/j.jcp.2018.10.045
  * Raissi et al. (2020) Raissi, M., Yazdani, A., & Karniadakis, G. E. 2020, Science, 367, 1026, doi: 10.1126/science.aaw4741
  * Rasmussen (2003) Rasmussen, C. E. 2003, in Summer school on machine learning, Springer, 63–71
  * Rudin (2019) Rudin, C. 2019, Nature Machine Intelligence, 1, 206
  * Sarkar et al. (2016) Sarkar, S., Weyde, T., d’Avila Garcez, A. S., et al. 2016, in CoCo@NIPS
  * Schwarzenberg-Czerny (1989) Schwarzenberg-Czerny, A. 1989, MNRAS, 241, 153, doi: 10.1093/mnras/241.2.153
  * Shappee et al. (2014) Shappee, B. J., Prieto, J. L., Grupe, D., et al. 2014, ApJ, 788, 48, doi: 10.1088/0004-637X/788/1/48
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Smolec (2005) Smolec, R. 2005, Acta Astron., 55, 59.  https://arxiv.org/abs/astro-ph/0503614
  * Szklenár et al. (2022) Szklenár, T., Bódi, A., Tarczay-Nehéz, D., et al. 2022, ApJ, 938, 37, doi: 10.3847/1538-4357/ac8df3
  * Torrealba et al. (2015) Torrealba, G., Catelan, M., Drake, A. J., et al. 2015, MNRAS, 446, 2251, doi: 10.1093/mnras/stu2274
  * Trevisan et al. (2020) Trevisan, P., Pasquato, M., Ballone, A., & Mapelli, M. 2020, MNRAS, 498, 5798, doi: 10.1093/mnras/staa2663
  * van der Veer et al. (2021) van der Veer, S. N., Riste, L., Cheraghi-Sohi, S., et al. 2021, Journal of the American Medical Informatics Association, 28, 2128, doi: 10.1093/jamia/ocab127
  * Welch & Stetson (1993) Welch, D. L., & Stetson, P. B. 1993, AJ, 105, 1813, doi: 10.1086/116556
  * Zhang & Bloom (2021) Zhang, K., & Bloom, J. S. 2021, MNRAS, 505, 515, doi: 10.1093/mnras/stab1248