A Machine Learning Approach to Integral Field Unit Spectroscopy Observations: III. Disentangling Multiple Components in H ii regions

The synthetic data is created using the ORB – Outil de Reduction Binoculaire – software package described in Martin et al. (2016) and Martin et al. (2012). The create$\_$cm1$\_$lines$\_$model is used to create the following emission lines: [S ii]$\lambda 6717$, [S ii]$\lambda 6731$, [N ii]$\lambda 6584$, H$\alpha$(6563Å), and [N ii]$\lambda 6548$. We use these lines since they are all present in the SN3 (651-685 nm) filter. ORB creates synthetic spectra by taking the spectral resolution, line amplitude, line function, line broadening, and velocity as user inputs. As in other papers, these spectra are representative of those expected from the SIGNALS (Star formation, Ionized Gas, and Nebular Abundances Legacy Survey; Rousseau-Nepton et al., 2019). We randomly sample the resolution from a uniform distribution varying between R=4800 and R=5000; this variation mimics expected resolution variations within the SITELLE field of view, that are observed in SITELLE data cubes.

The velocity value is randomly chosen from a uniform distribution between $-500$ $\mathrm{km\,s}^{-1}\,$ and $500$ $\mathrm{km\,s}^{-1}\,$; these values well encompass the typical rotation velocity range for galaxies found in the SIGNALS program. Additionally, we select the broadening from an uniform distribution between 5 $\mathrm{km\,s}^{-1}\,$ and 20 $\mathrm{km\,s}^{-1}\,$. We selected the lower limit since SITELLE cannot resolve broadening beneath this level (Rousseau-Nepton et al. 2019). We note that shocks in the ISM are not being considered here – they would require an increase in the maximum broadening. We verified that the chosen lower limit does not impact the classification process by testing different lower bounds (3 $\mathrm{km\,s}^{-1}\,$ and 1 $\mathrm{km\,s}^{-1}\,$). We selected 20 $\mathrm{km\,s}^{-1}\,$ since it is the expected upper limit of dispersion in HII regions based onon previous experience in SITELLE analysis (Rhea et al. 2020b). Moreover, we randomly vary the signal-to-noise ratio concerning the H$\alpha$ emission between 5 and 30. We propagate the noise to the other emission lines by adding it to each spectral channel. The noise factor itself is sampled from Gaussian with a sigma of 1. Therefore, each spectrum has a different signal-to-noise ratio described by a Gaussian centered between 5 and 30 with a sigma of 1. The last feature required to create synthetic spectra is the relative amplitude of the strong emission lines. Following the methodology described in detail in Rhea et al. (2020b), relative line amplitude of the HII regions were sampled from the Mexican Million Models database (3MdB; Morisset et al. 2015) BOND simulations (Asari et al. 2016).

We created a set of 10,000 spectra with a single component; therefore, there are five emission lines present ([S ii]$\lambda 6717$, [S ii]$\lambda 6731$, [N ii]$\lambda 6584$, H$\alpha$, and [N ii]$\lambda 6548$) with the same velocity and broadening values. An additional 10,000 spectra with two components was also created; therefore, there are two sets of five emission lines present with the same velocity and broadening values within each set. Each set has a different velocity and broadening value. The amplitudes for each set come from different instances of the BOND simulations from 3MdB. Thus, we have 20,000 spectra in total. Following standard procedure, we use 70% of the synthetic data for the training set, 20% for the validation set, and 10% for the test set (e.g., Breiman 2001). We note that individual lines may have different kinematics; however, we do not consider this in our analysis for simplicity’s sake.

In the past decade, Bayesian statistics have become house-hold tools for astronomers; Bayesian techniques are used to study eclipsing exoplanet signals (Taaki et al. 2020; Ruffio et al. 2018), fitting spectra in order to extract model parameters (Sereno 2016; Sharma 2017), and model comparison (Jenkins & Peacock 2011; Trotta 2007). In this section, we outline the mathematics behind Bayesian inference, how it can be used to compare models, and its practical implementations in python.

Bayesian inference is the process of uncovering underlying parameter distributions through the use of Bayes’ theorem. Bayes’ theorem states

where $\Theta_{M}$ describes a set of model parameters for model $M$, $Y$ is the set of observed data, and $I$ represents any assumed information. The left-hand side of the equation is known as the posterior distribution, and it describes the probability distribution of model parameters given the existing data. On the right-hand side, we have the likelihood distribution, $P(Y|\Theta_{M},I)$, which describes the distribution of the data given a set of model parameter values, the prior distribution, $P(\Theta_{M},I)$, describing our prior knowledge or believes regarding the distribution of the model parameters, and the evidence, $P(Y|I)$ or $Z$. In many scenarios, Bayes’ theorem is rewritten as

where the evidence is ignored. However, in the case of model comparison, we wish to calculate the Bayes factor, which is defined as the ratio of the evidences calculated using two distinct models. Therefore, the evidence must be calculated; note that the evidence is defined as the integral over all the likelihood multiplied by the prior model parameters.

We write the Bayes’ factor, $P(M|D)$, as

We adopt a standard Gaussian likelihood function of the following form:

where $Y_{M}(M_{\Theta})$ is the assumed model evaluated given the model parameters $M_{\Theta}$ and $\sigma_{0}$ are the measurement errors. For the sake of the calculations simplicity (discussed in §4.1), we employ a standard Gaussian to model the emission line. The Gaussian has the following form:

where $A$ is the amplitude of the Gaussian, $x$ is the wavelength channel, $x_{0}$ is the position of the line, and $\sigma$ is the width of the line. We emphasize that we calculate the posteriors for the parameters $M_{\Theta}=[A,x_{0},\sigma]$ in the case of a single Gaussian model. We adopt linearly uniform priors for all the fit parameters.

Unfortunately, the calculation of the evidence is exceptionally costly in high-dimensional spaces, and alternative methods to calculate the evidence are required. Nested sampling is one of the most popular methods for accurately and rapidly calculating the evidence (e.g. Skilling 2006; Chopin & Robert 2010). We, therefore, use the nested sampling code dynesty for our calculations (Speagle 2020). Since the calculations are computationally expensive, we explore a faster algorithm – a convolutional neural network.

Neural networks, and their extension convolutional neural networks, have been used extensively to solve astrophysical problems (e.g. Baron 2019; Davies et al. 2019; Aniyan & Thorat 2017; Kim & Brunner 2017). Recently, convolutional neural networks have been shown to be efficient at rapidly and accurately extracting parameters from spectra (e.g. Fabbro et al. 2018; O’Briain et al. 2021;Keown et al. 2019; Rhea et al. 2020b). For a detailed review of convolutional neural networks, we direct the reader towards the following works: Kuo (2016), Liu et al. (2016), Khan et al. (2020). In particular, by interpreting model selection as a classification problem, neural networks with a softmax activation on the last layer can be used as a substitute to the reference Bayesian model selection approach, with a substantial gain in computational complexity.

For this work, we adopt the structure of the network STARNET developed in Fabbro et al. (2018) and used in our previous work (Rhea et al., 2020b). We note that the hyperparameters (i.e. filters) are tuned separately and are thus different from those used in Fabbro et al. (2018). The network is as follows:

1. 1.

   Convolutional layer with 4 filters of 16 elements activated with the relu function

2. 2.

   Convolutional layer with 4 filters of 8 elements activated with the relu function

3. 3.

   Pooling layer with a size of 8 elements

4. 4.

   Flattening Layer

5. 5.

   20% Dropout

6. 6.

   Fully connected layer with 1000 nodes activated with the relu function and regularized using $\ell_{2}$ regularization

7. 7.

   20% Dropout

8. 8.

   Fully connected layer with 1000 nodes activated with the relu function and regularized using $\ell_{2}$ regularization

9. 9.

   Fully connected output layer with two nodes activated with the softmax function

We note that the use of the softmax activation function in the output layer not only generates a binary classification but also assigns a probability of accuracy to the classification.

Although we adopt the convolutional neural network, STARNET, described in detail in Fabbro et al. 2018, we apply a complete suite of hyper-parameter tuning. Using the sklearn grid search implementation, we optimized the number of convolutional kernels in each layer (4 in both layers) and the length of each filter (16 elements in the first layer and 8 elements in the second layer. We also optimized the number of spectra fed into the network at a time (i.e., the optimal batch size; we found this to be 2 spectra), the best optimization algorithm, Adam, the length of the convolutional pooling in each layer, 8. These results are summarized in table 2.3.

Since the resolution is free to vary uniformly from 4800 to 5000 and the network requires inputs of a constant size, we interpolate the spectra using linear cardinal splines using a grid taken from a reference spectrum with a resolution set to 5000. This does not affect the fidelity of the spectrum.

Following the procedure described in $\S$2.2, we fit a single Gaussian model and a double Gaussian model using the dynesty implementation of static nested sampling to the same test set used in the CNN calculations. We fit only the H$\alpha$ line in order to reduce the computational cost (see §4.1 for a more detailed discussion on the repercussions of this decision). We then used the reported outputs to calculate the Bayes’ factor. The double component model is considered favorable if the Bayes’ factor is greater than 1.08 (this value was experimentally determined to return the best classifications¹¹1In order to determine the cutoff Bayes factor we varied the cutoff from 1 to 5 in increments of 0.05 and calculated the f1-score given the chosen factor. We then took the Bayes factor that returned the best f1-score. ); otherwise, a single component model is favorable. We used the static nested sampling implementation over the more accurate dynamic nested sampling implementation since the latter increased the computational costs by more than an order of magnitude – this is further discussed in $\S$4. The results are shown in figure 1. We investigated the cases in which a double component model was mistaken for a single component model (approximately 15% of the time) and discovered that the majority of misclassifications occurred when the absolute velocity difference between the two components was less than 250 $\mathrm{km\,s}^{-1}\,$ (see the next subsection for further discussion).

The convolutional neural network was trained and validated using 70% and 20% of the synthetic data. At the end of the process, the overall accuracy of the algorithm was 93%. The mse (mean squared error) on both the test and validation set is 0.44; this indicates that the network is not over-fitting. Over-fitting occurs when the network learns the test set too well and cannot generalize to other data (i.e., the validation or test set). The confusion matrix for the test set (see figure 2) reveals near-perfect precision in classifying single line-of-sight component spectra. Moreover, the matrix demonstrates that the network classifies double line-of-sight component spectra accurately 90% of the time; the remaining 10% are misclassified as single components. We note that this is less than the 95% accuracy of the Bayesian methodology.

In order to establish that the network is not over-fitting, we applied the standard k-fold cross-validation with $k=10$. Additionally, we applied a modified $k$-fold cross-validation algorithm in which only the training and validation set are varied while the test set remains consistent across all folds. The reported accuracy of the model for each fold regardless of implementation remained constant (less than 5% variation). This is a strong indication that the network is not over-fitting, and, thus, it can be generalized to unseen data.

Moreover, we examined the regimes in which the network fails to predict the number of underlying components accurately. In figure 3, we plot a histogram of the absolute velocity difference for the test set spectra with double line-of-sight components. The stacked histogram shows that the incorrectly categorized spectra (rose) are primarily clustered around absolute velocity differences less than 250 $\mathrm{km\,s}^{-1}\,$. This trend is expected since, at these low absolute velocity differences, the components are blended considerably. We show four spectra illustrating this issue in appendix A.

In addition to using a convolutional neural network to estimate the number of emission-line components, we tested another common machine learning algorithm for binary classification problems: the random forest. While random forests are themselves complex machine learning algorithms, we do not discuss them in detail here since the purpose of this section is to demonstrate the efficacy of other machine learning algorithms to solve our problem; instead we point the reader to (Breiman 2001) and the references therein. Although we test several configurations for the random forest classifier, our final classifier has a total of 50 independent estimators (decision trees) with a maximum recursion depth of 5 levels. The split criterion is based on entropy calculations instead of gini calculations. As evidenced by figure 4, the random forest algorithm does not achieve the same fidelity in categorization as the neural network.

As evidenced by the confusion matrices shown in figures 1 and 2, the convolutional neural network outperforms the Bayesian inference classifications in the single component case. Notably, the CNN slightly (3%) under-performs the Bayesian Inference analysis for the classification of double line-of-sight components. Since the two methods were tested on the same set of data (the CNN’s test set), the methodologies’ comparison is unbiased with regards to the data they are tested on. Although the Bayesian inference methodology returns similar results to that of the CNN, the method takes two orders of magnitude more computational time, rendering it impractical. The Bayesian method takes approximately 200 hours to fit the test set, while the network takes 45 minutes to train and approximately 1 second to predict the test set. Additionally, the creation of the synthetic data set takes approximately 20 minutes. Therefore, the total time for the CNN is slightly over 1 hour.

Generally, Bayesian Inference sets a baseline for fitting results; however, our results indicate that the CNN is outperforming the Bayesian inference model estimates for the single line-of-sight component scenario since it is more accurate for single line-of-sight component classifications. We attribute this discrepancy to the chosen calculation of the Bayesian evidence and the model used. In order to reduce the number of parameters to a more computationally manageable size, we only fit a Gaussian (or two Gaussians) to the H$\alpha$ component. In doing so, we must marginalize over a 3 or 6 dimensional space, respectively (there are 3 components used to describe each Gaussian in the Bayesian approach). If we were to fit all five lines, we would be required to marginalize over a 15 or 30 dimensional space thus further increasing the required computational time. Therefore, we emphasize that the CNN is outperforming the Bayesian approach only because the Bayesian approach is incomplete since it does not fully treat the ILS of SITELLE nor does it take into account every line present in the SN3 spectra. We further note that a Bayesian approach taking into account all lines and the proper ILS is not implemented since it is too computationally expensive at this time.

We also note the work of González-Gaitán et al. (2019-01-21) in which the authors use spatial information and prior estimations on the spatial correlations to drastically reduce the required computational time using Bayesian methods. However, a direct comparison between our works is beyond the scope of this paper.

In order to port this methodology to observations that are not part of the SIGNALS collaboration, we create an alternative test set at a different resolution. We select a resolution of 3000 to test the algorithm’s efficacy on lower-resolution spectra. After creating a set of synthetic data following the same procedure described in section $\S$2.1 with $R\sim 3000$, we train, validate, and test the algorithm on this lower-resolution set. Figure 5 demonstrates that the categorizations’ fidelity does not change for single line-of-sight component spectra. At the same time, the accuracy is reduced by 10% for double line-of-sight component spectra. Additionally, we followed the same procedure for $R\sim 1000$ spectra; however, the network fails to achieve the same level of efficacy due to the shallow resolution. We, therefore, do not suggest using the network on spectra with a resolution considerably lower than $R\sim 3000$.

Although we have only considered the problem as a binary classification (i.e. either single or double component) due to the constraints of the SIGNALS observing campaign, it is possible that a spectrum can contain more than two components. Therefore, we constructed a set of 1,000 synthetic spectra comprised of three underlying emission components following the same methodology as outlined in $\S$2.1. We then apply our network to the spectra. Figure 6 demonstrates that the network classifies the majority ($\sim$95%) of the three-component spectra as having two components. Therefore, we suggest the following interpretation if a region is suspected to have more than two components: a classification as two components should be considered as at least two components. These results are similar to those found in Rhea et al., 2020a.

In order to demonstrate the utility of this methodology, it is applied to actual SITELLE observations at R$\sim$3000 of the interacting pair of spiral galaxies NGC 2207 and IC 2163 (P.I. R. Pierre Martin). This system (D = 35 Mpc) is undergoing a grazing collision and has been studied from X-rays to cm wavelengths (e.g. Kaufman et al. 2012; Elmegreen et al. 2017). The closest interaction of both galaxies is estimated to have happened 200–400 Myr ago (Kaufman et al. 2012) and, according to models, both galaxies are expected to merge in the future (Struck 2007). SITELLE observations cover the entire interaction, including zones not severely affected by the collision, tidal tails, and the zone of closest proximity between both galaxies. Therefore, this target presents an ideal test-bed for our algorithm. We note that no $R\sim 5000$ SITELLE observations of merging galaxies exist at the time of writing this paper.

In addition to applying the neural network to the SITELLE data cube, we obtain a deep image of the SITELLE cube using ORCS. The deep image is a 2d, stacked representation of the 3d data cube (i.e., the spectral information is compressed to a single pixel). The neural network returns a map of the galaxy where each pixel is designated as having either a single or double component. A comparison of these maps (the component map and the deep image) is shown in figure 7. The component map reveals several interesting features: double components in the bulges of each spiral galaxy, the double component emission in the current closest proximity zone (green circle), and single component emission in the spiral arms. Additionally, the component map classifies diffuse regions (indicated in magenta circles) as single component emission regions. We note that the HII regions in the center of each galaxy are also categorized by double component emission. These findings are consistent with the literature (i.e. Soto et al. 2012). Since we use a softmax activation function in the final layer of the network, we also obtain probability maps. These maps in general can be used to make probability cuts reflecting regions for which the network is sure of its classification. We discuss this in more detail in the appendix C. We hypothesize that since the bulges are being classified with a high probability as having two components these regions have at least two components; this is consistent with our findings in $\S$6. These results indicate that the network’s classifications should be accepted cautiously for regions which are heavily contaminated by astrophysical processes such as AGN or a bright stellar component. In regions that correspond to the deep image background (i.e., not part of either galaxy in the system), the component map is noisy. This indicates that the component map is reliable only for regions of sufficient flux – which is expected. We stress that we do not expect to be able to resolve more than two emission-line components given the resolution of SIGNALS SN3 observations (R$\sim$5000) and the distance to the objects. Expanding the network to be a multi-component classifier is a viable avenue for future research.

Determining the number of resolvable line-of-sight components in galaxy spectra carries several important scientific implications. In situations where two line of sight components are present and resolvable, it is important to fit both lines separately instead of treating the emission as a single component for several reasons. Distinguishing a region of multiple component line-of-sight emission from a large blended region will allow for the correct characterization of thermal broadening in H ii regions which will help untangle the reaction of the gas to the ionizing source. Moreover, if the two components have considerably distinct kinematic parameters (velocity and broadening), treating the emission as a single component will lead to a misclassification of the kinematics of the region which can lead to incorrect interpretations of the star forming mechanisms present. In effect, a misunderstanding of the number of resolvable line-of-sight components present will lead to a misrepresentation of the underlying physics.

Additionally, mapping out regions in which two line-of-sight components is crucial in disentangling galaxy systems undergoing a merger. In doing so, the kinematics of the component galaxies can be studied; moreover, treating each galaxy’s emission separately leads to more accurate calculations of the emission line fluxes. This is crucial for understanding the underlying ionization mechanisms at play in mixing regions of merger systems (e.g., Baldwin et al. 1981; Kewley et al. 2001; Kewley et al. 2019; Rich et al. 2015). Furthermore, by applying the networks developed in Rhea et al. (2020b) and Rhea et al. (2021) to the appropriate regions, kinematic parameters and emission-line ratios can be recovered. The authors note that those networks function only for a single component spectrum; future work will demonstrate their extension to double component spectra.

Although we focus on SITELLE data cubes in this article, the methodology described herein can be readily ported to any other high-resolution integral field unit that has access to the [N ii] doublet, the [S ii] doublet, and H$\alpha$ emission lines such as the Multi Unit Spectroscopic Explorer, MUSE (e.g., Bacon et al. 2010; Kreckel et al. 2019; Kreckel et al. 2020; Foster et al. 2020). When applying this methodology to other instruments, it is essential to repeat the creation of synthetic data and network training, validation, and testing if the instrumental response function is not a sincgauss.