Bayesian decomposition of the Galactic multi-frequency sky using probabilistic autoencoders

All-sky observations of the Milky Way show both Galactic and non-Galactic diffuse emission, for example from interstellar matter or the cosmic microwave background (CMB). The different emitters are partly superimposed in the measurements, partly they obscure each other, and sometimes they dominate within a certain spectral range. The decomposition of the underlying radiative components from spectral data is a signal reconstruction problem and often associated with detailed physical modeling and substantial computational effort. We aim to build an effective and self-instructing algorithm detecting the essential spectral information contained Galactic all-sky data covering spectral bands from $\gamma$-ray to radio waves. Utilizing principles from information theory, we develop a state-of-the-art variational autoencoder specialized on the adaption to Gaussian noise statistics. We first derive a generic generative process that leads from a low-dimensional set of emission features to the observed high-dimensional data. We formulate a posterior distribution of these features using Bayesian methods and approximate this posterior with variational inference. The algorithm efficiently encodes the information of 35 Galactic emission data sets in ten latent feature maps. These contain the essential information required to reconstruct the initial data with high fidelity and are ranked by the algorithm according to their significance for data regeneration. The three most significant feature maps encode astrophysical components: (1) The dense interstellar medium (ISM), (2) the hot and dilute regions of the ISM and (3) the CMB. The machine-assisted and data-driven dimensionality reduction of spectral data is able to uncover the physical features encoding the input data. Our algorithm is able to extract the dense and dilute Galactic regions, as well as the CMB, from the sky brightness values only.


Introduction
The interstellar medium (ISM) is a key element of the Milky Way and subject to both astrophysical and cosmological studies. It consists of localized components such as molecules and interstellar dust in cold clouds, atomic and ionized hydrogen, and hot plasma in the Galactic halo, as well as components which pervade the entire ISM like cosmic rays and magnetic fields. Our present knowledge about the existence of these components is based on the fact that they all contribute to the interstellar radiation field (e.g., Draine 2011). Each component, or the interplay of several components, generates radiation of a specific spectrum and can be reconstructed by component separation algorithms of varying complexity. Some components show very characteristic emission lines such as CO or neutral (HI) and ionized (HII) atomic hydrogen, which permits their Galactic distribution to be determined very precisely. Other components, however, can generate radiation distributed over completely opposite areas of the electromagnetic spectrum, making the component separation task more sophisticated: For example, the interaction of cosmic rays and magnetic fields generates synchrotron radiation and can be observed in the radio regime, as in the 408 MHz radio map (Haslam et al. 1982), while the interaction of cosmic rays with interstellar matter imprints in the γ-ray regime due to hadronnucleon collisions producing pions (e.g., Mannheim & Schlick-eiser 1994). Another example is hot ionized plasma, which radiates in the X-ray regime when generated by supernovae, whereas hot molecules ionized by collisions emit in the UV regime (e.g., Ferriere 2001). In this specific case, both the UV and soft X-ray photons are again absorbed by interstellar dust, which prevents observations in regions of high dust density like the Galactic plane. This interplay also shows the high complexity of radiative extinction, which needs to be taken into account when reconstructing single emission components.
A famous component that is not subject to dust extinction, but still difficult to measure is the non-Galactic cosmic microwave background (CMB). The reason is that the CMB is superimposed by Galactic foregrounds. Cosmological studies aim to extract the CMB radiation from multiple frequency channels by identifying and systematically removing these Galactic foregrounds. In frequencies below 100 GHz, Galactic synchrotron and free-free emission contaminate the CMB, while above 100 GHz, thermal dust emission and the cosmic infrared background dominate (Gold et al. 2011;Planck Collaboration 2016b). To distinguish between different sources of emission, members of the Planck Collaboration (2018b) developed several component separation algorithms, for example the Commander (Eriksen et al. 2008), Sevem (Leach et al. 2008), or Smica code (Delabrouille et al. 2003;Cardoso et al. 2008). The Article number, page 1 of 25 arXiv:2009.06608v1 [astro-ph.IM] 14 Sep 2020 A&A proofs: manuscript no. NEAT_VAE results obtained from those algorithms contain, among others, all-sky maps of the CMB, synchrotron, free-free, thermal and spinning dust, and CO line emission (Planck Collaboration 2016b). These approaches are however based on cosmological, astrophysical and instrumental parameters, and require preprocessed spectral templates or explicit knowledge about physical correlations. To verify these results in a manner independent of such assumptions, approaches that allow automated component identification, like machine learning techniques, have been increasingly pursued in recent years (Longo et al. 2019).
A broad range of machine learning algorithms was applied to cosmological and astrophysical problems (Fluke & Jacobs 2020). Here we just name a few: Beaumont et al. (2011) employed the Support Vector Machine (SVM) algorithm to classify structures in the ISM. The algorithm was able to identify a supernovae remnant behind a molecular cloud based on a sample of manually classified data. Ucci et al. (2018a) examined the composition of the ISM of dwarf galaxies by processing the available spectral information in a machine learning code. The so-called Game algorithm was trained on a large library of synthetic data (Ucci et al. 2018b) and recovered ISM properties such as metallicity and gas densities and their respective correlations on the basis of spectral emission lines. By training a convolutional neural network on synthetic spectra, Murray & Peek (2019) were able to decompose the thermal phases of neutral hydrogen HI. This selection of algorithms represents supervised learning approaches, which require labeled or pre-classified data in order to be trained. We want to investigate to what extend unsupervised approaches can be applied to Galactic observations. One unsupervised machine learning approach which automatically identifies relevant features within some input data is called representation learning (RL) (Hinton 1990;Bengio et al. 2013;Goodfellow et al. 2016). Given an observation, RL methods aim to extract the underlying causes which generated the data, in other words they learn the most informative representation of the observation (Bengio et al. 2013;Goodfellow et al. 2016). This can for example be achieved by reducing the dimension of the input data using a neural network to the so-called underlying, explanatory factors of variation (Hinton & Salakhutdinov 2006;Bengio et al. 2013;Goodfellow et al. 2016). This concept can be translated to the task of astrophysical component separation by investigating which underlying, data-generating components can be extracted from Galactic all-sky data. These components are constructed to encode mutually independent features of the data and are often used as the input for subsequent analyses. Especially, considering the large quantities of astrophysical data produced every day by surveys such as the Sloan Digital Sky Survey (SDSS), an effective preprocessing strategy needs to be developed to be able to analyze the vast amount of data (Kremer et al. 2017;Reis et al. 2019).
In the present study, we apply RL to a data set of 35 Galactic all-sky maps recorded in multiple frequencies provided by Müller et al. (2018). On their data set, the authors learned spectral pixel-wise relations using Gaussian Mixture Models (GMMs), which permitted to augment pixels with missing measurement data. They verified the pre-stated hadronic and leptonic component of the γ−ray sky (Selig et al. 2015) and presented a higher resolved hadronic component as well as a completion of non-observed information. However, the main component maps of the GMM did not not capture different astrophysical environments. This motivated the approach in the present study to explore other latent variable models, like autoencoders (Hinton & Salakhutdinov 2006), which are able to encode useful representations of the data in their latent space.
Based on the provided data set, we combine generative modeling with variational inference to learn a lower-dimensional representation of the data, which we call features. Our model approximates the posterior probability on these features, and we efficiently optimize the approximation using a state-of-the-art variational autoencoder (VAE) (Kingma & Welling 2013). Such an approach was, for example, used in the context of the previously mentioned Sloan Digital Sky Survey, as Portillo et al. (2020) successfully applied variational autoencoders to the SDDS data. The authors efficiently reduced the dimensionality of 1000 input pixels to six latent components, while the VAE was able to outperform principal component analysis considering the spectral dimensionality reduction. Their latent space separates types of galaxies or detects outliers, making it a very useful preprocessing step for large astrophysical data. However, the authors claim that the uncertainty quantification could be improved and suggest to include the pixel-level noise as a separate feature to improve latent variance representations.
Our variational autoencoder is based on the principles of information theory, meaning that we follow a Bayesian approach to track all relevant uncertainties and we enable our model to adapt to the introduced model noise. Our mathematical derivation of the loss function is, although based on the specific signal reconstruction problem, a general approach: Starting with a generative data model, we are able to explain all terms occurring in a classical VAE's loss function, but in addition we provide further terms which clearly result from our calculations in Sect. 2 and, simultaneously, deliver robust results in Sect. 3.

Data
Observations of the Milky Way can be visualized by all-sky maps showing the sky brightness in a certain frequency range. When we combine data from Galactic all-sky records in multiple frequencies, we can obtain a more complete picture of our Galaxy, but we also gain redundant information. Our goal is to determine a reduced representation of the observed sky, which contains only non-redundant or essential information. For this analysis, we use the aforementioned data set consisting of 39 Galactic all-sky maps distributed over the entire electromagnetic spectrum compiled by Müller et al. (2018). To generate the data set, the authors assembled information from all-sky surveys ranging from γ-ray to radio frequencies (Müller et al. 2018, and references therein). They incorporated all-sky data with at least 80 % spatial coverage in HEALPix format (Gorski et al. 2005) and used the highest resolution map for each frequency. The UV regime is not part of their data set, since the respective GALEX survey (Bianchi et al. 2017) shows too many unobserved regions. To homogenize and unify the data, they (1) converted the sky brightness values in the original maps to flux magnitude values, (2) reduced the noise level in X-ray data by smoothing with a Gaussian kernel, (3) removed extra-Galactic sources and calibration artifacts from all maps except the already cleaned γ-ray data from (Selig et al. 2015), and (4) unified the resolution of all maps to nside = 128.
Based on the resulting 39 Galactic all-sky maps, we build a data set D (see Table A.1) for the present study with three modifications: (1) we discard the hadronic and leptonic component maps by Selig et al. (2015), since these are derived from the Fermi data and thus contain redundant information, (2) we remove the AKARI far infrared map recorded at 65 µm, which exhibits a poor spatial coverage (we can only work with pixels that are covered by all maps), and (3) we neglect the CO line emission map, since this information is partly contained in the Planck 100, 217, and 353 GHz frequency channels (Planck Collaboration 2016b). We then define our data set consisting of k = 35 all-sky maps as an indexed set D = d 1 , . . . , d p , where p is the number of pixels and the magnitude vectors d i ∈ R k represent the magnitude flux values of all frequency maps for the i th pixel 1 .

Model design
The determination of a non-redundant and lower-dimensional representation of the frequency information in our data set is an inverse problem: We are looking for unknown quantities, our so-called features, and procedures applied to them that could have generated the full information in D. To describe this problem, we use generative modeling to define a data model with corresponding parameters Θ. The solution to the inverse problem is then given by the posterior probability distribution P(Θ | D) of the model parameters Θ, which we specify in the next paragraph. Using Bayes' theorem, this posterior can be expressed by the data likelihood and the prior distribution up to a Θ-independent factor, reading P(Θ | D) ∝ P(D | Θ) × P(Θ). We perform a pixel-based analysis of image data by assuming magnitude vectors of different pixels to be independent. This allows us to factorize the data likelihood and prior distributions. In the following, we will calculate the respective distributions per pixel and finally combine all derivations to one solution for the entire data set D.
We start by defining a generative process that leads from an abstract source S to the observations in the data set D. We do not know S at this point, but we aim to associate S a posteriori with some relevant physical quantities. We assume that each d i ∈ D is generated from a source vector s i ∈ R l with l ≤ k, and the collection of these vectors is building up a source set S = (s 1 , . . . , s p ). This relation is expressed in a data model where we define the observed data D to be composed of D ( d 1 , . . . , d p ), which is the output of a generative process G : R l → R k , and some model noise N model . The generative process G( · ) mapping the variables of S pixel by pixel to D, i.e. D = G(S), is unknown.

Prior distributions
We factorize the prior distribution of the set of source vectors S as P(S) = p i=1 P(s i ). The prior on P(s i ) can be arbitrarily complex, but without loss of generality we can find a transformation s i = T (z i ) of the source vectors to a set of latent vectors z i ∈ R l that provides a standardized prior distribution P(z i ) = G(z i , 1). This coordinate transformation into the eigenspace of the 1 In other words: i is the magnitude flux value of the λ th sky map at pixel i.
Approximate inference function with neural network In our mathematical derivation, we first define a generative data model and afterwards build an inference process using information theoretical methods. The combination of both processes allows the resulting algorithm to perform in the opposite direction: It uses its input d i to infer a latent representation z i of the data, which is (among other factors like priors) led by the generative process aiming to regenerate the initial data vector d i from the latent space again. The minimization objective, or loss function, directing the algorithm's learning process is Eq. (8).
prior, also known as random variate generation using inverse transform sampling (e.g., Devroye 1986), or reparametrization trick (Kingma & Welling 2013;Rezende et al. 2014;Titsias & Lázaro-Gredilla 2014), allows us to absorb all complex and unknown structure of the source space in the transformation T ( · ), and provides us with an easy to calculate unit Gaussian prior distribution. Using this definition, we can rewrite the generative process as a parametrized function of the latent variables d i = G(s i ) = G(T (z i )) f θ (z i ), with θ denoting the parameters of the transformed generative forward model.
We further define the model noise N model = (n 1 , . . . , n p ) as a set of p noise vectors n i ∈ R k and assume the pixel-wise noise is independent and identically distributed with a Gaussian distribution of zero mean and noise covariance N ∈ R k×k , meaning N). Later in the discussion we will see that the noise covariance induces a metric in data space, which is between two magnitude vectors d i and d i . N thus indicates how accurately the data set d i can be reconstructed using the latent variables z i . Since this accuracy is unknown, we introduce N as an inference parameter. We again perform a transformation to the prior eigenspace of the form N = t ψ (ξ N ) with the latent noise parameter ξ N ∈ R, being distributed as P(ξ N ) = G(ξ N , 1) with the transformation parameter ψ. For our specific application, we choose a log-normal mapping with mean µ N and standard deviation σ N . In other words, this is a diagonal noise covariance matrix N with transformation parameters ψ = (µ N , σ N ) T and a single global noise parameter ξ N , which we aim to learn.
Combining the definitions of the noise and the generative model, we can rewrite the data model from Eq. (1) as a pixel-wise expression: At this point, we can define the data model parameter vector Θ described in Sect. 2.2. It consists of the latent vectors z i , the latent noise parameter ξ N , which indirectly defines the model noise n i , and the parameters of the generative model θ. The parameters θ are the parameters of the neural network used to approximate the generative process, and are specified in further detail in Sect. 2.6. Since we do not have prior knowledge on these forward model parameters, we assume a uniform prior distribution on θ (see Table 1). Thus we have Θ = (Z, θ, ξ N ) with the latent space set Z = (z 1 , . . . , z p ).

Data likelihood
We can include our forward model in the pixel-wise likelihood by marginalizing over the model noise: where we assumed the noise n i to be a priori independent of Z and θ.

Approximating the posterior distribution
Combining the prior distributions listed in Table 1 and the data likelihood in Eq. (4), the posterior distribution P(Θ | D) reads: where we used P(θ) = const. We are not able to calculate expectation values from this high-dimensional probability distribution, but we can use variational inference (e.g., Blei et al. 2017) to approximate the posterior distribution P(Θ | D) with an easier to integrate distribution Q Φ (Θ | D) with variational parameters Φ. In the following, we define a suitable approximate distribution Q Φ (Θ | D) and use the Kullback-Leibler Divergence 2 as a measure to evaluate the dissimilarity of P(Θ | D) and Q Φ (Θ | D).
Assuming that Z, θ and ξ N are a posteriori independent, the approximate distribution Q Φ can be written as (Jaynes 1982). The principle states that if only the mean µ and covariance Σ of some data are known, then the knowledge contained in that data set can best be expressed by a Gaussian distribution with exactly this mean µ and covariance Σ (e.g., Enßlin 2019). Later we will see that by construction µ i and Σ i are functions of the input data, making the choice of a Gaussian distribution for Inserting the derived expressions from the previous sections and using Monte Carlo Methods to approximate integrals with finite sums, we arrive at where we absorb all constant terms into H 0 . The full calculations with all intermediate steps are carried out in Appendix B.

NEAT-VAE
Our goal is to infer a lower-dimensional representation of the data, which in our case is expressed by the latent source space Z. We can achieve this by minimizing Eq. (8), which describes the spurious amount of artificial information introduced by the approximation of P with Q Φ . In the following, we will explain how this objective function can be translated into the framework of a latent variable model called variational autoencoder (Kingma & Welling 2013).
A basic autoencoder learns a low-dimensional, latent representation of higher-dimensional input data. This is achieved by training the autoencoder to reconstruct the original input as accurately as possible from the reduced, latent representation (e.g., Rumelhart et al. 1985;Hinton & Salakhutdinov 2006;Goodfellow et al. 2016). In this context, training describes the minimization of an objective or loss function, for example the mean squared error between input data and the reconstructed output data. The dimensionality reduction of the input to the latent space occurs in the so-called encoder, the latent space itself is called bottleneck layer, and in the decoder, the latent space gets translated back to data space. Variational autoencoders (VAEs) offer a probabilistic framework to jointly optimize latent variable (or generative) models and inference models (Kingma & Welling 2019). Eq. (8) can be transformed to such a VAE framework (see Fig. 1) as follows: -Generative Process: Neural networks are generalized function approximators. Since the exact form of the data generating process f θ (z i ) = d i is unknown, we can use a neural network with input z i and output d i to approximate f θ and use back propagation to optimize the parameters θ of the network, which corresponds to the generative decoder.
We can build a function delivering posterior samples following this variational distribution, expressed by a neural network with input d i and output z i . Let the pixel-wise mean µ i ∈ R l and covariance Σ i ∈ R l×l be determined by a parametrized function of the input data e φ (d i ) = µ i , log Σ i , where φ contains the parameters of the function e φ and the matrix logarithm of Σ i is calculated. This parametrization ensures the variance to be positive and the calculation to be numerically stable, since the logarithmic function maps the small values of the variance to a larger space. By inverse transform sampling, we can then define the posterior latent space variables as z i = µ i + exp 1 2 log (Σ i ) · i = µ i + √ Σ i i with an auxiliary variable i and P( i ) = G( i , 1). In practice we approximate the function e φ (d i ) by the variational inference encoder. We take √ Σ i to be diagonal and describe it by its diagonal vector diag √ Σ i , allowing us to calculate µ i and diag √ Σ i as two distinct outputs of the encoder network.
-Independent representation: Based on the input data, the encoder network delivers latent space variables z i with mean µ i and covariance vector By using this definition we find the optimal independent approximation to the posterior P(Θ | D). This leads to a disentanglement of the input information, meaning each dimension of z i , which is equivalent to the hidden neurons in the bottleneck layer, encodes (approximately) mutually independent features of the data.
The minimization objective in Eq. (8) contains the loss function of a classic VAE with three modifications: (1) We include the size of the latent space l and its corresponding weight (see Appendix B for the derivation) to be able to compare different latent spaces with each other, (2) besides the network weights, we aim to optimize the noise covariance and thus the latent noise variable ξ N . The prior on this latent noise adds an extra term proportional to ξ N 2 to the objective function, and (3) by including noise in the data model, the likelihood contains a factor 1/t( ξ N ) and contributes an additional term tr ln t( ξ N ) from the normalization. Since we expand the VAE framework by the adaption to noise, we name our algorithm NEAT-VAE (NoisE AdapTing Variational AutoEncoder). Using the transformations illustrated before, the final objective function depends only on the generative decoder parameters θ, the variational encoder parameters φ and the latent noise parameters ξ N . We implement an autoencoder architecture (encoder network, bottleneck layer, decoder network) with Eq. (8) as the respective loss function in the PyTorch framework 3 , which is publicly available at https://gitlab.mpcdf.mpg. de/msara/neat_vae. The framework allows us to calculate derivatives of Eq. (8) with respect to θ, φ and ξ N using automated back propagation, and to minimize the loss with a build-in optimizer. The conducted experiments are described in Appendix C.

Results and discussion
Applied to our set of Galactic full-sky observations, the NEAT-VAE framework yields a posterior probability distribution of the latent space variables that capture the essential information in our data. The derived loss function forces the latent variables to be statistically independent of each other, and thereby to represent individual physical components. We obtained the posterior by simultaneously optimizing two processes: an inference mechanism reducing the observed data D to a lower-dimensional latent space Z, and a generative process mapping these latent space variables back to the higher-dimensional data space D. These two processes, described by artificial encoder and decoder networks, support each other during training: the decoder reconstructs the data space based on the latent variables the encoder delivers. This reconstruction is constantly compared to the input data by the likelihood term in Eq. (8), ensuring the encoder to adapt the inference function and thus to provide improved latent space variables to the decoder. In addition, we aim to find latent space variables that encode mutually independent features of the input data. In the following, we analyze the resulting latent space variables z i ∈ Z and their full-sky representations. The dimension of the latent space Z corresponds to the number of hidden neurons in the bottleneck layer. From here on, we call the subsequent hidden neurons 'features' and the resulting full-sky representations 'feature maps'.

Dimensionality reduction
We first analyze how the dimension of the latent space in the NEAT-VAE correlates with the reconstruction quality of the generative process. We quantify the reconstruction by the mean squared error of input maps and reconstructed maps with the number of pixels p, the 35-dimensional data vectors d i ∈ D and the corresponding reconstruction vectors d i ∈ D.
We observe that only three features are required to achieve an MSE of below 0.1, which describes an average deviation of the reconstructed values compared to the input values of a natural logarithm-based flux magnitude. For small values, the absolute uncertainty on logarithmic scale equals the relative uncertainty on linear scale, meaning MSE = 0.1 corresponds to a relative uncertainty of ≈ 10%. The MSE decreases for a further growing number of features, and stagnates around a value of 0.02, that is a relative uncertainty of ≈ 2% at approximately ten features. We interpret this as an indication for a high redundancy of the information contained in the 35 Galactic all-sky maps, since increasing the number of features that are able to encode the input data beyond this point do not increase the quality of the reconstruction any further.

Morphology of features
Based on the experiments with different latent space dimensions, we examine the spatial structures in the feature maps in more detail. We recognize some spatially correlated structures of the input data in the feature maps. We assume this to be a meaningful result, since the autoencoder only learns correlations among the magnitude flux values d i within one data pixel and is not informed about spatial structures. We also observe feature maps with the same morphology to occur in latent spaces of different dimensions. To investigate whether there is a pattern or an order among the feature maps with respect to information content, we calculate the significance of each feature for the reconstruction of the data by using the following measure where z feature map, i is the intensity value at the i th pixel, z feature map is the mean intensity value of a single feature map, and σ 2 feature, i denotes the posterior variance as calculated by the algorithm at the i th pixel. This means we calculate the feature map variance, which describes the fluctuations within the posterior mean map, weighted by the posterior feature variance averaged over all pixels p. The significance thus expresses the ratio of the magnitude of fluctuations within a map compared to the uncertainties of the map. In this context, a high significance marks the features that the autoencoder is most certain about to be required for the reconstruction, while a feature significance below 1 corresponds to an insignificant feature, as the posterior uncertainty is larger compared to the posterior mean values of the map. Averaged over all the experiments we conducted (see Appendix C), we call the three most significant features A, B and C.
For the configuration shown in Fig. 3, the features have significance values S featureA = 1.67 × 10 3 , S featureB = 1.41 × 10 2 , and S featureC = 4.68 × 10 1 . The remaining features have significance values ranging from 8.75 to 3.82 × 10 1 and encode artifacts and other morphologies of the input data, as displayed in Appendix D. A similar behavior was also observed by Müller et al. (2018): The authors build a Gaussian mixture model (GMM) to reconstruct observations from a certain frequency range based on data of complementary frequencies. The GMM components used in their study also encoded artifacts of the input data when the number of components was increased. However, posterior samples of the latent space of our algorithm show that the NEAT-VAE does not always assign information to each and every feature: Starting from twelve features and adding further neurons to the latent space, the significance of the added features drops below one. This means that the posterior variance of a feature map is greater than the fluctuations within the map, and the resulting posterior samples show white noise statistics. On average, ten to eleven features are significant throughout our experiments with varying latent space size. We assume that from this point on, our algorithm has identified all mutually independent features of the input data. When further dimensions are added to the latent space, the algorithm 'tunes out' those degrees of freedom by making them insignificant.
In the next sections, we focus on the three most significant features of the configuration with in total ten latent space features (since from this point on the reconstruction does not change significantly, see Fig. 2), and their physical interpretation. Visual inspections of the other features (Appendix D) indicate that the separation into independent components is not fully reached, possibly due to the finite amount of optimization.
Morphological structures similar to the North Galactic Spur or the Fermi bubbles imprint onto several features simultaneously, as well as measurement artifacts.
The posterior mean values in Fig. 3 reflect the internal representation of the three most significant latent space features by our algorithm. We observe that the overall sign of the maps changes for varying hyperparameter configurations. We did not observe the change of sign to follow a specific rule or to depend on the sign of other features from the same hyperparameter configuration. Hence we assume the overall sign to not have any physical interpretation and only to depend on the initialization of the network parameters. We however observe that the relation of positive and negative values within single maps is constant throughout different choices of hyperparameters, that is the overall feature map structure. Since these signs do not change the morphology of the map, we chose the sign of the color code in the map display of Fig. 3 according to that of the input map which most closely resembles it. This is positive for stronger emitting regions and negative for weaker emitting regions.
We can investigate which data information a feature map encodes by considering the generative process of the NEAT-VAE and using the back propagation algorithm: Since the amount of parameters in the latent space is much smaller compared to the data space, the autoencoder is forced to extract the essential information of the input data in order to generate it again. The generative process therefore offers, at least approximately, an explanation for the information flow within the autoencoder, and we can visualize to which extend specific features are generating the individual Galactic all-sky maps using sensitivity analysis (e.g., Zurada et al. 1994). Here, we compute the gradients of the reconstructed output maps with respect to the feature maps using the back propagation algorithm. The values of the resulting decoder Jacobian are displayed in HEALPix pixelization in Appendix E.
Another measure for quantifying which input information is encoded by the features is the mutual information. The mutual information I(X; Y) = x,y P(x, y) ln P(x, y) P(x) P(y) can be calculated from two-dimensional histograms of feature maps (x) and output maps (y). For a given number of bins, the joint probability distribution P(x, y) is represented in an (x, y)-shaped matrix by counts per bin, while the marginalized distributions P(x) and P(y) are obtained by summing over the respective y and x-axes of this matrix. In our interpretations, we will use both measures to evaluate the encoded information in the latent space (mutual information of feature and input maps), as well as the generative process from the latent space (decoder Jacobian maps), in order to determine the physical content of the features A, B and C.

Identifying the information encoded by feature A
Feature A, which is the most significant feature in 98% of the examined hyperparameter configurations (see Appendix C), is displayed in Fig. 3a. From a visual analysis, we recognize a positive color-coded Galactic plane and negative color-coded Galactic poles in the posterior mean. In the eastern part of the Galactic plane, we see a bright, circular structure in the Cygnus region. Further in the east, south of the Galactic plane, structures similar to the Perseus region occur. The circular shaped structure north of the Galactic center resembles the Ophiuchus region and southwestern of the Galactic center structures similar to the Small and Large Magellanic clouds can be recognized. In the western part of the Galactic plane, the bright structures look like the Orion region. The posterior variance shows a high certainty in the region of the Galactic plane and the structures of the surrounding latitudes, while at the southern and northern Galactic poles, the uncertainty increases.
Interpretation. We find compelling evidence that feature A traces the dense and dusty parts of the ISM: Based on calculations of the mutual information with 512 bins, the top three data sets contributing to feature A are the AKARI far-infrared 140 µm input data with mutual information I(X; Y) = 1.91, the IRIS infrared 100 µm input data with I(X; Y) = 1.84, and the AKARI farinfrared 160 µm input data, also with I(X; Y) = 1.84. This frequency band of the interstellar radiation field is dominated by the infrared emission of dust (e.g., Draine 2011, p.121). Feature A shares the least mutual information with the ROSAT X-ray 1.545 keV input data with I(X; Y) = 0.20. Next, we analyze the decoder Jacobian maps (see Appendix E), which display the gradients of the reconstructed Galactic all-sky maps with respect to latent space features in HEALPix format. We observe the following: The lower and mid-latitudes of the γ-ray regime strongly depend on feature A, and going from higher to lower energies, the overall-dependence on feature A grows. Especially, the lowenergy regime of the Fermi data is dominated by hadronic interactions of cosmic rays with the interstellar medium (ISM) and thus shows the gas distribution (Selig et al. 2015). In the X-ray regime, we have small to no dependence; however, the mid-and high latitudes of the soft X-ray data (0.212 keV and 0.197 keV) show negative gradients toward feature A. Here, negative gradients state that when increasing values in feature A, the corresponding values in the reconstructed X-ray maps will decrease. Such an inverse-proportional behavior is very plausible by the physical context: Radiation from X-rays is extincted by cold interstellar gas (Ferriere 2001), thus the absence of X-ray emission reveals regions where interstellar matter is present. The H α map at 656.3 nm, which displays emission due to hydrogen transitions occurring from the second excited state n = 3 to the first excited state n = 2 (e.g., Finkbeiner 2003), positively depends on feature A. The physical interpretation is that neutral hydrogen preferably resides in the dense ISM, as traced by feature A. We observe a very strong dependence of the Galactic planes of the 12 and 25 µm infrared data on feature A, while with further decreasing energies, the infrared data down to 545 GHz show an overall, positive dependence. In this regime, emission from thermal dust is observed (Klessen & Glover 2014; Planck Collaboration 2018a). With further decreasing energy, the Planck microwave data depends strongly on feature A, in particular the Galactic plane. The 21 cm line emission shows an overall dependence on feature A and describes the total neutral atomic hydrogen column density, since the displayed emission occurs due to the transition between two levels of the ground state of atomic hydrogen (Ewen & Purcell 1951). In the synchrotron radio regime at 1420 MHz and 408 MHz, both feature A and feature B play an important role in determining the emission structures, which we will address in Sect. 3.6. The positive dependencies of sky maps on feature A describe areas of our Galaxy with high density of interstellar matter,  In this case, the algorithm was trained to reduce its input of 35 Galactic all-sky maps to ten features in the latent space with the configuration described in Fig. 2. The mean feature maps (left panels) are displayed in order of significance according to Eq. (10), meaning these features show the highest ratio of feature fluctuations to feature uncertainties within the latent space (see Sect. 3.2 for details). The colors in the posterior mean of feature C (left panel in (c)) are inverted for illustration purposes. The grey pixels correspond to missing values in the input data.
while the negative gradients correlate with extinction of emissions by interstellar matter. We assume that positive gradients mark the pixels in latent space that are used to generate the corresponding pixels in data space, while negative gradients mark an anti-correlation of feature and data pixels. On the basis of this specific combination of gradients and the mutual information, we infer that feature A encodes dense regions of the ISM. Discussion. We can test our hypothesis by investigating the correlation of feature A with dust as a tracer for the ISM (Kennicutt & Evans 2012). Fig. 4b shows the relationship of feature A and the thermal dust emission (Fig. 4a) as calculated by the Planck Commander code (Eriksen et al. 2008;Planck Collaboration 2016b). The posterior mean of feature A has a strong, positive, and linear correlation with the thermal dust component. Both the Commander and the NEAT-VAE algorithm perform a Bayesian, pixel-wise analysis based on a data model containing the linear sum of a signal function and noise, but with two main differences: First, we only employ statistical information in our prior knowledge, while the Commander priors include detailed physical models for the various emission processes contributing to the radio to far-infrared sky, as well as calibration and correction factors, and a prior for the CMB. Second, the Commander algorithm has 11 million free parameters to tune (Planck Collaboration 2016b), while our model has just about 5, 000. The algorithms are not directly comparable, since the Commander code was especially developed to separate the Galactic foregrounds to reconstruct the CMB, while the NEAT-VAE only seeks to find an essential representation of the input data. But in this context, the NEAT-VAE summarized the input data into categories, one of which already contains dust emission. With this result, we assume it is possible to derive the dust component based on feature A with little computational effort. We investigate another tracer for interstellar matter, namely the hadronic γ-ray component derived by Selig et al. (2015), shown in Fig. 4c. It represents the γ-ray emission due to the interaction of cosmic ray protons with interstellar matter, and is positively correlated with feature A, see Fig. 4d. This correlation is reasonable, since Selig et al. (2015) composed the hadronic component of the low-energy γ-ray maps, on which feature A also depends on. These results meet our initial aim of finding a reduced representation of the input data that combines redundant information in one feature, in this case tracers for dense regions in the ISM. The high significance of feature A is likely to result from the choice of input data, since most of the maps in our data set D represent emission from the dense interstellar matter. We did not include the CO line emission data, but from the positive gradients of the Planck 100 -353 GHz channels with respect to feature A, which include the information of the CO line emission (Planck Collaboration 2016b), we assume that the Galactic plane of feature A most likely would generate the CO data. This would support our interpretation of feature A, since CO is a tracer for molecular interstellar gas (Scoville & Sanders 1987), and thus, for regions of high density.

Identifying the information encoded by feature B
Feature B, the second most significant feature in 98% of our experiments, is displayed in Fig. 3b. We see a negative colorcoded equator which resembles the Galactic plane, with positive color-coded bulges north and south of the plane. Especially the northern bulge structure looks similar to the morphology of the North Polar Spur. The positive color-coded, circular structure in the western part of the Galactic plane lies in the Vela region, whereas in the east, the location of the Cygnus region appears negative color-coded with a positive, ring-like structure surrounding it. The uncertainty of this map appears to be low in all latitudes, but is an order of magnitude higher compared to the variance map of feature A.
Interpretation. Feature B shows highest mutual information with the X-ray input data of ROSAT at 0.885 keV (I(X; Y) = 0.82), 0.725 keV (I(X; Y) = 0.78) and 1.145 keV (I(X; Y) = 0.67), and lowest mutual information with the Planck 70 GHz data (I(X; Y) = 0.18). In general, X-ray emission is assumed to arise from hot, ionized gas in regions of the ISM with low density (e.g., Ferriere 2001). Considering the gradient maps in Appendix E, there is little to no dependence of the reconstructed γ-ray data on feature B. Starting from X-ray data at 1.545 keV, the positive gradients get stronger with decreasing energy, and especially the bulge-like area of the reconstructed X-ray data strongly depends on feature B. The soft X-ray regime around 0.25 keV, which coincides with hydrogen cavities (Sanders et al. 1977;Snowden et al. 1990, e.g.,), also shows a weak, positive dependence on feature B. The infrared and microwave regime again show little to no dependence, we however observe small negative gradients to be present where in the corresponding gradients maps of feature A, strong, positive gradients occur. The reconstructed 1420 MHz and 408 MHz radio maps show positive dependence on feature B, especially in the bulge-like areas. This data set detects the synchrotron radiation generated by the interaction of cosmic ray electrons with magnetic fields in the ISM (Ginzburg & Syrovatskii 1965), and the intensity of this radiation depends on the density of relativistic particles and the magnetic field strength. The former one is predominantly found in the hot areas of the ISM, if only for the much larger volume occupation of this phase within the Milky Way (Cox 2005). Thus we assume that feature B encodes the enhanced presence of such cosmic ray electrons. With exception of the low energy X-ray and radio synchrotron regime, we observe the data predominantly generated by feature A to have little to no dependence on feature B, which we interpret as feature B being complementary to dense regions of the ISM. In combination with the indications from the mutual information with X-ray input data and positive gradients, we assume that feature B encodes tracers for dilute and hot regions of the ISM.
Discussion. Fig. 5, showing the correlation between Feature A and B, supports our interpretation: The features, in our case the dense and dilute regions of the ISM, are basically uncorrelated. This relationship is reasonable in the sense that the features describe two distinct categories of the ISM: The data generated by feature A are based on emissions from cold and warm ionized gas (as observed in the 21 cm and H α line emissions), dust, and interactions with interstellar matter (cosmic rays). The ionization of the warm gas occurs mainly due to photo-ionization by O and B stars, the hottest and most massive stars of the Milky Way. Due to their their ratio of mass and luminosity, these stars have a short main sequence life cycle and are thus found near their initial birthplaces, the dense ISM (e.g., Blome et al. 1997, p.60). Feature B, however, encodes data which represent radiation of an even hotter medium, the hot ionized plasma, which is generated by supernovae explosions (e.g., Kahn 1980). The radiative processes encoded by features A and B can thus be traced back to two fundamentally different origins, namely emission from interstellar matter and from stellar explosions. One regime of the electromagnetic spectrum missing in our input data is the ultraviolet (UV) frequency band. In this regime, the emission of very hot gas which gets ionized by collisions can be observed, but UV radiation does not penetrate dense regions of the ISM and is thus mostly absorbed in low and mid-latitudes (Ferriere 2001). Due to this extinction, it is likely that UV radiation would be interpreted by the NEAT-VAE as redundant with the soft X-ray data, which shows similar dust absorption patterns. This would support our interpretation of feature B to encode the hot and dilute ISM.

Identifying the information encoded by feature C
The third most significant feature in 74% of the investigated configurations, Feature C, is displayed in Fig 3c. Here, we inverted the colors of the posterior mean such that the fluctuations resemble the color-coding of the CMB, see Interpretation. The mutual information of feature C is highest with the Planck data of 70, 100 and 143 GHz with values of I(X; Y) = {1.01, 0.95, 0.82}, respectively, and the least mutual information is observed with the ROSAT X-ray 1.545 keV data with I(X; Y) = 0.12. According to the decoder Jacobian maps in Appendix E, feature C mostly generates the Planck 30 − 217 GHz channels, especially the 70 and 100 GHz data, in which the CMB can be observed (Planck Collaboration 2018a). All other reconstructed maps show little to no dependence on feature C, with exception of a weak, positive dependence of the soft ROSAT X-ray data on feature C that we will address in Sect. 3.6. Based on the mutual information and the strong positive gradients in the microwave regime around 100 GHz, we assume feature C to encode the CMB.
Discussion. Considering our physical interpretation of features A and B, it is reasonable that the CMB emission is encoded in a separate feature, since the underlying physical processes in the Early Universe shaping the CMB fluctuations are independent of the processes in the ISM. To test our hypothesis we compare feature C to the CMB all-sky map derived by the Planck Commander code in Fig. 6a. Feature C shows a positive, linear correlation (I(X; Y) = 0.88) with the CMB (Fig. 6b). A main difference between feature C and the CMB can be seen in the Galactic plane, which feature C encodes (in addition to the Cygnus region) with high intensities, but also with high uncertainty (see Fig. 3c). We assume the high intensities in the Galactic plane of feature C to originate from the model's architecture: Since we do not take spatial correlations into account, the algorithm learns spectral relations of each pixel independently. The Galactic planes of the input maps are dominated by high intensity values, while for other latitudes, more low-intensity values are present. We assume that this generates two distinct clusters of intensity ranges, leading to two different spectral relations learned by the algorithm. We can also observe this distinction in intensity by analyzing the contributions of pixels for data generation (see Appendix E): The structure of the Galactic plane is recognizable in most of the decoder Jacobian maps, even though the algorithm has no information about spatial correlations. This again indicates that the relations learned in different intensity regimes represent different processes in the NEAT-VAE structure, which however are rather to be associated with internal computations than with physical properties.

Non-physical interpretation of the NEAT-VAE
We finally examine the decoder Jacobian maps which do not have a clear physical interpretation. For example, we observe feature A (the dense ISM) contributing to the 1420 and 408 MHz radio data, as well as feature C (the CMB) to supposedly generate the soft X-ray data. The NEAT-VAE algorithm recombines all features of the latent space in a highly non-linear way to generate the output maps, meaning that the gradients in Appendix E cannot always be considered independently or in a linear way. Especially when the absolute gradient values of the reconstructed maps are large for more than one feature, a holistic analysis is required. In most of the displayed cases in Appendix E, there is one feature dominantly generating the reconstruction of input data. One exception, however, is the radio data, where we observe positive gradients with respect to both feature A and B. A speculative physical explanation is that feature A marks regions of enhanced magnetic field strength and feature B of enhanced relativistic electron densities, both quantities that in combination determine the synchrotron emission. However, we rather assume that these two features have no physical meaning for synchrotron emission, but the internal computations of our algorithm run most efficiently when latent values of features A and B are used for radio data generation. One always has to consider the fact that our algorithm has no information about physical relations and primarily optimizes a statistical function. Another example for such a behavior can be found in the soft X-ray data, where we observe strong negative gradients with respect to feature A and weak positive gradients toward both features B and C. This mixing of features in a partially contrasting manner indicates the nonlinearity of the function mapping from the latent space values to the reconstructed output data: We hypothesize that some values might be used by the algorithm just to tune some others out, since the decoder has to regenerate the data based on all available features. This discussion shows that the decoder Jacobians give a valid first interpretation of the strongest correlations, but have to be analyzed with care. Not each contribution to the gradient represents a meaningful relationship, especially since all gradients are entangled and cannot be considered independently.

Conclusions
In summary, we derived a probabilistically motivated machine learning framework, the NEAT-VAE, which successfully identifies the most significant sky emission components according to pixel-wise sky brightness values across the full electromagnetic spectrum. At this task, the algorithm performs computationally efficient and in a fully unsupervised manner. The three most significant resulting sky components express physical relationships in the considered data set and can be assigned to emission processes of the dense ISM, the hot and dilute ISM, and the CMB.
We achieved this performance by developing a Bayesian formulation of the component separation problem for astrophysical data, which serves as the minimization objective, or loss function, for a state-of-the-art variational autoencoder. Using Bayes' theorem, we combined a generative data model with a variational inference process, incorporated a very limited amount of prior knowledge of generic nature (e.g., independence of features, Gaussian statistics of typical noise processes), and approximated unknown functions by neural networks. The resulting algorithm is able to group essential information contained in a set of Galactic all-sky observations into mutually independent features. This property results from the algorithm's architecture and the diagonal covariance approximation we chose: First, an autoencoder maps its input to a lower dimen-sional latent space, from which it aims to reconstruct its input again. By reducing the dimension of the data in this so-called information bottleneck, the autoencoder is forced to learn a useful representation of the input data. By additionally stating a diagonal covariance matrix in the variational approximation of the latent space posterior distribution (see Sect. 2.5), we find the latent spaces' optimal independent approximation. We observe that there is an order among the latent space features based on their significance to encode the data set. This computational significance correlates with physical significance, with the most significant features having a clear physical interpretation and representing emissions from the dense ISM, the hot and dilute ISM, and the non-Galactic CMB, respectively. Our interpretations are based on the analysis of (1) the mutual information of the features and the input data, (2) the generative properties of the features to reconstruct the input data, and (3) the features' correlations with other astrophysical quantities. Thus we were able to verify that the NEAT-VAE algorithm detects the most informative components of emission into which the Galactic input data can be decomposed. However, as discussed in Sect. 3.6, the analysis of our generative process is not always straightforward and has to be performed with care, since the examined gradients might represent internal computation strategies of the algorithm rather than physically meaningful generative properties.
The relevance of our work lies in the insights we obtained in the context of representation learning (RL). In RL, the reduced representation of high-dimensional data is often used as a preprocessing step in order to perform the actual analysis more efficient. For example, the task of deriving the thermal dust emission is computationally very expensive. Our feature A, however, which encodes the dense ISM, is highly correlated with the thermal dust component. We hypothesize that a component separation algorithm only analyzing the information contained in feature A, or using feature A as a starting point, is likely to perform faster and more efficient compared to analyzing a larger data set with redundant and entangled information.
The loss function of the NEAT-VAE and the associated hyperparameter tuning mainly determine the performance of the algorithm. Although the loss function used to direct the learning is based on very general assumptions, the Bayesian framework in which it was derived describes the initial problem probabilistically. This allows us to track uncertainties and evaluate the significance of each feature. Hyperparameter tuning can be reduced by including those parameters in the Bayesian inference process, as in our case the noise covariance of the model noise. In general, the joint optimization of an inference and a generative process is a state-of-the-art approach to disentangle data, which we showed to be very successful for unsupervised learning applied on astrophysical observations. Especially the fact that astrophysical data are a superposition of several radiative processes, and thus are rich in relationships between emitting as well as absorbing components, makes them a very suitable target for decomposing algorithms.
The fact that we neither incorporated spatial correlations nor physical priors in the learning process shows how much information can be retrieved from data-driven approaches: Alone the sky brightness values of single pixels in Galactic all-sky maps are sufficient to detect emission from the dense ISM, the dilute ISM and the CMB. By additionally combining Bayesian methods with neural networks, we were able to track uncertainties and efficiently approximate unknown functions. This work shows that machine-assisted, Bayesian signal inference can be performed with neural networks, and that the minimization of a loss function based on information theoretic principles reveals meaningful relations in astrophysical data.

Inserting the definitions for
where H 0 collects all terms independent of the parameters. We rename the subsequent terms of the Kullback-Leibler Divergence as follows and calculate each expression respectively.
First energy term T1: The trace of a scalar is the scalar itself, resulting in ( Finally, the trace is a linear operator and therefore commutes with the expectation: In the upper calculation, we carry the term l (1 + ln (2π)) through our computations, since we change the number of latent space features l in our different experiments.
Second energy term T2: Third energy term T3: The expression z i z T i can be rewritten as z i z T i = (z i − µ)(z i − µ) T − µµ T + z i µ T + µz T i . We insert this in the line above: Fourth energy term T4: +tr(ln t ψ ( ξ N )) tr (ln t ψ ( ξ N )) + C 4 (B.5)