{'fraction_non_alphanumeric': 0.029952021932830707, 'fraction_numerical': 0.017563399588759424, 'mean_word_length': 4.5944210122699385, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Search for high-amplitude Delta Scuti and RR Lyrae stars in Sloan Digital Sky Survey Stripe 82 using principal component analysis

We propose a robust principal component analysis (PCA) framework for the exploitation of multi-band photometric measurements in large surveys. Period search results are improved using the time series of the first principal component due to its optimized signal-to-noise ratio.The presence of correlated excess variations in the multivariate time series enables the detection of weaker variability. Furthermore, the direction of the largest variance differs for certain types of variable stars. This can be used as an efficient attribute for classification. The application of the method to a subsample of Sloan Digital Sky Survey Stripe 82 data yielded 132 high-amplitude Delta Scuti variables. We found also 129 new RR Lyrae variables, complementary to the catalogue of Sesar et al., 2010, extending the halo area mapped by Stripe 82 RR Lyrae stars towards the Galactic bulge. The sample comprises also 25 multiperiodic or Blazhko RR Lyrae stars.


INTRODUCTION
During the last decade, wide-area and multi-epoch surveys have started to play a major role in astronomical research. Developments in astronomical instrumentation and in space observation techniques, together with the rapidly growing data storage facilities and the broadly available software for combining these data provide an enormous wealth of information. The traditional manual procedures must be replaced by quick automated methods for pre-processing, selection and analysis.
Analysis of variable stars has benefited greatly from data collected by large-scale surveys such as ASAS (Pojmanski 2002(Pojmanski , 2003, OGLE (Udalski, Kubiak & Szymanski 1997), MACHO (Alcock et al. 1997), or EROS (Aubourg et al. 1993;Spano et al. 2011). Studies of different, sometimes rare classes of objects become possible with unprecedentedly large sets of objects. These studies require specific pre-processing: a preliminary classification of the observed E-mail: Maria.Suveges@unige.ch objects, in order to enable an efficient extraction of homogeneous samples. The first step in this procedure is variability detection, providing a set of candidate variable stars. The next is characterization of the observed time series by the calculation of some numerical summaries of the observed time series. These are usually statistical moments (e.g. mean magnitude, skewness, kurtosis), derived quantities from period search (e.g. amplitudes and relative phases of harmonic components), and some astrophysical parameters such as colours. These parameters, called generally attributes in the context of classification, are then used to estimate the types of the objects. The volume of data requires fast automated data mining techniques like, for instance, those proposed by Eyer & Blake (2002. Most recently, Random Forest (Breiman 2001;Dubath et al. 2011;Richards et al. 2011;Rimoldini et al. 2012), multistage Bayesian networks (Debosscher et al. 2007;Sarro et al. 2009) and gradient boosting methods (Richards et al. 2011) were tested for this purpose with promising results.
Unfortunately, automatically distributed class labels can be less reliable than the results of a careful human in-spection. A more efficient use of the information contained in the data can improve on this. The goal of this study is to consider a new way of including colour information into automated classification procedures. Although variable stars with different origin of variability and different physical parameters show different amplitudes and light curve patterns depending on wavelength during their cycle, this variation is not easy to summarize in a concise numerical form, as it is a function of phase. We attempt to find summaries of these colour variations, and to use it as complementary information to the usual attributes in automated classification of multi-band survey data.
Principal component analysis (PCA, Jolliffe 2009;Hastie, Tibshirani & Friedman 2009) gives precisely such summaries. It considers the vector of observations in M bands at one time as a point in an M -dimensional space, and looks for a decomposition of the space of observations into directions with maximal variance. The first direction is called the first principal component or PC1. The projection of any point to the PC1 direction is a linear combination of the simultaneous measurements with constant coefficients, which is called PC1 score. It provides several potential advantages. It improves signal-to-noise ratio when searching for periods on the time series of the projections on the direction of PC1. It yields a variability criterion based on the presence of correlated variations across filters. Also, the direction of PC1 is characteristic to the origin of variability and to the physical properties of the star, and is useful to separate eclipsing binaries from some pulsating variables with symmetric light curves such as RR Lyrae type c (hereafter RRc). We demonstrate these advantages on 5-band time series from a flux-limited sample of Stripe 82 objects of the Sloan Digital Sky Survey (SDSS).
However, applying PCA on astronomical time series is not straightforward. Usually, the data in different filters have different error levels, depending on the measured magnitudes too. Outliers and non-normality can also strongly influence PCA, which is built entirely on the assumption of normality. Since usually not all type of errors can be anticipated and accounted for in advance, the robustness of the methods is very desirable: a good automated method should be able to function acceptably well even in the presence of contaminations from various unknown sources. We propose a combination of variance stabilizing transformation with robust principal component analysis to minimize the impact of these issues.
After constructing the methodology, we use it to select samples of candidate RR Lyrae and high-amplitude δ Scuti (HADS) from SDSS Stripe 82. RR Lyrae are interesting as they are Population II halo structure tracers, and obey welldetermined period-luminosity relations (e.g. Smith 1995). The HADS stars are also pulsating variable stars with spectral type from late A to early F, occurring in the instability strip on and just above the main sequence below the RR Lyrae (e.g. Breger 1980;Clement et al. 2001;Pigulski et al. 2006). Their pulsation is in majority radial-mode, though for some, there are indications of non-radial modes (Mazur, Krzeminski & Thompson 2003;Poretti 2003). They also satisfy period-luminosity relations depending on metallicity and oscillation mode (Nemec, Nemec & Lutz 1994;McNamara 1995McNamara , 1997McNamara , 2000, but Population I (δ Scuti) and II (SX Phoenicis) objects are mixed in the group, and they are affected also by mode identification problems. In Stripe 82 data, the confirmed RR Lyrae stars published by Sesar et al. (2010) will be used as a performance test and a training set for the PCA-based methodology. For HADS stars, we build a training set to obtain a clean sample in Stripe 82. Finally, we present a list of 129 new RR Lyrae and 132 HADS candidate variables.
In Section 2, we briefly present SDSS Stripe 82 (York et al. 2000). Section 3 summarizes the statistical background: variance stabilizing transformation and principal component analysis with its robust version, the tested period search variants, and the Random Forest classifier. Section 4 describes the results obtained by PCA in variability detection, period search, characterisation and classification. A short summary and a discussion follows in Section 5.

SDSS STRIPE 82
The Sloan Digital Sky Survey Data Release 7 (Abazajian et al. 2009) provides five-band (u, g, r, i and z) photometry of more than 11000 deg 2 of the sky. The photometric errors are around 0.02 mag for g < 16 mag, and around 0.04 mag for g ∼ 21 mag (see Fig. 2 of Sesar et al. 2010). Equivalent values for the noisiest u band are 0.02 mag at the bright end and 0.05 mag around 20 mag. The 95% completeness limits of the images are u = 22.0, g = 22.2, r = 22.2, i = 21.3, z = 20.5 (Abazajian et al. 2004). One of the southern regions, Stripe 82, was observed repeatedly during the first phase of SDSS (SDSS-I) and the following SDSS-II Supernova Survey (Bramich et al. 2008;Frieman et al. 2008). The measurements in the five bands were taken quasi-simultaneously, with around 1.2 minute time difference between band records.
Our data set is a flux-limited subset of Sesar et al. (2007), separated into around 68,000 variable and 200,000 non-variable objects with median g magnitude brighter than 20.5 mag. The variable flag used there was based on two criteria: (1) for the root mean square scatter in g and r, σr > 0.05 and σg > 0.05; (2) for the reduced chi-square in r and g, ξr 3 and ξg 3, respactively. This selection will be referred as 'flagged as variable' or 'flagged non-variable' hereafter, and will serve in comparisons with the PCA-based detection method. Details about the construction and testing of the catalogue are described in Ivezić et al. (2007). Sesar et al. (2010) published 483 confirmed RR Lyrae variables in Stripe 82; our subset of data contains 450 of these. This makes Stripe 82 data particularly adapted to test our methods, as we can check their performance by direct comparison to previous analyses of the same data.

Principal Component Analysis
Principal component analysis is a statistical tool developed for finding orthogonal directions of maximal variance in high-dimensional data sets. It is assumed that directions of large variance are of particular interest: in signal processing, they may contain the bulk of the information transmitted by a signal; in image analysis, they may concisely summarise the most characteristic pattern forms; or in classification, they might coincide with the directions that best separate some distinct groups of objects. In order to find them, we consider the point cloud of multi-filter observations plotted against each other, regardless of times. The goal is to find an orthogonal coordinate system adapted to this point cloud: the first axis is fixed so that it points to the direction of maximal variance. Next, the points are decomposed into projections to this first axis and to its orthogonal subspace. Then the second axis is chosen to point into the direction of largest variance in the orthogonal subspace. These steps are repeated until an orthogonal basis spanning the original space is found. Mathematically, we seek for a successive decomposition of the space of the observables into orthogonal directions to which the projections of the points have the highest residual variance.
Suppose that we observed M sequences X11, X21, . . . , XN1; . . . ; X1M , X2M , . . . , XNM at N times, with the values Xi1, Xi2, . . . , XiM observed simultaneously. This corresponds to a sequence of observations made in an M -dimensional space, where a joint observation Xi1, Xi2, . . . , XiM is represented by a point in the Mdimensional state space. If we define the matrix X so that its columns are the M sequences, and the rows correspond to the different M -dimensional observations, then we can write the sample covariance matrix as X T X/N , with the superscript T denoting transposition. Finding orthogonal directions of maximal variance turns out to be equivalent to finding the eigendecomposition Here the columns of V are the eigenvectors, or in other terms, the principal directions of X, and D 2 is a diagonal matrix with non-negative values d 2 1 , d 2 2 , . . . , d 2 M , representing the variances of the projections of the points on the principal directions. By convention, the order of the directions is such that d 2 1 > d 2 2 > · · · > d 2 M , and therefore the first principal component corresponds to the direction of the maximal variation. The vector of projections of the points on the first principal direction can be written as the linear combination z1 = Xv1, where v1 is the first column of the matrix V.
How to apply this for variable star analysis? For M time series of a star consisting of N simultaneously taken points, we consider the M -dimensional point cloud of the observations, as shown in Fig. 1. We intend to apply principal component analysis in order to find the direction where variability is maximal, so that we can find more easily the period of the variable star. However, the application immediately hits an obstacle visible in Fig. 1, the different average noise in the different bands. Even if a star shows coherent deterministic variations across the bands, this can be easily masked by the noise in one of the bands (in SDSS, the noisiest bands are generally the u and the z band). PCA will pick the direction of the noisiest band as shown by the blue ellipsoids in the upper middle and right panels of Fig. 1, and not the direction of coherent variations. The remedy is to use the estimated errors to scale the observations so that we have unit variance in every band. In the case of a non-variable star and near-independent errors, this implies that the point cloud appears as a sphere. For a variable star, the noise is added to a deterministically varying light curve, which causes the centre of the sphere to move in the M -dimensional space, and instead of a ball, we observe an elongated ellipsoid-like 20.4 20.8 21.2 u magnitude  shape. The first principal direction is the longest axis of this shape, providing the largest variability amplitude that can be obtained by a linear transformation. Assume that the distribution of the observations is Xim ∼ N (µm, σ 2 im ), that is, the star is constant with mean magnitude µm in band m, the standard deviation of the measurement Xim around the mean is σim, and the error distribution is normal. We have several options to obtain unit variance of the noise. whereμm andσ 2 m can be either the mean (non-robust scaling) or the median (robust scaling) of Xim and of σ 2 im , respectively. In the case of correlated errors that can be assumed normal, an analogous matrix transformation can be based on the covariance matrix. In this paper, we assume that the correlation between the errors is weak compared to the correlation between bands for an RR Lyrae or HADS star (Scranton et al. 2005).
(iii) Considering that the errors depend on true magnitudes, and this effect can be strong for some classes of variable stars or very faint objects, we tested also the variance stabilizing transformation (Everitt 2002). Suppose we have a random variable X from a distribution with mean µ and variance σ 2 for which σ 2 = g(µ) .
Define the function f as the solution to f (x) = [g(x)] −1/2 . Then the transformation Y = f (X) leads to a variable with unit variance: Var(Y ) = 1. In our case, we supposed a functional form σ 2 = g(x) = a exp(bx) between the variances σ 2 im and centred observed magnitudes xim with different coefficients for each band. This corresponds to a rough approximation with Poisson noise, which leads to an easily tractable closedform transformation. Fitting this function bandwise to observed magnitudes and squared standard errors of a sample of variable and non-variable stars resulted in the bandwise estimatesâm andbm and a transformation Yim = 2â −1/2 mb −1 m exp(−bmxim/2) for band m. Then the point cloud of Yim is centred by removing either the mean or the median from it.
The assumption of normality can be checked by quantile-quantile (Q-Q) plots (Chambers et al. 1983;Cleveland 1994). After scaling and centring, the ordered sample of transformed magnitudes of a star is plotted against the corresponding theoretical quantiles of the standard normal distribution. If the sample follows indeed the standard normal distribution, then the points should be aligned along a straight line with intersect 0 and slope 1. If the assumption of normality is true, but the mean and the variance are not 0 and 1, the points still fit on a straight line, but the intersect and the slope of the line change to the mean and the standard deviation, respectively. If even normality is not satisfied, the points deviate systematically from the straight line. Figure 2 shows the Q-Q plots of a constant star for all variants of scaling, with the expected standard normal line pictured in red, and the best-fitting normal distribution in blue. The sample is visibly not normal in any of the bands: the tails deviate from both lines, indicating heavier tails (the absolute values of the observations are larger than what is expected from a normal distribution). The blue and the red lines have only slightly different slopes in most of the plots, suggesting that the scaling did obtain approximately unit variance, though the assumption of normality is not valid.
This implies that the non-robust PCA, built on the assumption of normality, is not appropriate for our purposes. Another important issue that violates this distributional condition is the presence of outliers. We can see their effect in Fig. 1: a single outlier distorts completely the estimated normal distribution, outlined by the red ellipsoid. Robust versions of PCA are given in the statistical literature, among which we chose the minimum covariance determinant method (Rousseeuw 1985). This involves repeated random subsampling of the data, leaving out a fixed percent of the observations each time. Then the covariance matrix is computed using each subsample, and the one having minimal determinant is chosen as the robust estimate. The blue ellipsoids in Fig. 1 are the result of the application of this method. The percent of left-out data is fixed so that the effect of a single outlier is decreased, but that of a few consistently located outlying observations is preserved. This choice is motivated by two conflicting aims. One is to get rid off of the effects of true erroneous observations. The other is to preserve that of those that are true representatives of the light curve, though they are rare; for instance scarce measurements of the dips in detached eclipsing binaries. Indeed, the blue ellipsoids in Fig. 1 represent better the apparent distribution of the bulk of the data.
Robust PCA yields some directly useful quantities. The first are the coefficients v1 of the linear combination, which characterize the direction of the first principal component. It depends on the type of variability, and thus can be used in classification. The second are the scores of the observations on the first principal component, that is, the time series of the linear combinations z1 = Xv1. By convention, the linear coefficients v1 are fixed so that the variance of the noise after transformation remains one, so the improvement on the signal-to-noise ratio is M m=1 vm1Am, where Am is the signal amplitude in band m after scaling. The third is the variance d 2 1 of z1, which is related to the elongation of the point cloud. Its ratio to the total variance M m=1 d 2 m (called hereafter the variance proportion) gives indication about the presence of variability. Finally, the distances of the observations from the centre with respect to the robust variance-covariance matrix fit indicate outlyingness of the observations, and can be used to weight the observations in period search.

Period search
We used the generalised Lomb-Scargle method (Zechmeister & Kürster 2009) to obtain the periodogram of the first principal component, and tested various trimming and weighting options on simulated sinusoids with errors imitating SDSS error patterns.
The usual weighting with the inverse squared errors cannot be applied for the z1 time series, because the scaling and the construction of the principal components lead to approximate unit variance of the noise on z1. Instead, we can base a measure of outlyingness on the robust distances, and trim or weight the elements of the time series according to this measure. The squared distances, r 2 i , in an M -dimensional space under the hypothesis of approximate standard normality should follow a χ 2 M distribution. Similarly to the principle of Q-Q plots, we can determine what value is expected for the i-th largest distance, r 2 (i) , among a sample of size N according to the χ 2 M distribution: this will be the i/(N + 1) quantile, c M,i/(N +1) , of the χ 2 M distribution (N + 1 is taken instead of N to avoid exactly 0 or 1 values). Such a χ 2 5 Q-Q plot of the robust distances is shown in Fig. 3 for a variable star. We checked various alternatives of trimming thresholds and weight definitions based on the difference between the ith largest observed distance with simulations.
The goal of the simulations is to reproduce the SDSS error levels, error dependence on magnitude, outliers and observational cadences. In order to obtain this, we used the following procedure: 1. From 2000 randomly selected objects from the 68,000 Stripe 82 objects flagged as variable, we extracted the sets of five amplitudes based on the excess root mean square variability, the sets of five median magnitudes and the u-band observational times. 2. We generated 200 sine functions with random frequencies in the interval [0, 20] day −1 , with random phases, and randomly selected sets of five amplitudes and of median Std. normal quantiles Figure 2. Quantile-quantile plots of a non-variable star from Stripe 82. The columns are the five bands u,g,r,i,z. The rows correspond to different scaling methods: first row from top, local scaling; second row, non-robust scaling; third row, robust scaling; fourth row, variance stabilizing transformation. The red line corresponds to a standard normal distribution, the blue line to a robustly fitted best-fit normal distribution to the data. magnitudes from those extracted from our sample. Each sine function was multiplied by all five elements of one amplitude set, and one set of median magnitudes was added to it. We obtained thus 200 pure five-band light curves, each with coherent brightness variations in the bands. We sampled these sine waves with time cadences extracted from Stripe 82 in the previous step. 3. We added noise to the light curves. We divided the full magnitude range of the 2000 stars in each filter into 0.5 mag wide bins, and grouped the error bars according to the magnitude value with which it occurred. We obtained thus a sample of all error bars, in magnitude bins, for all bands. Then for each simulated magnitude values of the sinusoidal light curves in each band, we randomly selected an error bar value from the bin corresponding to the simulated magnitude value in that band. We generated a random number from the N (0, 2 ) distribution, and added this value as the realised noise to the simulated magnitude. Finally, we joined as the error bar to the simulated series. 4. To simulate outliers, we changed some observations to fainter values. The differences were random values from N (0, (6 ) 2 ), in a random number of filters simultaneously. The number of outliers followed a Poisson distribution with mean equal to 0.005.
Visual inspection and Q-Q plots showed that the observed and simulated light curves are very similar. The simulations were then used to investigate the performance of different combinations of scaling, weighting or trimming and period search. We found that the combination leading to the highest frequency recovery rate was the generalised Lomb- Figure 3. Robust distances from the center of the point cloud for the 5-dimensional observations of a star, on a square root scale for better visibility. The red line indicates the exact χ 2 5 distribution. Scargle method with weights defined by −1/2 , applied to either the z1 time series derived from variance stabilized observations or to observations scaled with robust estimators of the mean and standard deviation. For comparison, the generalised Lomb-Scargle method with the usual weights based on error bars was also applied to detect periodicity on the single g band, which generally had the best signal-to-noise ratio. The results of the best z1 and the single-band analysis are compared in Fig. 4. The logarithm of the difference between the found and the true frequency is shown versus the number of observations in the time series. Period search on z1 outperforms the singleband analysis. While the single-band analysis finds a yearly or daily alias in 31 cases, the analysis of z1 reduces this number to only 8.
Based on these simulation results, in the analysis of the variable stars in Stripe 82 we applied the generalised Lomb-Scargle method with weights defined by (1). For some variable stars, we also performed a multifrequency analysis. The folded z1 light curves were fitted with B-splines, restricting the smoothing parameter to provide a reasonable degree of smoothing (R procedure smooth.spline with the constraint 0.5 spar 1 (R Development Core Team 2010), and selecting its optimal value by leave-one-out cross-validation. The smoothed light curvesẑ1 were inspected visually, and the rare over-or undersmoothed cases were corrected by a manual selection of the smoothing parameter. The frequency analysis of the residuals was not performed on light curves where the equivalent degrees of freedom of the smoothing was high with respect to the number of data points, since this left too few residual degrees of freedom to obtain meaningful harmonic fits in the residual frequency spectrum. For the cases with enough residual degrees of freedom, the residuals r = z1 −ẑ1 were extracted from the B-spline model, and the same period search algorithm with the same weights as for z1 was performed on their time series. The significance of eventual peaks in the residual periodogram was assessed using a combination of non-parametric bootstrap (Davison & Hinkley 2009) and extreme-value methods (Coles 2001).

Random Forest
Supervised classification methods estimate the class (usually a discrete variable like 1, 2, . . . , L or RRAB, RRC, EA, EB, ...) for an object of unknown class based on some attributes (for example period, amplitude, colour, etc.), by fitting a model to a set of objects with measured attributes and known classes (the training set), and then using this for prediction. Random Forest (Breiman 2001;Hastie, Tibshirani & Friedman 2009) is a popular method which works excellently in a very wide range of data mining problems. It consists of building a collection ('forest') of classification trees on the training set, then passing any new instance down on all trees and obtaining the class estimate by the majority vote of the forest.
The training ('growing of the forest') begins with a selection of a large number of bootstrap samples from the training set. Then a classification tree is built separately for each sample with binary splits. First a relatively small subset of attributes is selected randomly. For each of these, the split-point is computed that separates the given bootstrap sample with the least mixing of classes by some criterion (for example, one of the subsets contains only classes A and B, whereas the other mostly C and D). Then the tree is split according to the attribute for which this separation is the cleanest. For the next step, both subsets (called nodes) are further split in the same way as was done in the first step: from a small random subset of attributes, the one is selected that splits best the node in question. The procedure is repeated until each final node is either perfectly clean or has a pre-defined minimal size. Then, using another bootstrap sample from the training set, a new tree is built. As a result, a large forest of many trees emerges. To classify a new object, the prediction of the class by each tree in the forest is calculated, and the class with majority vote is accepted as the estimate. An estimate of the probabilities to belong to each of the classes can be obtained by calculating the proportion of votes for each class.
The main advantage of this procedure is that it obtains the class estimate as an average of many individual estimates by the trees. These estimates are unbiased, but have high variance. Averaging preserves unbiasedness and reduces variance, and this variance reduction is larger, if the trees are less correlated (Hastie, Tibshirani & Friedman 2009). Random Forest achieves de-correlation by applying two tricks: first, it uses only bootstrap samples from the training sets, so the basic data set driving the learning process is not identical for each tree; and second, it uses the best of only a random subset of attributes, not of all attributes. The consequent higher variance of the individual trees is more than compensated by the decrease of the variance of the average because of the nearly non-correlated trees.
There are a few tuning parameters in the Random Forest procedure: the number of trees in the forest, the number of random attributes at the nodes, and the final node size. For most classification problems, growing several hundred trees is enough, and adding more trees does not improve the prediction accuracy. With respect to the final node size, optimal results are achieved by maximally growing the trees (that is, the minimal final node size is 1), but it is customary in larger problems with many instances and attributes to grow the trees only to a somewhat larger final node size to speed up the process without great loss of accuracy. The most influential tuning parameter is in general the number of attributes from which the best split is chosen at the nodes, though Random Forest is only weakly sensitive to this as well. Breiman (2001) proposes to choose the largest integer k such that k < √ K, where K is the number of attributes.

Principal Component Analysis
We considered only stars that have at least 10 complete ugriz observations from our data set. As discussed in Section 3.1, scaling and centring is necessary before applying robust PCA. We tested all four methods (local scaling, robust and non-robust average scaling and variance stabilization) outlined there for variability detection, characterization and period search. The requirements of these three tasks are not the same, and thus different scaling methods perform better in each one, but inspection of the results suggested that variance stabilization yields the best overall performance. We present here only the results based on this transformation.

Variability detection based on z1 variance
After scaling and centring the data, we expect to see a spherical point cloud, similar to a multivariate standard normal sample, if the errors are nearly non-correlated and the star is constant. Detection of variability is thus equivalent to detecting excess variance in the point cloud of the scaled observations as compared to the variance of the first principal component of a 5-dimensional standard normal point cloud. Samples from even a spherically symmetric distribution show stochastic distortions from the perfect sphere, and the smaller the number N of the points in the cloud is, the stronger this distortion is. As PCA selects the direction of maximal spread in the data, we cannot expect d 2 1 to be exactly 1. The distribution of d 2 1 as a function of N is the easiest to obtain with simulations. The zero hypothesis of non-variability of a source may be tested at a given confidence level α by comparing the observed d 2 1 of the source to the quantile c1−α of the simulated distribution: d 2 1 > c1−α rejects the hypothesis of non-variability of the source at the level α.
As a PCA-based selection criterion, we use the proportion of d 2 1 to the total variance. Figure 5 shows its comparison to the variable selection of Sesar et al. (2010). Five hundred objects that were flagged variable (black dots) and 500 flagged non-variable (filled grey triangles) are plotted, together with the RR Lyrae sample of Sesar et al. (2010) (empty blue circles). The 0.9999 quantile derived from the simulations is denoted by the red line. For sources that are above this line, non-variability is rejected at the level of 0.0001. The left panel shows the proportion of d 2 1 to the total variance. Most of the objects flagged variable (black dots) are above the red line, meaning that the variance proportion criterion and the classical cuts select approximately the same variable sample. Some difference nevertheless can be observed. Among the stars flagged non-variable, the PCA-based criterion found some variables: the grey triangles above the red line. The significant variance excess seems to be due to two reasons. When the bands are weakly correlated, the excess might be the consequence of either microvariability, or correlated and underestimated errors across bands (Scranton et al. 2005), so to detect microvariability in Stripe 82, a de-correlated version of scaling should be used. The other is that the variance stabilization did not succeed to obtain unit variance in one band. In this case, the point cloud is elongated along one axis. Correspondingly, the coefficients v1 of the first principal component contain usually one value close to 1 at the band which is the origin of the excess variance and a value close to 0 for all others. The characteristic profile of v1 helps to recognize these cases.
Conversely, the black dots below the red line in the left panel of Fig. 5 represent points that are flagged variable but found non-variable by PCA. These seem to belong to two groups: one for which the flag seems to originate from an outlier, and another that appears ball-like on the pairwise scatterplots like Fig. 1, but not with unit variance. The variable flag from the classical analysis in the latter case can be due to under-estimated but nearly non-correlated errors in all bands Scranton et al. (2005). If the under-estimation is of similar order in every band, then the PCA-based criterion is less biased than traditional cuts, as it is defined as a proportion. Moreover, v1 would again contain one value close to 1 and all the others near-zero, so it is possible to filter out these cases based on the v1 profile. In general, this form of v1 can be used to recognize and filter the cases where the excess scatter is due to underestimated errors in one or a few bands rather than true variability of the source.

Period search on PC1 time series
For the majority of the investigated variable objects, period search on the SDSS time series suffers from aliasing problems. The daily and yearly observational patterns (see e.g. Sesar et al. 2007Sesar et al. , 2010 create aliases with complex structures, displaying combinations of the variability frequency with the 1/day and 1/year frequencies with comparable amplitudes. According to the simulation results of Section 3.2, using the time series of z1 with weights defined by (1) improves on period search on g-band. We applied this procedure, complemented by a non-linear optimization to find the exact value of the frequency, to a small random sample of 2000 stars flagged variable, which contained 36 known RR Lyraes. The published periods of these RR Lyraes, as described in Sesar et al. (2010), were determined by a visual inspection of light curves folded with the five best periods returned by the SuperSmoother algorithm (Reimann 1994) restricted to the [0.2, 1] day interval. For 30 out of the 36 RR Lyrae, the result from our algorithm coincided with the true frequency. In the 6 other cases, daily alias periods were found. This suggests that we can expect good period search results from a simple automatic procedure even in the absence of visual inspection.
We applied the weighted period search method (without the nonlinear optimization) on the full variable star set. The dominant frequency was used as an attribute in the classification procedures, presented in Section 4.2. The classification selected a sample of 317 candidate RR Lyrae and HADS stars. For these, we modelled the folded z1 light curve, computed the residuals r = z1 −ẑ1, and performed another period search on r as described in Section 3.2. The found likely multi-period objects are discussed in Sections 4.2.2 and 4.2.3.

Characterisation of variability with PC1 spectrum
RR Lyrae variables show a characteristic wavelengthdependent amplitude pattern: their variability is stronger at shorter optical wavelengths than at the red end of the spectrum (Smith 1995). If we could observe RR Lyrae stars with SDSS filters with equal errors in all bands, we would see the first principal direction tilted towards the blue rather than the red bands, and correspondingly, larger coefficients v1 for blue than for red wavelengths. However, the error patterns distort this simple picture. We must scale the bands in order to get rid off of the effect of different average error levels, and the downscaling will be stronger in bands with larger errors, most notably in u and z. As a result, variability amplitudes get downscaled as well. Thus, the typical RR Lyrae v1 pattern (the PC1 spectrum) takes a characteristic shape. It is composed of, on the one hand, the amplitude pattern of pulsating variables governed by the physical parameters and pulsation mode, and on the other, of the scaling patterns of the survey. Typical PC1 spectra of identified RRab and RRc stars are shown with black lines in the middle panel of Fig. 6. They exhibit small coefficients vu1 and vz1 for the u-and z-band contributions, and the highest value is vg1 at the g-band, which has high variability amplitude and low errors.
HADS variables are in some aspects similar to the RR Lyraes (see e.g. Breger (1980); McNamara (1995); Petersen & Christensen-Dalgaard (1996); McNamara (1997); a light curve is shown in the bottom left panel of Fig. 6). Their effective temperature is roughly in the same range, and although their pulsation may be more complex than that of RR Lyraes, they are mainly radial-mode pulsators. They show similar patterns across the wavelengths with decreasing amplitudes from blue to red wavelengths (Pigulski et al. 2006), so we can expect their PC1 spectrum to be similar to that of the RR Lyraes. The middle panel of Fig. 6 presents the PC1 spectrum of several variables that were found to have periods, colours and folded light curves characteristic to HADS variables. These are plotted in blue, superposed on the lines of the RR Lyrae stars.
For some other types of variability, we can expect different PC1 spectra. Eclipsing binaries that are composed of two components of the same age and similar masses have similar colours, and therefore show only very weak colour changes and an almost-equal contribution of all bands to the light variation. This results in a flatter PC1 spectrum. Combined with the specific error pattern of SDSS, this yields a profile that is low at the noisy u and z bands and have a higher plateau at g, r and i. Several PC1 spectra of this type are also shown in the middle panel of Fig. 6; these objects have clear eclipsing binary-type folded light curves. Figure 6 shows the use of the PC1 spectrum for discriminating certain types. In the two panels on the righthand side, we show the light curve of an RRc star from the confirmed sample of Sesar et al. (2010) (upper right panel) and an eclipsing binary candidate with only slightly different minima and near-sinusoidal light variation (bottom right). This eclipsing binary, when folded by half of the period, shows a folded light curve very similar to that of an RRc variable, as the tiny difference in the depths of the minima is masked by the error bars. If such eclipsing binaries fall furthermore close to the colour-colour region of RRc stars, the automated classification is difficult. However, the PC1 spectrum is approximately flat-topped for a geometric origin, whereas it has a peak at g for a pulsating RRc. The    inclusion of the first principal direction v1 can thus be helpful in an automated classification procedure. For the visual selection of HADS training sample, we have expressed all possible selection criteria in the framework of PCA. The extraction of variable objects was performed with the aid of the variance proportion of the first principal component. To reflect the astrophysical properties of the sought stars, we require them also to show the characteristic colour changes during the cycle. In the language of PCA, this corresponds to have broadly the PC1 spectrum demonstrated by the blue lines in the middle panel of Fig.  6. Furthermore, they must have admissible u − g and g − r colours. Also, the dominant frequency Fz 1 of the time series of the first principal component z1 must correspond to the HADS frequency range. In summary, the selection criteria are as follows.

Preliminary selection
(i) PC1 has a variance proportion higher than the 0.9999-quantile of that of a 5-variate standard normal point cloud with the same number of observations. (ii) The PC1 spectrum v1 = (vu1, vg1, vr1, vi1, vz1) T satisfies the following cuts: 0 < vu1 < 0.8, 0.45 < vg1 1, 0.3 < vr1 < 0.75, 0.15 < vi1 < 0.6, 0 < vz1 < 0.5. (iii) The (extinction-corrected) median colours are in the region 0.7 < u − g < 1.4 and −0.2 < g − r < 0.4 of the u − g, g − r diagram. (iv) Fz 1 > 3/day. All cuts are given very broadly, so the resulting selection contains a majority of contaminating objects such as eclipsing binaries with flat PC1 spectra, quasars, mainsequence variables and even RRc stars. Visual inspection of the g-band light curves resulted in 117 plausible HADS candidates. The training set for the first cleaning iteration consisted of these as class HADS; the RRab and RRc variables of Sesar et al. (2010) as classes RRAB (379 objects) and RRC (104 objects), respectively; and a random selection of 2400 other (O) variables from the remaining variable set, yielding a training set of 3000 objects in total. For all objects we calculated the following attribute list: • the median (apparent) magnitudes in all bands, corrected for interstellar extinction using the dust map of Schlegel Figure 7. Attribute ranking by mean accuracy decrease when classifying without the attribute in question. The components v 1g and v 1i of the PC1 spectrum are highly ranked, comparable to the frequency, r − i and g − r colours.
Selection of a small but sufficient attribute set is nevertheless necessary: the performance of most data mining procedures is sensitive to the presence of noise-like attributes. The importance of any attribute can be measured by Random Forest by for example calculating the average accuracy loss on trees generated without the use of the attribute in question. The importance plot for all attributes in Fig. 7 indicates that the best two attributes are the frequency and the g-band component vg1 of the first principal direction. They are followed by the r − i and g − r median colours, then another two PCA-based attributes, vi1 and the IQR of the z1 time series. This latter proves to be more important than the g-band IQR. The relatively low accuracy loss due to omission of colours is a consequence of our preselection of colour range in the u − g, g − r plane and of the correlation between colours. According to a procedure similar to that applied in Sections 4.2-4.4 of Dubath et al. (2011), the most performant attribute set is the 11 top-ranking attributes: Fz 1 , vg1, r − i, g − r, vi1, IQR z 1 , u − g, Sg, i − z, Sz 1 and IQR g . This is what we used for the following cleaning of the training set and the extraction of our set of candidate RR Lyrae and HADS stars.
The next stage is the cleaning of the training set, to obtain clearly separated HADS and O classes. In the first iterative step of this, Random Forest was trained on the 3000 objects using the attribute list presented above. At this stage, the set O may contain unrecognized RR Lyrae and HADS stars. We assume that this contamination is low enough not to bias strongly the automated classification results for class O in the first run, as these are relatively rare compared to all other types taken together. An overall error rate less than 1% in the first iteration confirmed this assumption. After each run of Random Forest, we removed the objects from type O and HADS that were misclassified (keeping all RR Lyrae regardless of predicted type, since these are confirmed variables), and repeated classification. After 3 iterations, we ended up with 108 HADS and 2372 O stars, which were separated perfectly from each other and from the RR Lyrae classes. The only confusion arose from 10 RRAB and RRC stars that were classified as type O, and one RRAB classified as RRC. This procedure biases the subsequent extraction of HADS, RRAB and RRC towards purity against completeness, as we dropped all cases that might represent unusual specimens of classes, and the final selection does not fully take into account eventual real similarities between stars belonging to distinct types. The choice between purity and completeness is certainly subjective, and depends on the purpose of the study.
Re-processing the whole variable collection by Random Forest using the resulting training set gave 317 candidates, 163 HADS, 82 RRC and 72 RRAB stars. Multifrequency analysis as described in Section 3.2 was performed on them. The final visual check was done with the help of 'portraits' of the stars: summary plots ( Figures A1-A7) that contained the most important information that could have been extracted from the data, namely, z1 and residual frequency spectra, folded z1, g, g − i and residual light curves, the raw observed time series, colour-colour, colour-period and period-amplitude diagrams, and PC1 spectrum. These plots are presented in the Appendix. The visually selected 132 HADS, 68 RRab, 36 RRc and 25 multiperiodic or peculiar RR candidates are discussed in the next two sections.  Fig. 9; V : mean dereddened V -band magnitude, determined from a best-fit V -band template synthetized from the best-fitting g and g templates, A u , ...: amplitudes; φ u 0 : epoch of maximum brightness; u 0 : the magnitude at the epoch of maximum brightness corrected for the ISM extinction (the last three determined from the best-fit template); Tu, ...: best-fit template identifier number as given in Sesar et al. (2010). The asterisks indicate slightly underestimated values. d and V were computed only for those RRab variables that are plotted on Fig. 9, and are therefore missing for some stars. The complete table is available online.  Figure 8. g-band template fit for star 3874813, using the templates of Sesar et al. (2010), an example that cannot be fitted well by these. For comparison, see the second panel from top in the middle column of Fig. A1, which shows the same star fitted with a B-spline.

Candidate RR Lyrae variables
The list of the 68 RRab and 36 RRc type objects is given in Table 1, containing the identifier, coordinates, median mag- nitudes corrected for interstellar extinction, the frequency and the bandwise amplitude estimates derived from a Bsplines fit to the folded light curve. A summary plot for one of them is presented in Fig. A1   were several RRab stars that did not have acceptable template fits for all or part of the bands, as the template set is not all-encompassing (Sesar et al. 2010). Such light curves showed a sharper rising branch, a narrower peak and a decreasing branch composed of a steep initial fading from the brightest state and a more prolonged near-minimum state than the templates. One example, the same star that is presented in Fig. A1, is shown in Fig. 8. The template amplitudes of these stars are underestimated on average by around 10% (up to 20%); such cases are denoted in Table 2 by an asterisk following the amplitude value. For the RRab stars, we estimated heliocentric distance following Section 4.1.1 of Sesar et al. (2010). We assumed the same average halo metallicity of [Fe/H] = -1.5 and a mean value MV = 0.6 for their absolute magnitude. The extended heliocentric map of the halo is presented in Fig. 9. The majority of the new variables is situated close to the Galactic bulge, continuing smoothly the already known structures towards the bulge. The multifrequency analysis, performed on all candidates with sufficient residual degrees of freedom, yielded 23 further objects with multiple modes whose significance was confirmed by a combination of bootstrap and extreme-value procedures (though two with only very few residual degrees of freedom after smoothing), and two more showing clear light curve changes without indication of a second frequency. These are listed in Tables 3, 4 and 5.  Table 3 contains the found double-mode RR Lyrae candidates showing the theoretically expected fundamentalovertone frequency ratio close to 0.745 (12 stars, around 2% of the total Stripe 82 RR Lyrae sample found in Sesar et al. 2010 and here). The summary information plot for one of them is presented in Fig. A2. These stars fit very well the expected relation between the fundamental period and the fundamental-overtone period ratio, as shown in Fig. 10. The span of periods indicates a broad range of metallicities and/or masses to account adequately for the observed ranges of periods and ratios (see e.g. Alcock et al. 2000a). In addition, Table 3 lists two more objects for which the period search gave probably aliased results, and therefore the ratio is off for the frequencies corresponding to the maximum of the periodogram. However, supposing another alias to be the true frequency (Fr − 1 day −1 for the fundamental mode of object 1149344 and Fz 1 + 1 day −1 for the overtone in the case of 1283137), the stars are in the regions admissible for double-mode RR Lyrae variables on the colour-colour, colour-period and period-amplitude diagrams. With these frequency choices, they fit well also the relationship between the fundamental period and the period ratio, as shown by the empty circles in Fig. 10. Table 4 contains RR Lyrae candidates with two closely separated frequencies and/or showing multiple thread structure in their folded light curves, together with two others showing multiple threads without the presence of a resolved second frequency (an example is given in Fig. A3). Such properties may be due to the Blazhko effect, slow periodic or irregular modulations of the amplitude or of the period. The folded light curve panels in Fig. A3 indicate the Julian date of the observation with colours (the colour legend is given under the folded light curves), and the different threads visibly belong to different years of observations. For this object, the separation of the two frequencies is very small, suggesting possibly a long Blazhko period, the presence of period changes or phase shifts. We list two more objects (6058913 and 6086465) with very small number of data where the smoothing left very few residual degrees of freedom, but the bootstrap indicated presence of secondary modes. More ob- servations of these objects are needed before their type can be resolved. Several other RR Lyrae candidates showed two significant peaks at neither closely spaced nor fundamental-first overtone frequencies. Objects with similar frequency ratios were recently found by Soszyński et al. (2010) in OGLE data in the LMC and by Olech & Moskalik (2009) in ω Centauri. We present three of them in Table 5 for which the Monte Carlo assessment showed these secondary peaks significant. The first one, with its period ratio around 0.8, may be a candidate double-mode RR Lyrae pulsating in the first and second overtone (its summary information plots are given in the appendix, Fig. A4); the other two show period ratios around 0.6. These objects need further data for confirmation.

HADS candidates
Using the training set described in Section 4.2.1, the Random Forest procedure selected 163 candidates. This was reduced by a subsequent visual inspection to 132 accepted HADS candidates, 109 monoperiodic and 23 multiperiodic stars. These variables often have frequency spectra with complex structure, usually several secondary peaks beside a highly dominant frequency. The significance of the secondary peaks was much less prononunced in the multiperiodic HADS sample than among the double-mode RR Lyraes. Table 6 presents a summary of the basic properties of the monoperiodic sample, similarly to Table 1. Two examples, one with symmetric and one with asymmetric folded light curve are presented in Figs. A5 and A6. All frequencies given in Table 6 can be affected by aliasing, though this does not modify the classification of the stars, as the typical frequencies of the HADS are high, and the aliasing would hardly lead to RR Lyrae frequencies and a consequent misclassification. In the absence of metallicity estimates or spectroscopic information, their further division into Population I (δ Scuti) and Population II (SX Phoenicis) objects is not possible.
In this sample, we also found several suspected multiperiodic objects. We included a star into our multiperiodic set only if two conditions held: (i) the combined bootstrap and extreme-value methods furnished a P -value smaller than 0.05 for the peak in the residual spectrum, (ii) the visual check of the folded residual light curve showed perceptible systematic variation rather than several outliers grouped by the folding. The properties of these stars are summarized in Table 7. We divided the table into two parts, the top half containing 12 objects with residual peak more significant than 0.01, and the bottom half another 11 with secondary peaks with P -value between 0.01 and 0.05. Among our multiperiodic candidates, we observed a broad variety of ratios between the primary and the residual periods at various levels of significance: they range from around 0.65 to 0.86, in addition to a group at very high ratios (> 0.92), and one peculiar object with a ratio 0.19. The scatter of these values suggest a wide variety of objects of diverging types and with large differences in their characteristics such as the evolutionary state, mass, metal content or rotation. The majority of the most significant secondary peaks seem to form close doublets with the main frequency. Such behaviour is quite frequent among low-amplitude δ Scuti and SX Phoenicis stars (see e.g. Poretti 2003;Breger & Bishof 2002;Mazur, Krzeminski & Thompson 2003). However, aliasing makes it difficult to find the true frequency in both principal component and residual spectra, and the presence of a strong secondary peak may be also the consequence of a not sufficiently precise primary frequency identification due to the relatively scarce number of observations per time series. A few stars of our sample exhibits signs of amplitude or phase/period changes, light curve threads separated by observational year and showing a high dispersion of the peak values (the clearest example is 225151, shown in Fig. A7).

DISCUSSION
Multi-filter observations comprise much information about the characteristics of the observed star through the colour variations. We tested robust PCA as a way to extract information about these variations. We found that PCA produces several useful quantities: the proportion of the variance of the first principal component to the total variance, the time series of the first principal component, a robust notion of outlyingness, and the direction of the first principal component. We found that the proportion of PC1 variance to the total variance is a useful variability indicator. The PC1 spectrum can help detect cases when excess variance is due to underestimated errors, by having one dominant and four near-zero elements. The time series of the PC1 yields the best attainable signal-to-noise ratio. Also, weights for period search can be defined based on a robust notion of outlyingness. Simulations and a test on real data containing several known RR Lyrae stars confirmed that period search on the PC1 time series with robust weights decreases the impact of aliasing on the results.
The direction of the largest variation in the 5dimensional point cloud, which is termed the PC1 spectrum, proved to be useful in classification. For our sought sample of radially pulsating variables between 6500-7500 K, combined with bandwise SDSS error patterns, this direction points towards the g band, with decreasing contribution from r, i, u and z bands. It was used in the Random Forest classification procedure as a new attribute helping to separate pulsating variables from eclipsing binaries. The coefficients of the g and i bands ranked among the best five attributes, beside the period and the r − i and g − r colours.
The study produced 132 HADS, 68 RRab, 36 RRc and 25 multiperiodic or peculiar RR Lyrae variable stars. The RR Lyrae are new addition to the confirmed RR Lyrae list of Sesar et al. (2010). The vast majority of these new RR Lyrae stars were found in the regions not originally considered by Sesar et al. (2010). Only a few type-ab RR Lyrae stars were missed by their work in the range considered by that study (308 deg < RA < 60 deg and |Dec| < 1.25), making their sample about 99% complete. Among the new multiperiodic RR Lyrae, there are 14 double-mode candidates, several others showing signs of Blazhko effect, period or phase change, and we found candidates with unusual period ratios, of which one may be a double-mode star pulsating in the first and second overtone. Among the HADS candidates, the multiperiodicity seems to show a broad variety of ratios between periods. The HADS candidates, similarly to the multi-mode RR Lyrae, need eventual observational follow-up to clarify and confirm their type and pulsation modes.
This work yielded promising results about the utility of PCA, and opens the way to further improvements. Missing data reduce the number of observations in the time series of the first principal component to the number of data in the most scarcely observed band. Thus, though signal-to-noise ratio improves with the application of PCA, the decrease in the number of observations in the time series can imply worse performance of period search. Optimization of the period search performance can be achieved either by restricting the analysis to the combination of only the best bands, or by a statistical methodology that is able to deal with the missing values. Further investigations will explore these interesting possibilities.