Identification of 4FGL Uncertain Sources at Higher Resolutions with Inverse Discrete Wavelet Transform

The Fermi-LAT collaboration has recently released the incremental version of the 12 yr Fermi-LAT Gamma-ray Source Catalog (4FGL_DR3; Abdollahi et al. 2022), which contains 6659 sources, among which 3809 sources are AGNs, 559 sources are associated with non-AGNs (including pulsars, high-mass binaries, and supernova remnants, etc.), and 2291 uncertain sources (134 sources are associated with counterparts of unknown nature and 2157 unassociated sources). Within the AGNs, there are 3743 sources confirmed with blazars, among which 1456 blazars are associated with BL Lacs, 794 blazars are associated with FSRQs, and the other 1493 Blazars are BCUs that have not been tagged as BL Lacs or FSRQ.

To accomplish Mission A and B, we need to select features that can distinguish one from another. Variability is one of the widely known characteristics of AGNs with respect to others that are detected by Fermi-LAT (Abdollahi et al. 2020, 2022). Therefore, those features that contain variability information (Flux1000, Flux_Band, Variability_Index, and Frac_Variability) should be included for accomplishing Mission A. While FSRQs and BL Lacs demonstrate significant γ-ray spectra, the GeV γ-ray regime situates at the different regimes of the higher hump of blazar spectral energy distribution (Fan et al. 2016; Yang et al. 2022, 2023). Those features that contain spectral information (Flux_Band, Pivot_Energy, PL_Index) should be included for Mission B.

In this case, we compile the data of 13 attributes from 4FGL_DR3, as listed in Table 1.

For Mission A and Mission B, we randomly split the data into three subsets with a ratio of 8:1:1, as shown in Figure 1. The training set is used to fit the ML methods or DNN models, and the validation set aids in fine-tuning hyper-parameters for better performance. Finally, the generalization capability of the model is independently evaluated on the test set. Moreover, to explore the robustness of the model, this split policy was repeated 10 times and in the end 10 data sets were prepared for each mission. In this section, we analyze the attributes of the training set with the methods shown in Figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As mentioned earlier, it is found that real-world data often contain highly redundant and unimportant attributes. The decision tree (DT) is a popular technique to estimate the importance of attributes. The attributes that appear in the tree are considered important, with the frequency of their appearance being their attribute importance. An attribute that appears less frequently is assumed to be less important. RF is composed of multiple DTs and reduces the bias in estimating attribute importance (Breiman 2001). It is always the practice to remove the unimportant attributes and perform the so-called attribute selection (AS) for dimension reduction.

We employ 50,000 DTs to build the RF and take the entropy as the criterion. Other hyper-parameters align to the default values of scikit-learn (Pedregosa et al. 2011). The results of attribute importance on the training sets are averaged over 10 data sets and are shown in Figure 3. Interestingly, from the results, we find that the four most important attributes are Pivot_Energy, PL_Index, Variability_Index, and Frac_Variability in both Data Set A and B.

Zoom In Zoom Out Reset image size

Another method to study attributes is to project the data from the attribute space to a new space. In this new space, the direction with maximum variance is considered to contribute most to the data, while the directions with small variances can be removed without sacrificing crucial information. This method is called principal component analysis (PCA) (Jolliffe & Cadima 2016). For each split data set, we normalize the samples and perform PCA on training sets. The variance ratio is averaged and depicted in Figure 4. We find that most of the information of the training data is concentrated in the first three components (variance ratio larger than 0.1) in both Data Set A and B.

Zoom In Zoom Out Reset image size

Even though the RF model aggregates the attribute importance values from all the DTs (see Section 2.2), it may result in redundancy because it does not consider the correlations among important attributes. To address this issue, the FD theory is applied to estimate the potential contribution of each attribute to the data set and measure the correlations among attributes (Belussi & Faloutsos 1995, 1998). This method is based on the observation that independent attributes contribute more to the data set, while correlated attributes contribute less.

Mathematically, for a data set A with E attributes: ${\mathbb{A}}=\{{a}_{1},{a}_{2},\ldots ,{a}_{E}\}$, the FD theory supposes that the potential existence of correlated attributes leads the set of points in the original E-dimensional space to describe one spatial object in a dimension that is lower than or equal to E. The dimension of the object represented by the data set is called the intrinsic dimension (ID), denoted by D, $D\in {{\mathbb{R}}}^{+}$. The ceiling of the ID ⌈D⌉ is the minimum number of attributes that must be retained to keep the essential characteristics of a data set (de Sousa et al. 2007).

However, the ID of a data set is difficult to obtain. Alternatively, we consider the ID projecting the data set on a subset C ⊂ A , where C is defined by an attribute subspace ${\mathbb{C}}\subset {\mathbb{A}}$. It is named partial intrinsic dimension (PID) on ${\mathbb{C}}$: ${pD}({\mathbb{C}})$. Based on the definitions, the attribute ${a}_{i}\in ({\mathbb{A}}-{\mathbb{C}})$ increases ${pD}({\mathbb{C}})$ by at most its individual contribution (IC) according to the degree of the correlation between a_i and the attributes in ${\mathbb{C}}$. The IC of a_i , i.e., iC(a_i ), is the maximum potential contribution of the attribute a_i to ${pD}({\mathbb{C}})$. The greater the correlation between a_i and the attributes in ${\mathbb{C}}$, the lower its contribution to ${pD}({\mathbb{C}})$.

In addition, iC(a_i ) can be measured by pD({a_i }) and it ranges in [0, 1]. A more independent distribution of the values of a_i leads to iC(a_i ) closer to one, while a more structured distribution brings iC(a_i ) closer to zero (de Sousa et al. 2007). Thus, the E-dimensional data set can be seen as formed by adding the attributes with different contributions to the D-dimensional sub-data set.

Moreover, the degree of correlations among attributes can be measured by a threshold ξ: a sub-data set B ⊂ A is said to be ξ-correlated to another sub-data set C ⊂ A (their attribute spaces ${\mathbb{B}}\cap {\mathbb{C}}=\varnothing$) if every attribute ${a}_{i}\in {\mathbb{B}}$ does not contribute more than ξ ∗ iC(a_i ) to ${pD}({\mathbb{C}})$. The threshold ξ ∈ [0, 1) tunes how strong the correlation between attributes in ${\mathbb{B}}$ and attributes in ${\mathbb{C}}$ should be to be detected.

A greedy algorithm FDASE was developed by de Sousa et al. (2007) to find a subset of attributes whose PID approaches the ID of the whole data set. The resulting subset is called the attribute set core ξ C (ASC), with the given correlation threshold ξ and scale range n. Therefore, the ratio pD(ξ C)/ID normalizes the contribution of ASC to the whole data set. Similarly, for each split data set, we utilize the FDASE algorithm on the training sets with scanning correlation threshold ξ and compare the value of pD(ξ C)/ID with the number of ASC in Figure 5. We fix the scale range n = 50 and plot the IC for each attribute in Figure 6 under a suitable ξ. From the figures, we find that when ξ = 0.5 two attributes (PL_Index and Pivot_Energy) contribute around 98% information to Data Set A. While for Data Set B, four attributes (PL_Index, Pivot_Energy, Flux_Band1, and Flux_Band7) contribute around 92% information when ξ = 0.5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The attribute analysis methods introduced in Section 2 are used to find the important attributes or components that can be retained for the learning process, while other attributes or components are removed for dimension reduction. However, the crude removal may cause the loss of correlation features that degrade the performance.

We propose to retain all of the attributes but perform IDWT on the original samples for representations at higher resolutions. After IDWT, the correlation features are supposed to be highlighted, while the dimension is possibly reduced, which is important when dealing with big data.

The well-known discrete wavelet transform (DWT) is usually implemented with filtering operations for high efficiency. The filters are designed based on the standpoint of multiresolution: the difference of information between the approximation of a signal at the resolutions 2^m+1 and 2^m (where m is an integer) can be extracted by decomposing this signal on an orthonormal basis of wavelets (Mallat 1989). The pyramidal structure of the wavelet filter bank makes it possible to infer information at a low resolution from information at a high resolution. IDWT is the converse process of DWT and provides representations at high resolutions for DL.

However, in practical applications, a finite signal should be considered. The length of the signals varies at different resolutions, which is due to the operations of downsampling, upsampling, and filtering. In the IDWT process, if we denote p and s as the length of the signal at low resolution and high resolution, respectively, then Rajmic & Prusa (2014) give

$$\begin{eqnarray}s=\left\{\begin{array}{cc}2p-u+1, & \mathrm{for}\ \mathrm{even}\ s+u-1,\\ 2p-u+2, & \mathrm{for}\ \mathrm{odd}\ s+u-1,\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where u is the length of the filters.

To perform IDWT on the original data set, the attribute data of one sample will be seen as the signal at a lower resolution. As a result, the attributes should be rearranged in some order based on wavelet theory. This can be achieved with the alignment of the information in the wavelet domain and attribute space, which will be discussed next.

The process of DWT is to analyze the original signal from a fine scale to a coarse scale. The representation of a finite signal f(t) in the wavelet domain after DWT is a collection of vectors:

$$\begin{eqnarray}&&\{{{\boldsymbol{c}}}_{J}(k):k\in {\mathbb{Z}}\}\cup \{{{\boldsymbol{d}}}_{j}(k):1\leqslant j\leqslant J;k\in {\mathbb{Z}}\},\end{eqnarray} \tag{ 2 }$$

where $J\in {\mathbb{N}}$ is called decomposition level; c _J contains the approximation coefficients at level J, i.e., the lowest resolution; and d _j contains the detail coefficients at level j, i.e., the high resolutions. An example of the three-level pyramid transform is illustrated in Figure 7. We find that the original signal f(t) can be seen as the approximation coefficients at the highest resolution, i.e., c ₀.

Zoom In Zoom Out Reset image size

However, based on prior knowledge, it is often observed that at a specific level, the majority of information in a natural signal is typically presented in the approximation coefficients. Additionally, the information is typically reduced by half during the decomposition of the signal from level j to j + 1 in DWT, due to the downsampling operation. Hence, it can be found that the information of coefficients at different levels roughly respect

$$\begin{eqnarray}&&{I}_{c}({{\boldsymbol{d}}}_{J})\lt \cdots \lt {I}_{c}({{\boldsymbol{d}}}_{2})\lt {I}_{c}({{\boldsymbol{d}}}_{1})\lt {I}_{c}({{\boldsymbol{c}}}_{J}).\end{eqnarray} \tag{ 3 }$$

As an example, the three-level DWT is performed on three images to intuitively illustrate the relationship, as shown in Figure 8. These figures suggest that the approximation coefficients at each level retain the image contour and contain a significant amount of information, while the objects in the image cannot be easily identified through the detail coefficients. However, to some extent, the detail coefficients at a low level appear to provide more information than those at a high level, i.e., I_c ( d _j+1) < I_c ( d _j ).

Zoom In Zoom Out Reset image size

To transform the coefficients from lower resolutions to higher resolutions in the IDWT process, the wavelet decomposition vector (WDV) c and bookkeeping vector (BV) l are required according to multiresolution analysis. The WDV includes the coefficients shown in Equation (2):

$$\begin{eqnarray}&&{\boldsymbol{c}}=({{\boldsymbol{c}}}_{J},{{\boldsymbol{d}}}_{J},{{\boldsymbol{d}}}_{J-1},\,\ldots ,\,{{\boldsymbol{d}}}_{1}),\end{eqnarray} \tag{ 4 }$$

while the BV is made up of the number of coefficients in c :

$$\begin{eqnarray}&&{\boldsymbol{l}}=(\mathrm{len}({{\boldsymbol{c}}}_{J}),\mathrm{len}({{\boldsymbol{d}}}_{J}),\mathrm{len}({{\boldsymbol{d}}}_{J-1}),\,\ldots ,\,\mathrm{len}({{\boldsymbol{d}}}_{1}),\mathrm{len}(f(t))),\end{eqnarray} \tag{ 5 }$$

where len(·) indicates the number of the elements in a vector.

In the field of big data, an attribute's information reflects its contribution to the data set. As explained in Section 2.4, the IC estimates the potential contribution of an individual attribute to the data set. However, the presence of correlations among attributes means that the actual contribution of an attribute a_i , i.e., its information I_a ({a_i }), cannot be precisely determined if a_i involves other attributes.

Moreover, as discussed, ASC is the smallest subset of attributes that can fully characterize the entire data set. Given a correlation threshold ξ, the ASC ξ C can be found with the FDASE algorithm (de Sousa et al. 2007). The attributes in ASC are not ξ-correlated with each other and they contain the most information in the data set.

Based on the analysis, if ${\mathbb{A}}=\{{a}_{1},{a}_{2},\ldots ,{a}_{E}\}$ is denoted as the universal attribute set of the data set that has E attributes, then the real contribution of attribute ${a}_{i}\in ({\mathbb{A}}-\xi C)$ to the data set can be seen as the degree of correlation between a_i and the attributes in ξ C. The weaker the correlation between a_i and ξ C, the higher the contribution of a_i to the data set, i.e., attribute a_i contains more information. Therefore, the information of attribute a_i can be estimated by

$$\begin{eqnarray}&&{I}_{a}(\{{a}_{i}\})=| {pD}(\{{a}_{i}\}\cup \xi C)-{pD}(\xi C)| ,{a}_{i}\in ({\mathbb{A}}-\xi C),\end{eqnarray} \tag{ 6 }$$

where pD(·) is the PID of a sub-data set.

By estimating the information for each attribute ${a}_{i}\in ({\mathbb{A}}-\xi C)$, we can obtain a rough ordering of attributes based on their information. In this case, the data under the attributes can be treated as wavelet coefficients and arranged to form the WDV c and BV l for IDWT. The data set is then transformed into representations at a higher resolution where the correlation features are presented. This is the main idea of our method, which will be detailed next.

The J-level IDWT aims to transform the original data into representations at the J-higher resolution. To achieve this, the WDV and BV must be constructed with the attribute data of each sample. However, the positions of the data in the WDV are determined by the estimation of the attribute information with Equation (6). As a result, this approach is referred to as FDIDWT.

As discussed earlier, the attributes in ASC $\xi C\subset {\mathbb{A}}$ contain the most information of the data set defined in ${\mathbb{A}}=\{{a}_{1},{a}_{2},\ldots ,{a}_{E}\}$. It is assumed that the information of attributes in ${\mathbb{A}}-\xi C$ respect an ascending order:

$$\begin{eqnarray}&&{I}_{a}(\{{a}_{1}\})\lt {I}_{a}(\{{a}_{2}\})\lt \cdots \lt {I}_{a}(\{{a}_{i}\})\lt \cdots \lt {I}_{a}(\{{a}_{Q}\})\lt {I}_{a}(\xi C),\end{eqnarray} \tag{ 7 }$$

where ${a}_{i}\in ({\mathbb{A}}-\xi C)$ and 1 ≤ i ≤ Q. If the number of attributes in ξ C is P, then P + Q = E. Contrasting with the order of the information in the wavelet domain shown in Equation (3), we group the attributes in ${\mathbb{A}}-\xi C$ and obtain a similar order of the information in the attribute space:

$$\begin{eqnarray}&&{I}_{a}({D}_{J}^{M})\lt {I}_{a}({D}_{J-1}^{N})\lt \cdots \lt {I}_{a}({D}_{j}^{L})\lt \cdots \lt {I}_{a}({D}_{1}^{K})\lt {I}_{a}(\xi C),\end{eqnarray} \tag{ 8 }$$

where D^L _j is the attribute group containing L attributes whose information respects to Equation (7), i.e.,

$$\begin{eqnarray}&&\begin{array}{l}{D}_{J}^{M}=\{{a}_{1},{a}_{2},\,\ldots ,\,{a}_{M}\},\\ {D}_{J-1}^{N}=\{{a}_{M+1},{a}_{M+2},\,\ldots ,\,{a}_{M+N}\},\\ ...\\ {D}_{1}^{K}=\{{a}_{Q-K+1},{a}_{Q-K+2},\,\ldots ,\,{a}_{Q}\},\\ M+N+\cdots +K=Q.\end{array}\end{eqnarray} \tag{ 9 }$$

Thus, the data under the attributes of D^L _j and ξ C can be taken as the detail coefficients d _j and approximation coefficients c _J , respectively. Specifically, if denote x _j and x _{ξ C} as the data under the attributes in D^L _j and ξ C, respectively, for one sample of the original data set, then the WDV c and BV l for the IDWT process can be constructed as

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{c}} & = & \left({{\boldsymbol{x}}}_{\xi C},{{\boldsymbol{x}}}_{J},{{\boldsymbol{x}}}_{J-1},\,\ldots ,\,{{\boldsymbol{x}}}_{1}\right),\\ {\boldsymbol{l}} & = & \left(P,M,N,\,\ldots ,\,K,O\right),\end{array}\end{eqnarray} \tag{ 10 }$$

where O is the number of attributes in the transformed data set, i.e., the transformed dimension in the new attribute space, which can be derived with Equation (1). Table 2 shows the procedure of the proposed FDIDWT method. To better describe how the algorithm works, we fix the decomposition level J = 3. Moreover, as in the common case, the length of filters u for the pyramid transform is set to be even. Then:

• 1.  
  Steps 1–3: Given the correlation threshold ξ and scale range n, the FDASE algorithm (de Sousa et al. 2007) is used to find out the attribute set core ξ C. The information I_a ({a_i }) of the remaining attributes ${a}_{i}\in ({\mathbb{A}}-\xi C)$ is then calculated. The attributes of ${\mathbb{A}}$ are arranged into a set ${\mathbb{A}}^{\prime}$ according to the ascending order of I_a ({a_i }). ⁷ For instance, if the relationship in Equation (7) holds, then

  $$\begin{eqnarray}&&{\mathbb{A}}^{\prime} =\xi C\cup \{{a}_{1},{a}_{2},\,\ldots ,\,{a}_{i},\,\ldots ,\,{a}_{Q}\}.\end{eqnarray} \tag{ 11 }$$

  In our experiment, we scanned the correlation threshold ξ in the range [0.05, 1.0] with step 0.05, and set scale range n = 50. We got the ASC {a₃, a₁} and {a₃, a₁, a₅} for Data Set A and B, respectively, as shown in Figures 5 and 6. Then, the arranged attribute set ${\mathbb{A}}^{\prime}$ for the two data sets will be {a₃, a₁, a₅, a₁₃, a₁₀, a₂, a₉, a₁₂, a₄, a₈, a₁₁, a₇, a₆} and {a₃, a₁, a₅, a₁₃, a₄, a₁₂, a₁₁, a₇, a₈, a₆, a₉, a₂, a₁₀}.
• 2.  
  Steps 4–5: P is the number of attributes in ξ C. In the proposed method, the length of the approximation coefficients at level J is assumed to be not shorter than P. If this length equals P, then Figure 9 shows other groups divided from ${\mathbb{A}}^{\prime}$. ${C}_{3}^{P}$ denotes the approximation coefficient group of attributes in ξ C. Note that the original dimension E, i.e., the number of attributes of the original data set, may be smaller than the length of WDV c due to the given decomposition level J and filter length u. In this case, some necessary placeholder attributes "∗" are inserted following the ASC group ${C}_{3}^{P}$. The role of the placeholder attributes is to ensure that IDWT works while not changing the information of the data. Therefore, the values under the placeholder attributes are set to zeros so that they contribute nothing to the data set. In the signal processing field, signals are always transformed into the wavelet domain with DWT for processing, and the lengths of the wavelet coefficients at each level are recorded for the reconstruction with IDWT. However, there is no DWT in the proposed method. Hence, many cases exist for the lengths of the divided groups based on Equation (1):

  $$\begin{eqnarray}&&\begin{array}{l}M=P,\\ N=\left\{\begin{array}{cc}2M-u+1, & \mathrm{for}\ \mathrm{even}\ N+u-1\\ 2M-u+2, & \mathrm{for}\ \mathrm{odd}\ N+u-1\end{array}\right.,\\ K=\left\{\begin{array}{cc}2N-u+1, & \mathrm{for}\ \mathrm{even}\ K+u-1\\ 2N-u+2, & \mathrm{for}\ \mathrm{odd}\ K+u-1.\end{array}\right.\end{array}\end{eqnarray} \tag{ 12 }$$

  The length of the placeholder attributes L can be derived based on the equality of the length of WDV, i.e., P + M + N + K = E + L. Then, the WDV and BV are constituted as follows:

  $$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{c}} & = & \left({{\boldsymbol{x}}}_{\xi C},{{\boldsymbol{x}}}_{3},{{\boldsymbol{x}}}_{2},{{\boldsymbol{x}}}_{1}\right),\\ {\boldsymbol{l}} & = & \left(P,M,N,K,O\right),\end{array}\end{eqnarray} \tag{ 13 }$$

  where O is the transformed dimension computed with:

  $$\begin{eqnarray}O=\left\{\begin{array}{cc}2K-u+1, & \mathrm{for}\ \mathrm{even}\ O+u-1\\ 2K-u+2, & \mathrm{for}\ \mathrm{odd}\ O+u-1.\end{array}\right.\end{eqnarray} \tag{ 14 }$$

  We used the Daubechies filter, which can be indicated by its length u. For example, the filter length u = 8 indicates that the db4 wavelet filter was used. Under the assumption of the maximum decomposition level J = 3, we built the WDV and BV for two missions in Table 3.
• 3.  
  Steps 6–7: Finally, IDWT is performed on each sample of the normalized data set, and a new data set defined on ${\mathbb{B}}=\{{b}_{1},{b}_{2},\,\ldots ,\,{b}_{i},\,\ldots ,\,{b}_{O}\}$ is generated. b_i is the ith new attribute.

Zoom In Zoom Out Reset image size

Table 2. The Proposed FDIDWT Method

Input:	original data set defined on attribute space ${\mathbb{A}}$
	decomposition level J
	correlation threshold ξ
	scale range n
Output:	new data set defined on attribute space ${\mathbb{B}}$
1:	run FDASE algorithm to find out the ASC ξ C, denote its length as P;
2:	calculate Ia (ai ), ${a}_{i}\in ({\mathbb{A}}-\xi C)$, with Equation (6);
3:	arrange the attributes of ${\mathbb{A}}$ into ${\mathbb{A}}^{\prime}$ according to the order of Ia (ai ) as shown in Equation (7);
4:	based on P, calculate the length of other groups, insert placeholder attributes if needed;
5:	construct the WDV c and BV l ;
6:	perform IDWT on each normalized sample with c and l ;
7:	output the new data set defined on ${\mathbb{B}};$

Download table as:  ASCIITypeset image

Since the new attribute space ${\mathbb{B}}$ is considered to respect some natural order after FDIDWT, we believe that the convolution operation with a kernel will learn more features and achieve better performance. The kernel of CNN determines a visual field on the adjacent attributes in DL. Thus, we utilized a lightweight CNN model modified from the Matchbox net (Majumdar & Ginsburg 2020) but with a smaller size of only 94k parameters for high efficiency, which we refer to as MatchboxConv1D. Its structure is depicted in Figure 10.

Zoom In Zoom Out Reset image size

Extracting correlation features among the attributes of Fermi sources is a promising research field. The commonly used attribute analysis methods, such as the estimation of attribute importance (see Section 2.2) or principal components (see Section 2.3), aim to find the significant attributes or components, and remove the unimportant ones for dimension reduction. However, the correlation features will be removed with these methods, which results in inferior classification performance (see the results of Tables 6 and 7).

The FDASE gives a global perspective of the whole data set and estimates the correlation features among all of the attributes based on FD theory. The ASC is the resulting attribute subset of the FDASE algorithm and its data contains most of the information of the data set under some correlation threshold and scale range. Nevertheless, the retention of ASC also suffers from the loss of correlation features (see the results of FDASE attribute selection in Tables 6 and 7). Thus, we propose to rearrange the original attributes to highlight correlation features at higher resolutions based on FDASE and IDWT, which is referred to as the proposed FDIDWT method.

The resulting attribute space is also considered to respect some natural order after FDIDWT. We believe that the convolution operation will extract more correlation features from the ordered attributes. Meanwhile, the structure of the CNN gives more potentiality for better classification performance. Therefore, by combining FDIDWT with the MatchboxConv1D model, we obtain the best test accuracy for both Mission A and B in case #6. In addition, the dimension is transferred to be 8, which is significant for reducing the computing burden in big astronomical data processing. It may also be concluded that most features of Data Set A and B exist at one higher resolution because case #6 corresponds to decomposition level 1 (see Table 3).

The prediction results have been carried out with an accuracy of 96.65% ± 1.32% and 92.03% ± 2.2% for Mission A and Mission B, respectively, using the proposed method. The results are listed in Table 9. A general comparison between the predicted AGNs of Mission A and the originally confirmed AGNs is shown in Figure 13. A comparison between the predicted BL Lacs of Mission B and the originally confirmed BL Lacs is shown in Figure 14. Finally, a comparison between the predicted FSRQs of Mission B and the originally confirmed FSRQs is shown in Figure 15. From the comparison, we can find that the distribution shapes of the 13 attributes of the predicted sources resemble the shapes of the original 4FGL_DR3 sources with the corresponding classification, in general. The results indicate that our predicted sources are correctly classified.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We notice that the histograms of the attribute Variability_Index show a longer high-variation-index tail for the original 4FGL_DR3 sources than for the predicted sources (including predicted AGNs, predicted BL Lacs, and predicted FSRQs in Figures 13–15, respectively), and the predicted sources have more contribution on the high-variation-index head than the original 4FGL_DR3 sources. This result suggests that our method has the advantage of finding less variable sources from the uncertain sources. We also notice that the original 4FGL_DR3 sources give more contribution to the histogram tails than the predicted sources, while the predicted sources give more contribution to the histogram head for the attributes of multi-band intensities (Flux1000, Flux_Band1, Flux_Band2, Flux_Band3, Flux_Band4, Flux_Band5, Flux_Band6, Flux_Band7, and Flux_Band8). This suggests that our method has the advantage of finding relatively faint γ-ray sources from the uncertain sources. It also encourages us to believe that our method should be able to make efforts to identify less variable and faint sources in the era of survey telescopes, e.g., the Large Synoptic Survey Telescope (LSST; LSST Science Collaboration et al. 2009), China Space Station Telescope (CSST; Zhan 2011), etc.

By employing the proposed method, we managed to predict 2291 uncertain sources into 1731 AGNs and 560 non-AGNs, and predict 1493 BCUs into 948 BL Lacs and 545 FSRQs. The likelihood probabilities of these predicted uncertain sources and BCUs are shown in Table 9 and the distribution is displayed in Figure 16.

Zoom In Zoom Out Reset image size

We set a boundary of likelihood probability greater than 95%, which is shown as the dashed-red lines in Figure 16, to claim this source as a candidate of the corresponding class. In this case, we get further constrained results into 1354 AGN candidates in Mission A, 482 BL Lacs candidates and 128 FSRQ candidates in Mission B, as shown in the last column of Table 9.

Several methods have been employed to classify 4FGL BCUs as FSRQ candidates or BL Lac candidates, as Mission B of the present work, in the previous works. Kang et al. (2019) utilized three supervised ML methods (RF, SVM, and ANN) to make a classification for 1312 4FGL_DR1 BCUs, and carried out a combined classification result of 724 BL Lac candidates and 332 FSRQ candidates. By crossing matching the results from Kang et al. (2019) and our results, we found that there are 419 overlapping BCUs in these two works. Among the overlapping 419 BCUs, 324 BCUs are predicted as BL Lac candidates and 90 BCUs are predicted as FSRQ candidates in both works, which gives our results a consistency of 98.8% with Kang et al.'s work.

Similarly, five different supervised ML algorithms (RF, LR, XGBoost, CatBoost, and neural network) were applied to the 4LAC_DR3 BCUs in Agarwal (2023), and 610 BCUs were able to be classified as BL Lac candidates and 333 BCUs were classified as FSRQ candidates. A comparison between our work and Agarwal's work suggests that there are 481 overlapping BCUs in the two works, 392 BCUs are classified as BL Lac candidates and 87 BCUs are classified as FSRQ candidates in both works, which gives our results a consistency of 99.6% with their work.

Fan et al. (2022) employed three physical parameters and built diagrams among them (i.e., the photon spectrum index against the photon flux diagram, the photon spectrum index against the variability index diagram, and the variability index against the photon flux) to separate the known BL Lacs from known FSRQs. They then used the boundary to divide BCUs into BL Lacs and FSRQs. In their work, 751 BCUs were classified as BL Lac candidates and 210 BCUs were classified as FSRQ candidates. There are 492 overlapping BCUs in the two works, 409 BCUs are classified as BL Lac candidates and 83 BCUs are classified as FSRQ candidates in both works, which gives our results a consistency of 100% with their work.