Machine-learning Application to Fermi-LAT Data: Sharpening All-sky Map and Emphasizing Variable Sources

In this study, ML-based algorithms were applied to a gamma-ray ASM observed by Fermi-LAT. The first objective was to fabricate a high-quality ASM with sufficient photon statistics such that the ASM uncertainties are small compared to the counting uncertainties of shorter timescale maps when only a low-quality ASM is known and/or provided as an input. As discussed in Section 3, dictionary learning and U-net require the image pairs of the short-observation map and long-observation map for image prediction. Moreover, because the Noise2Noise demands image pairs with different noise patterns, two types of short-observation maps, corresponding to the first and second week in each year, were prepared. Therefore, a long-term observation map and two types of short-term observation maps in the energy range from 100 MeV to 300 GeV were created at the beginning of each year for years 2009–2019. In this study, all energy values were incorporated into a single map. We produced the ASMs using gamma rays in the P8R3_SOURCE event selection and the associated P8R3_SOURCE_V2 instrument response functions (IRFs). In addition, to correct the ASM for exposure, we calculated an exposure map with one energy bin using gtexpcube2, which is the standard analysis tool provided with Fermi-LAT ScienceTools. We also divided the ASM by the generated exposure map. In addition, the zenith angle cut was set to 90° to remove gamma rays from Earth's limb in the analysis. In this study, ASMs were produced in the galactic coordinates in the Hammer–Aitoff projection. We generated count maps with a resolution of 0.1 deg² per pixel.

The ML algorithms assume that all sources are persistent over 1 yr, keeping their fluxes exactly comparable to what has been observed in the first week. Thus, a large deviation of the ML maps and the GT map may suggest that variable sources are included in regions of interest (ROIs). In such cases, a light curve was generated by using binned analysis to investigate the flux variability for specific gamma-ray sources. For each source, a 15° × 15° ROI centered at its 4FGL catalog (Abdollahi et al. 2020) was chosen, and a single power-law model was applied to the data. When creating the light curve, the prefactor and the index of the central source were free, while the parameters of the other sources within the ROI were fixed. As the background model, the isotropic and diffuse backgrounds were modeled using ISO_P8R3_SOURCE_V 2_v1.txt and gll_iem_v07.fits, respectively. Moreover, only the flux bins with test static (TS) > 9 (corresponding to ∼3σしぐま) were used for each light curve.

The following four processing steps were applied to each ASM obtained in Section 2.1 before training the ML algorithms:

• 1.  
  applying Gaussian blurring;
• 2.  
  applying the natural logarithm function;
• 3.  
  standardization of the pixel values in each map;
• 4.  
  data augmentation.

First, Gaussian blurring was performed for each ASM to reduce the number of pixels whose value equals zero. The blurring process is necessary to prevent ML algorithms from learning only map smoothing. Subsequently, the natural logarithm function was applied to these maps owing to the high variance of pixel values; then, the pixel values were normalized from 0 to 1. Finally, random sampling was performed for data augmentation. In this process, the same regions were extracted from the ASM pairs with short-observation and long-observation maps, and the extracted images were randomly flipped. The size of the extracted images was 64 × 64 pixels, which corresponds to 6° × 6° in the galactic center. The random sampling was repeated 2000 times for each ASM pair. Detailed information about the components of the data sets is presented in Table 1.

Hereinafter, "input map" and "GT map" refer to the ASMs, which cover the time intervals MJD 58486.034–58493.026 (2009 January 3–9) and MJD 58486.034–58857.048 (2009 January 2–2020 January 9), respectively. The input map was blurred with a 11 pixel × 11 pixel Gaussian kernel. In the input map, standard deviation, which defines the smoothing effect, equals 3.0. Additionally, we will use "blurred map" to refer to blurred versions of the ASM, which were blurred with a 15 pixel × 15 pixel Gaussian blur kernel with a standard deviation equal to 100. "ML maps" refers to an ASM generated from the input map using three ML-based algorithms (dictionary learning, U-net, and Noise2Noise). We also note that the ML algorithm applied in this study only suggests the source variability in the corresponding image pixels. Therefore, detailed follow-up analysis such as drawing light curves of candidate sources was necessary to check the variability, as described in detail in Section 5.2.

Dictionary learning is generally used for dimensional reduction in the field of signal processing (Tosic & Frossard 2011). Unlike the deep-learning algorithms, dictionary learning is based on linear algebra; hence, one can sequentially trace the obtained results using this technique. Dictionary learning involves three steps: preparing a pair of databases, creating dictionaries, and generating images. First, five patches with 16 × 16 pixels were extracted randomly from each training sample of a short-observation map. The patches were arranged in columns and aligned to create a database of short-observation data (A_short). In the same manner as A_short, a database of long-observation data (A_long) was also created. Then, the arrays, called the "dictionary," were calculated. In this study, the dictionaries from short-observation data (D_short) and long-observation data (D_long) were created from A_short and A_long by sparse representation. After the dictionaries were calculated, the test samples from the input map were represented by D_short and sparse matrix C. Finally, D_short was replaced with D_long to generate images corresponding to the GT map. The details of the calculation in dictionary learning are provided elsewhere (Sato et al. 2020). In this study, scikit-learn (version 0.21.2), which is a library of ML algorithms used in Python, was used for dictionary learning and sparse coding.

Originally, a U-net was developed for masking the medical images (Dong et al. 2017; Han & Ye 2018; Zhou et al. 2018). An example of the algorithm architecture is illustrated in Figure 1. U-net consists of two paths: the contracting path and expanding path. Two convolution layers and a max pooling layer are repeated four times in the contracting path. In this study, 3 × 3 convolution was used with leaky rectified linear unit (ReLU) activation and batch normalization. The number of channels was doubled for each downsampling. The expanding path consists of two convolution layers and a deconvolution layer with 3 × 3 filter sizes. The ReLU function was used as the activation function after the convolution layers except the last one. After the last convolutional layer, a hyperbolic tangent function was used. In this study, the "GeForce 1080Ti" graphics processing unit (GPU) processer was employed, and the U-net was implemented with TensorFlow (version 2.1.0) and Keras (version 2.3.1), which are libraries of deep-learning algorithms used in Python.

Zoom In Zoom Out Reset image size

Noise2Noise, which is a semisupervised learning algorithm, is used to reconstruct clean signals from corrupted ones (Wu et al. 2019; Hasan et al. 2021; Díaz Baso et al. 2019). This algorithm can perform signal de-noising without the need for clean signal data sets. Therefore, Noise2Noise is applicable to cases where clean signals are difficult to obtain. The superresolution residual network (SRResNet; Ledig et al. 2017) used in this study consists of a convolution layer with parametric rectified linear unit (PReLU) activation, six layers of residual block, two layers of convolution layer with batch normalization, and a concatenated layer, as depicted in the left panel of Figure 2. The right panel of Figure 2 depicts the details of the residual block. The residual block consists of a convolution layer, batch normalization, and PReLU activation. In this study, the filter number and size in the convolution layer were set as 64 and 3 × 3, respectively. The software and GPU processer were used in the same way as for U-net in Section 3.2.

Zoom In Zoom Out Reset image size

To confirm the improvement of maps in detail, close-ups of the various regions in Figure 3 were extracted, as shown in Figure 4. Image noise resulting from low photon statistics is reduced by ML algorithms; thus, the parts of the point source are enhanced. Moreover, the feature of the diffuse gas emission around the galactic plane is also well reproduced.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To evaluate the image improvement in ML maps, we applied a root mean square error (RMSE), which is the standard evaluation index in computer vision. The RMSE represents the squared error of each pixel value in the two images and was calculated using the following equation:

$$\begin{eqnarray}&&\mathrm{RMSE}=\sqrt{\displaystyle \frac{1}{n}{{\rm{\Sigma }}}_{{\rm{i}}=1}^{{\rm{n}}}{({V}_{{\rm{i}}}^{\mathrm{image}1}-{V}_{{\rm{i}}}^{\mathrm{image}2})}^{2}},\end{eqnarray} \tag{ 1 }$$

where n and V_i denote the number of image pixels and the ith pixel value in each input image, respectively. An RMSE that is closer to zero represents a smaller error between the input images. Similarly, we also computed the structural similarity (SSIM) between the produced maps and the GT map. The SSIM, which represents visual similarity, was calculated using

$$\begin{eqnarray}&&\mathrm{SSIM}=\displaystyle \frac{(2{\mu }_{x}{\mu }_{y}+{{\rm{C}}}_{1})(2{\sigma }_{{xy}}+{{\rm{C}}}_{2})}{({\mu }_{x}^{2}+{\mu }_{y}^{2}+{{\rm{C}}}_{1})({\sigma }_{x}^{2}+{\sigma }_{y}^{2}+{{\rm{C}}}_{2}),},\end{eqnarray} \tag{ 2 }$$

where μみゅー and σしぐま denote the average value and standard deviation of the image, respectively. C₁ and C₂ are small constants for avoiding instability in the case in which ${\mu }_{x}^{2}+{\mu }_{y}^{2}$ or ${\sigma }_{x}^{2}+{\sigma }_{y}^{2}$ is close to zero. Details of the constants are described in Wang et al. (2004). Figure 5 depicts the results of RMSE and SSIM evaluation. In Figure 5, the horizontal axis represents the observation time, and the vertical axis represents the RMSE or SSIM values. We discuss the evaluation results in Section 5.3.

Zoom In Zoom Out Reset image size

For source detection, we compared the excess of pixel values in the ML maps corresponding to the source position, as already listed/known in the 4FGL catalog. The bottom panel of Figure 6 shows an example of a projection histogram, the region of which is shown in the top panel. In Figure 6, 15 representative sources with the highest average flux are indicated. The signal-to-noise ratio (S/N) in the projection was improved by ML algorithms. However, unlike input and GT maps based on photon statistics, pixel values in the ML maps lose any physical meaning; thus, we could not evaluate TS, which is often applied in the analysis of Fermi-LAT to evaluate the significance of source detection. We will develop such advanced algorithms for ML maps in the future; however, it is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

To find the correspondence between the ML maps and the GT map, the 1D projections along the same galactic latitude l were compared for various regions. An example of a projection histogram obtained by the GT map, U-net map, and Noise2Noise map is depicted in the top panel of Figure 7. The extracted part is the same as in Figure 6. The projection was created by integrating the pixel values for 15° in the direction of galactic latitude l. For a quantitative comparison, the values in the ML maps were divided by those in the GT map. The bottom panel of Figure 7 indicates the ratio of the ML maps to the GT map. Three regions with large differences are marked by orange circles. Figure 8 shows close-ups around region 2 to illustrate such a difference. While no evident sources exist in the input map and ML maps, a bright point source, which is identified as PKS 2052−47, is found only in the GT map. In the opposite case shown in Figure 7, there exist regions in which the ratio is smaller than unity, as represented by region 4 in Figure 9. Figure 10 shows close-ups around region 4. The details of regions with apparent "mismatch" are listed in Table 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this study, three types of ML algorithms were applied: dictionary learning, U-net, and Noise2Noise. Dictionary learning is based on linear algebra and can be used to trace the calculation process sequentially. However, dictionary learning has a low representation capability compared with other deep-learning algorithms and reflects a relatively low quality image compared with other algorithms presented in Figure 4. Thus, dictionary learning is less effective than other deep-learning algorithms, for the purpose of this study. In the U-net architecture, the input and output layers are connected by two paths: an extended path and a shortcut connection. Since these two paths are considered to capture detailed features and context, U-net can successfully predict the long-observation map from the input map with high precision (see Figure 4). However, in Noise2Noise, the algorithm captures features of statistical noise by comparing two types of short-observation maps. Therefore, the algorithm is used to generate the ASM by neglecting the statistical noise, which can sometimes overestimate the significance of each point source. Thus, U-net may be superior to Noise2Noise for predicting a long-observation map. In this study, U-net was considered to be the most reliable/precise for predicting the long-observation map. Furthermore, because these three ML algorithms are based on completely different theories, ML maps are more reliable if the outputs from each ML algorithm are similar to each other.

In this section, a possible reason for the mismatch between the GT map and ML maps, as shown in Figures 7 and 9, is discussed. First, the regions of the largest mismatch, indicated as orange circles in Figures 7 and 9, are discussed. The ML maps are generated based on the hypothesis that the sources are persistent over 1 yr, which is, however, not the case for variable and/or transient sources. Thus, the observed mismatch between the ML maps and GT map is likely associated with the variable and transient sources, such as the flaring activity of AGNs. To confirm this assumption, the sources listed in Table 2 are analyzed in more detail with the Fermi analysis tool, as described in Section 2.1. The light curves around each source are shown in Figure 11. For example, a variable source in region 1 corresponds to a flat spectrum radio quasar (FSRQ) 3C 454.3 with a redshift of z = 0.859 (Table 2). As shown in Figure 11, the flux of 3C 454.3 in the first week was much smaller than the average flux, shown as an orange line, integrated over a year. Similarly, regions 2 and 3 were identified as two FSRQs: PKS 2052−47 and PKS 0454−234, respectively. Contrary to the sources in which flux substantially increased after the first week, the average flux over the year is smaller than the first week's flux in the case of region 4, which corresponds to FSRQ PKS 1824−582. Obviously, such temporal variability in fluxes explains the observed mismatch between the ML maps and the GT map, as shown in Figures 7 and 9. As can be seen in this figure, point sources with increasing photon flux appear to be brighter in the bottom panel of Figure 3. On the contrary, point sources with decreasing flux appear faint.

Zoom In Zoom Out Reset image size

However, there is a systematic difference between the GT map and ML maps around the galactic longitude 270° < l < 330° in the bottom panel of Figure 7, which cannot be explained by either the statistical fluctuation or a combination of variable sources. We also confirm the difference in the bottom panel of Figure 3. Such a systematic gradient is thought to be caused by the nonuniformity of the exposures inherent in the 1-week observation, which were used as input for ML algorithms. Although exposure map correction is performed based on the position and angle of the satellite, such nonuniformity affects the reconstructed counting map, especially in the case of short observation times. In fact, we confirmed that such a gradient is reduced if 2-week ASM maps are used as input, but the optimization of the time duration for the input maps is beyond the scope of this study. Additionally, ASMs in this paper were produced with no energy binning; hence, the exposure map might not correct the actual exposure. This is because the accuracy of integrating the exposure over the energy band depends on the smoothness of the source spectrum across the energy band, and these are quite large energy bands.

As shown in Figure 5, the ML-based algorithms were effective in improving the ASM quality. In particular, U-net was the best among the three ML algorithms. The RMSE and SSIM assessments of the U-net map were equivalent to approximately 10 and 12 weeks, respectively. Considering that the exposure map and variable source were incorrect, the image quality was improved by more than 10 weeks. Despite the significant improvement, faint source detection was considered difficult for our algorithms owing to image blur. Some peaks around the source location were blurred, as shown in Figure 6. In general, for ML processing, the loss function must be defined. "Training ML algorithms" refer to parameter optimization, which could be used for minimizing the loss function. The RMSE, which is one of the standard loss functions, was used in this study. However, in the case of using MSE or RMSE, it is known that high-frequency components are difficult to reproduce. Due to the difficulty of point-source detection, some schemes are necessary to improve the detection efficiency. For example, adopting the sum of point-source detectability and RMSE as the loss function may improve the performance. Other approaches are mentioned in Section 5.5.

Compared with other techniques described in Section 1, our ML-based method has the following three advantages:

• 1.  
  The ML algorithm can be easily implemented by end users. Moreover, any assumption or modeling of foreground diffuse emission is unnecessary. Owing to its simplicity, the ML ASM can be obtained within less than 10 minutes.
• 2.  
  S/N on the image is improved in terms of both RMSE and SSIM, although the optimum algorithms to detect faint point sources remain to be determined in a future work.
• 3.  
  By comparing the ML and GT maps, variable sources can be detected automatically in ASM in either case (increasing and/or decreasing) of fluxes.

Conversely, the ML method generally loses the physical properties of the source images, such as the energy of each photon and the statistical significance. Thus, one should consider such an ML-based image as a starting point for finding possible interesting sources. These sources would be chosen for subsequent detailed analysis based on the standard analysis tool released by the Fermi-LAT collaboration. Detailed modeling of DGE, as well as evaluation of the TS, is essential even after the application of ML methods presented in this paper.

In this context, the D3PO inference algorithm, which is one of the latest modeling techniques, demonstrated excellent noise reduction. However, analysis with complex models is required, which is labor-intensive and computationally expensive. Our ML-based method, on the contrary, is model independent and can be completed by simple implementation, as mentioned above. While the method is cost-effective, ML maps lose physical properties, such as energy spectrum.

Finally, we discussed the variable source detection in comparison to existing methods. For source variability in the 4FGL catalog, the variability index, such as TS_var (Nolan et al. 2012), was calculated from light curves using 1 yr bin data. For bright and variable sources, light curves are automatically updated daily for cataloged or well-known sources. In contrast, our ML methods have the potential to automatically detect serendipitous variable sources. However, as mentioned earlier, our ML-based technique loses any physical properties of the sources; hence, a detailed analysis is necessary. In this context, FAVA, which is a simple and effective method for detecting variable sources, was employed. However, as FAVA detects the deviation between the map expected from long-term observations and the one from short-term observations, the S/N does not change. Our algorithms improved the S/N around the point source; hence, the variable source detection has the potential to achieve higher accuracy than FAVA.

In this study, an initial 1-week map obtained with Fermi-LAT in 2019 was prepared as a test sample. However, the prediction method is applicable not only to gamma-ray sky maps but also to other maps obtained in different energy bands, such as radio, optical, and X-ray ASMs. Moreover, one can set any time interval and duration depending on the research purpose, astrophysical phenomenon, and source. For example, one can use a timescale that is as short as minutes to detect various transient sources, such as GRBs, without human bias at any point by applying ML algorithms. Longer timescales of weeks to months are more favorable for searching tidal disruption events (Burrows et al. 2011; Gezari et al. 2012) and/or supernovae.

In this research, a 1-week map was used as an input map to demonstrate the concept and performance of various ML algorithms as applied to Fermi-LAT data. In planned future work, an attempt will be made to find transient sources by setting ASMs generated at different timescales, to widely cover scales of minutes to days as input. Moreover, another interesting possibility is that one can predict longer observation maps if the observations are continued, even after the mission has been completed. For example, if a 1 yr ASM is used as the input, it may be possible to predict a 100 yr ASM from existing data. Clearly, such predictions may or may not be correct, as shown in Section 5; however, it can provide strong motivation for future missions and important scientific aspects to be explored in the future. Similar ML approaches have been proposed to confirm the giant but weak emission structure of diffuse gas, such as the Fermi bubble (Su et al. 2010) and NPS/Loop I (Kataoka et al. 2018), which was proposed using Gaussian mixture algorithms (Müller et al. 2018).

For future research directions, three different approaches are promising. The first is to train the ML algorithm separately for the point-source and DGE regions. In this study, although we used ASM data for training, the DGE and point-source regions were separated to train the ML algorithms. The second approach is to apply other ML algorithms, especially deep-learning algorithms. In this study, U-net and Noise2Noise were considered to be effective for improving the ASM obtained by Fermi-LAT. However, it is necessary to search other ML algorithms, such as the generative adversarial network (GAN) type or Fourier convolutional neural network (FCNN) type, for more practical applications. By using the FCNN-type algorithms, it may be possible to obtain the features represented in the frequency domain. Another approach is to use a large amount of data for training. In this study, only actual observational data were prepared. The amount of training data can be easily increased using simulation data.