Feature Selection for Classification of Blazars Based on Optical Photometric and Polarimetric Time-Series Data

We used the data obtained with the Kanata telescope which was published in Itoh et al. (2016). The data includes $V$-band time-series photometric and polarimetric data for 45 blazars from 2008–2014. Figure 1 shows examples of light curves and variations in the degree of polarization (PD).

Panels (a), (b), and (c) of figure 1 show examples of FSRQs: PKS 1510$-$089 and 3C 454.3 in 2010 and in 2008, respectively. The peak-to-peak amplitudes in the light curves are large, over 2 mag in all cases, while the light curve profiles are diverse: a solitary, short flare appears in panel (a), while a number of short flares superimposed on long outbursts appear in panels (b) and (c). The light curves change their apparent characteristics year by year even for the same object, as shown in panels (b) and (c). Panels (d), (e), and (f) of figure 1 show examples of BL Lac objects: BL Lac (LBL), 3C 66A (IBL), and Mrk 501 (HBL), respectively. The peak-to-peak amplitude of the light curve in panel (d) is comparable to that of FSRQs, while the light curve profile looks different. The characterization and classification of these variations are the main subjects of this paper. The variation in polarization could give rise to some interesting features. For example, PD flares are associated with FSRQ flares, though no clear correlation can be seen in the light curve and PD variations in panels (d) and (e). Panels (e) and (f) show that the variation amplitude apparently decreases from LBL to HBL, as mentioned in the previous section.

We use the OU process for our time-series analysis. The OU process is a stochastic model based on the multivariate normal distribution whose covariance between the data at time $t_{i}$ and $t_{j}$, $S_{ij}$, is given as:

$\displaystyle S_{ij}=A_{\rm exp}\exp\left(-\frac{|t_{i}-t_{j}|}{\tau}\right),$

(1)

where $\tau$ and $A_{\rm exp}$ represent the characteristic variation timescale and amplitude at $\tau$, respectively (Uhlenbeck & Ornstein (1930)). For the time-series data followed by the OU process, $f(t)$, the observed data, $m(t)$, is given by $m(t)=f(t)+\mathcal{N}(0,\sigma_{\rm OU}^{2})$, where the second term is the noise defined by the normal distribution having zero mean and variance $\sigma_{\rm OU}^{2}$. We can extract the characteristic features, $A_{\rm exp}$, $\tau$, and $\sigma_{\rm OU}^{2}$, from the observed time-series data using the OU process regression.

The time-series data introduced in § 2.1 have different observation periods for each object. The time-series data of each object was divided into one-year segments, each of which is regarded as a sample in this paper. For modeling the light curves with the OU process, the magnitude values were translated to fluxes on a logarithmic scale, simply dividing by $-2.5$. We assumes that the light curves are approximated with the OU process with a characteristic time-scale less than a few tens of days for our sample. The short time-scale is supported by the data in which erratic variations are detected over measurement errors in all samples. If our assumption is true, the power spectrum should be flat for frequencies ($f$) lower than the characteristic frequency, and decays as $f^{-2}$ for higher frequencies (Kelly et al. (2011)). A strong linear trend in the time-series data breaks this assumption because the power becomes larger in lower frequencies. We consider that the linear trend has an origin different from the short-term variations governed by the OU process. The presence of such distinct short- and long-term variations are reported in AGN (Arévalo et al. (2006); McHardy et al. (2007); Kelly et al. (2011)) and also in blazars (Sobolewska et al. (2014)). Hence, we first subtracted the linear trend from the samples, and then performed the OU process regression. The slope value of the linear trend can be considered an indicator of the power at the lowest frequencies, and we use it as a feature for the classification in the next section.

The OU process is identical to the Gaussian process with an exponential kernel. We used the python package for the Gaussian process, GPy, which includes a package for the Markov chain Monte Carlo (MCMC) method for the estimation of the posterior probability distributions of the parameters. In the present work, we estimated the posterior distribution of $A_{\rm exp}$ with a flat prior probability and that of $\tau$ with a positive flat prior. We fixed $\sigma_{\rm OU}^{2}$ with a typical measurement error of the data. We estimated the posterior probability distributions of $A_{\rm exp}$ and $\tau$ using the MCMC method for each light-curve sample. We set $\sigma_{\rm OU}^{2}=10^{-5}$. Figure 2 shows trace plots of $A_{\rm exp}$ and $\tau$, their posterior distributions, and the observed and model light curves for the sample S5 0716$+$714 between MJD 55050 and 55389. The MCMC samples converge to a stationary distribution and the posterior distributions have single-peaked profiles. In this case, we successfully obtained unique solutions of $A_{\rm exp}$ and $\tau$.

On the other hand, we found that $A_{\rm exp}$ and $\tau$ were not uniquely determined for several of the samples, mainly because the data size is not large enough to make a meaningful estimate of the parameters. A significant number of the samples from Itoh et al. (2016) have only $<30$ data points. Even in the samples with larger data size, $\tau$ is not uniquely determined if it is too long. Figure 3 shows an example, AO 0235$+$16 between MJD 54617 and 54946. The MCMC samples do not converge to a stationary distribution and $\tau$ can be very large, reaching over 300 d. Kozłowski (2017) reports that the OU process model is degenerate when the baseline of the light-curve sample is shorter than ten times $\tau$. The result in figure 3 is probably an example of such a case..

In this paper, we used only samples for which $A_{\rm exp}$ and $\tau$ were uniquely determined, as in figure 2. This selection reduces the number of samples to 38 for 18 objects. The selected samples include 12 samples for 8 FSRQs and 26 samples for 10 BL Lac objects. The samples are listed in table A in the Appendix. The designation of the objects to the subtypes FSRQ, LBL, IBL, and HBL is taken from Itoh et al. (2016).

Blazars occasionally exhibit large prominent flares, as shown in panel (a) of figure 1, which definitely arise from a non-stationary process, whereas the OU process model assumes a stationary stochastic process. In order to characterize the non-stationarity, we calculate the cross-validation error (CVE) using the OU process regression, as follows: First, the sample is divided into 25-d bins, a sub-sample of which is for validation while the others are for training. Then, the OU process regression is performed with the training subsets. Then, the log-likelihood is calculated from the validation subset and the optimized model. Using the other sub-samples as validation data, we obtained about $10$ log-likelihoods for each sample. The CVE is defined as their mean. A large CVE means that the validation data has a large deviation from the prediction of the model constructed from the training data. Hence, a large value of CVE indicates a high degree of non-stationarity.

The analysis of the time-series PD data was performed in the same manner as for the light curves, that is, dividing it into one-year-segments, converting to a logarithmic scale whose linear trends were subtracted. The OU process parameter $\sigma_{\rm OU}^{2}$ was fixed to $10^{-4}$. $A_{\rm exp}$ and $\tau$ were uniquely determined for the PD variations in all the samples except for seven. An example of one of the seven samples is shown in figure 4. Although, $A_{\rm exp}$ is uniquely determined, the MCMC samples of $\tau$ do not converge, and $\tau$ can be quite small (note that the scale of $\tau$ is logarithmic in figure 4). As a result, the model of PD is simply the mean of the data, as shown in panel (e). These results suggest that the characteristic timescale is too short to be properly determined with our data. We set $\tau=0.0$ for the seven samples. While a value of zero for $\tau$ is physically undefined, it works for training and evaluating the classifier representing very short timescales. CVEs were also calculated for the PD data. In addition, we also used the median of the PD as a feature parameter.

In total, we obtained nine features from the light curve and PD data: the four features from the light curve samples, $A_{\rm exp}$, $\tau$, the slope of the linear trend, and CVE, and five features from the PD samples, $A_{\rm exp}$, $\tau$, linear slope, CVE, and PD median. The values of the features are listed in table A in the Appendix.

We construct a classifier for FSRQs and BL Lac objects based on the nine features described in § 2.2. We use sparse multinomial logistic regression (SMLR) to determine the classifier (Krishnapuram et al. (2005)). We consider the problem of defining an $M$-class classifier with $N$ labeled samples, each of which has a $K$-dimensional feature vector, $\bm{\theta}_{i}=\{\theta_{i,1},\theta_{i,2},\cdots,\theta_{i,K}\}$ $(i=1,2,\cdots,N)$. A sample that belongs to the $j$-th class can be expressed with a vector $\bm{y}=\{y^{(1)},y^{(2)},\cdots,y^{(M)}\}$ such that $y^{(j)}=1$ and the other elements are $0$. Multinomial logistic regression gives the probability that a sample belongs to the $j$-th class, as follows:

$\displaystyle P\left(y^{(j)}=1|\bm{\theta},\bm{w}\right)=\frac{\exp\left({\bm{w}^{(j)}}^{T}\bm{\theta}\right)}{\sum_{j=1}^{M}\exp\left({\bm{w}^{(j)}}^{T}\bm{\theta}\right)},$

(2)

where $\bm{w}^{(j)}$ is the weight vector for the $j$-th class. The log-likelihood function is given by the data $\bm{\theta}$ as

$\displaystyle\ell(\bm{w})=\sum_{i=1}^{N}\log P\left(\bm{y}_{i}|\bm{\theta}_{i},\bm{w}\right).$

(3)

Then, the solution of SMLR is expressed as

$\displaystyle\hat{\bm{w}}=\mathop{\rm arg~{}max}_{\bm{w}}\limits\{\ell(\bm{w})-\lambda\|\bm{w}\|_{1}\},$

(4)

where $\|\bm{w}\|_{1}$ is the $\ell_{1}$ norm, $\|\bm{w}\|_{1}=\sum_{i}|w_{i}|$, and $\lambda$ is a sparsity parameter that controls the complexity of the model.

SMLR gives a linear classifier against the observed features if it is used as $\bm{\theta}$. In this case, SMLR can select the important features because the $\ell_{1}$ term makes $\bm{w}$ sparse. On the other hand, a non-linear classifier can be obtained if the observed features are transformed with non-linear kernel functions. Then, we can avoid over-fitting due to the $\ell_{1}$ term. In the present study, the features listed in table A were normalized and the feature vector of the $i$-th sample, $\bm{x}_{i}$, was obtained. The $j$-th element of $\bm{\theta}_{i}$ was obtained from $\bm{x}_{i}$ and $\bm{x}_{j}$ with the RBF kernel as follows:

$\displaystyle\theta_{i,j}=\exp\left\{-\frac{|\bm{x}_{j}-\bm{x}_{i}|^{2}}{2\sigma_{\rm RBF}^{2}}\right\},$

(5)

where $\sigma_{\rm RBF}^{2}$ is the bandwidth.

As mentioned in § 2.2, the number of samples, $N$, is 38. The number of classes, $M$, is two: FSRQs and BL Lac objects. Because of the small sample size, the three subtypes LBL, IBL, and HBL are combined as one BL Lac type, while the characteristics of the subtypes are discussed in § 4.2. The classifier is evaluated from the so-called “area under the curve” (AUC), which is defined by the receiver operating characteristic (ROC) curve. The simple accuracy of the classifier is inadequate because the number of BL Lac objects is larger than that of FSRQs in our sample (12 FSRQ samples and 26 BL Lac samples). The AUC is calculated by leave-one-out cross-validation for estimating the generalization error of the classifier. Optimization of the model and the calculation of the cross-validated AUC were performed with the Java-based application SMLR¹¹1http://www.cs.duke.edu/~amink/software/smlr.

A good classifier could incidentally be obtained in high dimensional problems even if all the features are not related to the real characteristics of the samples. We tested the significance of the obtained classifier described in the previous section using artificial data sets. The sets of artificial data consist of 12 FSRQs and 26 BL Lac objects, as for the case in table A, with random values for the nine features. The random numbers were uniformly distributed between $0$ and $1$. We made 100 sets of data and obtained $511\times 100=51100$ AUC values, in the same manner as described in the previous section. Figure 7 shows histograms of the AUC values from the real samples (red) and from the artificial samples (black). The distribution from the real samples exhibited systematically higher AUCs than that of the random data sets. The AUC values obtained from the random data sets are concentrated in the area AUC $<0.8$. Thus, it is unlikely that the obtained best AUC values from the real data ($\sim 0.9$) are incidentally obtained.

Here, we discuss the implications of the results of § 3. Figure 8 shows scatter plots of the four features. The correctly classified FSRQs ($P_{\rm BL}<0.5$) are indicated by the filled blue circles, while the misclassified FSRQs are indicated by the open circles. As can be seen from the top panel, the correctly classified FSRQs have high CVE and/or high $A_{\rm exp}$, while the misclassified FSRQs have values for these features comparable to those of the BL Lac objects. The high value of CVE indicates the presence of prominent non-stationary flares which deviate from the stationary OU process. The high value of $A_{\rm exp}$ indicates a large amplitude of variation at the characteristic timescale, $\tau$. We propose that FSRQs are characterized by rare and large flares which have a time-series structure distinct from ordinary variations. If the frequency of the flares is relatively high, a few times a year say, then the light curve can be reproduced by the OU process with a high $A_{\rm exp}$. If the frequency is low, such as once a year, then the light curve can be divided into two distinct periods, that is, the stationary state and the non-stationary flare, which causes a high CVE. The misclassified FSRQs may be objects in which the flare frequency was so low that no flare was detected in the year.

It is not evident that the characteristics of the light-curve CVE and $A_{\rm exp}$ originate from a different structure and/or physical condition for the jets between the blazar subtypes (e.g. Itoh et al. (2016)). It is possible that it is simply due to the $\nu_{\rm peak}$ effect. In order to investigate this point, we analyzed the X-ray data of the HBL Mrk 421 using the X-ray light curve presented in Yamada et al. (2020). The data was obtained with XRT/Swift from 2009 to 2014. The time-interval of the X-ray light curve is 1 d. We analyzed the data in the same manner as for the optical light curves, that is, dividing it into one-year-segments, calculating the flux density in a logarithmic scale whose linear trends were subtracted, and performing the OU process regression for each segment. Table A shows the estimated CVE, $A_{\rm exp}$, and $\tau$ for each sample. We successfully obtained values for the four segments listed in the table. We could not obtain those values for the segment MJD 55939–56078 mainly because of the small sample size ($N=34$). The estimated values are indicated by crosses in the top and middle panels of figure 8. They are definitely in the regime of FSRQs, especially regarding the large $A_{\rm exp}$. This result suggests that the large $A_{\rm exp}$ does not originate from different jet properties in the blazar subtypes, but from the $\nu_{\rm peak}$ effect.

The variation timescale of the light-curve, $\tau$, was also selected as an important feature for classification. However, as shown in the middle panel of figure 8, we cannot find any clear differences between the $\tau$ distributions of FSRQs and BL Lac objects. This feature was selected mainly because it is useful for the classification of only one FSRQ sample, 3C 454.3 in MJD 54542–54930. This sample has a small CVE ($=0.60$) and a not very high $A_{\rm exp}$ ($=3.15$), from which the object cannot be distinguished from BL Lac objects, but has an exceptionally large $\tau$ ($=21.30$). We consider that our analysis does not provide enough evidence to determine the importance of $\tau$.

It is proposed that the beaming factor of FSRQs is systematically larger than that of BL Lac objects. Giommi et al. (2012) proposed that the two subtypes have a common nature, except for the beaming factor. Itoh et al. (2016) propose that the jet volume fraction occupied by the fast “spine” should be larger in FSRQs than BL Lac objects. The difference in the beaming factor would also change the characteristics of the variability. For example, a shorter variation timescale is expected with a higher beaming factor. However, our analysis provides no strong evidence that $\tau$ of FSRQs is systematically smaller than that of BL Lac objects, while the uncertainty of $\tau$ is large.

In the bottom panel of figure 8, we can see a trend that BL Lac objects have high PD medians compared with FSRQs. This characteristic is stronger when HBLs (the black triangles in the figure) are neglected. The low PD medians of HBLs are probably due to a large contamination of the unpolarized emission from their host galaxies (Shaw et al. (2013)). Figure 9 shows examples of polarization variations in the Stokes $Q/I$–$U/I$ plane. The left and right panels show the data of the FSRQ PKS 1510$-$089 having a low PD median and that of the LBL OJ 287 having a high PD median, respectively. An increase in PD is occasionally associated with the flares of blazars, and in general the PD remains relatively low when the object is faint. In the left panel of figure 9, most of the data points exhibit low PDs, except for a few data points with high PDs over $>20$%, which are associated with the prominent flare shown in panel (a) of figure 1. The high value of the PD median in BL Lac objects indicates that PD is relatively high even in the faint state. The right panel of figure 9 shows an example: the object had a relatively high PD throughout the year. It has been proposed that the LBL object BL Lac has two polarization components: short-term variations superimposed on a stable or semi-stable component (Hagen-Thorn et al. (2002); Sakimoto et al. (2013)). The fact that the PD median was selected in our analysis suggests that the presence of a stable polarization component is a characteristic feature of BL Lac objects.

The origin of the PD median characteristic is unclear. The values of the PD median apparently correlate with $\nu_{\rm peak}$ in BL Lac objects, being lowest in HBLs and highest in LBLs. However, FSRQs have low PD medians, although they have the highest $\nu_{\rm peak}$. Hence, the characteristic is not due to the $\nu_{\rm peak}$ effect, but possibly due to a difference in the jet structure between the blazar subtypes. In this case, the stable polarization component should have a different emitting site or physical condition from the short-term variations since the characteristics of the short-term flux variability can be interpreted as the $\nu_{\rm peak}$ effect, as discussed above. The presence of the stable component may suggest that the accelerated electrons have a long life-time with a long cooling timescale. According to Kaspi et al. (2005), the size of the broad line region (BLR) has a positive correlation with the AGN luminosity, and the AGNs with the highest luminosity have large BLRs up to sub-pc scale. FSRQs form a sub-group of blazars with the highest luminosity, not only of the jets, but also of the AGNs (Fossati et al. (1998); Shaw et al. (2013); Ghisellini et al. (2017)). The lack of a stable component in FSRQs may be reconciled with the presence of a strong radiation field induced by a large BLR causing strong Compton cooling of the electrons even in the sub-pc region, which is the source of the stable polarization component in BL Lac objects.

Giommi et al. (2012) reported that LBLs include both low luminosity FR I objects and high luminosity FR II objects. If this is the case, there may be LBLs with a low PD median. Our samples included only two LBL objects, BL Lac and OJ 287. The number of samples is so small that we cannot make conclusions about the population of LBLs. Further studies are required to understand the relationship between the presence/absence of the stable polarization component and FR types or AGN luminosity.

In this paper, we used features derived from both the light curves and PD variations, while the only PD feature selected as being useful for classifying FSRQs and BL Lac objects was the PD median. Figure 10 shows a scatter plot of $\tau$ of the light curve and of the PD variations. In this figure, we can see that the timescale of the PD variation tends to be shorter than that of the light curve. Most of the objects have a PD $\tau$ shorter than 5 d. As mentioned in § 2.2, the PD $\tau$ was too short to be uniquely determined in seven cases. These results imply that the real $\tau$ could be too short to be correctly estimated from our data. If this is the case, the PD features were not selected in our analysis possibly because they were not good indicators for the nature of the PD variability. The fact that the PD $\tau$ is significantly shorter than the light-curve $\tau$ suggests that the relaxation timescale of the ordered magnetic field is shorter than the cooling timescale of the accelerated electrons.

On the other hand, the presence of a stable polarization component in BL Lac objects can also cause a lack of PD variation features in the selected features. The observed Stokes parameters are a sum of those of multiple components. If the contamination of the stable component is strong, the PD variation of a flare is diluted by the stable component. For example, the increase in PD associated with a flare is canceled if the direction of polarization of the flare component is perpendicular to that of the stable component. This effect also causes the PD features to be poor indicators of the real PD variability. In future work, we will extract the features of the PD variation for both the short-term flares and the stable, or long-term, variation component by separating these components (Uemura et al. (2010)).