Analysing the 21cm signal from the Epoch of Reionization with artificial neural networks

TL;DR

This paper investigates the use of artificial neural networks to invert 21cm signal data from the Epoch of Reionization, aiming to efficiently constrain astrophysical parameters from upcoming observations.

Contribution

It demonstrates that ANNs can effectively reconstruct astrophysical parameters from 21cm power spectra, offering a computationally efficient alternative to traditional methods.

Findings

ANN performance depends on network architecture and data sensitivity.

Including multiple redshifts improves parameter reconstruction.

Thermal noise and sample variance reduce reconstruction accuracy.

Abstract

The 21 cm signal from the Epoch of Reionization should be observed within the next decade. While a simple statistical detection is expected with SKA pathfinders, the SKA will hopefully produce a full 3D mapping of the signal. To extract from the observed data constraints on the parameters describing the underlying astrophysical processes, inversion methods must be developed. For example, the Markov Chain Monte Carlo method has been successfully applied. Here we test another possible inversion method: artificial neural networks (ANN). We produce a training set which consists of 70 individual sample. Each sample is made of the 21 cm power spectrum at different redshifts produced with the 21cmFast code plus the value of three parameters used in the semi-numerical simulations that describe astrophysical processes. Using this set we train the network to minimize the error between the parameter values it produces as an output and the true values. We explore the impact of the architecture of the network on the quality of the training. Then we test the trained network on the new set of 54 test samples with different values of the parameters. We find that the quality of the parameter reconstruction depends on the sensitivity of the power spectrum to the different parameters at a given redshift, that including thermal noise and sample variance decreases the quality of the reconstruction and that using the power spectrum at several redshifts as an input to the ANN improves the quality of the reconstruction. We conclude that ANNs are a viable inversion method whose main strength is that they require a sparse exploration of the parameter space and thus should be usable with full numerical simulations.

Tables

Table 1. Table 1: The mean chi-square values for the EoR model parameter values computed by the ANN from the PS as functions of the values used in the simulations. We show mean chi-square values for four types of data. The first column of this table is the chi-square values when we do not include any noise, in second column we take thermal noise and sample variance into account. We use the 21cm PS at z = 9 𝑧 9 z=9 for both cases. The third column is when we take redshift evolution( z = 9 , 10 , 11 ) z=9,10,11) into account and the fourth column is when we take both thermal noise and sample variance into account and we reduce the number of training data samples.

	${\chi }_{\mathrm{wo}/\mathrm{noise}}^2$	${\chi }_{w/\mathrm{noise}}^2$	${\chi }_{w/\mathrm{noise},\mathrm{zevolution}}^2$	${\chi }_{w/\mathrm{noise},\mathrm{reduced}}^2$
$R_{\mathrm{mfp}}$	0.228	0.258	0.172	0.262
$\zeta$	0.271	0.288	0.168	0.290
$\log \left(T_{\mathrm{vir}}\right)$	0.027	0.038	0.019	0.029

Table 2. Table 2: We compare the mean chi-square value obtained by the ANN with the 1 σ 𝜎 \sigma fractional error obtained by 21CMMC for HERA and SKA observations at z = 9 𝑧 9 z=9 . For comparison, we adopt the noise used in 21CMMC for HERA.

	${\chi }_{\mathrm{SKA}}^2$	${\chi }_{\mathrm{HERA}}^2$	1${\sigma }_{\mathrm{SKA}}$	1${\sigma }_{\mathrm{HERA}}$
$R_{\mathrm{mfp}}$	0.258	0.278	0.178	0.184
$\zeta$	0.288	0.354	0.167	0.220
$\log \left(T_{\mathrm{vir}}\right)$	0.038	0.040	0.024	0.033

Equations

$$\delta T_b(\nu )$$

$$\displaystyle\quad\sim 27x_{{\rm H}}(1+\delta_{m})\bigg{(}\frac{H}{dv_{r}/dr+H}\bigg{)}\bigg{(}1-\frac{T_{\gamma}}{T_{{\rm S}}}\bigg{)}$$

$$\displaystyle\quad\times\bigg{(}\frac{1+z}{10}\frac{0.15}{\Omega_{m}h^{2}}\bigg{)}^{1/2}\bigg{(}\frac{\Omega_{b}h^{2}}{0.023}\bigg{)}\bigg{(}\frac{\Omega_{b}h}{0.031}\bigg{)}[{\rm mK}].$$

$$⟨\delta T_b(\mathbb{k})\delta T_b({\mathbb{k}}^{{}^{\prime }})⟩=(2\pi {)}^3\delta (\mathbb{k}+{\mathbb{k}}^{{}^{\prime }})P_{21}(\mathbb{k}).$$

$$s_i=\sum _{j=1}^n{w}_{ij}^{(1)}{x}_j.$$

$$y_i=\sum _j^k{w}_{ij}^{(2)}{t}_j$$

$$E=\sum _{n=1}^{N_{\mathrm{train}}}{E}_n=\sum _{n=1}^{N_{\mathrm{train}}}\left[\frac{1}{2}\sum _{i=1}^m(y_{i,n}-d_{i,n}{)}^2\right],$$

$$\Delta w_{ij}^{(l)}=-\eta \frac{\partial E}{\partial w_{ij}^{(l)}}=-\eta \sum _{n=1}^{N_{\mathrm{train}}}\frac{\partial E_n}{\partial w_{ij}^{(l)}}$$

$$\frac{\partial E_n}{\partial w_{ij}^{(2)}}$$

$$\frac{\partial E_n}{\partial w_{ij}^{(1)}}$$

$$\mathrm{RMSE}=\sqrt{\frac{1}{N_{\mathrm{train}}}\sum _{i=1}^{N_{\mathrm{train}}}{X}^2}$$

$$\delta P_{\Delta T}^{21}(k,\theta )\approx \sqrt{\frac{1}{N_c}}\frac{A_e{x}^{2}y}{{\lambda }^2{B}^2}{C}^N(k,\theta ),$$

$$C^N(k,\theta )={\left(\frac{{\lambda }^{2}B{T}_{\mathrm{sys}}}{A_e}\right)}^2\frac{1}{B{t}_{\mathbf{k}}},$$

$$\delta P_{\Delta T}^{21}(k)={\left\{\sum _{\theta }{\left[\frac{1}{\delta P_{\Delta T}^{21}(k,\theta )}\right]}^2\right\}}^{-1/2}.$$

$$P_N^{ij}(k)=P^i(k)+N^j(k)$$

$${\chi }^2=\sqrt{\frac{1}{N_{\mathrm{test}}}\sum _{i=1}^{N_{\mathrm{test}}}\frac{({\theta }_{\mathrm{ANN}}-{\theta }_{\mathrm{true}}{)}^2}{{\theta }_{\mathrm{true}}^2}}$$

$$\frac{\partial E_n}{\partial w_{ij}^{(l)}}=\frac{\partial E_n}{\partial s_i^{(l)}}\frac{\partial s_i^{(l)}}{\partial w_{ij}^{(l)}},$$

$$\frac{\partial E_n}{\partial s_i^{(l)}}=\sum _k\frac{\partial E_n}{\partial s_k^{(l+1)}}\frac{\partial s_k^{(l+1)}}{\partial s_i^{(l)}},$$

$$\frac{\partial s_k^{(l+1)}}{\partial s_i^{(l)}}=w_{ki}^{(l+1)}{\varphi }^{\prime }(s_i^{(l)}).$$

$${\delta }_i^{(l)}=\sum _k{\delta }_k^{(l+1)}{w}_{ki}^{(l+1)}{\varphi }^{\prime }(s_i^{(l)}).$$

$$\frac{\partial E_n}{\partial w_{ij}^{(l)}}={\delta }_i^{(l)}{t}_j^{(l-1)}.$$

Full text

††pagerange: Analysing the 21cm signal from the Epoch of Reionization with artificial neural networks–Analysing the 21cm signal from the Epoch of Reionization with artificial neural networks

1 Introduction

According to the standard hierarchical structure formation theory based on the $\Lambda$CDM cosmology model, density perturbations grow into small haloes where the first stars and galaxies form (e.g. Abel et al., 2002; Bromm et al., 2002; Yoshida et al., 2006; Bromm, 2013). These astrophysical sources produce ultraviolet (UV) ionizing photons that escape into the intergalactic medium (IGM) which is eventually reionized. The epoch when this process occurs is called the “Epoch of Reionization” (EoR). The detailed astrophysics of the EoR is currently poorly understood because no observation has yet probed the early and middle stages of the EoR although current observations have provided fruitful information on late stage of the EoR. For example, the absorption spectra of high-$z$ quasars indicate that reionization was complete by $z\sim 6$ (e.g. Fan et al., 2006) and the number density of Ly$\alpha$ emitter galaxies at $z=6-7$ implies that the neutral hydrogen fraction increases at $z>6$ (e.g. Konno et al., 2014). On cosmological scales, reionization induces Thomson scattering of CMB photons off free electrons. The optical depth of Thomson scattering measured by Planck is $\tau_{e}=0.058\pm 0.012$ which corresponds, using an instantaneous reionization toy model, to a reionization at $z=8.8\pm 0.9$ Planck Collaboration et al. (2016).

To further improve our understanding of the EoR, observations targeted on the cosmological redshifted 21cm signal from the EoR are on-going. The 21cm signal emission is due to the hyperfine structure of neutral hydrogen atoms and is expected to be a powerful tool to probe the neutral IGM, yielding both astrophysical and cosmological information such as matter density fluctuations and the ionization state and thermal history of the IGM at high redshift Furlanetto et al. (2006); Pritchard & Loeb (2012).

Recently, some first-generation radio interferometers have been attempting to detect statistically the 21cm signal from the EoR, such as the Murchison Wide field Array (MWA)  Tingay et al. (2012), the LOw Frequency ARray (LOFAR) Rottgering (2003) and the Precision Array for Probing the Epoch of Reionization (PAPER) Jacobs et al. (2015). These observational efforts have resulted in upper limits of the 21cm power spectrum and on the ionized state of the IGM at $z\sim 8$  (e.g. Parsons et al., 2014; Ali et al., 2015; Pober et al. 2015, ; Greig et al., 2016; Beardsley et al., 2016). Furthermore, future instruments such as the Square Kilometre Array (SKA) Mellema et al. (2013) and Hydrogen Epoch of Reionization Array (HERA) DeBoer et al. (2016) are designed to detect the 21cm signal power spectrum with higher signal to noise ratio and at higher redshift, during the Cosmic Dawn. The SKA should also be able to image the signal in 3D, which requires high sensitivity on small enough scale (that is sufficient collecting area in a large enough core).

So we expect a wealth of 21cm data in the near future. Then we face the fundamental question: what we can learn from the data? From theoretical and numerical works, we have already some insights on the process of reionization. For instance, considering galaxies in relatively massive host haloes results in larger and more uniform ionized bubbles and imprints a peak at larger scales in the 21cm power spectrum (e.g. McQuinn et al., 2007; Iliev et al., 2012). On the other hand, abundant and small minihaloes, which serve as absorption systems by self-shielding from ionizing photons, result in small and disjointed ionized regions (e.g. Ciardi et al., 2006). To extract information on the EoR from the observed 21cm signal, we need to be able compute the 21cm signal from the basic physics of reionization (e.g. Semelin & Iliev, 2015). The most self-consistent method is to run numerical simulations. There are currently two approaches. The first is to implement the full radiation hydrodynamics (RHD). This type of simulations is relevant if small scales, where feedback effects such as photo-heating play an important role, are resolved (e.g. Ciardi et al., 2000; Ricotti et al., 2002; Hasegawa & Semelin, 2013). The second approach is to run large scale simulations where radiative transfer (RT) is computed in post-processing (e.g. Iliev et al., 2006; Trac & Cen, 2007; Iliev et al., 2014). This type of simulations are able to account for large scales fluctuations in the signal but they require sub-grid modelling to treat processes on galaxy scales (feedback, ionizing photons escape fraction, etc.).

Although numerical study based on simulations is the most consistent method, it implies a large computational cost. This disadvantage is somewhat alleviated by using a semi-numerical approach instead of RHD or RT simulations. In most semi-numerical simulations, the production of ionizing photons is calculated based on the excursion set formalism, using analytically derived halo mass functions, and density fluctuations are obtained with linear perturbation theory Mesinger & Furnaletto (2007); Santos et al. (2010); Mesinger et al. (2011). However, the growth of HII regions is simply evaluated by considering the balance between the production of ionizing photons and the number of neutral hydrogen atoms Furlanetto et al. (2004). Recently, the effect of recombination is also taken into account (e.g. Sobacchi & Mesinger., 2014). The results produced by the semi-numerical approach shows good agreement with those by RT simulations on large scale ($\gtrsim$ Mpc) Zahn et al. (2007, 2011); Mesinger et al. (2011).

To maximize the scientific return of the upcoming observations it is important to establish systematic procedures to derive constraints on the EoR modeling parameters from the observed 21-cm data. Recently, several works have studied how such constraints can be obtained by exploring the EoR parameter space with semi-numerical simulations. For example, Fisher analysis Pober et al. (2014); Kubota et al. (2016); Shimabukuro et al. 2016b ; Liu & Parsons (2016) and Bayesian parameter inference such as the Markov Chain Monte Carlo (MCMC) approach Greig & Mesinger (2015); Bernardi et al. (2016) have been applied.

In this work, we suggest a new approach for parameter reconstruction based on a machine learning method. Machine learning is one of the hot topics in data science as a method to deal with big data. It is currently applied to many fields such as pattern recognition or search engine. The main purpose of machine learning is to find approximate functions that, given the input produce the desired outputs. This is achieved by “learning" from training datasets with known inputs and outputs (e.g. Ball & Brunner, 2010; Bloom & Richards, 2012).

Recently, machine learning methods have been applied in the field of astronomy. For example, learning from huge galaxy image catalogs helps with morphological classification of galaxies Folkes et al. (1996); Lahav et al. (1996); Banerji et al. (2010). Machine learning can also help with the selection and classification of transients  Morii et al. (2016); Mukund et al. (2016). Applying machine learning to a large sample of spectroscopic and photometric galaxies data can improve the accuracy of estimates of photometric redshifts Collister & Lahav (2004); Vanzella et al. (2004); Samui & Samui Pal (2017). Using simulated gravitational wave data with noise as learning sample, machine learning can help us search for gravitational wave signals from noisy real data Kim et al. (2015). In the context of cosmology, machine learning is used to model galaxy formation Kamdar et al. 2016a ; Kamdar et al. 2016b or to make templates of nonlinear matter power spectrum Agarwal et al. (2012, 2014). Closer to our field of interest, there is a study applying a machine learning method to estimate the escape fraction of ionizing photons during the EoR Jensen et al. (2016). In this study, they show how machine learning can estimate the Lyman continuum escape fraction by using mock spectroscopic simulation data. In our work, using a simple astrophysical parameterization of the EoR, we apply artificial neural network (ANN), which is one of the machine learning methods, to reconstruct the parameter values from the 21cm power spectrum data.

This paper is organised as follows. In section 2, we introduce the cosmological 21cm signal and the EoR parameterization we focus on. In section 3, we describe our artificial neural network and test the impact of its chosen architecture. In section 4 we show our main results, and we give a summary and discussion in section 5. Throughout this paper, we employ the best fit cosmological parameters obtained by Planck Collaboration et al. (2016).

2 Cosmological 21cm signal

2.1 Introduction of 21cm power spectrum

The brightness temperature for the $21$ cm signal is given by (e.g. Mellema et al., 2013):

[TABLE]

Here, $T_{\rm S}$ and $T_{\gamma}$ represent the local spin temperature of the IGM and the CMB temperature, respectively. $\tau_{\nu_{0}}$ is the local optical depth in the 21cm rest frame frequency $\nu_{0}=1420.4~{}{\rm MHz}$, $x_{\rm H}$ is the local neutral fraction of the hydrogen gas, $\delta_{m}({\mathbb{x}},z)\equiv\rho/\bar{\rho}-1$ is the evolved matter overdensity, $dv_{r}/dr$ is the local gradient of the gas velocity along the line of sight and $H(z)$ is the Hubble parameter. All quantities are evaluated at redshift $z=\nu_{0}/\nu-1$. As we can see, the 21cm signal includes both astrophysical and cosmological information. Thus, we can hope to use the 21cm signal to disentangle and quantify them Pritchard et al. (2015).

Let us now introduce the power spectrum of the 21cm fluctuations. We define the 21cm power spectrum as

[TABLE]

In our context, we use the dimensional 21cm power spectrum, $k^{3}P(k)/2\pi^{2}$.The 21cm PS would describe the statistical properties of the 21cm fluctuations perfectly if they were a gaussian random field. However, the 21cm fluctuations are expected to deviate from a gaussian behaviour because of astrophysical effects such as ionization and X-ray heating. Thus it is useful to compute higher order statistics and one-point statistics such as the bispectrum, the variance and the skewness to estimate non-gaussian features in the 21cm fluctuations (e.g. Barkana & Loeb, 2008; Harker et al., 2009; Baek et al., 2010; Patil et al., 2014; Shimabukuro et al., 2015; Shimabukuro et al. 2016a, ).

In order to generate the 21cm PS for a given set of astrophysical parameters, we use the publicly available code 21cmFAST Mesinger et al. (2011). This code is based on a semi-numerical model of cosmic reionization and thermal history of the IGM. It quickly generates maps and PS of the brightness temperature, matter density, velocity, spin temperature and ionization fraction at designated redshifts.

We performed simulations in a $200$ ${\rm Mpc}^{3}$ comoving box with $300^{3}$ grid cells for a wide range of EoR parameters described in the next section. In our calculation, we use the 21cm power spectrum in the range $0.06~{}\rm Mpc^{-1}\leq{\it k}\leq 1.4~{}\rm Mpc^{-1}$ divided into 14 bins.

2.2 EoR model parameters

It is common to characterize EoR models with parameters and then examine the effect of changing the parameters on the 21cm signal. We employ three key parameters which are often used. Let us briefly define these three parameters:

1. $\zeta$, the ionizing efficiency: $\zeta$ is the combination of several parameters related to ionizing photons escaping from high redshift galaxies and is defined as $\zeta=f_{\rm esc}f_{*}N_{\gamma}/(1+\overline{n}_{\rm rec})$ Furlanetto et al. (2004, 2006). Here, $f_{\rm esc}$ is the fraction of ionizing photons escaping from galaxies into the IGM and $f_{*}$ is the fraction of baryons locked into stars. These parameters are extremely uncertain at high redshift Gnedin et al. (2008); Wise & Cen (2009). $N_{\gamma}$ is the number of ionizing photons produced per baryon in stars and $\overline{n}_{\rm rec}$ is the mean recombination rate per baryon. In our calculation, we explore a range of $10\leq\zeta\leq 60$.

2. $T_{{\rm vir}}$, the minimum virial temperature of haloes producing ionizing photons: $T_{{\rm vir}}$ parameterizes the minimum mass of haloes producing ionizing photons during the EoR. Typically, $T_{\rm vir}$ is chosen to be $10^{4}{\rm K}$, corresponding to the temperature above which atomic cooling becomes effective. $T_{\rm vir}$ parameterizes the physics of star formation in high redshift galaxies. In haloes with virial temperature $>10^{4}{\rm K}$ atomic cooling is sufficient to trigger gravothermal instability and thus star formation. However, star formation is quenched if AGN or Supernovae feedback is effective and the IGM is heated up. This leads to effective minimum virial temperature larger than $10^{4}$[K]. In haloes with viral temperature $<10^{4}{\rm K}$, hydrogen molecule cooling is necessary. However, if stars begin to form in a halo, radiative feedback such as the photodissociation of $\rm H_{2}$ by Lyman-Werner photons may become effective and prevent the gas from cooling Sobacchi & Mesinger. (2013); Fialkov et al. (2013). Conversely, positive feedback, such as the enhancement of $\rm H_{2}$ molecules formation due to an increase in the free electrons density, tends to push the minimum virial temperature to lower value because cooling becomes more effective. Thus $T_{\rm vir}$ parameterizes the uncertainty in the efficiency of radiative feedback. In our work, we explore $T_{\rm vir}$ ranging from $10^{3}{\rm K}$ to $5\times 10^{5}{\rm K}$.

3. $R_{\rm mfp}$, the mean free path of ionizing photons : The propagation of ionizing photons through the ionized IGM strongly depends on the presence of absorption systems and the sizes of ionized regions are determined by the balance between sinks and sources of ionizing photons (e.g. McQuinn et al., 2011). This process is modelled by the maximum mean free path of ionizing photons, $R_{\rm mfp}$Sobacchi & Mesinger. (2014). Physically, the mean free path of ionizing photons corresponds to the typical distance traveled by photons within ionized regions before they are abosorbed and is determined by the number density and the optical depth of Lyman-limit systems. In our calculation, we explore $R_{\rm mfp}$ from 10 Mpc to 60 Mpc.

3 Artificial neural networks

3.1 Architecture of ANN

In this section, we introduce artificial neural networks (ANN). ANNs are one of the machine learning methods and are a mathematical model inspired by the natural neuron network in our brain. The main purpose of ANNs is to construct approximate functions which associate input data with output data. In order to construct such a function, the ANN has to learn from “ training data”. The architecture of a simple class of ANN consists of three layers: the input layer, the hidden layer and the output layer. Each of them has a number of neurons as shown in Fig.1. In a more general case, we could choose the number of hidden layers and the number of neurons at each layer arbitrarily. In our study, we use 1 hidden layer. Note that it is mathematically proven that neural networks with only 1 hidden layer can approximate any function with any accuracy if we use a large enough number of neurons Cybenko. (1989); Hornik et al. (1989).

Let us briefly describe the architecture of our ANN. The input data $x_{j}$ is fed to the $j$-th neurons in the input layer. Each neuron in the input layer is connected to the $i$-th neuron in the hidden layer and a weight $w^{(1)}_{ij}$ is associated with the connection. The input to the $i$-th neuron in the hidden layer $s_{i}$ is a linear combination of all the input neurons with weight $w^{(1)}_{ij}$ :

[TABLE]

Here, $n$ is the number of input data. In the hidden layer, the $i$-th neuron is activated by an activation function $\phi$ such as its output is $t_{i}=\phi(s_{i})$. Generally, the activation function is a nonlinear function. We use the sigmoid function $\phi(s)=1/(1+e^{-s})$. The properties of the sigmoid function are such that it saturates and returns a constant output when the absolute value of the input is large and that it is a smooth and differentiable function. Thanks to the nonlinear activation function, a trained ANN can express any function.

In the output layer, we compute linear combinations of the activated outputs of the neurons in the hidden layer with weights $w^{(2)}_{ij}$ and obtain the output vector:

[TABLE]

Here $k$ is the number of neurons in the hidden layer. Note that we do not activate the output value. The aim of training the ANN is to find a set of weights that ensures the output vectors produced by the ANN for a set of input vectors is sufficiently close to the desired output vectors. Once we adjust the weights to reach this goal on a training sample, we can make predictions for output vectors for arbitrary input vectors outside of the training sample (for example, new observational data).

A popular algorithm to compute the trained weights is the “Back propagation algorithm” Rumelhart et al. (1986). We will describe this algorithm briefly in the following section.

3.2 Back propagation algorithm

In this section, we present the back propagation algorithm for the 1 hidden layer case. We also show the back propagation algorithm for multiple hidden layers case in Appendix.

In order to quantify how well the output obtained by the ANN approximates the desired output for the training data set, we define the (total) cost function as:

[TABLE]

where $N_{\rm train}$ is the number of training input vectors and $m$ is the number of neurons at the output layer. $y$ and $d$ are outputs of the ANN and the (desired) training output data, respectively. Our purpose is to find the weight set that minimises the cost function. In order to find this weights set, we need to compute the partial derivative of $E$ with respect to the individual weights $w^{(l)}_{ij}(l=1,2)$ and find the local minimum of $E$ using gradient descent.

The weights are updated by gradient descent following the formula:

[TABLE]

Here, $\eta$ is a learning coefficient which controls how fast the weights are updated. We used $\eta=0.01$. We only need to calculate the derivative of the cost function for each training input vector and then sum over all input vectors as shown in eq.6. First, let us consider the derivative with respect to the weights between output layer and hidden layer. In this case ($l$=2), we can simply calculate the derivative of $E$ as

[TABLE]

Next, we calculate the derivative of $E$ with respect to the weights between the hidden layer and the input layer. In this case ($l$=1), the derivative of $E$ is

[TABLE]

In the second line, we use the chain rule for derivative because $E$ depends on the activated neuron $\phi(s)$ in the hidden layer only through the output neuron $y$. Here, $\phi^{\prime}$ denotes the derivative of the activation function with respect to $s$. Using eqs.(6), (7) and (8), we can iterate on the gradient descent until the outputs obtained by the ANN converge to the desired outputs (minimum of the cost function).

The back-propagation algorithm can be summarised as follows:

Starting with random weights, compute the output of the ANN using eq.(3) and eq.(4) for all input vectors in the training set (forward propagation) 2. 2.

Compute the cost function 3. 3.

Compute the derivative of the cost function with respect to the weights between output layer and hidden layer with eq.(7) and then the derivative with respect to the weights between input layer and hidden layer with eq.(8) (back-propagation). 4. 4.

Update the weights with eq.(6). 5. 5.

Go back to (i) and iterate until the cost function converges to a minimum.

3.3 Training data sets and architecture of the ANN in our work

In our case, we prepared 70 training data sets. Each set consist of the 21cm PS $P(k)$ obtained with 21cmFast as input data and the corresponding EoR parameters used in the simulation, $\theta~{}(=R_{\rm mfp},T_{\rm vir},\zeta)$ as output data. The architecture of the ANN is the following: (i) 14 neurons in the input layer, (ii) 14 neurons in the hidden layer, (iii) 3 neurons in the output layer. As we mentioned in section 2.1, we use 14 bins for the 21cm PS, and we use 3 EoR parameters. This is why the number of neurons in the input layer and in the output layers are 14 and 3, respectively. Note that the number of redshifts used to train will always match the number of redshifts being fit.

To analyse in details how the ANN performs we look at partial cost functions: we use the normalised root mean square error (RMSE) defined by

[TABLE]

Here, $X=(\theta_{\rm ANN}-\theta_{\rm data})/\theta_{\rm data}$ ( with $\theta$ one of $R_{\rm mfp},T_{\rm vir}$ or, $\zeta)$. $\theta_{\rm ANN}$ and $\theta_{\rm data}$ are the EoR parameters evaluated by the ANN and from the training data, respectively.

First, we study how the RMSE depends on the number of iteration for each EoR parameter. The result is plotted in Fig.2. In this figure, we fix the number of neurons in the hidden layer at 14 and perform calculation with $10^{6}$ iterations. After initial fluctuations, the RMSE decreases and converges for $>100000$ iterations. Next, we show how the RMSE depends on the number of neurons in the hidden layer after $10^{5}$ iterations. The RMSE is only weakly dependent on the number of neurons in the explored range. In the following section, we perform $10^{5}$ iterations for the back-propagation algorithm with 14 neurons in the hidden layer unless specifically mentioned.

4 Result

We apply the trained ANN to 54 test datasets. As for the training datasets, these test datasets consist of the 21cm PS and the corresponding EoR parameters. We use the 21cm PS data as input for the ANN and obtain estimated values for the EoR parameters. Then, we compare the estimated values with the true values that were used to compute the PS in the 21cmFast simulation, and evaluate how the ANN performs.

First, we apply the ANN to the 21cm PS at a single redshift without including the effect thermal noise or sample variance. In Figs.4 and 5, we plot the EoR parameters estimated with the ANN from the 21cm PS data against the true values. In Fig.4, we show the result at $z=12$. As we can see, $\zeta$ and $T_{\rm vir}$ estimated with the ANN are in relatively good agreement with the true values. On the other hand, for $R_{\rm mfp}$, the correct values are not recovered. This is because the 21cm PS is not sensitive at $z=12$ to $R_{\rm mfp}$. On the other hand, the 21cm PS becomes sensitive to $R_{\rm mfp}$ at lower redshift when reionization advances much more (since $R_{\rm mfp}$ physically expresses the maximum bubble size, the effect of $R_{\rm mfp}$ on the 21cm PS is remarkable mostly at lower redshift). Thus the cost function is insensitive to changes in $R_{\rm mfp}$ values; the learning process of the ANN is incomplete and systematic deviations remain. Note that the reason why the estimated value of $R_{\rm mfp}$ seems to be constant at $\sim 30$ is that $R_{\rm mfp}$=30 occurred more often than other values in the training datasets.

In Fig.5 we show the same plots obtained from the PS at $z=9$. The agreement between recovered and expected values extends over a larger range than at $z$=12, in particular for $R_{\rm mfp}$. Indeed the 21cm PS at $z=9$ is sensitive to $R_{\rm mfp}$ and thus the ANN learning process works better.

In the previous two cases, we considered the 21cm PS without any source of noise. We will now consider both the contribution of thermal noise and sample variance. We model the thermal noise as a gaussian random field characterized by its power spectrum. In an annulus in Fourier space with $N_{c}$ cells, the thermal noise power spectrum is (e.g. Morales, 2005; McQuinn et al., 2006; Yoshiura et al., 2015):

[TABLE]

where $\theta$ is the angle between k and line of sight, $A_{e}$ is the effective area of antenna, $\lambda$ is the observed wavelength and $x,y$ are factors converting $uv$ space units into comoving wavenumber units and are determined by cosmology. $C^{N}$ is expressed by

[TABLE]

where $T_{\rm sys}$ is the total system temperature, $B$ is the bandwidth and $t_{\bf k}$ is the effective observing time. The spherically averaged noise can be obtained by summing over $\theta$ in a shell with radius $k$ shown in eq.12.

[TABLE]

We recommend to read McQuinn et al. (2006). In our case, we assume SKA specification Dewdney. (2013). We estimate the sample variance from 10 simulations using different realizations of the initial conditions.

In the general case, we would include the noise at the level of the visibilities and compute the resulting noisy power spectrum. Assuming that the noise is a random gaussian field, we simply generate a noisy power spectrum by using the following formula:

[TABLE]

where $P^{ij}_{N}(k)$ is the noisy power spectrum, $P^{j}(k)$ is the power spectrum produced by the simulation for a given set of parameters generically labeled by the index j, and $N^{i}(k)$ is a random draw from a Gaussian probability distribution with variance equal to the thermal noise power spectrum or the sample variance. A different and independent draw is required for each value of $k$. For the learning sample we include only the sample variance. For each parameter set (70 possible values for $i$), we generate 50 realizations of the noise (labeled by $j=1\dots 50$). This means that we use 3500 different PS for the learning sample. This is necessary as we are not aware of a standard technique for directly including an uncertainty in the inputs of an ANN. For the test datasets, we add both thermal noise and sample variance to the 21cm PS, using the same procedure.

In Fig.6, we show the EoR parameters found by the ANN as functions of the values used in the simulations, using the 21cm signal PS at $z=9$, including both thermal noise and sample variance. The difference between Fig.5 and Fig.6 is not obvious at a glance. In order to quantify the difference, we compute the mean chi-square value, $\chi^{2}$:

[TABLE]

$\theta_{\rm ANN}$ is the $R_{\rm mfp}$, $T_{\rm vir}$, $\zeta$ reconstructed by the ANN and $\theta_{\rm true}$ is the value of the corresponding parameter used in the simulation. In table.1, we compare the $\chi^{2}$ values for each of the parameters, with and without noise. As we can see, the $\chi^{2}$ values for the 21cm PS with noises are worse than those without noises. Noise alters the efficiency of the learning process for the ANN.

Until now we have used as input data for the ANN the 21cm PS at a single redshift. We now expand our input datas to take the redshift evolution of the 21cm PS into account. We use the 21cm PS, which includes thermal noise and sample variance, at three redshifts $z=9,10$ and $11$ as input datas for the training and test datasets. Thus each set as 14$\times$3 input values and the number of neurons in the input layer is also 14$\times$3=42. We also increase the number of neurons from 14 to 42 in the hidden layer. Note that the number of training datasets is unchanged. In Fig.7, we show parameter values computed by the ANN as functions of the values used in the simulations when the ANN is using the 21cm PS at three different redshifts. Compared with the single redshift case, it seems that the scatter is smaller. We quantify this in table.1, and check that the scatter is indeed smaller for all three parameters. The reason why using the 21cm PS at multiple redshifts improves the accuracy of the ANN results is simply that we increase the information the ANN has to work with. At each redshift, the dependence and sensitivity of the 21cm PS on the parameters is different. Therefore, the trained network can learn much more if we consider the redshift evolution of the 21cm PS. We also notice that the chi-square values obtained by using multiple redshifts and including noise are smaller than those for a single redshift without any noise.

We next check the influence of the size of the training sample. If we reduce the number of training datasets, we expect that the learning efficiency of the ANN will decrease. To check this, we train the ANN with a subsample of the training datasets ($N_{\rm train}$=20). We show the results in Fig.8. As expected, we can see that the scatter of reconstructed parameter values is larger when using fewer training datasets. This is quantitified in table.1. In particular, the scatter of $R_{\rm mfp}$ increases more than for the other parameters. This is connected to the fact already established that the PS is not very sensitive to $R_{\rm mfp}$. Thus the size of the learning sample is crucial to make accurate EoR parameter predictions from observed data.

5 Summary & Discussion

In this paper, we applied artificial neural network (ANN) to estimate EoR parameters from the 21cm power spectrum (PS) . We used 21cmFAST to produce 21cm PS for different values of the following parameters: $R_{\rm mfp},\zeta$ and $T_{\rm vir}$. We ran $70$ simulations, that provided us with $70$ training datasets made of the 21cm PS and the corresponding EoR parameter values. With these datasets, we trained the ANN, which consists of an input layer, one hidden layer and an output layer. We achieved the training by applying the back-propagation algorithm.

In order to estimate the quality of the learning process of the ANN, we first show the root mean square error (RMSE) for the training datasets depends on the number of iterations of the back-propagation algorithm and on the number of neurons in the hidden layer. We found that the back-propagation algorithm converges as the number of iteration increases and weakly depends on the number of neurons in the hidden layer.

We next applied the trained network to test datasets made of the 21cm PS at a single redshift including or not thermal noise and sample variance. As expected, we found that the chi-square values for the ANN-recovered values of the parameters as functions of the true (simulation) values are larger when including thermal noise and sample variance. We also found that the $\zeta$ and $T_{\rm vir}$ values obtained by the ANN were a better match to the true values than for $R_{\rm mfp}$ due the different sensitivity of the 21cm PS to the different parameters.

Then we increased the information in the inputs given to the ANN by including the PS at three different redshifts. We found that the accuracy of parameters computed by the ANN was improved compared with those obtained based on the PS at a single redshift. This is simply that the neural network takes advantage in its learning process of the fact that the input hold more information and that the 21cm PS is sensitive to different parameters at different redshifts.

One of the remarkable points of our study is that we can recover the EoR parameters from the power spectrum using the ANN with good accuracy although we only use 70 training datasets. This number is quite small compared with that required for MCMC (e.g. $10^{5}$ realizations of the PS are used in 21CMMC Greig & Mesinger (2015)). In table.1, we show the mean chi-square value for each parameter calculated by the ANN. To enable a quantitative comparison with 21CMMC, we perform a parameter estimation with the ANN for HERA and SKA observations and compare to the $1\sigma$ uncertainties obtained with 21CMMC (for HERA, we actually use the noise data provided in 21CMMC package). We show the result in table. 2. Using a single redshift, the accuracy of our reconstruction is of the same order as that of 21CMMC, although somewhat worse. Remember however that it does improve when using several redshifts, as does the accuracy of 21CMMC. However, note that the MCMC method finds the local shape of the likelihood space, while the ANN method only finds the average error over a swath of parameter space. So, while the ANN method requires a far sparser exploration of the parameter space, 21CMMC does provide more detailed informations. In addition, we need to mention the effect of the foregrounds. The foregrounds wedge may hinder the use of part of the power spectrum. Thus taking this into account would likely affect the constraints on the $R_{\rm mfp}$ because it depends heavily on the available range of $k$-space. We defer an ANN analysis incorporating the effect of foregrounds to a future work. Both the training process for the ANN and the parameter estimation from the observed data using the trained ANN have a very small computing cost compared to that required for producing training datasets.

In this work, we used the 21cm PS generated by 21cmFAST. However, it is worth applying the ANN to other EoR modelling methods including full simulations in order to check its ability to reconstruct the astrophysical parameters. In our model, we focused on the EoR and used the 21cm PS. Thus, a possible development is to expand the model beyond the EoR and include the Cosmic Dawn. During the Cosmic Dawn, X-ray heating of the IGM plays an important role. Thus, it is important to incorporate X-ray heating models and study their impact on the 21 cm Baek et al. (2010); Mesinger et al. (2013); Ewall-Wice et al. (2016); Fialkov et al. (2017). Recently, Cohen et al. (2016) explored the behaviour of the 21cm global signal in a large parameter space from $z=40$ to $z=6$. We believe it is worth applying the ANN method using the 21cm global signal as an input and try to reconstruct parameter characteristic of the Cosmic Dawn and EoR both.

Another possible application of ANNs is to use the 21cm tomography as an input. Then the learning sample is made of light cone data from 21cm simulations instead of the 21cm PS. Although the 21cm PS is a powerful tool, the $21$ cm signal is not gaussian, and thus the PS does not contain all information. The 21cm light cone data on the other hand holds all the information that will ever be accessible trough observations. It also contains irrelevant information that the ANN will have to filter out: the random phases of the initial density field.

If we consider a large survey field of view and depth, that is a large box size for the learning sample simulation, we will need to treat large amount of data with the ANN. Machine learning method have been developed to treat such big data. One of the methods that have received attention is deep learning (e.g. Schmidhuber, 2014; LeCun et al., 2015). In our case, we only used 1 hidden layer in the ANN. A deep learning neural network consists of one input layer, multiple hidden layers and one output layer. Deep learning has been applied successfully to image recognition. Therefore, one of the our future works will be to apply deep learning to 21 cm light cones from full simulations data.

Appendix

Here, we show the back propagation algorithm in the case of multiple hidden layers. By analogy with eq.(8), we can express the derivative of $E_{n}$ with respect to the weight in the $l$-th layer, $w_{ij}^{(l)}$ $(l=1,2,...L-1)$ as

[TABLE]

where $s_{i}^{(l)}=\sum w_{ij}^{(l)}t^{(l-1)}_{j}$. Since the changes in $s_{i}^{(l)}$ are transmitted to $E_{n}$ through each $s_{k}^{(l+1)}$ in the ($l$+1)-th layer, the derivative of $E_{n}$ with respect to $s_{i}^{(l)}$ can be expressed as

[TABLE]

Remember that $s_{k}^{(l+1)}$ can be expressed with the activation function $\phi$ as $s_{k}^{(l+1)}=\sum w_{kj}^{(l+1)}t^{(l)}_{j}=\sum w_{kj}^{(l+1)}\phi(s_{j}^{(l)})$,

[TABLE]

Here, we define $\delta_{i}^{(l)}\equiv\partial E_{n}/\partial s^{(l)}_{i}$, then we can re-write eq.(16) as

[TABLE]

Combining eq.(18) with $\partial s^{(l)}_{i}/\partial w^{(l)}_{ij}=t_{j}^{(l-1)}$, the eq.(15) can be re-written simply as

[TABLE]

This form tells us that we can easily obtain the derivative of the cost function with respect to $w_{ij}^{(l)}$, which connects the neuron $j$ in the ($l$-1) th layer to the neuron $i$ in the $l$ th layer, as the product of $\delta_{i}^{(l)}$ and $t_{j}^{(l-1)}$. As shown in eq.(18), we start to calculate $\delta_{i}^{(l)}$ from the output layer ($l$=L) to the input layer. This is why this algorithm is called “back propagation”. If we use eq.(5) as the cost function, we easily calculate $\delta_{i}^{(L)}=y_{i}-d_{i}$.

Acknowledgement

This work is benefited from a grant from the French ANR funded project ORAGE (ANR-14- CE33-0016). We thank to G.Mellema, A. Fialkov, S. Majumdar, S.Giri and K. Hasegawa for their useful comments and thank to S. Yoshiura for providing the thermal noise data.