A Machine Learning Based Morphological Classification of 14,245 Radio AGNs Selected from the Best–Heckman Sample

We have designed an autoencoder, namely MCRGNet (Morphological Classification of Radio Galaxy Network), which consists of an encoder block and a decoder block by convention (Yu & Deng 2014; Goodfellow et al. 2016), based on a deep CNN to learn the morphological features (i.e., the intrinsic representations) of radio AGN images and to precisely reconstruct the original input images by applying the learned image features.

Many of the state-of-the-art deep neural networks, e.g., AlexNet (Krizhevsky et al. 2012), VGG (Simonyan & Zisserman 2014), and GoogLeNet (Szegedy et al. 2015), have achieved good performance in image classification and object detection competitions (e.g., ILSVRC; Russakovsky et al. 2015). They have also been applied to different but related tasks through transfer learning (e.g., Pan & Yang 2010; Esteva et al. 2017). Note that although nearly all of the above networks have been trained on image sets drawn from ImageNet (Deng et al. 2009), which is composed of 14 million images of 1000 generic object classes spanning a wide range of subjects that cover many aspects of either human activity or life and natural phenomena on Earth (Esteva et al. 2017), no generic object class is directly related to celestial objects (Lukic et al. 2018). Therefore, negative transfer problems (Mahmud & Ray 2008; Pan & Yang 2010), e.g., dissimilar risks (see Rosenstein et al. 2005 for details), may be caused when these networks are directly applied to the task of radio AGN classification for transfer learning. In addition to this, all of these networks are composed of tens of layers, which contain millions of parameters; thus they must be run on workstations with multiple graphics processing units (GPUs) to be trained. In order to consume fewer device resources and enable relevant studies even on portable computers, we have decided to design a small-scale neural network instead of using the above large networks (e.g., Diamantis et al. 2018; Lukic et al. 2018).

The schematic architecture of the autoencoder is shown in Figure 2(a), which basically consists of several convolutional layers assigned with multiple 3 × 3 pixel kernels according to state-of-the-art neural network theories (e.g., AlexNet and VGG; Krizhevsky et al. 2012; Simonyan & Zisserman 2014; Ma et al. 2017). The encoder block functions via a feature extractor, which focuses on local image features from lower to higher orders through successive convolution operations performed in multidimensional (multichannel) convolutional layers and employs a classifier consisting of one fully connected layer to transform the derived feature maps into feature vectors to implement the classification task (see Figure 2(b)).

Zoom In Zoom Out Reset image size

The decoder block is a mirror of the encoder block. It is constructed by reversing the positions of the convolutional and fully connected layers in the encoder, and the layers in the decoder share the same structures of weights as their mirror counterparts in the encoder. The image reconstruction function of the decoder not only assists the unsupervised learning and training of the network parameters through analyzing the residuals of the model but also acts as a generator of radio AGN images when fed with proper feature vectors sampled from the probability distribution of the extracted features by the encoder. Similar architectures have been widely applied in a variety of areas, especially when simulating or generating data is important, such as in bottleneck training in speech recognition (e.g., Yu & Deng 2014), in the optimization of discriminator and generator in image simulation (e.g., Goodfellow et al. 2014), and in artistic-style learning (e.g., Radford et al. 2015).

Based on a series of tests, we embed five convolutional layers into the encoder block (Figure 2(a)) to hierarchically extract image features, so that the sizes of the input images reduce reasonably from the original 140 × 140 pixels in the first layer to 5 × 5 pixels in the fifth layer. After the fifth convolutional layer we use one fully connected layer, which is named the encode layer and is used to conduct the classifications, to further compress the derived features into a 64-dimension vector; including more fully connected layers will not significantly increase the performance of classification in this case (see Section 5.3 for a discussion). Immediately after each convolutional layer we apply a rectified linear unit to generate the activation function, which is nonsaturating and can help increase the convergence speed by truncating negative output into zero (Glorot et al. 2011). The number of kernels (i.e., the number of channels) of each convolutional layer and the corresponding kernel sizes are listed in Table 5.

Because the decoder block is designed as the reverse of the encoder block, we feed the output of the encoder (i.e., the extracted feature vector) into the decoder and compare the output of the latter, i.e., the reconstructed images, with the original images using the mean squared error to evaluate the performance of the autoencoder. Obviously, such an autoencoder can also be used to generate radio AGN images as a simulator by feeding the randomly obtained 64-length feature vectors sampled from the distributions of the extracted morphological features by the encoder into the decoder block (see, e.g., our recent work in Ma et al. 2018).

Given that the size of the labeled sample is relatively small, we apply an effective training strategy that integrates the "pre-training" with the "fine-tuning" process (Hinton et al. 2012; Girshick et al. 2014; Christodoulidis et al. 2016) by adopting the following steps (Figure 3): (1) Pre-training—Apply the self-taught autoencoder to the unLRG sample (14,245 unlabeled radio AGNs) to pre-train a large number of weight and bias parameters of the deep CNN model, in order to deeply learn the representations of the radio AGN images and to avoid as much as possible the overfitting problem usually associated with learning on a small sample. (2) Fine-tuning—Append a softmax layer (e.g., Hinton et al. 2012) to the end of the encoder that has been pre-trained in the first step, and use it to classify the LRG sample (1442 radio AGNs) with an emphasis on the training of the fully connected layer, which is used to carry out classification. (3) Classification—Apply the fine-tuned encoder to classify the unLRG sample.

Zoom In Zoom Out Reset image size

Our approach of training the neural network with the above three steps (PFC steps hereafter) has several advantages. First, when there is a lack of sufficiently large labeled samples (a sample containing ≳10⁴ labeled radio AGNs may be necessary in our case, given the level of image complexity; see a similar case in LeCun et al. 2015 and the theoretical discussion in Esteva et al. 2017) for the training of a complex neural network, a straightforward training with small labeled samples does not help lower the risk of model overfitting and poor generalization capability, i.e., low predictive precision on new images. For example, the neural network is likely to overreact to minor, unimportant fluctuations in the training data and learn inessential random features. This problem is overcome by applying the PFC steps. Second, our approach enables one to perform the classification of a sample of about 15,000 radio galaxies efficiently and accurately in several seconds even on a portable computer equipped with a single GPU of moderate memory (≲10 GB), not just on a workstation with multiple GPUs (≳100 GB), which other approaches usually require. This is because the resource consumption of the deep neural network in the encoder block of our autoencoder is significantly lower than that of other complex networks, e.g., GoogLeNet (Szegedy et al. 2015; Aniyan & Thorat 2017). This allows astronomers to inspect a newly obtained sample conveniently. In addition, our autoencoder can be used as a random-sample generator for morphology-based studies of radio galaxies.

3.3.1. Convolution

A convolution operation is usually a linear transformation that multiplies an input image with a kernel that slides across the input image to produce a matrix called an output feature map (Dumoulin & Visin 2016). For a CNN containing a set of L layers, each convolutional layer is assigned a set of kernels W^l_n and biases b^l to form ${N}_{\mathrm{ch}}^{l}$ channels, where $l=1,\,2,\,...,\,L$ and $n=1,\,2,\,...,\,{N}_{\mathrm{ch}}^{l}$ denote the ordinal numbers of the layers and channels, respectively. Thus for the nth channel in the lth layer of the network, the output feature map ${O}_{n}^{l}$ can be written as

$$\begin{eqnarray}{O}_{n}^{l}=\left\{\begin{array}{ll}{W}_{n}^{l}\ast {I}_{0}+{b}_{n}^{l}, & \mathrm{if}\,l=1,\\ \displaystyle \sum _{m=1}^{{N}_{\mathrm{ch}}^{l-1}}{W}_{n}^{l}\ast {I}_{m}^{l}+{b}_{n}^{l}, & \mathrm{if}\,l=2,\,3,\,...,\,L,\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where I₀ is the original image fed into the network, ${I}_{m}^{l}\equiv {O}_{m}^{l-1}$ (l ≥ 2) is the mth input image ($m=1,\,2,\,...,\,{N}_{\mathrm{ch}}^{l-1}$), the notation "*" is the convolution operator, and b^l_n is the bias.

Once an input is given, the corresponding output of a convolutional layer is uniquely determined by its kernels and biases and the choices of stride and zero-padding modes. For simplicity, in this work we always scale the images, as well as the kernels, as square matrices. The initial values assigned to the elements in the kernels are randomly generated with truncated normal Gaussian distribution, and corresponding biases are initialized with constants (Abadi et al. 2015). We select half-padding (i.e., same padding) $p=\lfloor k/2\rfloor$, where k is the size of the squared kernel, to keep the shape and set the size of the input. We also set the stride, the step length to slide the kernel, to s = 2, which leads to a linear pooling (i.e., subsampling or downsampling), to decrease the number of free parameters and meanwhile extract the most remarkable features.

3.3.2. Deconvolution (Transposed Convolution)

The decoder block of the autoencoder can be regarded as a mirror of the encoder block, which can be derived by reversing the inputs and outputs of the CNN described in Section 3.3.1 in turn. Similar to the convolution operation (Equation (1)), deconvolution can be expressed as

$$\begin{eqnarray}\left\{\begin{array}{ll}{O}_{n^{\prime} }^{l^{\prime} }=\displaystyle \sum _{m=1}^{{N}_{\mathrm{ch}}^{l^{\prime} }}{\left[{W}_{m}^{l^{\prime} }\right]}^{T}\,{* }^{T}{{I}}_{m}^{l^{\prime} }+{b}_{n^{\prime} }^{l^{\prime} }, & \mathrm{if}\,l^{\prime} =1,\,2,\,...,L-1,\\ {O}^{l^{\prime} }=\displaystyle \sum _{m=1}^{{N}_{\mathrm{ch}}^{l^{\prime} }}{\left[{W}_{m}^{l^{\prime} }\right]}^{T}\,{* }^{T}{{I}}_{m}^{l^{\prime} }+{b}^{l^{\prime} }, & \mathrm{if}\,l^{\prime} =L,\end{array}\right.\end{eqnarray} \tag{ 2 }$$

for the n'th kernel channel ($n^{\prime} =1,\,2,\,...,\,{N}_{\mathrm{ch}}^{L-l^{\prime} }$, and $n^{\prime} \equiv 1$ when l' = L) in the l'th layer, where ${I}_{m}^{1}$ ($m=1,\,2,\,...,\,{N}_{\mathrm{ch}}^{L}$) denotes the images of the first layer in the decoder block (i.e., the feature maps recovered from the Lth encode layer), ${I}_{m}^{l^{\prime} }\equiv {O}_{m}^{l^{\prime} -1}$ (l' ≥ 2) is the mth input image of layer l' ($m=1,2,\,\ldots ,\,{N}_{\mathrm{ch}}^{L-l^{\prime} }$), ${b}_{n^{\prime} }^{l^{\prime} }$ is the bias, and the superscript T denotes the transpose arithmetic (Dumoulin & Visin 2016). A detailed explanation of the transpose convolution is illustrated in Appendix A with examples.

In order to evaluate the performance of our autoencoder, we visualize the classification of the LRG sample, which is achieved by analyzing the multidimensional feature vectors, in two-dimensional space (the visualization space) by applying the t-SNE algorithm, a dimension reduction algorithm that can preserve the similarities between vector points (i.e., distances) when mapping data from high-dimensional space onto lower dimensions (e.g., Van Der Maaten & Hinton 2008; Esteva et al. 2017).

To be specific, we use a set of multidimensional vectors $\{{{\boldsymbol{x}}}_{1},\,{{\boldsymbol{x}}}_{2},\,...,\,{{\boldsymbol{x}}}_{{N}_{\mathrm{LRG}}}\}$, where N_LRG is the number of the LRG sample. These vectors form a distribution ${ \mathcal P }$, in which p_ij (i, j = 1, 2, ..., N_LRG) denotes the probability of similarity between any two vectors ${{\boldsymbol{x}}}_{i}$ and ${{\boldsymbol{x}}}_{j}$. If we further use ${p}_{j| i}$ to define the conditional probability of ${{\boldsymbol{x}}}_{j}$ when ${{\boldsymbol{x}}}_{i}$ is given, by assuming the Gaussian distribution, we have

$$\begin{eqnarray}&&{p}_{j| i}=\displaystyle \frac{\exp (-\parallel {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{j}{\parallel }^{2}/2{\sigma }_{i}^{2})}{{\displaystyle \sum }_{k\ne i}\exp (-\parallel {{\boldsymbol{x}}}_{i}-{{\boldsymbol{x}}}_{k}{\parallel }^{2}/2{\sigma }_{i}^{2})},\end{eqnarray} \tag{ 3 }$$

where ${\sigma }_{i}^{2}$ is the variance of ${{\boldsymbol{x}}}_{i}$ (see Hinton & Roweis 2002 for details), and

$$\begin{eqnarray}&&{p}_{{ij}}=\displaystyle \frac{{p}_{j| i}+{p}_{i| j}}{2{N}_{\mathrm{LRG}}}\end{eqnarray} \tag{ 4 }$$

by symmetry.

Similarly, in the visualization space we use a vector set $\{{{\boldsymbol{y}}}_{1},\,{{\boldsymbol{y}}}_{2},\,...,\,{{\boldsymbol{y}}}_{{N}_{\mathrm{LRG}}}\}$ to describe the distribution ${ \mathcal Q }$ for the classified radio AGNs and q_ij to represent the probabilities of similarities between any two vectors. As shown by Van Der Maaten & Hinton (2008), q_ij appropriately preserves the similarities between the vectors in high-dimensional space, i.e., p_ij. Unlike the case of the distribution ${ \mathcal P }$, for which Gaussian distribution is assumed, Van Der Maaten & Hinton (2008) employed a long-tailed Student's t-distribution with one degree of freedom to describe ${ \mathcal Q }$ to avoid abnormal points, which is written as

$$\begin{eqnarray}&&{q}_{{ij}}=\displaystyle \frac{{\left(1+\parallel {{\boldsymbol{y}}}_{i}-{{\boldsymbol{y}}}_{j}{\parallel }^{2}\right)}^{-1}}{{\displaystyle \sum }_{k\ne m}{\left(1+\parallel {{\boldsymbol{y}}}_{k}-{{\boldsymbol{y}}}_{m}{\parallel }^{2}\right)}^{-1}}.\end{eqnarray} \tag{ 5 }$$

By applying the Kullback–Leibler divergence (Kullback & Leibler 1951)

$$\begin{eqnarray}&&\mathrm{KL}({ \mathcal P }| | { \mathcal Q })=\displaystyle \sum _{i\ne j}{p}_{{ij}}\mathrm{log}\displaystyle \frac{{p}_{{ij}}}{{q}_{{ij}}}\end{eqnarray} \tag{ 6 }$$

as the objective function, we can evaluate the similarity between distributions ${ \mathcal P }$ and ${ \mathcal Q }$ and locate ${{\boldsymbol{y}}}_{i}$ by minimizing Equation (6) with respect to ${{\boldsymbol{y}}}_{i}$ with

$$\begin{eqnarray}&&\displaystyle \frac{\partial \mathrm{KL}}{\partial {{\boldsymbol{y}}}_{i}}=4\displaystyle \sum _{j}({p}_{{ij}}-{q}_{{ij}})({{\boldsymbol{y}}}_{i}-{{\boldsymbol{y}}}_{j})(1+{\left(\parallel {{\boldsymbol{y}}}_{i}-{{\boldsymbol{y}}}_{j}{\parallel }^{2}\right)}^{-1}.\end{eqnarray} \tag{ 7 }$$

We apply an approach that consists of the following three steps to prepare both the labeled and unlabeled images before we feed them into the autoencoder. First, we define a region of interest (ROI) in each image by cropping a box region with a size of 140 × 140 pixels (42 × 42) from the original 150 × 150 pixel (45 × 45) image. Next, in order to improve the contrast of the ROI of each image, we carry out sigma clipping to limit the image dynamic range by setting the parameter σ to 3 and the number of iterations to 50, and discard pixels dominated by background noise. Finally, in order to further avoid overfitting of the network, we perform an augmentation operation to increase the number of images in both the LRG and unLRG samples by applying the following operations to each image:

• 1.  
  Randomly flip the sigma-clipped image. The flipping mode is either left-to-right, top-to-bottom, or diagonal.
• 2.  
  Rotate the flipped image with an angle θ generated randomly between 0° and 360°.

For each unLRG image, the number of augmentation operations is 10, which results in a total of around 1.43 × 10⁵ images for the pre-training. For each LRG image, however, the number of augmentation operations is determined according to the morphological type of the radio AGN, in order to balance the numbers of different morphological types in the dichotomous classification tree. Given that on the classification tree (Figure 1) the two subbranches of any branch should have nearly the same AGN numbers, the ratios of the numbers of compact, FRI, FRII, BT, XRG, and RRG types are set to 8:2:2:2:1:1 after the augmentation operation, which results in about 3.9 × 10⁴ images in total (38,882, 19,554, 9731, 9823, 4900, and 4923 images for the binary classification branches compact–extended, typical–atypical FR, typical FRI–FRII, BT–irregular, FRI-like–FRII-like, and XRG–RRG, respectively; Table 6). One example of the flipping and augmentation operations is displayed in Figure 4.

Zoom In Zoom Out Reset image size

Table 6.  Augmentations of the LRG Sample

AGN Type	Label	a ${N}_{\mathrm{Trn}+\mathrm{Val}}$	b ${N}_{\mathrm{Tst}}$	c ${R}_{\mathrm{Aug}}$	d ${N}_{\mathrm{Aug}}$	eBranch
compact	1	302	75	64	19,328	1
FRI	2	169	42	29	4901	1, 2, 3
FRII	3	345	86	14	4830	1, 2, 3
BT	4	245	62	20	4900	1, 2, 4, 5
XRG	5	67	17	37	2479	1, 2, 4, 6
RRG	6	26	6	94	2444	1, 2, 4, 6
Total	⋯	1154	288	⋯	38,882	⋯

Notes.

^aNumber of sources in the training (Trn) + validation (Val) subsets. ^bNumber of sources in the test (Tst) subsets. ^cAugmentation rates for the corresponding radio galaxy types. ^dNumber of augmented images for the corresponding radio galaxy types. ^eBranches on the dichotomous tree: (1) compact–extended, (2) typical–atypical FR, (3) typical FRI–FRII, (4) BT–irregular, (5) FRI-like–FRII-like, and (6) XRG–RRG.

Download table as:  ASCIITypeset image

The augmented unLRG sample is then divided into independent training and validation subsets, which are used to optimize the parameters of the neural network and to evaluate the learning performance (i.e., the ability to reconstruct the images from which features are learned) to avoid the overfitting risk for each epoch, respectively. The training and validation subsets are divided at the ratio of 4:1; therefore the two subsets have about 1.14 × 10⁵ and 2.85 × 10⁴ augmented images, respectively.

In the original LRG sample, the numbers of fine-tuning images (i.e., training and validation images) and testing images are set to 1154 (80%) and 288 (20%), respectively. Data augmentation is performed only on the images used for fine-tuning, which are divided with a ratio of 4:1 to obtain 31,106 training images and 7776 validation images (Table 6). Note that the subset of testing images is used for evaluating the classification performance (i.e., the ability to classify radio AGNs correctly) of the trained network; therefore it contains only the original LRG images that are not augmented.

Following the strategy described in Section 3.2, we first feed the unLRG training and validation subsets into the autoencoder to pre-train the weight and bias parameters, and then input the LRG training, validation, and test subsets into the neural network, which is complemented with a softmax layer at the end of the encoder block, to fine-tune the network parameters with an emphasis on the fully connected encode layer.

4.2.1. Approach to Training and Evaluation

In this work we pre-train and fine-tune the neural network using the strategy of batch iteration (Abadi et al. 2015). Instead of using the whole training data set, we separate it into batches (100 images for each) and input them into the network batch by batch to adjust the network parameters iteratively. One round of training with all the batches is called an epoch.

For the pre-training stage the augmented images in the unLRG training or validation subset are fed into the encoder block in batches to generate 64-dimension feature vectors, which are successively fed as the input of the decoder block to reconstruct the images. We compare the reconstructed images with the corresponding original input images by calculating the mean squared residual error between them as the batch loss of pre-training, which is

$$\begin{eqnarray}&&{R}_{{\rm{P}}}=\displaystyle \frac{1}{{N}_{{\rm{B}}}}\displaystyle \sum _{i=1}^{{N}_{{\rm{B}}}}\displaystyle \sum _{j=1}^{{N}_{\mathrm{row}}}\displaystyle \sum _{k=1}^{{N}_{\mathrm{col}}}| {I}_{i,j,k}^{\mathrm{CAE}}-{O}_{i,j,k}^{\mathrm{CAE}}{| }^{2},\end{eqnarray} \tag{ 8 }$$

where N_B is the number of images in one training batch, N_row = N_col = 140 are the number of rows and columns of the images, I^CAE represents the input images of the encoder block, and O^CAE represents the reconstructed images of the decoder block. According to Goodfellow et al. (2016), the loss of each batch is back-propagated to the network along the gradient decreasing direction to adjust network parameters until the loss converges to a constant value. This is achieved after about 500 epochs of pre-training are completed.

For the fine-tuning stage, a softmax layer is appended to the end of the encoder block, and the batch loss of fine-tuning R_F is defined by averaging the cross-entropy of the augmented LRG images (Goodfellow et al. 2016) as

$$\begin{eqnarray}&&{R}_{{\rm{F}}}=-\displaystyle \frac{1}{{N}_{{\rm{B}}}}\displaystyle \sum _{i=1}^{{N}_{{\rm{B}}}}\displaystyle \sum _{j=1}^{{N}_{c}}{y}_{i,j}\cdot {\mathrm{log}}_{2}{p}_{i,j}^{S},\end{eqnarray} \tag{ 9 }$$

where N_c = 2 is the number of subbranches for each branch on the dichotomous classification tree and y_i,j (i = 1, 2, ..., N_B, and j = 1, 2, ..., N_c; y_i,j = 1 if the source belongs to type j, and y_i,j = 0 if otherwise) is the one-hot probability of the source. ${p}_{i,j}^{S}$ represents an output of the softmax layer, which is the normalized probability of the source being classified as type j in a certain batch. The batch loss converges to a constant after about 200 epochs of the fine-tuning are completed.

In the pre-training and fine-tuning stages we have employed the following standard rules.

• 1.  
  Dropout of weights with values approaching zero (Srivastava et al. 2014), which helps avoid overfitting, is performed each time after the calculation in a convolutional or a fully connected layer is completed. The dropout rate is fixed at 50% for each layer.
• 2.  
  In the pre-training stage the number of epochs is set to 500, which is determined by the convergence of training curves. A high initial learning rate of 0.0005 together with a decay rate of 0.95 to form the adaptive learning rate (Zeiler 2012; Abadi et al. 2015) and an adaptive moment optimization function (Kingma & Ba 2014) are both adopted. In the fine-tuning stage the number of epochs of each subclassifier is set to 200 also according to the convergence of the training curves, and a lower initial learning rate of 0.0001 is adopted to carefully fine-tune the parameters.
• 3.  
  In either the pre-training or the fine-tuning stage, the corresponding validation loss R_P or R_F, which is obtained by feeding the unLRG or LRG validation subset into the autoencoder or CNN, is used to evaluate the overfitting risks for each epoch. If the validation loss does not decrease or even increase while the training loss decreases steadily to convergence, the network parameters may be overfitted. Note that the validation loss is not used for parameter optimizations.

4.2.2. Results of Training and Evaluation

In Figures 5 and 6, we illustrate the training and validation loss curves obtained in the pre-training stage and fine-tuning stage, respectively. For comparison, in the fine-tuning stage we also plot the loss curves obtained when the pre-training process is switched off for comparison as well as the classification accuracy R_acc, which is defined as the ratio of the number of correctly classified sources to the actual source number (see also Equation (11)). To evaluate the robustness of the network, each training process has been conducted 10 times, and the means and 68% confidence limits of the losses and accuracies are calculated for every 25 epochs (see Figure 6). In either the pre-training stage or the fine-tuning of all the subclassifiers, both the training loss and the validation loss decrease monotonically, and the training loss becomes lower than the corresponding validation loss after about 10 epochs, and this remains unchanged until the pre-training or fine-tuning is completed. These indicate that the trained parameters in the neural network are not overfitted.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the fine-tuning stage, except for the typical–atypical FR branch, the average fine-tuning loss R_F decreases significantly faster, and the neural network achieves better classification accuracies after pre-training is applied. As for the typical–atypical FR branch, because (1) the sources therein are comparatively easier to classify due to their distinctive morphological features and (2) the number of images in the training subset is 19,554, the pre-training becomes more or less insignificant.

In order to further evaluate the performance of the neural network, we have employed three indices, i.e., sensitivity ${R}_{\mathrm{sen}}^{t}$, accuracy ${R}_{\mathrm{acc}}^{t}$, and precision ${R}_{\mathrm{pre}}^{t}$, to describe the classifier performance for the LRG samples (Khalaf et al. 2016; Ma et al. 2017). They are defined as follows:

$$\begin{eqnarray}&&\,{R}_{\mathrm{sen}}^{t}={N}_{\mathrm{TP}}^{t}/({N}_{\mathrm{TP}}^{t}+{N}_{\mathrm{FN}}^{t}),\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\,{R}_{\mathrm{acc}}^{t}=({N}_{\mathrm{TP}}^{t}+{N}_{\mathrm{TN}}^{t})/{N}_{\mathrm{LRG}},\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{R}_{\mathrm{pre}}^{t}={N}_{\mathrm{TP}}^{t}/({N}_{\mathrm{TP}}^{t}+{N}_{\mathrm{FP}}^{t}),\end{eqnarray} \tag{ 12 }$$

where the superscript t denotes the morphological type (t = 1, 2, ⋯, 6, corresponding to compact, typical FRI, typical FRII, BT, XRG, and RRG, respectively), N_LRG represents the total number of radio AGNs in the original LRG sample, ${N}_{\mathrm{TP}}^{t}$ (${N}_{\mathrm{TN}}^{t}$) is the number of type t (non–type t) radio AGNs that are correctly classified (T), and ${N}_{\mathrm{FN}}^{t}$ (${N}_{\mathrm{FP}}^{t}$) is the number of type t (non–type t) radio AGNs that are misclassified (F). We list the calculated indices in Table 7 with 68% error bars and find that our neural network can achieve a satisfactory classification. The classifications of compact, typical FRI, typical FRII, and BT types are very accurate, although the classifications of XRG and RRG types are relatively poor, possibly due to insufficient training/learning caused by their small sample scale.

Table 7.  Evaluation of the CNN Based on the Classification of the LRG Sample

AGN Type	Pre-trained			Non-pre-trained			Gain
	Rsen	Racc	Rpre	Rsen	Racc	Rpre	Rsen	Racc	Rpre
	(%)			(%)			(%)
compact	98.41 ± 0.25	98.96 ± 0.16	97.88 ± 0.63	98.14 ± 0.27	98.82 ± 0.17	97.12 ± 0.67	0.27 ± 0.17	0.14 ± 0.13	0.76 ± 0.58
FRI	88.72 ± 1.69	96.94 ± 0.48	95.39 ± 2.73	83.46 ± 2.13	95.42 ± 0.56	93.94 ± 2.74	5.26 ± 1.32	1.52 ± 0.57	1.45 ± 2.80
FRII	91.07 ± 2.12	95.66 ± 0.39	92.24 ± 1.01	80.77 ± 1.97	88.49 ± 0.46	88.41 ± 1.48	10.90 ± 1.52	7.27 ± 0.33	3.83 ± 1.27
BT	94.27 ± 1.89	95.25 ± 0.46	86.56 ± 0.92	85.57 ± 2.61	89.88 ± 0.52	76.57 ± 1.88	8.70 ± 1.32	5.37 ± 0.43	9.99 ± 1.50
XRG	75.05 ± 3.54	96.72 ± 0.34	72.34 ± 2.70	64.29 ± 4.78	91.19 ± 0.53	55.76 ± 4.27	10.78 ± 2.44	5.61 ± 0.58	16.58 ± 4.32
RRG	74.62 ± 3.64	97.04 ± 0.30	65.89 ± 3.65	70.42 ± 6.47	95.15 ± 0.28	60.21 ± 5.88	4.20 ± 4.52	1.89 ± 0.45	5.68 ± 7.66
Total	⋯	⋯	92.44 ± 1.94	⋯	⋯	84.67 ± 2.82	⋯	⋯	7.77 ± 3.02

Note. R_sen, R_acc, and R_pre are the evaluation indices, namely the sensitivity, accuracy, and precision, of the corresponding type.

Download table as:  ASCIITypeset image

The results obtained by the no pre-training strategy are also listed in Table 7 for comparison. The averaged precision gain between pre-training and no pre-training is 7.77% ± 3.02%, which indicates that pre-training indeed improves the network performance. Besides, we show the confusion matrices of classifying the LRG sample with and without pre-training strategies in Figure 7 for visualizing the difference between their performance. Though the gains in precision and sensitivity after applying pre-training on the AGN types with large source numbers (e.g., compact, FRI, etc.) might be relatively insignificant, they are significant for the XRG and RRG sources of small numbers, especially the gains in sensitivity. The same conclusion could also be drawn in the classification of the unLRG sample (Figure 8; see Section 4.3 for details).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We apply MCRGNet pre-trained by the proposed PFC strategy on the LRG sample and list the results, including the corresponding 68% errors, in Table 7, along with the results obtained without applying pre-training. The corresponding averaged gains of the indices obtained after carrying out pre-training for each subtype are also shown in the same table. Compared with the results of Aniyan & Thorat (2017), who trained their neural network with a relatively small labeled sample (716 sources), the classification precisions of both our typical FRI and FRII sources are higher (95% versus 91% and 92% versus 75%), although that of BT sources is lower (87% versus 95%). Aniyan & Thorat (2017) also reported that the recall (i.e., the index R_sen in this paper), which reflected the ability to avoid misclassification due to the incorrect inclusion of sources of other types, of the BT sources obtained with their neural network was 79%. They argued that this was because some sources in the BT sample were misclassified as FRIIs, which was supported by the high recall but low precision of FRII sources. As can be seen in Table 7, the sensitivities R_sen for the BT and FRII sources derived by our approach are both >91%. This, combined with the corresponding high precisions, suggests an improved classification of BT and typical FRII radio AGNs.

4.2.3. Feature Visualization of the LRG Sample

We show the results of the dimension-reduced feature visualization of the LRG sample classification with t-SNE in Figure 9, in which the subsamples of different types are marked with different colors. We find that the points of each type gather into a cluster that overlaps only slightly with those of other types, except for the case between FRIIs and XRGs and the case between compacts and FRIs. By carrying out visual inspection, we find that the FRII–XRG confusion is caused when the minor pair of jets in the XRG is dominated by the major pair, and the compact–FRI misclassification is caused when an FRI shows small-scale and weak jets.

Zoom In Zoom Out Reset image size

To show the merit of using a dichotomous classification tree structure (Figure 1), we have trained a CNN which shares the same structure as the encoder block illustrated in Figure 2(a) but has a six-node softmax layer as a classifier and used it to directly classify the LRG images into the six types. As shown in Figure 9(c) the result is unacceptable, because the points of the different types overlap with one another severely. Based on this we may conclude that, compared with direct six-type (non-tree) classification, the application of a dichotomous tree leads to significantly better performance.

In Figure 7 we show the 6 × 6 confusion matrices of the LRG classification obtained with and without pre-training and that obtained when direct six-type classification is performed. Compared to the dichotomous tree based methods, the direct training leads to the worst results, especially for XRG and RRG sources, whose numbers are relatively small. The same conclusion could also be found in Figure 8 on the classification of the unLRG sample (see Section 4.3 for details).

After deeply training the neural network with both unlabeled and labeled samples, we classify 14,245 radio AGNs in the unLRG sample and list the results in Tables 4 and 8. In addition to the CNN classification, for each source we have also estimated the joint probability of the classification by cascadingly multiplying the probabilities of the nodes decided and given by the softmax layer at the corresponding branches (Goodfellow et al. 2016; see Appendix B for details). Furthermore, if a source is classified as a BT, a subtype label (1 for FRI-like sources or 2 for FRII-like sources) is further presented, which is obtained by the CNN trained at the typical FRI–FRII branch.

In order to cross-check the results, we have also visually inspected all the unLRG images and classified the sample manually. We divide the images into three parts, and each part is inspected by at least three of the authors independently. If the three authors arrive at the same result, the visual classification is determined. Otherwise, an additional two authors will join, and the visual classification is determined when at least three authors agree on the result. The visual classification is also listed in Table 4.

We find that in general the CNN classification agrees with the human classification. About 80% of the sources have the same CNN and manual labels, and this ratio is even higher for compact sources (89%) and typical FRI sources (80%). About 3% of the small-scale FRIIs are misclassified by the neural network as FRIs, due to the small distances between the two lobes and their weak fluxes. About 7% of the typical FRs are misclassified by the neural network as BTs, XRGs, and RRGs, because these FRs tend to have ambiguous morphologies and low PAPRs (≤10.0 dB). In most cases, a low or relatively low (≲0.90) classification probability is found for the sources misclassified by the neural network. A further discussion on this issue will be presented in Section 5.2.

Fanaroff & Riley (1974) proposed a clear radio power transition between FRIs and FRIIs based on their study of the 3C (Mackay 1971) radio sources. Ledlow & Owen (1996) suggested that the transition was a function depending on optical luminosity and was most likely a result of the different lifetimes for sources with different radio powers. Subsequent studies (Best 2009—1083 FRIs and FRIIs; Wing & Blanton 2011—1228 FRIs and 683 FRIIs; Capetti et al. 2017a—122 FRIIs), however, do not support the claim that FRIs and FRIIs reside in distinctly different regions in the radio luminosity–optical magnitude plane (L_radio−M_opt plane). With the large sample of radio AGNs classified in this work, which includes 5048 FRIs and 2083 FRIIs, it will be very interesting to investigate to what extent the distributions of FRIs and FRIIs differ in the L_radio−M_opt plane.

The 1.4 GHz radio luminosities of the sources in our sample are calculated with the flux densities measured at 1.4 GHz, as listed in Table 1 of Best & Heckman (2012). The optical r-band luminosities are calculated by utilizing photometric parameters and spectroscopic redshifts drawn from SDSS (York et al. 2000). As recommended by SDSS, we adopt the cmodel magnitude, one of the magnitude measurements provided in the SDSS database, to estimate the optical fluxes of radio AGNs. By using the SDSS photometric pipeline to combine the de Vaucouleurs and exponential models (into a linearly composite model or cmodel) to obtain best fits of the SDSS images, we obtain the composite flux of each source as follows:

$$\begin{eqnarray}&&{F}_{\mathrm{composite}}={f}_{\mathrm{deV}}\times {F}_{\mathrm{deV}}+(1-{f}_{\mathrm{deV}})\times {F}_{\exp },\end{eqnarray} \tag{ 13 }$$

where f_deV is the coefficient of the de Vaucouleurs model; meanwhile F_deV and F_exp are the best-fit fluxes of the de Vaucouleurs and exponential models, respectively. With the apparent r-band magnitude m_r, which is calculated from F_composite, we have the absolute r-band magnitude M_r:

$$\begin{eqnarray}&&{M}_{r}={m}_{r}-5{\rm{l}}{\rm{o}}{\rm{g}}\left(\displaystyle \frac{{D}_{L}}{{\rm{p}}{\rm{c}}}\right)+5-A-K,\end{eqnarray} \tag{ 14 }$$

where D_L is the luminosity distance, A is the correction for the Galactic extinction derived with the map of Schlegel et al. (1998), and K is the k-correction calculated following Blanton & Roweis (2007).

In Figure 10 we show the distributions of typical FRI (blue and cyan circles), typical FRII (red and coral circles), FRI-like (blue and cyan triangles), and FRII-like (red and coral triangles) sources in our sample on the L_radio−M_opt plane, where the classification based on the cross-check between the MCRGNet prediction and human inspection is marked in blue and red and the classification purely made by humans (i.e., for sources not correctly classified by MCRGNet, which are termed "human-supplemented") is marked in cyan and coral. For the cross-checked sources, a probability threshold of 0.90 is applied to filter those sources with smaller predicted probabilities (i.e., sources with high uncertainties; see Section 5.4 for a detailed discussion).

Zoom In Zoom Out Reset image size

The corresponding number distributions as a function of either L_{1.4 GHz} or M_r, as well as their best Gaussian fits, are plotted in the lower panel of Figure 11 for comparison. We find that at any given M_r the distribution of typical FRIIs completely overlaps that of typical FRIs, although the latter extends to lower radio luminosities by about one order of magnitude. Meanwhile the FRI-likes and FRII-likes tend to share the same distribution. In addition to this, we find that the division line between FRIs and FRIIs that was reported in Ledlow & Owen (1996) crosses either of these distributions approximately in the middle. These results, which are obtained with a sufficiently large sample, confirm those reported in previous studies (Best 2009; Lin et al. 2010; Wing & Blanton 2011; Capetti et al. 2017a) with good consistency, indicating that an abrupt separation between FRIs and FRIIs in the L_radio−M_opt plane cannot be reproduced. A further discussion on this issue with an additional criterion on the samples will be presented in Section 5.4.

Zoom In Zoom Out Reset image size

In order to search for the possible redshift evolution of the L_radio−M_opt distributions, we have divided the sources into z ≤ 0.3 and 0.3 < z ≤ 0.7 subgroups and plot their distributions in Figure 12, where the markers and colors of the sources are the same as in Figure 10. We find that roughly the same conclusion as above can be drawn for either of the two redshift ranges. Compared with FRIs located in z ≤ 0.3, FRIs in 0.3 < z ≤ 0.7 tend to be confined in relatively higher L_{1.4 GHz} intervals, which is likely to be caused by the selection effect as fainter sources are hard to observe at higher redshifts and/or by the 5 mJy cut of NVSS to the sample. Therefore, we conclude that no redshift evolution is observed with our sample.

Zoom In Zoom Out Reset image size

As can be seen in Section 4.2.2, the fully pre-trained and fine-tuned MCRGNet has achieved very satisfying classification precision on the LRG sample. Comparing the results obtained with MCRGNet and visual classification on the unLRG sample in Section 4.3, we find that the classification precisions for the compact and FRI sources are high (≳80%), but the precisions are low (≲70%) for the FRII and BT sources and are actually worse for the XRGs and RRGs. As for the FRII and BT sources, we note that about 14% of the former (308 FRIIs) and about 12% of the latter (56 FRII-like BTs) are misclassified into each other by MCRGNet, which reproduces the result of Aniyan & Thorat (2017). In addition, about 1% of the FRIIs are misclassified into the compact category, and this can be attributed to relativistic beaming, by which a bright nearer lobe dominates its fainter counterpart in the images. As explained in Section 4.2.2 for the LRG sample classification, the low precision for the XRG and RRG sources in the unLRG sample may also be due to their small numbers (78 XRGs and 48 RRGs) in the sample (14,245 sources) and their ambiguous morphologies.

By cross-checking the misclassified sources, we find that most of them require more of us authors to join to vote them into certain AGN types. Also for these sources, most are found with low probabilities (≲0.90) obtained by the CNN. Tentatively, we set a threshold to the probabilities of the final labels that are predicted by the dichotomous tree classifier, i.e., objects with probabilities higher than the threshold are kept to calculate the precisions. The evolution of the precision of the classification for the six AGN subtypes is illustrated in Figure 13 as a function of the probability thresholds. We find that for the compact and FRI sources, the classification precisions are not nearly affected by changes in the probability threshold, which is consistent with the corresponding high probabilities (≳0.95) of the labels predicted by our MCRGNet (see Table 4). However, for the XRGs and RRGs, which have relatively complex morphologies, the classification precision significantly increases as the threshold increases.

Zoom In Zoom Out Reset image size

Actually, the problems described above are caused by the fact that the completeness of the labeled sample used to train the neural network is still inadequate. Although the pre-training on a large unlabeled sample (the unLRG sample) may be able to extract most of the features necessary to describe the radio galaxies, the fine-tuning process tends to enhance the features of existing instances carried in the LRG sample while weakening the uncarried. For example, no typical relativistically beamed FRII source is correctly involved in the LRG sample, as well as in other labeled FRII samples available in literature, which prohibits the recognition of the corresponding sources. This may be called a data-unseen problem, which could introduce uncertainty to the network (see Section 5.4 for a detailed discussion) and has been discussed in many machine learning based works, e.g., on image classification (e.g., Chatzilari et al. 2015) and object detection (e.g., Zhu et al. 2016). The solution to this problem is straightforward, i.e., enlarging the sample. Actually, providing astronomers with a large morphologically classified sample which is as complete as possible is part of the motivation of our work. By repeating such a process to enlarge the labeled sample data scale, the data-unseen problem can be efficiently mitigated. This is called iteratively training thinking (e.g., Ma et al. 2017).

In our neural network (Figure 2(b)), there is only one fully connected layer, which encapsulates the extracted features of the radio galaxies in our sample. However, in typical CNNs, such as LeNet-5  (LeCun et al. 1995), AlexNet (Krizhevsky et al. 2012), and GoogLeNet (Szegedy et al. 2015), more than one fully connected layer is employed. The reason for using only one fully connected layer in this work is that in a series of tests we find that more than one is not effective in the morphological classification task in this work. In fact we find that by adding more than one fully connected layer, the classification precision decreases sharply instead of increasing. This phenomenon has been detected in many CNN-based works and is named degradation, i.e., the performance of the network saturates or even declines gradually with more layers than needed (e.g., He et al. 2015, 2016). Thus, based on a series of tests, only one fully connected layer is adopted.

The length of the extracted features, i.e., the number of nodes in the fully connected layer, also affects the precision of the radio galaxy classification task and the complexity of the network, which has been discussed in many related machine learning based problems (e.g., Krizhevsky et al. 2012; Szegedy et al. 2015; Ma et al. 2017). In Figure 9 we illustrate the results of tree structure classification using either 64-dimension or 32-dimension vectors. We find that the degree of separation of the six subtypes (i.e., compact, typical FRI, typical FRII, BT, XRG, and RRG) obtained with 64-dimension feature vectors dominates that obtained with 32-dimension vectors, especially for the types of typical FRIs/FRIIs and BTs. This is because the shorter feature vectors squeeze the learned representations of the radio galaxy morphologies, so that some weak but useful information is lost. On the other hand, we choose not to further increase the dimensions of the features because a larger fully connected layer will exponentially increase the scale of the neural network, which will lead to overfitting risk, affect the efficiency of the network, and thus consume more resources (Han et al. 2015; Li et al. 2016).

Two types of uncertainties, namely aleatoric and epistemic uncertainty, are commonly discussed in machine learning based problems (e.g., Kendall & Gal 2017; Perreault Levasseur et al. 2017). The former describes the intrinsic uncertainty related to the quality of the images (e.g., occlusions, lack of visual features, and over/underexposure; see Kendall & Gal 2017 for details), which is expected to be relieved after enough high-quality images are obtained in future astronomical observations. The latter measures the uncertainty caused by factors that the model does not take into account due to the lack of a complete training data set. We will focus on the latter by applying Bayesian inference as shown below.

In the usual way, to measure epistemic uncertainty, the prior distribution of the weight and bias parameters in the model is directly inferred, often referred to as the Bayesian neural network (Gal & Ghahramani 2015). However, such direct calculations will cost a lot of time and computing resource, especially for a large-scale network. Gal & Ghahramani (2015) and Kendall & Gal (2017) have proved that the optimization of the neural network with dropout, which zeros out neurons randomly as a Bernoulli distribution, is approximately equivalent to a form of Bayesian inference. Thus we utilize the dropout layers already appended after each of the convolutional and fully connected layers (Section 4.2.1) by switching them on using the approach of Gal & Ghahramani (2015), which is called Monte Carlo dropout. Note that the dropout layers are switched off at the final stage of the labeling of the unLRG sample (see Table 4). In Figure 14 we illustrate the test loss curves together with the 68% errors as a function of the keep probability of the dropouts for the five branches on the dichotomous tree. As expected in Section 5.2, for a given branch the uncertainty of the classification of the unLRG sample is in general larger than that of the LRG sample. When classifying the unLRG sample, we find again that the losses of the FRI–FRII (green line) and XRG–RRG (purple line) branches are significantly higher than those of the other branches by taking the uncertainty into consideration. This may be attributed to the unLRG sources with ambiguous morphologies (e.g., the FRIIs with relativistic beaming, which is, in contrast, a clear feature under human inspection). These results agree with what we have obtained in Section 5.2. Besides, from the behavior of the error bars and the small fluctuations on the curves as the dropout keep probability decreases, we are confident that the proposed MCRGNet with dropout layers could be robust.

Zoom In Zoom Out Reset image size

If we increase the threshold set to the probability, which is used to quantify how sure it is that an image belongs to a specific radio galaxy subtype (see Appendix B for details), the precision of the predictions tends to increase, especially for the XRGs and RRGs, the numbers of which are small in both the LRG and unLRG samples (see Figure 13 for details). To be on the safe side we suggest an empirical probability threshold of 0.90 for selecting sources from our unLRG sample (Table 4).

In addition, we have eventually propagated the uncertainty of our network (i.e., selecting sources with the probability threshold) to Figures 10–12, by which our conclusion remains the same as that in Section 5.1.