The Extremely Luminous Quasar Survey in the SDSS Footprint. I. Infrared-based Candidate Selection

For all SDSS sources, we use the point-spread function asinh magnitudes (Lupton et al. 1999) in the five optical band passes (ugriz; Fukugita et al. 1996). Throughout this paper, we only use AB magnitudes. We therefore converted the SDSS u-band and z-band magnitudes with ${u}_{\mathrm{AB}}=u^{\prime} -0.04\,\mathrm{mag}$ and ${z}_{\mathrm{AB}}\,=z^{\prime} +0.02\,\mathrm{mag}$. All magnitudes are corrected for galactic extinction using the extinction values from the Casjobs Data Release 13 (DR13) PhotoObjAll or PhotoPrimary tables (Schlafly & Finkbeiner 2011). The imaging data of the SDSS survey were completed in 2009 and cover a unique area of 14,555 deg². The magnitude limits (95% completeness for point sources) in the five optical bands are $21.6,\,22.2,22.2,\,21.3,{\rm{a}}{\rm{n}}{\rm{d}}\,20.7$ for u, g, r, i, and z, respectively.

WISE mapped the entire sky at 3.4, 4.6, 12, and 22 μm (W1, W2, W3, and W4). The recent AllWISE data release combines the data from the cryogenic and post-cryogenic (Mainzer et al. 2011) phases of the mission.⁴ The AllWISE source catalog achieved 95% photometric completeness for all sources with limiting magnitudes brighter than 19.8 and 19.0 (Vega: 17.1, 15.7) in W1 and W2. Saturation affects sources brighter than $8,7$ in the W1 and W2 bands. We restrict ourselves to photometry of the W1 ($3.4\,\mu {\rm{m}}$) and W2 ($4.6\,\mu {\rm{m}}$) infrared bands and convert them to AB magnitudes using ${\rm{W}}{1}_{\mathrm{AB}}={\rm{W}}1+2.699$ and ${\rm{W}}{2}_{\mathrm{AB}}={\rm{W}}2+3.339$. The WISE photometry is then extinction-corrected using the extinction coefficients ${A}_{{\rm{W}}1},{A}_{{\rm{W}}2}\,=0.189,0.146$ with the extinction values from the SDSS photometry.

2MASS has mapped the entire sky in the near-IR bands J ($1.25\,\mu {\rm{m}}$), H ($1.65\,\mu {\rm{m}}$), and ${K}_{{\rm{s}}}$ ($2.17\,\mu {\rm{m}}$). We use the 2MASS point source catalog (PSC), which was prematched to the WISE AllWise source catalog. The 2MASS PSC is generally complete at a level of $10\sigma$ photometric sensitivity for all sources brighter than 16.7, 16.4, and 16.1 (Vega: $15.8,15.0,14.3$) in the J, H, and ${K}_{{\rm{s}}}$ bands, respectively. Yet, the 2MASS PSC includes all sources with at least a signal-to-noise ratio of ${\rm{S}}/{\rm{N}}\geqslant 7$ in one band or ${\rm{S}}/{\rm{N}}\geqslant 5$ detections in all three bands. Furthermore, due to the confusion of sources close to the galactic plane, the photometric sensitivity is a strong function of galactic latitude. Based on the online documentation,⁵ we estimate the limiting magnitudes of the 2MASS PSC for higher latitudes to be $J=17.7$, $H=17.5$, and ${K}_{{\rm{s}}}=17.1$. We generally convert the 2MASS magnitudes to AB using ${J}_{\mathrm{AB}}=J+0.894$, ${H}_{\mathrm{AB}}=H+1.374$, and ${{K}_{s}}_{,\mathrm{AB}}={K}_{s}+1.84$ and correct for galactic extinction (${A}_{J},{A}_{H}$, ${A}_{{{\rm{K}}}_{{\rm{s}}}}=0.723,0.460,0.310$). While we test the machine-learning methods using 2MASS photometry, the final selection of the quasar candidate catalog uses only the 2MASS bands for the JKW2 color cut.

We selected all missed quasars with ${m}_{{\rm{i}}}\lt 17.0$ and $2.5\leqslant z\lt 5$ that fall into the well-surveyed spring sky of the SDSS footprint to understand why they were missed. We analyze whether the non-fatal and fatal image quality flags and the high-redshift inclusion regions of the original selection (Richards et al. 2002) play a role in them being overlooked. All of the quasars below with $z\geqslant 2.8$ are selected with our selection method and will be included in the ELQS spring sample (Paper II).

3.1.1. Q1208+1011

3.1.2. APM 08279+5255

3.1.3. SDSS J1622+0702A

3.1.4. B1422+231

3.1.5. CSO 167

3.1.6. CSO 1061

While we have shown the overlap between the missed quasars and the stellar locus in optical color space in Figure 2 panels (a)–(c), we also show the 2MASS and WISE J–K_s, K_s–W2 color space in panel (d). Here, the stellar locus clearly separates from the distribution of missed quasars. After examining the distribution of known quasars to known stars in this color–color space, we qualitatively determined the color cut that separates these distributions best. In Vega magnitudes, the color cut reads

$$\begin{eqnarray}&&{\mathrm{Ks}}_{\mathrm{Vega}}-{\rm{W}}{2}_{\mathrm{Vega}}\geqslant 1.8-0.848\cdot ({J}_{\mathrm{Vega}}-{\mathrm{Ks}}_{\mathrm{Vega}}),\end{eqnarray} \tag{ 1 }$$

while in AB magnitudes it changes to

$$\begin{eqnarray}&&\mathrm{Ks}-{\rm{W}}2\geqslant -0.501-0.848\cdot (J-\mathrm{Ks}).\end{eqnarray} \tag{ 2 }$$

The color cut separates quasars from stars because of their difference in the K–W2 flux ratio (or K–W2 color). This flux ratio is fundamentally different in stars (up to spectral type T) and quasars (up to redshifts of $z\approx 5$). For those stars and quasars, the stellar flux declines more strongly than the quasar flux between the K-band and the W2-band. Hence, quasars will have a redder K–W2 flux ratio (or color). Emission lines, like ${\rm{H}}\alpha$ in the K-band at $z\approx 2.1\mbox{--}2.50$ or ${\rm{H}}\beta$ in the K-band at $z\approx 3.2\mbox{--}3.8$, do affect the K–W2 flux ratio, but their influence is negligible. As a result, the JKW2 color cut cannot discriminate between quasars at different redshifts.

Six known quasars lie beneath the JKW2 color cut within the stellar locus in Figure 2, panel (d). No spectra were available for any of these objects in the literature. A closer examination of their identifications and photometry concluded that either the photometry or the identification reference seems to be unreliable.

We test the color cut in a $\sim 70\,{\deg }^{2}$ region (${\rm{R.A.}}\,=120-130\,\deg$, $\mathrm{decl}.\,=\,40-50\,\deg$), shown in Figure 4. For this purpose, we selected all sources in SDSS DR13 (PhotoObjAll) brighter than SDSS ${m}_{{\rm{i}}}=18.5$, matched the sources with the AllWISE source catalog, including matched 2MASS PSC photometry, and retrieved spectral identification from the SDSS DR13 (SpecObj) catalog, where possible.

Zoom In Zoom Out Reset image size

The full test field includes a total of 1,327,439 sources detected with full SDSS photometry, of which 1,155,203 sources have full 2MASS and WISE W1 and W2 photometry.

For our purposes, we estimate the efficiency of the color cut using the knowledge about the spectroscopically identified quasars (at all redshifts), stars, and galaxies, and the total number of sources with 2MASS and WISE W1 and W2 colors. The results are summarized in Table 2.

Table 2.  Selection Criterion on the $70\,{\deg }^{2}$ Test Region: JKW2 Color Cut

Data Sample	Quasars	Stars	Galaxies	Full Sample	No Spectral Id.
Full SDSS colors	526	3327	5685	1,327,439	1,317,901
SDSS+2MASS+W1W2	301	2138	5517	1,155,203	1,147,247
JKW2 color cut	298	8	2706	25,470	22,458
Fraction	99.00%	0.37%	49.05%	2.20%	1.96%
SDSS+2MASS+W1W2	207	2081	39	841,926	839,599
petroRad_i $\leqslant 2.0\,\mathrm{arcsec}$
JKW2 color cut	205	7	15	1296	1,069
petroRad_i $\leqslant 2.0\,\mathrm{arcsec}$
Fraction	99.03%	0.34%	38.46%	0.15%	0.13%

Download table as:  ASCIITypeset image

The color cut clearly separates stars and quasars. It reduces the total number of sources in the test region to 25,470 out of 1,155,203 ($2.20 \%$). The stellar contamination is greatly reduced. Only $0.37 \%$ (8/2138) of spectroscopically identified stars in SDSS DR13 make the color cut. Of all 526 quasars in this region, 301 have 2MASS and WISE W1 and W2 colors, and of these, 298 are included in the JKW2 color cut. Considering only the quasars with full near-IR photometry, the color cut has an estimated completeness of $\sim 99 \%$.

The quasar sample in the $70\,{\deg }^{2}$ test region is rather small. Therefore, we test our color cut with all quasars from the DR7 and DR12 quasar catalogs that have 2MASS photometry in the J and K_s bands and WISE W2 photometry. This sample includes 5945 quasars, out of which 5927 are included in the color cut, a fraction of more than 99%.

Galaxies, however, straddle the JKW2 color cut and only half of all spectroscopically identified galaxies (2706/5517) are excluded. Therefore, we believe the majority of the remaining 25,470 sources above the color cut to be galaxies.

In total, the JKW2 color cut manages to eliminate the majority of stellar sources while retaining a highly complete sample of bright quasars uniform in redshift and SDSS i-band magnitude. However, in order to efficiently use the color cut, we need, first, to reject galaxy contaminants and second, to estimate the photometric quasar redshifts.

The 2MASS and WISE surveys do not reach the depth of optical surveys, like SDSS. Therefore, fainter sources might only show spurious detections or no detections at all. Even for a bright quasar survey, this is a concern regarding the photometric completeness. We analyze the fraction of sources in the MQC without detections in the 2MASS and WISE bands necessary for the JKW2 color cut. More than $99.5 \%$ of ${m}_{{\rm{i}}}\lt 18.0$ quasars at $z\gt 2.5$ in the SDSS-matched MQC are detected with ${\rm{S}}/{\rm{N}}\gt 5$ in the W1 and W2 bands. 2MASS is shallower but still allows for about $\gt 80 \%$ of ${m}_{{\rm{i}}}\lt 18.0$ quasars at $z\gt 2.5$ in the SDSS-matched MQC to be detected in J and K_s. The photometric completeness of all surveys will be taken account when we calculate the overall survey completeness in Paper II.

The SDSS survey uses the type parameter to differentiate between different source types. It allows point sources (type = 6) and sources that are classified as extended (type = 3) to be distinguished from one another. However, we do not generally want to exclude quasar lenses in our selection, and thus do not use this parameter. Instead, we turn to a more quantitative measure, the Petrosian radius.

At redshifts above $z=2.5$, even lensed quasars will be reasonably compact (petroRad_i $\lt 2^{\prime\prime}$) if they are not strongly distorted, but might still be categorized as extended objects (SDSS flag type = 3; e.g., the lensed quasar B1422+231).

In Figure 5, we show that a Petrosian radius of $2^{\prime\prime}$ in the SDSS i-band strongly reduces the number of galaxies in the total test region while mainly reducing the number of quasars between $0\leqslant z\leqslant 1.0$, where its host is likely resolved in the SDSS photometry. In the targeted redshift range of $2.5\leqslant z\leqslant 4.0$, this cut on the Petrosian radius in the SDSS i-band does not reject any quasars.

Zoom In Zoom Out Reset image size

We summarize the effects of the cuts in the Petrosian radius on the $70\,{\deg }^{2}$ test region (Figure 4) in Table 3. The limit on the Petrosian radius reduces the number of spectroscopically identified quasars from 301 to 207, out of which 205 ($99 \%$) are within the JKW2 color cut. As expected, the Petrosian radius limit does not significantly alter the number of stars within the color cut. However, it reduces the total number of spectroscopically identified galaxies that make the JKW2 color cut from 2706 to 15, and the total number of sources from 25,470 to 1296. The majority of the rejected sources will be galaxies with a contribution from low-redshift ($z\leqslant 1.0$) quasars.

Table 3.  Selection Criterion on the $70\,{\deg }^{2}$ Test Region: Petrosian Radius

Criterion		petroRad_i $\leqslant 3.0\,\mathrm{arcsec}$	petroRad_i $\leqslant 2.0\,\mathrm{arcsec}$	petroRad_i $\leqslant 1.5\,\mathrm{arcsec}$
Full sample	1,155,203	965,887	841,925	544,643
Galaxies	5517	578	39	12
Quasars	301	257	207	166
Stars	2138	2112	2081	1781
No spectral id.	1,147,247	962,940	839,598	542,684

Download table as:  ASCIITypeset image

5.1.1. The Empirical Quasar Catalog

In order to devise training sets for the classification and photometric redshift estimation, we combined the SDSS quasar catalogs with 2MASS All-Sky and WISE AllWISE photometry.

We select all quasars from the SDSS DR7 and DR12 quasar catalogs (DR7Q and DR12Q, respectively) and matched them with SDSS DR13 (PhotoObjAll) sources in a $2^{\prime\prime}$ radius to obtain additional information about the sources. The resulting catalog is then matched to WISE AllWISE and 2MASS All-Sky photometry in a $1^{\prime\prime}$ radius to make sure that the matches are accurate. We exclude magnitude outliers ($2\lt i$ and $i\gt 29;$ $i\equiv$ apparent i-band magnitude) from the catalog. The magnitudes are then corrected for galactic extinction (${m}_{{\rm{i}}}\equiv$ extinction-corrected apparent i-band magnitude). Fluxes and flux errors (where possible) are then calculated from the extinction-corrected magnitudes. To ensure that the magnitude error distributions are free of outliers, objects with positive values in the SDSS flags INTERP and DEBLEND_AT_EDGE are excluded.

The resulting quasar catalog is termed the "empirical" quasar catalog and has a total of 215,087 objects with full SDSS photometry. An overview of the different sample sizes drawn from this catalog is given in Table 4.

Table 4.  Different Training Sets for the Photometric Redshift Estimation and Quasar Classification

Data Set	Photometry	Constraints	Size
Emp. QSOs	SDSS	⋯	215,087
Emp. QSOs	SDSS	${m}_{i}\lt 18.5$	12,408
Emp. QSOs	SDSS+W1W2	⋯	153,890
Emp. QSOs	SDSS+W1W2	${m}_{i}\lt 18.5$	12,388
Emp. QSOs	SDSS+2MASS+W1W2	⋯	4815
Emp. QSOs	SDSS+2MASS+W1W2	${m}_{i}\lt 18.5$	4021
Emp. QSOs	SDSS+2MASS+W1W2	${m}_{i}\lt 18.5$, JKW2	4015
DR13 Stars	SDSS	⋯	387,854
DR13 Stars	SDSS	${m}_{i}\lt 18.5$	219,375
DR13 Stars	SDSS+W1W2	⋯	245,326
DR13 Stars	SDSS+W1W2	${m}_{i}\lt 18.5$	197,798
DR13 Stars	SDSS+2MASS+W1W2	⋯	174,218
DR13 Stars	SDSS+2MASS+W1W2	${m}_{i}\lt 18.5$	159,211
DR13 Stars	SDSS+2MASS+W1W2	${m}_{i}\lt 18.5$, JKW2	209

Download table as:  ASCIITypeset image

The column "photometry" refers to the photometric information included in this training set. Not all sources in the full empirical quasar catalog have information in the 2MASS or WISE filter bands because the 2MASS and WISE surveys do not always reach the same depth as the SDSS photometry. If we mandate values in the WISE W1 and W2 bands or the 2MASS bands and their corresponding errors, the size of the training set decreases.

One additional question that naturally arises is: Can we find a bright population of quasars that does not exist in the training set? The distribution of the brightest quasars does not have inherently different spectral slopes. As a consequence, their flux ratios (or colors) do not differ substantially from less brighter ones. Since the main features used in the machine-learning methods are the flux ratios of adjacent photometric bands, bright quasars should be comprehensively selected without problems.

5.1.2. The Empirical Star Catalog

The random forest method (Breiman 2001) uses an ensemble of decision (or regression) trees to vote for the most popular class (value) regarding a classification (regression) problem. The method is a supervised machine-learning method. Therefore, it requires a training set to learn from. Furthermore, random forests are non-parametric, inherently allow for multiple ($\gt 2$) classes, and avoid the problem of overfitting.

Since random forests rely on decision (or regression) trees, we will shortly introduce their operation.

The features of the training set, which are fluxes and flux ratios for our purposes, generate the multidimensional input space for the classification or regression. A decision tree divides this input space into cuboid regions along the feature axis. Each of these regions contains one or more data points that determine the target class or target value for that cuboid.

This structure is built to correspond to a binary tree, where at each node a decision on one input feature is made on how to split the tree into two branches. This is called recursive binary partitioning. The decision is reached using a criterion that encourages the formation of regions, where the majority of the data points belong to only one class.

To predict the target class or target value for a new data point, the class or value of the cuboid region that the data point falls into is chosen.

Decision trees are invariant under scaling of the feature values and to the inclusion of irrelevant features. They also produce reviewable models of the data that are easy to visualize and inherently work with multiple ($\gt 2$) classes.

In the random forest method, each decision tree is fit on a bootstrap subsample of equal size to the original input sample. During the construction of the decision tree, a random subset of all available features is used to determine the best splits for all nodes. In our case the number of random features considered is taken to be the square root of the total number of features. A full decision tree is grown without pruning the tree during the construction. The forest contains a large number of these trees ($\gt 100$), each of which gives a classification (target value) for every object. The class with the most votes is then chosen to be the final class of the object. For random forest regression, the mean of all target values is calculated to be the final regression value.

While one decision tree is prone to overfitting the training sets, random forests overcome this problem through the averaging process of many randomized decision trees.

Breiman et al. (2003) first proposed this machine-learning technique as an astronomical classification tool to find quasars. In astronomy, the method has since been applied for the classification of variable stars (Dubath et al. 2011; Richards et al. 2011) and variable quasars (Pichara et al. 2012), photometric redshift estimation (mainly aimed at galaxies; Carliles et al. 2010; Carrasco Kind & Brunner 2013), and quasar classification (Carrasco et al. 2015).

We are using the implementation of the random forest classifier and regressor provided by the scikit-learn (Pedregosa et al. 2011) Python library with many of the default parameters.

For the construction of the binary tree, we use the Gini impurity to determine the best split at each step. While the original random forest method (Breiman 2001) lets each classifier vote for the final class, this implementation averages the probabilistic predictions of each classifier to find the final probabilities for each class.

We adopt the hyperparameters that control the size of the decision trees (min_samples_split, max_depth) as well as the size of the forest (n_estimators) to find the best classification/regression model. The best-fit values for these hyperparameters are found using a limited grid search.

SVMs (Vapnik 1995, 1998; Burges 1998) offer a sparse kernel method for classification, regression, and novelty detection. Similar to random forests, SVMs belong to supervised machine-learning methods and therefore rely on a well-constructed training set. Fundamentally a two-class classifier, the extension of SVMs to multiclass ($\gt 2$) classification is problematic.

In the case of classification, the algorithm calculates a decision boundary (hyperplane) to divide the multidimensional feature space into two regions according to the two classes. The decision boundary is constructed to maximize the smallest distance between itself and any of the data samples. In the end, only a small subset of original data points is necessary to define the decision boundary. These data points are the support vectors. One strength of the SVM lies in this reduction of information to only a few data samples that define how to split the feature space. As a result, the SVM is very fast in making predictions.

SVMs as kernel methods use an algorithm that allows for the use of a kernel function to implicitly transform the original feature space into a higher-dimensional feature space. In such an algorithm, the input vector enters only in the form of its scalar product and is then substituted by the kernel. Therefore, the coordinates in the higher-dimensional feature space do not have to be calculated explicitly, only their inner product. The kernel allows for complex decision boundaries that are nonlinear in the original feature space and allow for more complicated distributions of the two classes.

However, many data sets do not have fully separable classes even if nonlinear kernels are used. As a solution, the algorithm allows for training data to lie on the "wrong" side of the decision boundary. These data points are then misclassified but ignored by the algorithm. They are called slack variables and allow for a better generalization of the classification. In older formulations of the SVM algorithm, the amount of slack variables are controlled with the parameter C.

Similar to other regression problems, SVM regression also seeks to minimize a regularized error function. This error function incorporates the previously introduced slack variables as well as an -insensitive error function (Vapnik 1995). Therefore, data points within of the regression model have no associated error, and the number of data points outside of this region is controlled using the slack variable parameter C.

SVMs have been widely used in astronomy for galaxy (Wadadekar 2005; Wang et al. 2008) and quasar (Han et al.2016) photometric redshift estimation and source classification (Gao et al. 2008; Huertas-Company et al. 2008; Kim et al.2012; Peng et al. 2012; Kurcz et al. 2016).

We use the scikit-learn (Pedregosa et al. 2011) implementation of the support vector regression for a comparison of the photometric redshift estimation with the random forest method above. We use the radial basis function kernel (kernel = "rbf"), which adds the hyperparameter gamma. In order to estimate the optimal hyperparameters C, epsilon, and gamma, we carry out limited grid searches.

Furthermore, it is important to note that the features of the input and test vectors need to be normalized because SVMs are not scale invariant.

With the JKW2 color cut rejecting the majority of stars and the limit on the Petrosian radius filtering out unwanted galaxies, lower-redshift quasars become major contaminants for our survey targeted at $z\geqslant 2.8$ quasars. Hence, it is critical to estimate the quasar redshift and use this photometric redshift as a criterion for our candidate selection. Instead of relying on optical-color cuts (e.g., Richards et al. 2002), we aim to utilize the full photometric information given with SVM regression (SVR) and random forest (RF) regression.

The training sets for both algorithms are drawn from the empirical quasar catalog, which is based purely on the SDSS DR7 and DR12 quasar catalogs, as described above.

We use three different training sets built from the empirical quasar data to test the effects of different feature sets and the size of the training sets on the regression outcome. The first (SDSS; Table 6, rows 1 and 5) includes all sources with full SDSS photometry, the second (SDSS+W1W2; Table 6, rows 3 and 7) includes W1 and W2 information in addition to full SDSS photometry, whereas the last (SDSS+W1W2, ${m}_{{\rm{i}}}\lt 18.5;$ Table 6, rows 4 and 8) uses only sources with SDSS i-band magnitude brighter than ${m}_{{\rm{i}}}\lt 18.5$ of the SDSS+W1W2 subset. In general, it is advisable to always include the largest amount of information possible. However, because we want to evaluate the benefit of including the W1 and W2 bands in addition to the SDSS photometry, we limit ourselves to the SDSS photometry for a test case on the SDSS+W1W2 subset (Table 6, rows 2 and 6).

Table 6.  Results of the Photometric Redshift Estimation Methods

Data Set	Training/Test Size	Constraints	Features	Algorithm	${\delta }_{0.3}$	${\delta }_{0.2}$	${\delta }_{0.1}$	σ	R2
DR7DR12Q	172069/43018	SDSS fl.r.	SDSS	RF	0.87	0.81	0.65	0.483	0.654
DR7DR12Q	123112/30778	SDSS+W1W2 fl.r.	SDSS	RF	0.84	0.76	0.58	0.503	0.624
DR7DR12Q	123112/30778	SDSS+W1W2 fl.r.	SDSS+W1W2	RF	0.96	0.94	0.86	0.277	0.884	⋆
DR7DR12Q	9910/2478	SDSS+W1W2 fl.r., ${m}_{{\rm{i}}}\lt 18.5$	SDSS+W1W2	RF	0.98	0.97	0.93	0.189	0.937
DR7DR12Q	172069/43018	SDSS fl.r.	SDSS	SVR	0.87	0.82	0.66	0.491	0.642
DR7DR12Q	123112/30778	SDSS+W1W2 fl.r.	SDSS	SVR	0.85	0.77	0.60	0.512	0.614
DR7DR12Q	123112/30778	SDSS+W1W2 fl.r.	SDSS+W1W2	SVR	0.96	0.92	0.81	0.291	0.873
DR7DR12Q	9910/2478	SDSS+W1W2 fl.r., ${m}_{{\rm{i}}}\lt 18.5$	SDSS+W1W2	SVR	0.98	0.96	0.89	0.194	0.933

Download table as:  ASCIITypeset image

We do not build a training set based on full SDSS, 2MASS, and W1 and W2 photometry because the number of quasars with full 2MASS photometry is too small (∼1000) to allow for sufficient training in the large feature space. In addition, a cut on the i-band magnitude at ${m}_{{\rm{i}}}\lt 18.5$ strongly reduces the number of higher-redshift quasars in the training set. As a consequence, the feature space is not sufficiently populated at higher redshifts, and the method is biased against those quasars.

For training sets with only SDSS features, we use the four adjacent flux ratios (u/g, g/r, r/i, i/z) from the five photometric bands and the SDSS i-band magnitude. When we include WISE photometry (SDSS+W1W2), we expand the flux ratios accordingly (+z/W1, W1/W2) and add the W1 magnitude to the feature set.

For each subset, we use a grid search on the hyperparameters of the random forest and SVM regression to determine the best regression model. We calculate the regression results for all subsets using both machine-learning methods to compare them to each other.

The grid for the random forest regression has the following hyperparameters: n_estimators = [50, 100, 200, 300], min_samples_split = [2, 3, 4], and max_depth = [15, 20, 25] (36 combinations).

For the SVM regression, we use a grid of C = [10, 1.0, 0.1], gamma = [0.01, 0.1, 1.0], epsilon = [0.1, 0.2, 0.3] (27 combinations).

We use five-fold cross-validation on the full training set and always train on 80% of the full training set, while the remaining 20% are used to test the regression. Before we continue to discuss the results of the calculations, we will first introduce the typical regression metrics.

5.4.1. Regression Metrics

To evaluate the success of the photometric redshift estimation methods, we use the standard R² regression score as well as the residual between the photometric and spectroscopic redshifts, ${\rm{\Delta }}z={z}_{\mathrm{spec}}-{z}_{\mathrm{reg}}$.

The R² score, also called the coefficient of determination, gives a measure of the goodness of fit of a model. Let us assume that the true redshifts are denoted by z_i and the predicted redshift values are denoted by ${\hat{z}}_{i}$. The R² score is then calculated from the total sum of squares ${{SS}}_{\mathrm{tot}}$ and the residual sum of squares ${{SS}}_{\mathrm{res}}$, where $\bar{z}$ is the mean redshift, according to

$$\begin{eqnarray}&&{{SS}}_{\mathrm{tot}}=\displaystyle \sum _{i}{({z}_{i}-\bar{z})}^{2},{{SS}}_{\mathrm{res}}=\displaystyle \sum _{i}{({z}_{i}-{\hat{z}}_{i})}^{2},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{R}^{2}=1-\displaystyle \frac{{{SS}}_{\mathrm{res}}}{{{SS}}_{\mathrm{tot}}}.\end{eqnarray} \tag{ 4 }$$

In the best-case scenario, the residual sum of squares will go to 0 and the R² score will reach 1. If the model always predicts the mean value of the training data $\bar{z}$, then it results in an R² score of 0. For bad models, it is possible that the R² reaches negative values. We use R² as a measure to compare models with one another to find the model that best represents the data. Therefore, we focus on the relative values of the different models and do not interpret the absolute one.

In the literature, the goodness of the photometric redshift estimation is often measured by the fraction of test quasars N with absolute redshift residuals $| {\rm{\Delta }}z| =| {z}_{i}-{\hat{z}}_{i}|$ smaller than a given residual threshold e (Bovy et al. 2012; Peters et al. 2015; Richards et al. 2015),

$$\begin{eqnarray}&&{\rm{\Delta }}{z}_{e}=\displaystyle \frac{N(| {z}_{i}-{\hat{z}}_{i}| \lt e)}{{N}_{\mathrm{tot}}}.\end{eqnarray} \tag{ 5 }$$

Typical values chosen are $e=0.1,0.2,0.3$. However, in many cases the redshift normalized residuals are used instead,

$$\begin{eqnarray}&&{\delta }_{e}=\displaystyle \frac{N(| {z}_{i}-{\hat{z}}_{i}| \lt e\cdot (1+{z}_{i}))}{{N}_{\mathrm{tot}}}.\end{eqnarray} \tag{ 6 }$$

5.4.2. Results

The results of the regression calculations using RF regression and SVR are shown in Table 6.

The top half of the table details the results from the RF method, whereas the bottom half shows the corresponding results for the SVR. Both algorithms perform similarly well for all four training and feature sets used. If the WISE W1 and W2 photometry is included, the RF method performs slightly better.

A comparison between the qualitative results of the two methods is given in Figure 6. The left panel shows the results of the second and sixth rows of Table 6 while the right panel shows the third and seventh rows of the same table. The orange and blue curves are the histograms of the redshift residuals ${\rm{\Delta }}z$ for the SVR and RF regression, respectively. While the RF method shows a tighter distribution of the histogram around ${\rm{\Delta }}z=0$, the different colored curves show the same general qualitative behavior in both panels.

Zoom In Zoom Out Reset image size

The main difference between both methods lies in the computing time. For the subset referenced with the ⋆ (Table 6), the computing time for RF regression amounted to $281\,{\rm{s}}$, whereas the same subset calculated with SVR took $4682\,{\rm{s}}$ for similarly good results. This is an inherent drawback of the SVR method as the computing timescales with ${ \mathcal O }({N}^{3})$, where the number of objects in the training set is N. For comparison, the random forest method computed timescales with the number of trees in the Forest T and the depth of each tree D as ${ \mathcal O }(T\cdot D)$. Hence, for large training sets, the random forest regression should always be preferred if both regression methods perform equally well.

If one compares the R² score of the second and third rows of Table 6, it becomes obvious that the inclusion of more photometric features clearly improves the regression results. While the R² score summarized the quantitative effect in one number, we show the qualitative differences in Figure 6. From the left to the right panels, the distribution of test quasars tightens significantly around ${\rm{\Delta }}z=0$ and the second bump around ${\rm{\Delta }}z\approx -0.8$ disappears.

The price is the reduction of the training sample from 172,069 quasars to 123,112. However, a comparison between the first and second rows of Table 6 shows that a larger training set improves the regression results if the same number of features are considered. This suggests that an addition of training objects still leads to a better characterization of the target values in the feature space.

The best regression results are achieved on the SDSS+W1W2 training set with ${m}_{{\rm{i}}}\lt 18.5$. Because we ignore the distribution of fainter quasars here, the photometric errors on this training set are smaller. As a result, the regression achieves better results. However, we have to advise caution here. The ${m}_{{\rm{i}}}\lt 18.5$ cut severely reduces the number of training and test objects. Because the redshift distribution of the empirical training set is dependent on the i-band magnitude of the quasars, it biases the training set against higher-redshift sources.

We choose the RF method with the SDSS+W1W2 training set and all available features for our photometric redshift estimation (marked with a ⋆ in Table 6).

Figure 7 shows the density map of the random forest regression test results using this training set as a function of the measured spectroscopic redshift. The color of the map corresponds to the number of test objects per rectangular bin. The photometric redshift of the majority of test objects is well estimated. However, a distribution of outliers still persists at spectroscopic redshifts of $z\approx 0.8,1.6$, and 2.0–2.3.

Zoom In Zoom Out Reset image size

5.4.3. Comparison to Other Methods

The difficulty in comparing different methods for redshift estimation to another is that the same training set should be used in all cases to create the model. Otherwise, model comparisons, even using the same regression metrics, are not meaningful.

In a recent study on the photometric selection of quasars (Yang et al. 2017, ApJ in press), the authors incorporate asymmetries into a model of the relative flux distributions of quasars. As part of their work, they carry out a comparison (see their Table 1; Yang et al. 2017, ApJ in press) between their new method (Skew-t) and a range of frequently used algorithms in the literature (XDQSOz, Bovy et al. 2012; KDE, Silverman 1986; CZR, Weinstein et al. 2004) based on the same photometric test and training sample. Their method achieves slightly better results than the others on the same training set based only on SDSS photometry.

While their training/test sample includes a range of high-redshift quasars from the literature, the majority of the sample is also based on the DR7 and DR12 quasar catalogs. Therefore, we compare our results (Table 6) with theirs in Table 7. The RF and SVR are outperformed by the Skew-t method on the SDSS feature/training set (see $\delta z=0.1$); meanwhile, the RF method performs equally well, once the WISE W1 and W2 bands are included. Hence, random forests can keep up with other modern methods of photometric redshift estimation.

Table 7.  Comparison of the Photo-z Regression

Method	Features/Training Set	${\delta }_{0.2}$	${\delta }_{0.1}$	Total Data Set Size
Skew-t	SDSS	0.82	0.75	304,241
RF	SDSS	0.81	0.65	215,087
SVR	SDSS	0.82	0.66	215,087
Skew-t	SDSS+W1W2	0.93	0.87	229,653
RF	SDSS+W1W2	0.94	0.86	153,890
SVR	SDSS+W1W2	0.92	0.81	153,890

Download table as:  ASCIITypeset image

To classify our quasar candidates as stars or quasars, we only use the random forests, because the SVM algorithm is not suited for simultaneous classification with more than two classes.

The training set for the classification problem is built from the empirical quasar and star catalogs described above. Both catalogs are added to form the full empirical training set on which the random forest classification for different stellar spectral classes and quasar redshift classes will be trained.

However, we limit ourselves to stars with spectral classes A, F, G, K, and M since these classes have enough objects for sufficient training and the main contaminants of quasars with redshifts z = 2–5 are K and M stars. Most L and T dwarfs are too faint for our ELQS and conflict only with even higher-redshift quasars in color space. We also excluded O and B stars because they are very rare, and their blue colors are hardly confused with $z\gt 2$ quasars.

Since quasar colors change considerably over redshift, we divide quasars into four redshift classes similar to Richards et al. (2015). The classes are designed to split the quasars at redshifts, where dominating features change the u-band to g-band flux ratio. The redshift classes are "vlowz" with $0\lt z\leqslant 1.5$, "lowz" with $1.5\lt z\leqslant 2.2$, "midz" with $2.2\lt z\leqslant 3.5$, and "highz" with $3.5\lt z$. At z = 1.5, the Lyα emission line is just blueward of the u-band and the CIV emission line is still in the g-band. The second break at z = 2.2 marks the point where the Lyα line is leaving the u-band and the last break at z = 3.5 is marked by a strong flux decrease in the u-band as the Lyα forest absorbs flux blueward of the Lyα line. We use these four quasar classes and five stellar classes to form the nine main labels for our classification problem. For evaluation purposes, we introduce the binary labels "STAR" and "QSO," which encompass all stellar classes and all quasar classes, respectively. We will use them to calculate the completeness of the quasar selection against stars in general.

The RF algorithm is trained only on the photometric information. For classifications using only the SDSS color space, the features are the four adjacent flux ratios (u/g, g/r, r/i, i/z) from the five photometric bands and the SDSS i-band magnitude. When we include the 2MASS or WISE photometry, we expand the flux ratios accordingly and add the J-band magnitude or the W1 magnitude to the feature set.

We investigate the effects of two magnitude limits in the SDSS i-band, ${m}_{{\rm{i}}}\lt 21.5$ and ${m}_{{\rm{i}}}\lt 18.5$, and the use of different sets of photometric features (SDSS, SDSS+W1W2, and SDSS+2MASS+W1W2) on the results of the random forest classification. The six subsets of the full empirical training set have different sizes, since not every object always has the full photometric information required. If information in one photometric band or its associated error is missing, we reject the source from the training set.

For each of the six subsets, we find the best combination of the training hyperparameters by calculating the classification results on a grid of n_estimators = [50, 100, 200, 300], min_samples_split = [2, 3, 4], and max_depth = [15, 20, 25] (36 combinations). We use five-fold cross-validation and always train on 80% of the full training set, while the remaining 20% are used to test the classification.

5.5.1. Classification Metrics

In order to measure the performance of the classification, we will introduce three standard classification metrics: the precision, the recall, and the F₁ score (Bishop 2006).

Precision (p) is defined as the ratio of the true positives (t_p) to the sum of true and false positives (${t}_{p},+\,{f}_{p}$). Usually in quasar selection one speaks of the purity/efficiency of the selection synonymous to the precision of the selection.

Recall (r) is the ratio of true positives to the sum of true positives and false negatives (${t}_{p}+{f}_{n}$). Regarding quasar selections, the recall is equivalent to the completeness of the selection,

$$\begin{eqnarray}&&p=\displaystyle \frac{{t}_{p}}{{t}_{p}+{f}_{p}}\,r=\displaystyle \frac{{t}_{p}}{{t}_{p}+{f}_{n}}.\end{eqnarray} \tag{ 7 }$$

The harmonic mean of the precision and recall values is the traditional F-measure or balanced F-score. The F₁ score reaches its best value at 1 and the worst score at 0,

$$\begin{eqnarray}&&{F}_{1}=2\cdot \displaystyle \frac{\mathrm{precision}\cdot \mathrm{recall}}{\mathrm{precision}+\mathrm{recall}}.\end{eqnarray} \tag{ 8 }$$

For multiclass classification problems, one can define a precision, recall, and F₁ score for each class individually against all other classes.

Also, an average precision, recall, and F₁ score weighted by the number of occurrences in each true class can provide an idea of how well the classifier generally works for a problem with multiple classes.

A helpful visualization for the results of the classification problems is the confusion matrix C. Each entry ${C}_{i,j}$ is the number of objects known to be in class i, but predicted in class j. Therefore, the entries ${C}_{i,i}$ show the true positives (t_p) for each class i. All other values in row i show the number of false negatives (f_n), objects predicted to belong to other classes, although they truly belong to class i. All other values in column i are false positives (f_p), the values predicted to belong to class i, although they truly belong to other classes in the sample.

5.5.2. Results

We present the results of the random forest classifications in Table 8. It shows the six different subsets of the full empirical training set along with their respective constraints and features used. We show the results of the classification as the precision (p), recall (r), and F₁ measures for the "highz" quasar class as well as the grouped classes of all quasars ("QSO") and stars ("STAR").

Table 8.  Results of the Random Forest Classification on the Full Empirical Training Set

Training/Test size	Constraints	Features	p/r/F1 (highz)	p/r/F1 (QSO)	p/r/F1 (STAR)
183,108/45,778	SDSS fl.r., ${m}_{{\rm{i}}}\lt 18.5$	SDSS	0.92/0.71/0.80	0.88/0.96/0.92	1.00/0.99/1.00
167,076/41,770	SDSS+W1W2 fl.r., ${m}_{{\rm{i}}}\lt 18.5$	SDSS+W1W2	0.97/0.86/0.91	0.91/1.00/0.95	1.00/0.99/1.00
129,934/32,485	SDSS+2MASS+W1W2 fl.r., ${m}_{{\rm{i}}}\lt 18.5$	SDSS+2MASS+W1W2	0.88/0.88/0.88	0.93/0.99/0.96	1.00/1.00/1.00
442,529/110,634	SDSS fl.r., ${m}_{{\rm{i}}}\lt 21.5$	SDSS	0.87/0.87/0.87	0.77/0.94/0.85	0.96/0.85/0.90
313,453/78,364	SDSS+W1W2 fl.r., ${m}_{{\rm{i}}}\lt 21.5$	SDSS+W1W2	0.92/0.95/0.93	0.88 /1.00 /0.94	1.00/0.92/0.96	⋆
141,837/35,460	SDSS+2MASS+W1W2 fl.r., ${m}_{{\rm{i}}}\lt 21.5$	SDSS+2MASS+W1W2	1.00/0.69/0.82	0.90 /0.99 /0.94	1.00/1.00/1.00

Download table as:  ASCIITypeset image

As expected, it is evident from Table 8 that the inclusion of more photometric features always leads to better classification results at the expense of the training sample size.

The subsets with the same amount of features but different limits on the SDSS i-band show that the brighter samples (${m}_{{\rm{i}}}\lt 18.5$), the ones with the qualitatively better photometry, also show more accurate classifications. Again, the training sets for those three subsets are much smaller and will not be able to populate the entire feature space as well as a fainter limit of ${m}_{{\rm{i}}}\lt 21.5$ would allow.

Since the number of higher-redshift quasars is a strong function of the i-band magnitude, any limitation of the training set in this regard will bias the classification of those objects. For example, the training set with full SDSS+2MASS+W1W2 photometry and ${m}_{{\rm{i}}}\lt 18.5$ (${m}_{{\rm{i}}}\lt 21.5$) has only 38 (42) "highz" ($z\gt 3.5$) training objects and 8 (13) test objects. This demonstrates that these classification results cannot be fully trusted because of the low number of objects that the precision, recall, and F₁ score is based upon.

Based on these insights, we adopt the SDSS+W1W2 subset with ${m}_{{\rm{i}}}\lt 21.5$ as the best training and feature set for the quasar classification in our ELQS quasar selection (see the ⋆ in Table 8). It achieves the highest completeness (recall) of "highz" quasars and "QSO" in general of all subsets in Table 8 with the hyperparameters set to n_estimators = 300, min_samples_split = 3 and max_depth = 25.

For this particular training and feature set, we show the full confusion matrix in Figure 8. The rows correspond to the true class of the object, while the columns show the predicted labels. The values on the diagonal are the correctly classified objects. Each entry in the matrix shows the total number of objects and the percentage of the objects in that entry with respect to the true class (the full row). The entries are color-coded based on this percentage, with a darker blue color corresponding to a higher percentage.

Zoom In Zoom Out Reset image size

Although the stellar classification encounters difficulties for the A, F, and G stars, the later spectral types and the quasars are well classified with diagonal entries above 80%. Only a very small percentage of stars are classified as quasars (top-right corner), with the majority of stellar contaminants falling into the "midz" quasar class ($2.2\lt z\lt 3.5$). This is the redshift range in which the quasar distribution overlaps with the stellar locus the most in optical-color space.

Conversely, the bottom-left corner shows quasars classified as stellar sources. Again, the majority of these objects fall into the "midz" quasar class.

We can simplify this confusion matrix by grouping all stellar spectral classes and all quasar classes in the binary classes "QSO" and "STAR." The binary classification results are shown in Table 9. Here, the top part summarizes the results of the confusion matrix in Figure 8. In the bottom part, we used the same training set but reduced the feature set to only include the SDSS photometry. As a result, the off-diagonal values are much higher, reflecting a larger fraction of stellar contaminants and a lower quasar completeness. This demonstrates how the information gained by including the WISE W1 and W2 bands enhances the performance of the classifications.

Table 9.  Binary Label Classification Matrices of the SDSS+W1W2 Subset (${m}_{{\rm{i}}}\lt 21.5$) with All Features (Top) and Only the SDSS Features (Bottom)

SDSS+W1W2	pred. STAR	pred. QSO
STAR	48,015	309
QSO	61	29,979
SDSS	pred. STAR	pred. QSO
STAR	47,747	577
QSO	333	29,707

Download table as:  ASCIITypeset image

We demonstrated that random forests efficiently classify quasars and estimate their photometric redshifts using the SDSS optical bands in concert with WISE W1 and W2 photometry. In this context, we need to ask whether the JKW2 color cut is really necessary for our quasar candidate selection.

To evaluate these questions, we assume the RF classification contamination as shown in Table 9 and apply it to the spectroscopically unidentified sources in the $70\,{\deg }^{2}$ test regions (Table 2). In this case, the number of stellar contaminants would at least rise from the ∼1000 which make the color cut to ∼1700. This illustrates how important the information of the 2MASS J- and K-bands is to reject stars.

Unfortunately, only a small fraction of all known quasars and stars in the training sets have well-detected 2MASS J and K photometry. The resulting training samples are too small to properly populate the multidimensional feature space to allow for sufficient training of the RF model. In addition, the quasar training set would be strongly biased against higher-redshift ($z\gtrsim 3$) objects, which have less bright apparent magnitudes and are therefore less likely to be detected by 2MASS. Even though we are focusing on the brightest quasars, it is important for the redshift range above $z\approx 3$ to be well populated in the training set to achieve reliable results in the classification and regression.

As long as J- and K-band photometry are limited to the brightest objects in the training sets, the JKW2 color cut will play an important role in rejecting stellar sources.

The ELQS survey includes all SDSS photometry with galactic latitudes $b\lt -20$ or $b\gt 30$. To estimate the effective area of the ELQS survey, we use the Hierarchical Equal Area isoLatitude Pixelization (HEALPix; Górski et al. 2005). The process of the calculation and the general parameters used are identical to the description in Jiang et al. (2016).

The effective area of the full ELQS survey is ${\rm{11,838.5}}\,\pm 20.1\,{\deg }^{2}$, with a contribution of $7601.2\,\pm 7.2\,{\deg }^{2}$ from the spring ($90\,\deg \lt {\rm{R.A.}}\lt 270\,\deg$) sky and a contribution of $4,237.3\pm 12.9\,{\deg }^{2}$ from the fall (${\rm{R.A.}}\gt 270\,\deg$ and ${\rm{R.A.}}\lt 90\,\deg$) sky.

With all tools at hand, we begin the construction of the ELQS quasar candidate catalog in the SDSS footprint. An overview is given in Figure 3. The general source selection is based on the WISE AllWISE catalog matched with photometry from the 2MASS PSC. Both surveys are all-sky and therefore include the galactic plane, where the source density of stars is extremely high and leads to confusion in the near- and infrared surveys. Therefore, we restrict the quasar candidates to larger galactic latitudes and only include sources with either $b\lt -20$ or $b\gt 30$. We also require the WISE W1 and W2 photometry to have ${\rm{S}}/{\rm{N}}\gt 5$ and positive J-band magnitudes for all objects. All sources that pass these criteria and obey the JKW2 color cut $K-W2\geqslant 1.8-0.848\cdot (J-K)$ are selected in our WISE–2MASS–Allsky catalog. It comprises a total of 3,376,354 sources.

We then proceed to match all of these sources to the SDSS DR13 (PhotoPrimary) catalog in a $3\buildrel{\prime\prime}\over{.} 96$ aperture. We do not reject any objects based on their photometric flags. None of the fatal or non-fatal flags of Richards et al. (2002) are evaluated in our selection. We instead inspect the images of all quasar candidates (∼400 objects) in the very end to be as complete as possible. The matched SDSS–WISE–2MASS catalog has a total of 1,690,813 objects.

In the next step, we apply the criterion on the Petrosian radius (petroRad_i ${\leqslant }\,2.0;$ see Section 4.3) to reject the majority of galaxies. We also require all sources to satisfy the i-band magnitude cut of ${m}_{{\rm{i}}}\lt 18.5$.

For the remaining candidates, we calculate photometric redshifts and evaluate their quasar probabilities using RF regression and classification (see Section 5). This demands that all sources have quantified photometric errors for all SDSS and WISE W1 and W2 photometry. The classification calculates the most probable class of objects (rm_emp_mult_class_pred), the general quasar or star class (rf_emp_bin_class_pred), and the total probability of the object to belong to the quasar class (rf_emp_qso_prob). The regression delivers the best estimate for the photometric redshift (rf_emp_photoz).

With this information at hand, all objects that obey the photometric redshift cut of ${z}_{\mathrm{reg}}\gt 2.5$ and are generally classified as quasars (rf_emp_bin_class_pred = QSO) or fall into the high-redshift inclusion boxes defined in Richards et al. (2002) are considered quasar candidates.

In addition to these primary candidates, we allow all objects that pass the photometric redshift cut and also have a probability of $\gt 30 \%$ to belong to the quasar class (rf_emp_qso_prob $\geqslant 0.3$) to be included as additional candidates. This results in a candidate catalog of 2253/1735 total/primary sources, out of which 920/876 are known quasars from the literature and 1333/859 are valid candidates for spectroscopic follow-up.

From this sample, we prioritize bright high-redshift objects with ${z}_{\mathrm{reg}}\geqslant 2.8$ and ${m}_{{\rm{i}}}\leqslant 18.0$, which will make up the ELQS spectroscopic survey. These criteria leave 742/594 total/primary candidates out of which 341/327 are known and 401/267 need to be followed up.

In the last stages, every candidate's photometry will be inspected before observation to check for image defects or obviously extended sources. Of the total/primary sample, roughly 59%/69% are good candidates, 10%/8% are blended with other sources in WISE, 9%/6% are extended, and 22%/17% have bad image quality in at least one of the bands or are contaminated by bright sources close by.

The final ELQS quasar candidate catalog includes a total of 237 targets, out of which 184 are primary candidates.