\section{Introduction}

Blazars are a sub-class of active galactic nuclei (AGN) with
relativistic jets that point toward us. The jet emission is
amplified by the beaming effect and dominates the observed flux at
almost all wavelengths
(\cite{bla78bllac,bla79agn,urr95agn}). Synchrotron emission from the
jet is dominant in the radio--X-ray regime. The observed
X-ray--$\gamma$-ray emission is mostly due to inverse-Compton
scattering by relativistic electrons in the jet. Blazars exhibit 
violent variability, which provides a hint for understanding the
physical conditions and structure in AGN jets (e.g. \cite{ulr97var}).

Blazars consist of two subtypes: flat spectrum radio quasars (FSRQs)
and BL Lac type objects. The former was originally defined by strong
emission lines observed in the optical spectra
(equivalent width $> 5$~\AA; \cite{sti91bllac,sto91bllac}), while the
latter was defined by weaker lines or featureless spectra. In
addition, FSRQs have lower peak frequencies in the synchrotron
emission, $\nu_{\rm peak}\lesssim 10^{14}\;{\rm Hz}$, in their spectral
energy distribution (SED), while BL~Lac objects have a
wide range of $\nu_{\rm peak}$
($10^{14}\lesssim \nu_{\rm peak}\;[{\rm Hz}]\lesssim 10^{18}$)
(\cite{abd10sed}). The luminosity of blazars has a negative
correlation with $\nu_{\rm peak}$; FSRQs form the most luminous class
of blazars, while BL~Lac objects are less luminous. In SEDs, the
relative strength of the inverse-Compton scattering component to the
synchrotron component is larger in FSRQs than in BL~Lac objects. These 
regularities are known as the ``blazar sequence''
(\cite{ghi98seq}).

In addition to the classification based on the emission line strength,
a classification scheme based on $\nu_{\rm peak}$ is also used for
blazars, with low synchrotron peaked (LSP) blazars for objects with
$\nu_{\rm peak}\lesssim 10^{14}\;{\rm Hz}$, intermediate synchrotron
peaked (ISP) blazars with 
$10^{14}\lesssim\nu_{\rm peak}\;[{\rm Hz}]\lesssim 10^{15}\;{\rm Hz}$,
and high synchrotron peaked (HSP) blazars with
$\nu_{\rm peak}\gtrsim 10^{15}\;{\rm Hz}$ (\cite{abd10sed}).
Most FSRQs are LSP blazars. In this paper, we call LSP, ISP, and
HSP BL~Lac objects LBLs, IBLs, and HBLs, respectively.

The nature and links between blazar subtypes are still incompletely
understood. \citet{ghi08seq} proposed that FSRQs are AGN
having a radiatively efficient accretion disk (a ``standard'' disk; 
\cite{sha73disk}) with a high accretion rate, while BL~Lac objects
have a radiatively inefficient accretion flow (RIAF;
\cite{qua01riaf,nar95adaf}) with a low accretion rate. The accretion
rate is considered to be linked to the extended radio morphology of
radio galaxies, that is, the Fanaroff--Riley (FR) classification
(e.g., \cite{bau95fr}). It is proposed that FSRQs and all or some LBLs 
are beamed counterparts to FR type II radio galaxies with high
luminosity, and IBLs, HBLs, and possibly some LBLs are counterparts to FR
type I objects with low luminosity (\cite{mey11type,gio12type}).
\citet{gio12type} report that known LBLs are inhomogeneous and contain
both FR~I and II subtypes. 

The variability characteristics of the flux and polarization have also
been discussed for the different subtypes of blazars, particularly for
the optical waveband in which all subtypes have been frequently
monitored. It is well known that the optical activity apparently
depends on $\nu_{\rm peak}$; LSP blazars are more variable than HSP
blazars (e.g., \cite{bau09var,ike11blazar,hov14var}). A similar
$\nu_{\rm peak}$ dependence has also been reported in the polarization
variations, though the number of previous studies is limited
(\cite{ito16fermi,ang16pol}). The mechanism of the effect of the
$\nu_{\rm peak}$ on the observed flux and polarization variability is
unclear. High $\nu_{\rm peak}$ objects show less activity, possibly
because a large number or large area of emitting regions blur each
short flare (\cite{mar10multi,ito16fermi,ang16pol}), or possibly
because the jet volume fraction of a slower ``sheath'' component
increases (\cite{ito16fermi,ghi05spine}). 

In this paper, we focus on blazar variability. We have performed
photometric and polarimetric monitoring of blazars using the 1.5-m
Kanata telescope in Hiroshima since 2008 
(\cite{ike11blazar,ito16fermi}). The present study has two major
objectives: to establish the observational features of the flux and
polarization variability for characterizing the subtypes, and thereby
to investigate the nature of the subtypes, for example, whether FSRQs
and LBLs have a common origin and whether the jet structure of FSRQs
is different from that of BL Lac objects. Our analyses can be divided
into two parts, the extraction of features from the observed
time-series data and the selection of the features which are
discriminative for the two subtypes. For the feature extraction, 
in past studies the blazar variability was occasionally characterized
only by the features based on the variance of the whole data,
while the variation timescale was not considered. We use the
Ornstein--Uhlenbeck (OU) process to estimate both the timescale and
the amplitude from the data. The OU process and more advanced models
based on it have been used to characterize the variations observed in
AGN and also in blazars
(\cite{kel09agn,kel11mixou,rua12ou,sob14ou}). For the feature 
selection, we propose a data-driven approach to select the best set of 
features for classifying blazars by maximizing the
generalization error of a classifier. 

The structure of this paper is as follows: In \S~2, we describe
the data (\S~2.1) and methods used in this paper, namely, the OU
process for the feature extraction (\S~2.2) and sparse multinomial
logistic regression for the classifier (\S~2.3). In \S~3, we present
the results of the feature selection. In \S~4, we evaluate the
classifier and discuss the implications for the selected features. 

\section{Data and methods}

\subsection{Data}

\begin{figure*}
 \begin{center}
  \includegraphics[angle=90,width=17cm]{fig-samples.pdf} 
 \end{center}
 \caption{Examples of light curves and PD time-series data used in
 this paper. (a) PKS~1510$-$089 (FSRQ), (b) and (c) 3C~454.3 (FSRQ) in
 2010 and 2008, respectively, (d) BL~Lac (LBL) in 2013, (e) 3C~66A
 (IBL) in 2009, and (f) Mrk 501 (HBL). The upper and lower panels show
 the light curves and PD variations for each, respectively. The
 vertical and horizontal scales are common in all
 panels.}\label{fig:samples}
\end{figure*}

We used the data obtained with the Kanata telescope which was published in
\citet{ito16fermi}. The data includes $V$-band time-series photometric
and polarimetric data for 45 blazars from 2008--2014.
Figure~\ref{fig:samples} shows examples of light curves and variations
in the degree of polarization (PD).

Panels~(a), (b), and (c) of figure~\ref{fig:samples} 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~\ref{fig:samples} 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. 

\subsection{Feature extraction with the OU process}

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: 
\begin{eqnarray}
  S_{ij}=A_{\rm exp}\exp\left(-\frac{|t_i-t_j|}{\tau}\right),
\end{eqnarray}
where $\tau$ and $A_{\rm exp}$ represent the characteristic variation
timescale and amplitude at $\tau$, respectively (\cite{ou30}). 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 \S~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 (\cite{kel11mixou}). 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 (\cite{Arevalo2006, McHardy2007, kel11mixou}) and
also in blazars (\cite{sob14ou}). 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 {\tt python} package for the Gaussian
process, {\tt 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~\ref{fig:good} 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$. 

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{fig-Kanata-BL-12-55050-55389.pdf} 
 \end{center}
 \caption{Results of MCMC estimation of $A_{\rm exp}$ and $\tau$
   for the sample S5~0716$+$714 between MJD~55050 and
   55389. Panels~(a) and (c) are the trace plots of the MCMC samples
   of $A_{\rm exp}$ and $\tau$, respectively. Panels~(b) and (d) are
   their posterior probability distributions. Panel~(e) is the
   observed and model light curves. The filled circles are the
   data and the solid line and shaded region indicate the mean of the
   model prediction and its 95\% confidence interval,
   respectively. This is an example in which both $A_{\rm exp}$ and
   $\tau$ are uniquely determined.}\label{fig:good}
\end{figure}

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 \citet{ito16fermi}
have only $<30$ data points. Even in the samples with larger data
size, $\tau$ is not uniquely determined if it is too long.
Figure~\ref{fig:bad} 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. \citet{koz17drw} 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~\ref{fig:bad} 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~\ref{fig:good}. 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~\ref{tab:log} in the
Appendix. The designation of the objects to the subtypes FSRQ, LBL,
IBL, and HBL is taken from \citet{ito16fermi}. 

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{fig-Kanata-BL-06-54617-54946.pdf} 
 \end{center}
\caption{As for figure~\ref{fig:good} with the light curve
  data of AO~0235$+$16 between MJD~54617 and 54946. This is an example
  in which $\tau$ is too long.}\label{fig:bad}
\end{figure}

Blazars occasionally exhibit large prominent flares, as shown in
panel~(a) of figure~\ref{fig:samples}, 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~\ref{fig:pdbad}. 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~\ref{fig:pdbad}). 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. 

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{fig-Kanata-PD-12-55746-56125.pdf} 
 \end{center}
\caption{As for figure~\ref{fig:good} with the PD data of
  S5~0716$+$714 between MJD~55746 and 56125. Note that the scale of
  $\tau$ is logarithmic. This is an example of PD analysis in
  which $\tau$ is too short.}\label{fig:pdbad}
\end{figure}

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~\ref{tab:log} in the Appendix.

\subsection{Sparse multinomial logistic regression}

We construct a classifier for FSRQs and BL~Lac objects based on the
nine features described in \S~2.2. We use sparse multinomial
logistic regression (SMLR) to determine the classifier
(\cite{smlr}). 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:
\begin{eqnarray}
  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)},
\end{eqnarray}
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
\begin{eqnarray}
  \ell(\bm{w})=\sum_{i=1}^N \log P\left(\bm{y}_i|\bm{\theta}_i,\bm{w}\right).
\end{eqnarray}
Then, the solution of SMLR is expressed as
\begin{eqnarray}
  \hat{\bm{w}}= \argmax{\bm{w}} \{ \ell(\bm{w}) - \lambda\|\bm{w}\|_1 \},
\end{eqnarray}
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~\ref{tab:log} 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:
\begin{eqnarray}
  \theta_{i,j} = \exp\left\{ -\frac{|\bm{x}_j-\bm{x}_i|^2}{2\sigma_{\rm RBF}^2} \right\}, 
\end{eqnarray}
where $\sigma_{\rm RBF}^2$ is the bandwidth.

As mentioned in \S~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 \S~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
{\tt SMLR}\footnote{http://www.cs.duke.edu/\~{}amink/software/smlr}.

\section{Results}

We investigate the features that are discriminative for FSRQs and
BL~Lac objects based on SMLR and cross-validated AUC using the nine
features obtained from the data. SMLR has two hyper-parameters,
$\sigma_{\rm RBF}^2$ and $\lambda$. We first consider appropriate 
values for these two parameters for our study. 

A small $\sigma_{\rm RBF}^2$ leads to a complicated model with a
large number of samples retained in the classifier. Such a small
$\sigma_{\rm RBF}^2$ occasionally creates an island-like boundary.
Figure~\ref{fig:island} shows examples of the probability map of
BL~Lac type samples calculated with $\sigma_{\rm RBF}^2=1.0$ (left)
and $5.0$ (right) using two features, the light-curve CVE and PD
median. We set $\lambda=0.1$ in this case. As can be seen in the left
panel, the high probability region forms an ``island'' within the
surrounding low probability area. However, it is unlikely that the two
subtypes of blazars have such a complicated boundary. A linear or
slightly non-linear model, like that in the right panel, is more
reasonable. We confirmed that a classifier with large bandwidths
($\sigma_{\rm RBF}^2\gtrsim 5$) does not have an island-like boundary 
using our samples. In the following analysis, we set 
$\sigma_{\rm RBF}^2=5.0$. 

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{cv-med-island.pdf} 
 \end{center}
 \caption{Examples of complicated and simple boundaries. The color
   map indicates the probability map of a BL~Lac type sample calculated
   from the light-curve CVE and PD median with SMLR. The left and
   right panels show those obtained with bandwidth parameters of
   1.0 and 5.0, respectively. The blue and red circles indicate FSRQ
   and BL~Lac samples.}\label{fig:island} 
\end{figure}

The sparsity parameter, $\lambda$, also controls the complexity of the
model. We investigated the best AUC of all combinations of the nine
features against various values of $\lambda$. The result is given in
figure~\ref{fig:lambda}, showing that AUC becomes maximum around
$\lambda=1.0$. A model obtained with $\lambda>1$ is too simple to
appropriately classify the samples. On the other hand, a small
$\lambda$ ($<1$) leads to over-fitting. We set $\lambda=1.0$ in the
following analysis.

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{lambda_cv.pdf} 
 \end{center}
 \caption{Optimal AUC against $\lambda$.}\label{fig:lambda} 
\end{figure}

\begin{table*}
 
  \tbl{Variables used, AUC, and accuracy of top 20 models. A bullet symbol indicates the parameter was included in that model.}{%
    \begin{tabular}{ccccccccccrr}
      \hline
      \multicolumn{4}{c}{Light curve} &
      &\multicolumn{5}{c}{Polarization degree} & & \\
      \cline{1-4} \cline{6-10}
      CVE & Slope & $A_{\rm exp}$ & $\tau$ & & Median & CVE & Slope &
      $A_{\rm exp}$ & $\tau$ & AUC & Accuracy\\
      \hline
      $\bullet$ & --- &$\bullet$ & $\bullet$ & &
      $\bullet$ & --- &$\bullet$ & --- & --- & 0.907 & 0.842\\
      $\bullet$ & $\bullet$ &$\bullet$ & $\bullet$ & &
      $\bullet$ & --- & --- & --- & --- & 0.904 & 0.816\\
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & &
      $\bullet$ & --- & $\bullet$ & --- & --- & 0.888 & 0.842\\
      $\bullet$ & $\bullet$ &$\bullet$ & $\bullet$ & &
      $\bullet$ & --- & --- & --- & $\bullet$ & 0.885 & 0.789\\
      $\bullet$ & $\bullet$ &$\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & --- & --- & --- & 0.881 & 0.868\\
      $\bullet$ & --- & $\bullet$ & $\bullet$ & &
      $\bullet$ & --- & $\bullet$ & --- & $\bullet$ & 0.881 & 0.737\\
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & 0.878 & 0.816\\
      $\bullet$ & --- &$\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & $\bullet$ & --- & --- & 0.878 & 0.789\\
      $\bullet$ & --- & $\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & --- & --- & $\bullet$ & 0.875 & 0.789\\
      $\bullet$ & $\bullet$ &$\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & --- & $\bullet$ & $\bullet$ & 0.872 & 0.789\\
      $\bullet$ & $\bullet$ &$\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & --- & --- & $\bullet$ & 0.872 & 0.816\\
      $\bullet$ & $\bullet$ &$\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & $\bullet$ & --- & --- & 0.869 & 0.842\\
      $\bullet$ & --- & $\bullet$ & $\bullet$ & &
      --- & $\bullet$ & --- & --- & --- & 0.865 & 0.842\\
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ &$\bullet$ & $\bullet$ & --- & 0.862 & 0.763\\
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & &
      $\bullet$ & --- & $\bullet$ & --- & $\bullet$ & 0.862 & 0.737\\
      $\bullet$ & $\bullet$ &$\bullet$ & --- & &
      $\bullet$ & --- & $\bullet$ & --- &$\bullet$ & 0.859 & 0.789\\
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & &
      $\bullet$ & --- & $\bullet$ & $\bullet$ & $\bullet$ & 0.856 & 0.763\\
      $\bullet$ & --- & $\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & --- & 0.856 & 0.763\\
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & $\bullet$ & --- & $\bullet$ & 0.853 & 0.816\\
      $\bullet$ & $\bullet$ & $\bullet$ & $\bullet$ & &
      $\bullet$ & $\bullet$ & --- & $\bullet$ & --- & 0.853 & 0.789\\
      \hline
\end{tabular}}\label{tab:auc}
\end{table*}

We made an exhaustive test of all the parameters to find the most
important features (e.g., \cite{iga18ex}). The number of combinations
of the nine features is $2^9-1=511$. Using SMLR, we developed 511
classifiers using models with different combinations of parameters,
and calculated the AUCs for each. Table~\ref{tab:auc} lists the top 20
classifiers in the order of AUC values. For example, the classifier
with the highest AUC ($=0.923$) uses six features, that is, the CVE,
$A_{\rm exp}$, $\tau$ of the light curve, and the median, CVE, and
$\tau$ of the PD. It has an accuracy of 0.842. As can be seen in the
table, the correlation between the accuracy score and AUC is low in
the 20 models. This is probably due to the small sample size, and
indicates that a small difference in the AUC is not
important. We found that CVE and $A_{\rm exp}$ of the light
curve are used in all the top 20 models, and that $\tau$ of the light
curve and PD median are used in 19 models.  This result
suggests that these four features are essential to classify FSRQs and
BL~Lac objects. 

The probability of the BL~Lac type ($P_{\rm BL}$) for each sample
obtained with the classifier using the four parameters is listed in
table~\ref{tab:log} in the Appendix. We can determine the class of
each sample based on $P_{\rm BL}$. Table~\ref{tab:err} is the error
matrix for several different decision criteria: $P_{\rm BL}=0.5$,
$0.6$, and $0.7$. In the case of $P_{\rm BL}=0.5$, all the samples
classified as BL~Lac objects are indeed BL~Lac objects (Accuracy
$=1.0$). On the other hand, only six of the 12 FSRQs are correctly
classified as FSRQ, while the other six samples are misidentified. The
BL~Lac prediction accuracy improves with increased decision criterion
($P_{\rm BL}$), while the prediction accuracy of FSRQs decreases in
that case. The high rate of misidentified FSRQs suggests that a
significant portion of FSRQs cannot be distinguished from BL~Lac
objects based on the four features. 

\begin{table}
  \tbl{Error matrix and accuracy.}{%
    \begin{tabular}{crrr}
      \hline
      \multicolumn{4}{c}{$P_{\rm BL}=0.5$}\\
      \hline
             & \multicolumn{2}{c}{Reference} & Accuracy\\
      Classification & BL Lac & FSRQ & \\
      BL Lac         & 26     & 6    & 0.81\\
      FSRQ           & 0      & 6    & 1.00\\
      \hline
      \multicolumn{4}{c}{$P_{\rm BL}=0.6$}\\
      \hline
             & \multicolumn{2}{c}{Reference} & Accuracy\\
      Classification & BL Lac & FSRQ & \\
      BL Lac         & 25     & 4    & 0.86\\
      FSRQ           & 1      & 8    & 0.89\\
      \hline\multicolumn{4}{c}{$P_{\rm BL}=0.7$}\\
      \hline
             & \multicolumn{2}{c}{Reference} & Accuracy\\
      Classification & BL Lac & FSRQ & \\
      BL Lac         & 21     & 3    & 0.88\\
      FSRQ           & 5      & 9    & 0.64\\
      \hline
\end{tabular}}\label{tab:err}
\end{table}

\section{Discussion}
\subsection{Significance of the classifier}

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~\ref{tab:log}, 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~\ref{fig:auchist} 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.

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{auc_hist.pdf} 
 \end{center}
 \caption{Histograms of the AUC obtained from real samples (red)
   and artificial data generated from random numbers (black). See the
   text for details.}\label{fig:auchist}
\end{figure}

\subsection{Implications from the four features}

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{scatplot.pdf} 
 \end{center}
 \caption{Scatter plots of the selected features. The top, middle,
   and bottom panels show CVE, $\tau$, and PD median against
   $A_{\rm exp}$, respectively. The blue filled and open circles
   denote the correctly classified and misclassified FSRQ  samples,
   respectively. The red, orange, and black triangles are LBL, IBL,
   and HBL samples, respectively. The gray bars represent the
   95\% confidence intervals.}\label{fig:scat}
\end{figure}

Here, we discuss the implications of the results of \S~3.
Figure~\ref{fig:scat} 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. \cite{ito16fermi}). 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 \citet{yam19mrk421}. The data was obtained with
XRT/{\it 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~\ref{tab:xrt} 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~\ref{fig:scat}. 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~\ref{fig:scat}, 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. \citet{gio12type} proposed that
the two subtypes have a common nature, except for the beaming
factor. \citet{ito16fermi} 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~\ref{fig:scat}, 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 (\cite{sha13spec}). Figure~\ref{fig:qusample} 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~\ref{fig:qusample}, 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~\ref{fig:samples}. 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~\ref{fig:qusample} 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 (\cite{hag02bllac,sak13bllac}). 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. 

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{qu_sample.pdf} 
 \end{center}
 \caption{Polarization variations in the Stokes $Q/I$--$U/I$ plane of
   the FSRQ PKS~1510$-$089 in MJD~54759--55124 (left) and the LBL
   OJ~287 in MJD~55045--55451 (right).}\label{fig:qusample} 
\end{figure}

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 \citet{kas05blr}, 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
(\cite{fos98blseq,sha13spec,ghi17blseq}). 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. 

\citet{gio12type} 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. 
  
\subsection{Features of polarization variability}

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~\ref{fig:flxpd} 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 \S~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 (\cite{uem10bayes}). 

\begin{figure}
 \begin{center}
  \includegraphics[width=8cm]{flx-PD.pdf} 
 \end{center}
 \caption{Correlations of the flux and PD $\tau$. The blue
   circles and the red, orange, and black triangles denote the FSRQ, LBL,
   IBL, and HBL samples, respectively. The dashed line indicates
   $\tau\;{\rm (PD)}=\tau\;{\rm (Flux)}$.}\label{fig:flxpd} 
\end{figure}

\section{Summary}
We characterized the optical variability of blazars using the OU
process and investigated the features which are discriminative for the
two blazar subtypes, FSRQs and BL~Lac objects. Our summarized findings
are as follows:  
\begin{itemize}
\item Four features, namely, the variation amplitude, 
  $A_{\rm exp}$, characteristic timescale, $\tau$, and
  non-stationarity, CVE, from the light curve and the PD median 
  are essential to classify blazars into FSRQs and BL~Lac objects.
\item FSRQs are characterized by rare and large flares based on a
  large $A_{\rm exp}$ and/or CVE. We found that the X-ray variability
  of the HBL Mrk~421 also has large $A_{\rm exp}$ comparable to the
  optical variability of FSRQs. Hence, the characteristics of
  $A_{\rm exp}$ and CVE are governed not by the differences in the jet
  structure between the subtypes, but by the $\nu_{\rm peak}$ effect.
\item The high PD median of BL~Lac objects suggest that they tend to
  have a stable polarization component. The lack of such a component
  in FSRQs is possibly due to strong Compton cooling from a
  large BLR in sub-pc scale jets.
\item The variation timescale of PD is significantly shorter than
  that of the light curves. This may indicate that the relaxation
  timescale of the ordered magnetic field is shorter than the cooling
  timescale of the accelerated electrons.
\end{itemize}


\begin{ack}
This work was supported by a Kakenhi Grant-in-Aid (No. 25120007) from
the Japan Society for the Promotion of Science (JSPS).  
\end{ack}


