Where’s Swimmy?: Mining unique color features buried in galaxies by deep anomaly detection using Subaru Hyper Suprime-Cam data

This study employs multi-band ($grizy$) imaging data from the S20A internal release of the HSC-SSP (wide layer), which is presently available as the third public data release (HSC-SSP PDR3; Aihara et al. (2021)). All the images were reduced using a dedicated pipeline, hscPipe (Bosch et al., 2018). The HSC-SSP survey field is divided into $\sim 1.7\times 1.7$ square degree areas called tracts and further subdivided into $\sim 12\times 12$ square arcmin areas called patches. We excluded the patches with full-width-half-maximum (FWHM) larger than 1.0 arcsec or shallow imaging depth of the bottom 5th percentiles in PSF limiting magnitude in one or more filters. The thresholds of the $5\sigma$ PSF limiting magnitude were $g<25.82$, $r<25.49$, $i<25.23$, $z<24.47$, or $y<23.59$. After removing the inferior-quality patches, we obtained images of 815 square degrees.

We restrict our analysis to 49 319 luminous galaxies in the above area that were brighter than 20 mag in the $i$-band cmodel magnitude (composite model; Abazajian et al. (2004); Bosch et al. (2018)) and within a redshift range of $z=0.05$–0.2 (figure 1), where spectroscopic redshifts are obtained by looking for counterparts in the SDSS DR15 (Blanton et al., 2017; Aguado et al., 2019; Bolton et al., 2012) and the GAMA DR3 (Baldry et al., 2018; Hopkins et al., 2013; Baldry et al., 2014) with a search radius of 1 arcsec (Aihara et al., 2021). In addition, we removed sources that were prominently affected by nearby stars, cosmic rays, corrupted pixels, or saturated pixels by referring to various flags (is False) stored in the database: pixelflags_edge, pixelflags_interpolatedcenter, pixelflags_saturatedcenter, pixelflags_crcenter, pixelflags_bad, sdssshape_flag, pixelflags_bright_objectcenter, mask_brightstar_halo, mask_brightstar_ghost, and mask_brightstar_blooming (Coupon et al., 2018; Aihara et al., 2021; Bosch et al., 2018). Moreover, the point sources were excluded by imposing psfflux_mag$-$cmodel_mag$>$0.2 in the $i$ band (see also Strauss et al. (2002); Baldry et al. (2010)).

Subsequently, we categorized the sample into several classes to train and test our selection model. First, we cross-matched the sample with the MPA-JHU catalog (Brinchmann et al., 2004; Kauffmann et al., 2003; Tremonti et al., 2004), which summarizes detailed spectral measurements for the SDSS DR7 spectra (Abazajian et al., 2009) and is useful for constructing training data and investigating the physical properties of anomaly candidates. In the anomaly selection, it is essential to exclude uncommon sources from the training data as many as possible to ensure that the proposed model does not learn outlier features from the training sample. Thus, we established a data-set based on the MPA-JHU catalog, excluding known XELGs, quasars, and AGN candidates. Herein, the XELG galaxies were defined with the observed equivalent width (EW) of either [Oii]$\lambda\lambda$3726,3729, H$\beta$+[Oiii]$\lambda\lambda$4959,5007, or H$\alpha$+[Nii]$\lambda\lambda$6548,6584, being larger than 100 Å, and 701 sources satisfied this definition. In addition, we detected 23 sources that have a counterpart in the latest quasar catalog (DR16Q v4: Lyke et al. (2020)). The XELG and quasar samples were used to test the performance of the proposed anomaly selection model, as discussed in subsection 4.1. In addition, possible AGN and low-ionization nuclear emission-line regions (LINER) candidates were removed from the training data based on the flag (BPTCLASS $\geq 3$) in the SDSS library, which classifies galaxy spectra using the Baldwin, Phillips, & Terlevich (BPT) diagram (Baldwin et al., 1981; Kauffmann et al., 2003). Based on these selection criteria, we established the training data that amounted to $N=$ 21 151 (table 1). Again the training data should ideally not contain anomalies; however, the training data may have undetected anomaly objects or unknown unknowns because these selection criteria depended on only the SDSS 1-D spectral data. Additionally, although this study does not delve into the remaining 23 582 spec-$z$ sources, they were still adopted for use in anomaly detection as well as the MPA-JHU sample.

Moreover, we constructed color-selected samples without spectroscopic counterparts for a complementary analysis, as detailed in appendix C.

For each galaxy, we retrieved an image of $64\times 64$ pixels ($11\times 11$ arcsec², with a pixel scale of 0.168 arcsec per pixel: the default size in the HSC-SSP) for each of the five $grizy$ broad-bands. Thereafter, we scaled the pixel values to a range of [0, 1] for the proposed deep-learning algorithm. There are various ways of flux normalization such as logarithmic, power-law, or square-root stretches. Although the appropriate selection of a normalization formula for deep learning of galaxy imaging data is still debatable, this research followed the approach by Lupton et al. (1999, 2004) that adopts an inverse hyperbolic sine function (arcsinh). Moreover, the arcsinh stretch has been commonly applied to galaxy images in public engagement programs such as Galaxy Zoo¹¹1https://www.zooniverse.org/projects/zookeeper/galaxy-zoo/ (Lintott et al., 2008, 2011; Willett et al., 2013) and GALAXY CRUISE²²2https://galaxycruise.mtk.nao.ac.jp. More importantly, this technique can avoid color saturation as the pixel values were normalized across color channels, and thus allow the model to interpret local anomalous color features within the galaxy image.

Furthermore, we slightly modified the code for arcsinh normalization written by the HSC-SSP pipeline team³³3https://hsc-gitlab.mtk.nao.ac.jp/ssp-software/data-access-tools/tree/master/pdr3/colorPostage to create multicolored imagine cubes for the hscMap⁴⁴4https://hsc-release.mtk.nao.ac.jp/hscMap-pdr3/app/, which is accessible on the HSC-SSP database and citizen science program, GALAXY CRUISE. This algorithm was compiled to apply the arcsinh stretch to five broadband images acquired from the HSC-SSP, $x\in A=\{g,r,i,z^{\ast},y^{\ast}\}$, where $x$ is a cut-out image in each band with a zero-point magnitude of 19 mag. Nevertheless, we scaled the $z$ and the $y$ band data by 1/1.2 and 1/1.5, respectively, corresponding to the typical flux to the $i$ band flux ratios. Such additional scaling aided in improving the balance of normalized fluxes between the broadbands. Otherwise, the proposed anomaly detection model will assign a higher weight to the $z$ and the $y$ band data, which are shallower than the $i$ band data (Aihara et al., 2021). The normalized flux $f(x)$ can be calculated as

where $I\equiv{\displaystyle\Sigma_{x\in A}}/5$, and the original values were used in the hscMap for the scaling parameter, $\alpha=\exp(10)$, and the bias parameter, $\beta=0.05$. Moreover, the flux range was limited to [0, 1], defined as follows:

The conducted scaling process and color visualizations based on the three bands of five broadband images are illustrated in figure 2. The standardization using the mean flux over the entire broadband ($I$) allows us to avoid color saturation (Lupton et al., 1999). Note that the $izy$ bands trace the stellar continua of targets at the rest-frame $\geq 6000$ Å, and the $z$ and the $y$ band fluxes are scaled to the $i$ band flux across the board, which causes an unclear color gradient in these three images. For instance, in case there is a strong H$\alpha$ line contribution to the $i$ band flux, a significant color excess would be detected in the $i$ band (section 4). Therefore, all the images were smoothed with a Gaussian kernel prior to flux scaling to match the observed FWHM to 1.0 arcsec, based on the observed information stored in each patch (subsection 2.1).

The denoising convolutional autoencoder (DCAE), which was applied in this study to detect anomalous galaxies, is an autoencoder (AE; Hinton and Salakhutdinov (2006)) equipped with a layer that add noise (Vincent et al., 2008) and Convolutional Neural Networks (CNN; Lecun et al. (1998); Krizhevsky et al. (2012); Lecun et al. (2015)). The AE is a machine-learning model that learns a latent representation of the input data $x$ using two Neural Networks (NN): an encoder and a decoder. The encoder $f\left(x;\bm{\theta}\right)$ with $\bm{\theta}$ as a parameter vector reduces the dimension of the input image $x$, and represents $x$ in a low-dimensional data $z=f\left(x;\bm{\theta}\right)$ in the latent space. Subsequently, the decoder $g\left(z;\bm{\phi}\right)$with $\bm{\phi}$ as a parameter vector reconstructs the input data $x^{\prime}=g\left(z;\bm{\phi}\right)$ from the encoded data $z$. Thereafter, the two NNs are optimized to minimize the value of the loss function $L\left(x,g\left(f\left(x\right)\right)\right)$, which is a function that evaluates the success of the learning process using the original data $x$ and reproduction data $g\left(f\left(x\right)\right)$. Upon considering the mean–squared error (MSE) as the loss function, the optimization relation can be expressed as

Prior to encoding, we added noise to the input data using a GaussianNoise layer that added Gaussian noise with a mean value of 0. The standard deviation of the Gaussian noise, $r_{\rm gauss}$, is a hyperparameter, and we adopted $r_{\rm gauss}=0.02$ (i.e., Gaussian noise with a standard deviation of 2 percent to the max count) through a trial-and-error process as described in appendix A.

Here, we briefly explain how our DCAE model can serve as a platform for anomaly detection. They trained a reproduction model with enormous data by updating the model weight parameters to minimize the loss function (see e.g., Chen et al. (2018); Zimmerer et al. (2018); Baur et al. (2021)). Therefore, the trained model learns the overall structure representation of the training data, but not the uncommon features existing in them. Upon inputting a galaxy with anomalous features to the trained model, the model is unable to appropriately reconstruct those features and produce an image without them. Consequently, the anomaly candidates can be selected by searching objects that manifest large residuals between the original and reconstructed images.

In this study, we adopted TensorFlow (version 2.5.0; Abadi et al. (2016)) and Keras (version 2.5.0; Chollet et al. (2015)) to construct the DCAE model. In the training process, we leveraged the Google Colaboratoy (Bisong, 2019) to expedite the calculations with a GPU environment. The model architecture is illustrated in figure 3. The input images were five-color ($grizy$) images, each containing $64\times 64$ pixels, as described in (section 2). The encoder contained convolutional layers (Conv2D), batch normalization layers (BatchNormalization; Ioffe and Szegedy (2015)), activation function layers (activation), and max-pooling layers (MaxPooling2D). Convolutional layers spatially slide $3\times 3$ filters over the images and produce feature maps, and then max-pooling layers take the max values in $2\times 2$ windows. Batch normalization layers normalize the output of prior units (Conv2D in this study) by every batch, which increases training efficiency and prevents overfitting (Ioffe and Szegedy, 2015). Activation function layers convert outputs of prior units using an activation function that is usually nonlinear and enables complex representations throughout the NN (CNN). The following rectified linear unit (ReLU; Glorot et al. (2011)) was adopted as the nonlinear activation function except the last layer using the sigmoid function:

Through Conv2D and MaxPooling2D layers in the Encoder, the input data dimensions are reduced to $8\times 8\times d$ in the latent space (figures 3 and 4), wherein $d$ denotes the number of channels in the final Conv2D layer within the encoder. Thus, a smaller $d$ value provides stronger constraints on the latent representation. We selected $d=8$ as a trial-and-error process (appendix A).

On the other hand, the decoder consisted of Conv2D, BatchNormalization, Activation, and upsampling layers (UpSampling2D). Upsamling layers simply resize the input images by a scale of two. In addition, the ReLU was adopted as the activation function in the activation layers, except for the final layer in which a sigmoid function was applied:

The filter numbers in each Conv2D layer were 8, 32, 64, and 5, respectively (figure 4).

The deep-learning model was trained using the training data prepared in section refs2, and the Adam optimizer (Kingma and Ba, 2014) was employed to adjust the weight parameters of the model and minimize the loss function value. The optimizer updates the weight parameters using the backpropagation algorithm (Rumelhart, Hinton, & Williams, 1986). A learning rate is a scalar that determines how much weight parameters are updated in each step. We trained the model with a learning rate of $l_{r}=0.001$ until a learning epoch of 10, after which it was altered to $l_{r}=0.0002$ (referred to as the learning decay technique in the machine-learning training). Figure 4 explains how imaging data of spiral and early-type galaxies were encoded and decoded by our reproduction model. Additionally, the examples of reconstructed images and residuals in the five bands for normal galaxies from the training data, DR16Q quasars, and XELGs (as defined in subsection 2.2) are illustrated in figure 6. Although the proposed model appropriately reconstructed the normal galaxies, we observed significant residuals between the original and reconstructed images for most of the quasars and XELGs. This is because the trained model was not optimized for uncommon color features originating from sources such as nuclear UV radiation in quasars and strong emission lines, as seen in XELGs.

To select objects with large residuals, such as quasars and XELGs, we defined the anomaly score $S_{\rm anom}$ as:

where $A^{\prime}$ means a calculation area described below. Thus, $S_{\rm anom}$ represents the MSE standardized by the square of the pixel counts of the reconstructed image, which gives the magnitude of the deviation from the model prediction. In particular, the area $A^{\prime}$ over which $S_{\rm anom}$ is calculated can be varied according to the target type of anomaly detection. Because this study aimed to detect anomalies that have uncommon features in the nuclear regions such as quasars and some XELGs, we only used the central circular area with a diameter of 6 pixels as the calculation area $A^{\prime}$, which broadly corresponds to the seeing FWHM (1 arcsec; subsection 2.3). Note that we may not be able to detect XELGs with their emission-line dominant components in the outer regions based on this definition. The most suitable weight data was defined to enable the selection of maximum number of DR16Q quasars and XELGs based on a model with high selection completeness. The model selection procedure is detailed in appendix A. Sigma excess values of $S_{\rm anom}$ for quasars and XELGs are provided in figure 6. The significant improvements in the anomaly scores in most of the DR16Q quasars and XELGs were compared to those of the normal galaxies derived from the training sample.

From the entire sample, we selected 5955 anomaly candidates by imposing the criteria: $S_{\rm anom}$ is in the top-5th percentile corresponding to the $2\sigma$ or higher excess (refer to figure 5), and the residual flux in the central area is positive in at least one of the $g$, $r$, and $i$ bands. Their identification numbers (object_id), sky coordinates in the HSC-SSP PDR3 (Aihara et al., 2021), and anomaly scores are summarised in table 3. Furthermore, examples of the selected anomaly candidates in each filter are represented in figure 7. The most anomalous galaxies in the current sample are discussed in subsection 4.4.

We determined that 61% of the DR16Q quasars $\left(=14/23\right)$ and 63.2% of the XELGs $\left(=443/701\right)$ satisfy the above-mentioned criteria. In contrast, there exist sources that have sufficiently high anomaly scores but are excluded from selection owing to their negative residual values, e.g., the numbers of those galaxies are 8, 322, and 14 for the $g$, $r$, and $i$ bands, respectively. These objects are further discussed in appendix B. Notably, as we evaluate the anomalies using only the central regions (equation 6), the chance projections of the background galaxies such as red galaxies will not affect the selection till they are not present in the central regions.

Our method selected 61% of the DR16Q quasars (= 14/23) and 63.2% of the XELGs (= 443/701). As these fractions are much higher than the anomaly fraction over the entire sample, 12.07% (= 5955/49319), the proposed method can preferentially select quasars and XELGs. To evaluate the performance of our anomaly detection method, we compared known quasars and XELGs selected as anomaly sources by the proposed method with the unselected known sources. We focused on radial profiles for quasars and emission-line EWs for XELGs because they are, respectively, the vital features characterizing these two populations. This comparison will assist us in interpreting the kind of features that can be identified by the proposed machine-learning model.

First, the quasars were examined. For ease of explanation, we focused on the quasars that are selected and not selected in a specific band, the $g$ band ($N=7$ and $16$, respectively; see table 2). The original and residual radial profiles for the quasars selected and not selected in the $g$ band are presented in figure 8. As compared to the unselected objects, the selected objects exhibit significant nuclear excesses in the residual image. Accordingly, we confirmed that the model can select quasars in cases of sufficiently large $g$ band excesses present in the center. As the model was not optimized for the unique color features observed in type-I AGNs and quasars, model reconstructions were generated without the blue (the $g$ band) nuclear components, as identified in the residual profiles (figure 8). However, the proposed model could accurately reproduce the $g$ band images of the unselected quasars in the $g$ band, which exhibited inadequate $g$ band excesses, regardless of the clear nuclear components in the $i$ band as well as the selected quasars. Note that the quasars and AGNs could be variable objects, and the SDSS spectra were obtained a decade ago. Therefore, certain known quasars could be inactive during the operation of the HSC-SSP (March 2014–January 2020 for HSC-SSP PDR3; Aihara et al. (2021)).

Furthermore, the XELGs were examined by focusing on the EWs of strong emission lines, as their intense emission-line contributions toward the selection band would yield high anomaly scores (see also subsection 4.3 about the redshift dependence of the anomaly selection rate). In addition, the parameters $EW_{g}$, $EW_{r}$, and $EW_{i}$ ($EW_{\rm band}$) were defined to quantify the emission-line contributions to the $g$, $r$, and $i$ bands, respectively. Herein, $EW_{\rm band}$ denotes the sum of the observed EWs of the strong emission lines in each band. In this calculation, we considered the actual throughput of the corresponding HSC filters (Aihara et al., 2021). Moreover, depending on the spectroscopic redshift of each source, the following emission lines provided in the MPA-JHU DR8 catalog (Brinchmann et al., 2004; Kauffmann et al., 2003; Tremonti et al., 2004) were considered: [Oii]$\lambda\lambda$3726,3729, [Neiii]$\lambda$3869, H$\beta$, [Oiii]$\lambda\lambda$4959,5007, H$\alpha$, [Nii]$\lambda$6548,6584, and [Sii]$\lambda\lambda$6717,6731.

The cumulative distributions of $EW_{\rm band}$ for the selected and unselected XELGs are presented in figure 9. The selection rates of objects with $EW_{\rm band}>100$ Å was 90.7% for the $g$ band, 84.4% for the $r$ band, and 64.4% for the $i$ band. Thus, the selection method can appropriately select the XELGs, given that they have prominent emission-line contributions to the target broadband images. More importantly, the high selection rate in the $g$ band could be caused by the $g$ band excess, which is a relatively uncommon trend in the training data.

Based on these results, we concluded that the proposed model tends to select known rare objects with more anomalous features, specifically, quasars with a larger excess in the residual image and XELGs with intense emission lines. Note that the proposed model is not an identifier only for quasars and XELGs, and thus we do not discuss selection purity of these sources in this paper.

In this subsection, the performance of the proposed model is determined based on the AGNs and quasars that were not used in the model training or selection process. For this purpose, we cross-matched our anomaly source sample with the Milliquas (Million Quasars) Catalog v7.2 (Flesch, 2021) that contains all the public AGNs and quasars as on April 30, 2021. Because we used the DR16Q sample in the model selection process, we excluded the DR16Q sample from the sample before matching with the Milliquas catalog. AGNs and quasars in the Milliquas catalog matched with the test data are archived from the SDSS DR16 (Ahumada et al., 2020), the low-redshift AGN catalog by Liu et al. (2019), the 2dF galaxy survey (Colless et al., 2001), the 2dF QSO Redshift Survey (Croom et al., 2004), the Principal Galaxy Catalog (Paturel et al., 2003), the Large Sky Area Multi-Object Fiber Spectroscopic Telescope DR3/DR2 (Dong et al., 2018), the WISE AGN catalog (Secrest et al., 2015), and other papers (Ge et al., 2012; Stalin et al., 2010; Francis et al., 2004; Croom et al., 2009; Lacy et al, 2007; Zhou et al., 2006; Gosset et al., 1997; Appenzeller et al., 1998; He and Chen, 1993; Markelijn et al., 1970; Bauer et al., 2000; Massaro et al., 2009; MacAlpine et al., 1981; Trichas et al., 2012; Trump et al., 2007; Lacy et al., 2013; Akiyama et al., 2015; Esquej et al., 2013; Garilli et al., 2014; Rowan-Robinson et al., 2013; Foltz et al., 1989; Chang et al., 2019).

The number of cross-matched sources within 1 arcsec apertures is summarized in table. 4, where object categories follow the object types given in the Milliquas catalog: an optically core dominated object is labeled as “quasar” and a host dominated object is labeled as “AGN”. In addition, Flesch (2019) claimed that “Type-$\rm II$ AGNs” may be contaminated by objects such as normal star-forming galaxies, LINERs, or starburst galaxies.

We determined that the proposed model selected 70% $\left(=16/23\right)$ of the quasars and 45.5% $\left(=130/286\right)$ of type-$\rm\,I\,$ AGNs, which are much higher than the anomaly fraction over the entire sample, 12.07% $\left(=5955/49319\right)$. Thus, these significant fractions suggested that the proposed method preferentially selected the quasar and type-$\rm\,I\,$ AGN samples. On the contrary, the selection rate of type-$\rm II$ AGNs was 12.8% $\left(=74/580\right)$, which was as low as the anomaly fraction. Thus, the proposed method is not advantageous in searching type-$\rm II$ AGNs. The lower selection rate of type-$\rm II$ AGNs was potentially caused by their smaller impacts on host galaxies in comparison to type-$\rm II$ AGNs in the optical regime. Thus, we concluded that an anomaly detection method, once optimized with appropriate test data, would be applicable for a blind search of quasars and type-$\rm\,I\,$ AGNs.

Next, the bandpass dependence of the selection rate is discussed. For quasars, type-$\rm\,I\,$ AGNs, and type-$\rm II$ AGNs, the selections based on the $g$ and $i$ bands yield three to five times higher selection rates than those based on the $r$ band (table 4). Such high selection rates in the $g$ band would be attributed to luminous UV radiation from the accretion disk, while those in the $i$ band would be caused by strong H$\alpha$ and [Nii] emission lines existing in this band. For the remaining AGN class (“others”), the selection rate was similar among the three bands. In contrast, the relatively lower selection rates in the $r$ band could be potentially because of their low $EW_{r}$ ($EW_{r}<20$ Å). Notably, the Milliquas catalog is incomplete, and therefore, the above-mentioned AGN fractions may increase upon considering the quasars or AGNs missing in the previous surveys.

Until the previous subsection, we determined that only 14 (from the DR16Q sample) and 234 (from the Milliquas catalog) sources ($\sim 4\%$) among the 5955 anomaly candidates bear a counterpart in the currently available quasar and AGN catalogs. We performed complementary analyses to examine the characteristics of all the candidates using literature data in this subsection. This characterization will provide valuable insights toward future improvements of the proposed anomaly detection model.

Initially, we examined the redshift dependence of the selection rate shown in figure 10. The selection rate in a given band is defined as the ratio of the number of selected galaxies in that band, $N_{\rm sel}$, to the number of all galaxies in the test data, $N_{\rm all}$. As can be observed from figure 10a, the selection rate in the $g$ band decreases with the redshift. This trend can be explained based on the behavior of the throughput of [Oiii] and [Oii] emission lines with the redshift, i.e., the stronger [Oiii] line is captured by the $g$ band only at $z<0.09$, whereas the weaker [Oii] line is within the band in most of the target redshift range. Such a behavior can also explain an increase in the $r$ band selection rate at $z>0.09$. Similarly, a spike in the $r$ band selection rate at the lowest redshift and a relatively weak redshift bin dependence in the $i$ band selection rate are caused by the H$\alpha$ line. These results indicated that a significant fraction of the anomaly candidates were selected because of their strong emission line(s) and the selection rate in a given band is strongly dependent on what lines are redshifted into that band.

The influences of $EW_{\rm band}$ on the selection process were further investigated by directly comparing the $EW_{\rm band}$ with the anomaly score $S_{\rm anom}$ for 25 714 galaxies that have EW measurements in the MPA-JHU catalog in figure 11, where the XELG sample is included as well. Again, we found that galaxies with high $EW_{\rm band}$ tend to exhibit high $S_{\rm anom}$, especially at $EW_{\rm band}>50$ Å, which indicated that the proposed model preferentially selected the extreme emission-line objects with high anomaly scores. However, figure 11 depicts certain intriguing sources with high $S_{\rm anom}$ but low $EW_{\rm band}$ ($EW_{\rm band}<1$ Å). More than half of these outlier objects apparently exhibited artificial-like features, as seen in appendix B, although we cannot expunge the possibility that they have intense line-emission only at the center, which cannot be judged by one-dimensional SDSS spectra captured with the $3"$-aperture fibers owing to dilution. Moreover, figure 11 indicates certain galaxies with high $EW_{\rm band}$ ($EW_{\rm band}>100$ Å) but low $S_{\rm anom}$, especially in the $i$ band (see also figure 9). However, 70% of these galaxies were selected in other band(s) as well, which implied that their low $S_{\rm anom}$ values can be caused by the flux excess in other bands. For instance, XELGs with low $S_{\rm anom,i}$ at $z=$ 0.1–0.2 can be explained based on the intense [Oiii] emission in the $r$ band. In such a case, the model may confuse XELGs with red galaxies and thereafter, reproduce their $i$ band excesses better than expected.

Furthermore, the dependence of the anomaly score on apparent magnitude was examined. As $S_{\rm anom}$ denotes the square sum of the normalized residuals, it is expected not to be sensitively dependent on apparent magnitude, except for the faint sources bearing low $S/N$ ratios per pixel. The anomaly score $S_{\rm anom}$ is illustrated as a function of apparent magnitude for the entire sample in figure 12. The top 5% and 0.7% (correspond to $2\sigma$ and $3\sigma$) criteria of $S_{\rm anom}$ for the respective bands are depicted by solid and dashed lines, respectively. In addition, Spearman’s correlation test could not derive any correlation in the three bands, which assures that the proposed method can evaluate anomalies independently of the apparent magnitude in the sample. Note that the anomalies selected in the $r$ band tend to be biased toward fainter magnitudes at $m_{r}>16.5$ because the luminous objects with $m_{r}<16.5$ and large $S_{\rm anom}$ tend to exhibit negative residuals. These sources are further examined in appendix B.

Figure 13 indicates the distribution of the galaxies in the MPA-JHU catalog on the SFR–$M_{\rm star}$ plane. In figure 13, the upper and lower crowds of objects are known to be the star-forming (SF) main sequence and the passive sequence (Brinchmann et al., 2004), respectively. Note that the quasars and AGNs are not plotted because their SFR and $M_{\rm star}$ would exhibit large uncertainties (e.g., Toba et al. (2018); Toba et al (2021)). In all the bands, the proposed model preferentially selects the galaxies in or above the SF main sequence, as indicated by the adequately small $p$-values ($\ll 0.01$) in the Kolmogorov–Smirnov (KS) test. Thus, the result is broadly consistent with the trend that the current method preferentially selected the strong emission-line galaxies. We also confirmed a similar trend as in the GALEX-SDSS-WISE Legacy Catalog 2 (GSWLC-2: Salim et al. (2016); Salim, Boquien, & Lee (2018)), which offers more reliable SFR and $M_{\rm star}$ measurements than the MPA-JHU catalog because of its multi wavelength coverage from UV to mid-IR. However, significant fractions of anomaly candidates (e.g. 78% in the $r$ band selection) that do not appear in this case owing to the absence of counterparts and failures of the SED fit in the GSWLC-2 (see also subsection 4.4), and thus we decided not to use the measurements by the GSWLC-2 in this paper.

Ultimately, we examined the emission-line strengths, morphology, and star-formation activities of “the most anomalous sources”, which are defined as with the $S_{\rm anom}$ higher than the $3\sigma$ (top $0.3\%$). The dashed lines in figures 5, 11, and 12 indicate an excess of $3\sigma$. In addition, figures 14 and 15 depict the original multicolored images and the SDSS DR16 (Ahumada et al., 2020) spectra of the most anomalous sources, respectively.

First, the XELGs were determined as dominant in the most anomalous sources. In particular, 71%, 48%, and 65% of the objects over the 3$\sigma$ excess from the MPA-JHU catalog were XELGs in the $g$, $r$, and $i$ bands, respectively. These results indicated that the proposed anomaly selection method with the $3\sigma$ excess, especially in the $g$ and $i$ bands, facilitated the blind searches of XELGs. However, as certain objects over the 3$\sigma$ excess in the $r$ band exhibited a small $EW_{r}$ ($EW_{r}<1$ Å; figure 11), the XELG fraction in the $r$ band was relatively low.

As can be observed in figure 14, a number of the most anomalous sources exhibited merger-like features or close neighbors. In the former case, we confirmed that the residual images of such galaxies tend to have multiple components even in their residual images, suggesting that they are probably dual AGNs or quasars or have multi-emission-line dominant components. As observed in figure 14, some of the most anomalous sources selected in the $g$ band portraying a blue component, such as the blueberry galaxies, also known as XELGs (Yang et al., 2017). Similar trends can be observed in the $r$ or $i$ bands, as they selected galaxies with a green or red component, respectively. In addition, certain sources with a purple component were caused by the strong emission-line contributions from the [Oiii] (or [Oii]) and H$\alpha$ to the $g$ and $i$ bands, respectively. Such intense emission-line features can be confirmed in their optical spectra, as depicted in figure 15.

The colored diamond plots in figure 13 illustrate the relationship between the SFR and $M_{\rm star}$ for the most anomalous sources. In all the band, the majority of the most anomalous sources are located above the SF main sequence. Besides, the most anomalous sources tend to be less massive compared to the entire sample. Although some of the most anomalous sources in the $r$ band are located in the passive sequence, in appearance they are the artifact-like sources as discussed in appendix B.

We also comment about fractions of the invalid SFR–$M_{\rm star}$ measurements in the GSWLC-2 for the most anomalous sources (see subsection 4.3). We noticed that $\sim 16$% of the most anomalous sources cannot be appropriately fitted with the SED templates in the GSWLC-2 (Salim, Boquien, & Lee, 2018), which are remarkably higher than that of the complete test data, 0.67%. Such a high fraction suggests that the proposed method selected anomalies that exhibited too anomalous photometric colors to apply the SED fitting, owing to extreme emission lines or emission from SMBH, as observed in figure 15. Therefore, specialized templates are required to appropriately derive these physical parameters, which constitute as the in-depth analyses in future research.

Conclusively, the present anomaly candidate sample contained unique populations. Further investigations of such sources will help in the search for interesting objects such as dual quasars, XELGs, or galaxies with unusual merger-like features. Thus, we expect that the application of the proposed deep-anomaly detection method on unprocessed larger data-sets can result in the discovery of extremely unique objects, including unknown unknowns. This aspect was demonstrated by applying the current model on a magnitude-limited color-selected sample without spectroscopic counterparts, as described in appendix C.