A novel stellar spectrum denoising method based on deep Bayesian modeling

We now develop our proposed model for stellar spectrum denoising based on the basic Bayesian model (1)–(2). Suppose the signal spectrum of a star is s ∈ ^D . It is a D-dimensional unobserved vector, and we want to recover it from noisy observations. Modern astronomical surveys take multiple exposures to get several noisy observations y₁,⋅, y_n ∈ ^D of the clean signal spectrum s. We express the observations in a signal-plus-noise model,

$$\begin{eqnarray}&&{{\bf{y}}}_{i}={\bf{s}}+{\epsilon }_{i},\,\,{\rm{for}}\,\,i=1,\cdots,n,\end{eqnarray}$$

where _i is a D-dimensional noise vector. The noise could have complex correlation structure across pixels. Most of the time, only a single average spectrum is used, we can directly set n = 1 in the above. However, our method is more powerful and can exploit multiple-exposure data where only the red or blue part of the spectrum is recorded in each exposure. Our denoising framework is inspired by the following fundamental but powerful Bayesian model

$$\begin{eqnarray}\begin{array}{ll}{{\bf{y}}}_{1},\cdots,{{\bf{y}}}_{n}|{\bf{s}}\sim & p({\bf{y}}|{\bf{s}}),\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{ll}{\bf{s}}\sim & p({\bf{s}}).\end{array}\end{eqnarray} \tag{ 2 }$$

In the above, p(s) is a prior density encoding the likelihood of the signal s; p(y|s) is the probability density of the observed y_i given the signal vector s. Based on some state-of-the-art deep density estimation methods (see Sect. 3.2), we will construct an expressive prior p(s) and an observation model p(y|s) by neural networks. Given the trained model and the observations y₁,⋅, y_n , the true signal vector s can be inferred from the posterior distribution p(s| y₁,..., y_n ).

Several additional adjustments of the model are necessary. It is not straightforward to model a clean spectrum s ∈ ^D in a high dimensional space ^D . To effectively construct a model for the signal spectrum, we exploit that, for a collection of astronomical spectra, the signal vectors s of various stars typically reside over a low dimensional manifold. This intrinsically low-dimensional structure can be effectively employed for spectrum data analysis. For example, Lawlor et al. (2016) used a local non-linear dimension reduction technique to discover the manifold structure. The low dimensional nature of the data implies that all clean spectra can be represented by a variable in a low dimensional space. Denote this low-dimensional variable by z ∈ ^L , and we can construct a mapping ${\mathcal{G}}$ such that the signal $s={\mathcal{G}}({\bf{z}})$ is the mapped value of z. Based on this mapping, the prior over the signal s can be directly expressed as a prior over the latent space p(z). More specifically, we will take the stellar class C into account and construct the prior p(s|C) conditional on each stellar subclass. The final hierarchical model is summarized as below.

$$\begin{eqnarray}\begin{array}{ll}{\rm{y}}_{1},\cdots,{\rm{y}}_{n}|s\sim & p({\rm{y}}|{\rm{s}}),\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{ll}s= & {\mathcal{G}}({\bf{z}}),\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{ll}{\bf{z}}|C\sim & p({\bf{z}}|C),\end{array}\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{ll}C\sim & p(C).\end{array}\end{eqnarray} \tag{ 6 }$$

In this model, the class label variable C has a uniform prior distribution over all stellar subclasses. The construction and training of the latent priors p(z|C) and the generator function ${\mathcal{G}}$ will be addressed in Section 3.3. The observation model p(y|s) will be discussed in Section 3.4. The final model training and denoising workflow will be summarized in Section 3.5.

Given our trained model in Equations (3)–(6), the posterior distribution of the latent variable z and the class label C is proportional to the joint distribution

$$\begin{eqnarray}&&p({\bf{z}},C|{\rm{y}}_{1},\cdots,{\rm{y}}_{n})\propto \displaystyle \prod _{j=1}^{n}p({\rm{y}}_{j}|{\mathcal{G}}({\bf{z}}))\times p({\bf{z}}|C)\times p(C).\end{eqnarray}$$

Monte Carlo Markov chain (MCMC) methods can be applied to draw samples from the posterior distribution. However, the MCMC technique is known to be computationally expensive and not scalable to large datasets. For large-scale modern astronomical surveys, we can adopt the faster maximum-a-posterior (MAP) estimation. In this way, the latent variables z and the class label C are found by

$$\begin{eqnarray}\begin{array}{ll}\hat{{\bf{z}}},\hat{C}= & \displaystyle \mathop{{\rm{\arg }}\,{\rm{\max }}}\limits_{{\bf{z}},C}\{\displaystyle \sum _{j=1}^{n}\mathrm{log}p({\rm{y}}_{j}|{\mathcal{G}}({\bf{z}}))+\mathrm{log}p({\bf{z}}|C)\\ & +\mathrm{log}p(C)\}.\end{array}\end{eqnarray} \tag{ 7 }$$

With the computed MAP estimation $\hat{{\bf{z}}}$, the cleaned and denoised spectra are given by $\hat{s}={\mathcal{G}}(\hat{{\bf{z}}})$.

This subsection reviews some works on deep density estimation. These works serve as the basis for building our deep Bayesian denoising model. Suppose x ∈ ^D represent a D-dimensional observation, for which we want to estimate its distribution. In the deep learning literature, most density estimation methods are based on the intuition that we can transform a simple base density π_z(z) into a more complex one p(x) via an invertible differentiable transform x = f(z). The function f is parameterized by a neural network to achieve flexible and adaptive transformation. By the basic formula of density transformation, it holds that

$$\begin{eqnarray}&&p({\bf{x}})={\pi }_{{\bf{z}}}({f}^{-1}({\bf{x}}))\left|\det \left(\displaystyle \frac{\partial {f}^{-1}}{\partial {\bf{x}}}\right)\right|.\end{eqnarray} \tag{ 8 }$$

Several methods have been proposed for deep density estimation. The approach of normalizing flows (Dinh et al. 2014) chooses f as a sequence of composite functions, i.e., f = f₁° f₂ ° ⋅ ° f_K . In this way, p(x) can be regarded as an invertible and differentiable transformation f of a base density π_z(z). The base density can be simple multivariate Gaussian distribution. The relationship between the observation data x and latent variable z can be represented as $x\mathop{\to }\limits^{{f}_{1}}{{\bf{h}}}_{1}\mathop{\to }\limits^{{f}_{2}}{{\bf{h}}}_{2}\cdots \mathop{\to }\limits^{{f}_{K}}{\bf{z}}$ via a stack of hidden variables. Under the invertible (or bijective) assumption of f, z can be calculated as z = f⁻¹(x) given x. Based on Equation (8), the log probability density can be directly computed as

$$\begin{eqnarray}\begin{array}{ll}\mathrm{log}p(x) & =\mathrm{log}{\pi }_{{\bf{z}}}({\bf{z}})+\mathrm{log}\left|\det \left(\displaystyle \frac{d{\bf{z}}}{dx}\right)\right|\\ & =\mathrm{log}{\pi }_{{\bf{z}}}({\bf{z}})+\displaystyle \sum _{i=1}^{K}\mathrm{log}\left|\det \left(\displaystyle \frac{d{{\bf{h}}}_{i}}{d{{\bf{h}}}_{i-1}}\right)\right|,\end{array}\end{eqnarray} \tag{ 9 }$$

where the scalar value log|det(dh_i /dh_{i – 1})| is the logarithm of the absolute value of the determinant of the Jacobian matrix (dh_i /dh_{i – 1}), also called the log-determinant. The series of transformations are required to be easily invertible and the log-determinant should be easy to compute. Rezende & Mohamed (2016) devised planar and radial flow as basic blocks for f. NICE (Dinh et al. 2014) adapts additive coupling layers to form a normalizing flow, and its successor Real NVP (Dinh et al. 2017) extends the transformation by stacking affine coupling layers. They all have a tractable triangular Jacobian matrix for the bijective mapping.

Autoregressive flow is another popular and tractable approach to density estimation. It factorizes the joint density as a product of conditional densities p(x) = Π_dp (x_d |x_{1:d – 1}) via the chain rule of probability (Uria et al. 2016). This factorization makes the Jacobian tractable as the Jacobian becomes a triangular matrix. The inverse autoregressive flow (IAF, Kingma et al. 2017) and the masked autoregressive flow (MAF, Papamakarios et al. 2018) take this approach. Under this approach, the prediction of current value depends on all of its past values, which is referred to as the autoregressive property. They use independent standard Gaussian distributions as the base density. The mean and variance of x_d are functions of the preceding observation vector x_{1:d – 1} or the preceding random numbers. They both use MADE (Germain et al. 2015) as their basic building blocks for the function mapping.

We first detail our spectrum generator ${\mathcal{G}}$ and the prior density p(z|C) for the latent variable. The generator is obtained from the auto-encoder framework, but with an additional local isometry constraint. In the standard auto-encoder framework, an encoder ${\mathcal E}$ maps an observation y to the latent ${\bf{z}}= {\mathcal E} ({\rm{y}})$, and then maps z ∈ ^L back to the original high-dimensional space ^D by the generator (decoder) ${\mathcal{G}}$. The training target is to minimize the reconstruction loss

$$\begin{eqnarray}&&\displaystyle \mathop{\min }\limits_{{\mathcal{G}}, {\mathcal E} }{{\rm{E}}}_{{\rm{y}}}\Vert {\rm{y}}-{\mathcal{G}}({\mathcal E} ({\rm{y}})){\Vert }^{2}.\end{eqnarray}$$

The standard convolutional auto-encoder architecture can be specified for ${\mathcal E}$ and ${\mathcal{G}}$. However, the generator ${\mathcal{G}}$ obtained hereby could create distortion in the latent space ^L , affecting density estimation. The left subfigure of Figure 1 shows an example of distorted mapping, where some low-density points (on the bottom right) are mapped to a high-density region (on the top right). To avoid the issue, we construct a locally isometric mapping such that ${\mathcal{G}}$ preserves the distance between samples in the latent space ^L . In other words, it holds that

$$\begin{eqnarray}&&\Vert {\mathcal{G}}({\bf{z}})-{\mathcal{G}}({{\bf{z}}}^{{\boldsymbol{^{\prime} }}})\Vert \approx \Vert {\bf{z}}-{{\bf{z}}}^{{\boldsymbol{^{\prime} }}}\Vert,\end{eqnarray}$$

where z, z' ∈ ^L is a pair of latent variables satisfying ||z – z'|| ≤ δ for some δ. The local isometry property is illustrated in the right subfigure of Figure 1, where the distance between the red and blue points is preserved under the mapping ${\mathcal{G}}$. The isometry property allows us to obtain the spectrum density p(s) by directly accessing the low dimensional z. It approximately holds that $p(s)=p({\mathcal{G}}({\bf{z}}))=p({\bf{z}})$. This helps us to approach the density estimation and effectively avoid the computation of $\det \displaystyle \frac{\partial {f}^{-1}}{\partial {\rm{x}}}$ for a complex transformation function f in Equation (8).

Zoom In Zoom Out Reset image size

With the additional local isometry constraint, the objective function to train the generator becomes

$$\begin{eqnarray} &&\displaystyle \mathop{\min }\limits_{{\mathcal{G}}, {\mathcal E} }{{\rm{E}}}_{{\rm{y}},{\boldsymbol{\delta }}}\{{[\Vert {\mathcal{G}}({\bf{z}})-{\mathcal{G}}({\bf{z}}+{\boldsymbol{\delta }})\Vert -\delta]}^{2}+\Vert {\rm{y}}-{\mathcal{G}}({\mathcal E} ({\rm{y}})){\Vert }^{2}\}.\end{eqnarray} \tag{ 10 }$$

In the above, the latent value is $z= {\mathcal E} (\rm{y})$. The perturbation δ is a random variable, drawn from a uniform distribution over the sphere of radius δ in ^L . Similar loss functions had been employed in the works (Geng et al. 2020; Atzmon et al. 2020) to train the auto-encoder to learn the manifold structure.

After training ${\mathcal E},{\mathcal{G}}$ based on Equation (10), we can compute the latent variables z for all training samples. Then, based on this collection of latent variables, a kernel density estimator (KDE) is deployed over z for the training samples of each stellar subclass C. This helps us to obtain p(z|C) for each subclass C.

This subsection constructs the observation model p(y|s) for y_n given the true signal $s={\mathcal{G}}({\bf{z}})$. Recall that we have used the additive noise model _n = y_n – s. Our goal is equivalent to construct a noise density model p(_n ) and set p(y_n |s) = p(_n ). The density model is built upon the idea of Papamakarios et al. (2018) such that the observation model has an auto-regressive structure. Compared with their model, our model only depends on a local group of pixels within an observation, is more parsimonious in parameterization, and focuses on extracting the local correlation structure of the noise.

Suppose the noise vector is written as _n = (_n1,_n2,..., _nK ). For d = K + 1,..., D, the noise value _nd for the d-th pixel depends on its immediate K preceding pixels via

$$\begin{eqnarray}&&p({\epsilon }_{nd}|{\epsilon }_{n,(d-K):(d-1)})={\mathcal{N}}({\epsilon }_{nd}|{\mu }_{nd},{(\exp {\alpha }_{nd})}^{2}),\end{eqnarray} \tag{ 11 }$$

where μ_nd = f_μ (_{n,(d – K):(d – 1)}), α_nd = f_α (_{n,(d – K):(d – 1)}), and f_μ and f_α are two functions expressed by neural networks. The architecture of f_μ and f_α will be specified later. The two functions f_μ and f_α determine the mean and the standard deviation for the noise _nd at the d-th pixel. In particular, we have

$$\begin{eqnarray}&&{\epsilon }_{nd}={\xi }_{nd}\exp ({\alpha }_{nd})+{\mu }_{nd},\end{eqnarray} \tag{ 12 }$$

where ${\xi }_{nd}\sim {\mathcal{N}}(0,1)$ follows the standard Gaussian distribution. Equivalently, the random variable can be expressed by the following inverse expression

$$\begin{eqnarray}&&{\xi }_{nd}=({\epsilon }_{nd}-{\mu }_{nd})\exp (-{\alpha }_{nd}).\end{eqnarray} \tag{ 13 }$$

The above specifies the local dependence structure for d = K + 1,⋅, D. For the first K pixels, there are not enough preceding pixels for us to determine the conditional likelihood (11). Instead, for the first K pixels, the noise _nd is imposed to follow a univariate Gaussian distribution with a fixed mean μ_d and a fixed standard deviation σ_d = exp(–α_d ). In summary, our model parameters to be trained include: the scalar values μ_d , α_d for d = 1,⋅, K; and two neural network functions f_μ and f_α .

As μ_nd and α_nd depend on ε_{n,(d – K):(d – 1)} by design, the Jacobian in density transformation Equation (8) is an upper triangular matrix. As a result, the absolute value of the determinant can be easily computed as

$$\begin{eqnarray}&&\mathrm{log}\left|\det \displaystyle \frac{\partial {f}^{-1}}{\partial \epsilon }\right|=-\displaystyle \sum _{d=K+1}^{D}{\alpha }_{nd}.\end{eqnarray} \tag{ 14 }$$

It follows that the log density of the observation model becomes

$$\begin{eqnarray}\begin{array}{ll}\mathrm{log} & p({\rm{y}}_{n}|{\rm{s}})=p({\epsilon }_{n})\\ = & \displaystyle \sum _{d=1}^{K}\mathrm{log}p({\epsilon }_{nd})+\displaystyle \sum _{d=K+1}^{D}\mathrm{log}p({\epsilon }_{nd}|{\epsilon }_{n,(d-K):(d-1)})\\ = & \displaystyle \sum _{d=1}^{K}[-(1/2)\exp (-2{\alpha }_{d}){({\epsilon }_{nd}-{\mu }_{d})}^{2}-{\alpha }_{d}]\\ & \,+\displaystyle \sum _{d=K+1}^{D}[-(1/2)\exp (-2{\alpha }_{d}){({\epsilon }_{nd}-{\mu }_{nd})}^{2}\\ & \,-{\alpha }_{nd}]+{\rm{const}}.\end{array}\end{eqnarray} \tag{ 15 }$$

The last equality holds up to some irrelevant constant const.

For our observation model, we develop a neural network architecture for local noise feature extraction, as shown in Figure 2. In the first layer, local features are extracted by a one-dimensional convolution layer. This convolution layer has one input channel and C output channels, kernel size K, one stride and zero padding. Its output h is of dimension N × C × (D – K + 1), where N is the mini-batch sample size. In accordance of the aggressive structure, we drop one redundant column (the first column) in the third dimension of the hidden state h¹, such that its dimension becomes N × C × (D – K). We then take a transpose to exchange the second dimension with the third dimension. The resulting h is of dimension N × (D – K) × C. In this way, for the d-th pixel (d = K + 1,..., D), we get a C-dimensional feature vector for it.

Zoom In Zoom Out Reset image size

After that, h is used as the input of the following L – 1 linear hidden layers with RELU, to sequentially reduce the dimension of hidden variables from N × (D – K) × C to N × (D – K) × C' for some C' < C. In the sequel, there are two separate linear hidden layers mapping from N × (D – K) × C' to N × (D – K) × 1, one is for μ and the other is for α, and we transpose μ and α back to N × 1 × (D – K). At this moment, vector μ and α contain the mean and standard deviation information for _nd with d = K + 1,⋅, D. As for the first K pixels, their scalar mean and standard deviation values (μ_d , α_d for d = 1,..., K) get included via the final masked linear layers.

Moreover, the observation model developed hereby can naturally handle a spectrum with missing flux values due to bad pixels. We can create a mask vector m_n such that m_nd = 1 if all of y_{n,d – K} , ..., y_{n,d – 1},y_n,d are observed, and m_nd = 0 otherwise. The observation log-likelihood for y_n becomes

$$\begin{eqnarray}\begin{array}{ll}\mathrm{log}\,p({\rm{y}}_{n}|{\rm{s}})= & \displaystyle \sum _{d=1}^{K}{m}_{nd}\mathrm{log}\,p({\epsilon }_{nd})\\ & +\displaystyle \sum _{d=K+1}^{D}{m}_{nd}\mathrm{log}\,p({\epsilon }_{nd}|{\epsilon }_{n,(d-K):(d-1)}).\end{array}\end{eqnarray} \tag{ 16 }$$

In other words, the likelihood of the d-th pixel is taken into account if and only if the d-th pixel and its immediate K preceding pixels are observed, which is shown in Figure 3. The masked likelihood allows us to deal with partially observed spectrum without resorting to ad-hoc missing value imputation.

Zoom In Zoom Out Reset image size

The main workflow of our model is summarized in Figure 4. In order to train and apply our proposed model, basically four main steps are involved:

• Train an encoder ${\mathcal E}$ and a locally isometric generator ${\mathcal{G}}$ based on a collection of high-SNR optical stellar spectra.
• Use Gaussian kernel density estimation (Gaussian KDE) to estimate the prior distribution of the latent variable z for each stellar subclass obtained from the encoder trained in the first step.
• Train a NoiseFlow observation model with the observational noise extracted from low-SNR optical stellar spectra.
• For each test spectrum, use the Gaussian KDE, the NoiseFlow model and the generator ${\mathcal{G}}$ trained in the first three steps to obtain its corresponding latent variable and its clean and complete spectrum with iterative optimization of Equation (7).

Zoom In Zoom Out Reset image size

To test the model performance, we randomly select an independent group of spectra with high SNR for each stellar subclass, as described in Section 2. These selected spectra are regarded as the ground truth that the model endeavors to predict. Then we extract another noisy dataset from an independent group of stellar spectra with SNR ranging from 10 to 20. The noise is randomly sampled and added to the above clean test spectra. These constitute our benchmark test dataset for method comparison. Each stellar subclass is created with 300 test spectra.

Figure 7 and Figure 8 exhibit some examples of denoised spectra for ten stellar subclasses. We choose one representative result from each of the ten stellar subclasses to demonstrate the power of our method. The left column shows the spectrum before and after noise contamination, where the purple curve is the high-SNR spectrum and the grey curve is the clean spectrum with added noise. The noisy spectrum gets cleaned by the standard auto-encoder and our proposed method. The denoised spectrum is shown in the right column of Figure 7 and Figure 8. In each subfigure of the right column, the purple spectrum is the true spectrum that both denoising algorithms try to recover. Though the clean purple spectrum is unknown to the denoising algorithm, our proposed method shows promising ability to recover it from the noisy observation. The predicted spectrum from our model (yellow curve) is much closer to the purple clean spectrum than the standard auto-encoder result (blue curve). Our proposed method also has the capacity to remove strong noisy emission lines (see the second row of Fig. 8) and keeps the signal emission lines (see the third row of Fig. 8).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows the overall comparison between our method and the convolutional denoising auto-encoder in spectrum denoising for various stellar subclasses. Each yellow point in a subfigure represents one synthetic noisy stellar spectrum. The noisy spectrum gets cleaned and compared with the true high-SNR spectrum, and the reconstruction loss is computed. The reconstruction loss is computed as follows. Suppose that s is the true signal spectrum and $\hat{\rm{s}}$ is the denoised spectrum by a model. Then, the reconstruction loss is $\sum _{d=1}^{D}{({s}_{d}-{\hat{s}}_{d})}^{2}/D$. In each subfigure, the horizontal axis is the reconstruction loss for the convolutional denoising auto-encoder, and the vertical axis is the reconstruction loss of our method. The blue line indicates where the two methods have equal performance. We can see that most points in each subfigure fall below the blue line, indicating that our method has smaller reconstruction loss. The title of each subfigure reports the proportion of spectra for which our method has smaller reconstruction loss. Within each stellar subclass, our method produces higher-quality denoised spectra for more than 80% of the testing samples. For subclasses such as carbon class and WD class, our method demonstrates improved performance for almost 100% of the test samples.

Zoom In Zoom Out Reset image size

We further consider how the signal-to-noise ratio affects our model performance. For the synthetic spectra, we measure their signal-to-noise ratio by the formula given by Stoehr et al. (2008). The computed signal-to-noise ratio is denoted as DER_SNR. To decrease DER_SNR of the synthetic data to a specific level, we add additional gaussian noise with various noise levels to each spectrum. This procedure results in a new group of test dataset. The left panel of Figure 10 shows how the boxplot of the reconstruction loss varies across distinct DER_SNR levels. The reconstruction loss does not increase very quickly as DER_SNR decreases. When the DER_SNR is below 3, the median reconstruction loss is still about 10⁻³. To illustrate how the spectrum looks like at this level of reconstruction loss, we plot a few examples in Figures 11–12. The grey curve in the left panel of Figure 11 is one synthetic spectrum with DER_SNR equaling to 2.79. The reconstruction loss between the true spectrum and the denoised spectrum in the right panel is about 1.0 × 10⁻³. Figure 12 shows one more example where the DER_SNR of the synthetic spectrum is 2.12 and the reconstruction loss is about 1.3 × 10⁻³. The right panel of Figure 10 compares the DER_SNR before and after denoising for the above dataset. For most spectra, their DER_SNR after denoising is above 100. From the right panel of Figure 10, we can also find that the DER_SNR of our model output is stable regardless of the DER_SNR of the input spectrum. This is because our model prior component and the generator in Equations (4)–(6) are fixed after model training. The denoised spectra, which are generated by the prior and the generator, will always have the same level of SNR of the training dataset.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Besides the above test dataset, another test dataset is created to evaluate model performance for spectrum denoising with missing flux values. To construct this benchmark data, almost the same procedure of Section 4.1 is taken. Realistic noise is added to the clean spectrum with high SNR. In addition, we randomly remove the flux values over an interval range of wavelength. The interval of missing pixels is also randomly selected for each test spectrum.

Based on the idea of masked likelihood Equation (16), our trained model can be directly employed to denoise spectrum with partially missing flux values. In other words, our model does not require re-training to deal with this kind of data. However, the standard denoising convolutional auto-encoder requires re-training to adapt to this dataset. Its original training data also get randomly removed flux values over a random wavelength range. The missing values are imputed by zeros, and the denoising convolutional auto-encoder aims to reconstruct the original full spectrum from the zero-imputed spectrum. The full spectrum is employed for its loss computation and parameter update.

Figure 6 shows the performance of our model depending on different positions and widths of the missing range. The three panels correspond to the cases where the missing flux occurs at the left end (blue end), the middle and the right end (red end) of the spectrum, respectively. In each panel, the horizontal axis is the number of missing pixels out of the total D = 2048 pixels. The vertical axis is the reconstruction loss. As expected, the reconstruction loss increases with the growing number of missing pixels. When the length of missing pixels is 800 (i.e., 40% of the whole spectrum), the median reconstruction loss is still below 10⁻³. Generally, the model performance is not sensitive to the location of the missing pixels, but missing at the left end (blue end) incurs slightly higher loss. This is due to the fact that the blue end is feature-rich and contains more information for spectrum reconstruction for most spectra.

Figure 13 and Figure 14 plot a few examples for the ten stellar subclasses. The left column shows the spectrum before and after noise addition and flux value removal. The purple spectrum is the original spectrum with high SNR. The grey spectrum has noise added, but at the same time, a random interval of flux values is removed. The missing flux is plotted as an interval of zeros. The purple spectrum is the ground truth spectrum that both algorithms try to recover. The denoised and imputed spectrum is shown in the right column. Our predicted spectrum is shown in yellow, and the result of convolutional auto-encoder is shown in blue. Our resulting spectra are much closer to the true spectra in these cases. The overall result for this test data is summarized in Figure 15. The interpretation of the figure is similar to that of Figure 9. Although our model has a moderate lead over the convolutional auto-encoder in F-type and G-type classes, the performance difference margins are wider in the other stellar subclasses.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size