Please trust the notebook to view the visualizations.

CS109B: Advanced Topics in Data Science

Measuring the Shape and Brightness of Galaxies with Neural Networks¶

Group 75: Dmitry Vukolov, Ning Xu, Rohit Beri, Sunil Chomal¶

Harvard University, Spring 2020
Instructors: Pavlos Protopapas, Mark Glickman and Chris Tanner
Project Advisors: Douglas Finkbeiner and Jun Yin

Outline¶

1. Problem Statement
2. Data Generation
3. Exploratory Data Analysis
4. Model Building
5. Evaluation
6. Future Work
7. References

Problem Statement¶

Vast amounts of data: As astronomers collect more and more image data, there is a need for further development of automated, reliable, and fast analysis methods.

Traditional methods: For example, to describe a galaxy, software is used to fit a parametric or a non-parametric model and estimate the galaxy's brightness, shape, size, and orientation. A commonly used, though not a perfect model, is the Sérsic profile, developed in 1963.

New methods: In recent years, researchers have started exploring the application of deep learning methods to analyze astronomical data. For instance, classification of galaxies using CNNs has demonstrated high levels of accuracy. In this work, we study the possibility to use neural networks to estimate galaxy parameters from simulated data.

Our objective: Explore deep neural network architectures to estimate galaxy parameters from simulated data.

Data Generation¶

We use GalSim, an open-source library, to generate the images of galaxies, following these steps:

1. Sérsic profile: Define a galaxy's Sérsic profile parametrized by the Sérsic index and the radius that encloses half of the total flux
2. Flux & shear: Add flux and shear (defined by ellipticity and orientation) to complete the galaxy definition
3. PSF: Convolve the galaxy profile with the Point Spread Function, which is determined by the telescope optics and the atmosphere (for ground-based telescopes).
4. Noise: Add Poisson noise (e.g. detected photoelectrons) and Gaussian noise (e.g. read noise) to the generated the image
5. Signal-to-noise ratio: Compute the signal-to-noise ratio (SNR) based on the pixel values and the noise level (assuming faint galaxies):

We preserve only the images with SNR between 10 and 100.

Datasets¶

We generate three datasets with increasing complexity for parameter estimation:

Number of samples: We simulate 200,000 observations in total, using 180,000 samples for training and 20,000 for validation. Our choice to generate a large number of samples is driven by literature review, e.g. Tuccillo, 2017 and Madireddy, 2019, who are using 55,000 and 120,000 images respectively and suggest scaling up to 1 million. Our subsequent observations confirm that even on synthetic data, model performance improves as we add more data.

Discretization: One major challenge in simulating a large number of observations has to do with GalSim, which is relatively slow especially when generating images with higher Sérsic indices. The developers of GalSim recommend considering the use of only discrete values of the Sérsic index for large simulations. This takes advantage of internal caching mechanisms and speeds up generation considerably. We follow that advice, discretize Sérsic index into 1,000 equally spaced values, and increase GalSim's cache size. This cuts down CPU time for data generation from approximately 90 hours to 1.2 hours per dataset.

Target variables¶

Our goal is to estimate the following five parameters from image data, which are recorded alongside the generated images:

• Sérsic profile: determined by the Sérsic index and the Sérsic radius
• Galaxy flux
• Reparameterized ellipticity and orientation: $g_1$ and $g_2$

Source code: The code for data generation is based on the template provided by the module leaders and is available in the GitHub repo.

Exploratory Data Analysis¶

Milestone 2 EDA: The majority of EDA has been carried out for Project Milestone 2 and is available in the corresponding EDA.ipynb notebook. We avoid repeating the analysis in the current report and instead refer the reader to the original document. Here we present only the findings that informed our model building approach.

Interactive app: As part of the analysis we created an interactive web app that laid the foundation for understanding the data. We used it throughout our project to study the relationship between galaxy parameters and the resulting images. It also helped us identify which parameter combinations could present a particular issue for the model versus which would be relatively easier to estimate.

Key Findings¶

The noisy images of the galaxies are characterized by three different aspects:

• The galaxy itself, independent of an observer. Sérsic index, Sérsic radius, flux, and shear are some of its parameters.
• The observational distortion determined by the telescope specification and the state of the atmosphere. The Point spread function (PSF) is used to encode this information.
• Contributions of noise, which can include Poisson statistics of detected photoelectrons and Gaussian noise such as read noise.

These are independent of each other, however, the final image is affected by each of these. We need to train a network to "deconvolve" these aspects and predict the parameters.

Considerations for Model Building¶

The knowledge of the data-generating procedure, together with the EDA informed our subsequent model building as follows:

• The Sérsic index is quite difficult to estimate visually making it an important parameter when doing model selection.
• The target variables need to be rescaled so that the loss is not driven by the variables with the largest magnitude.
• The simulated images form a perfect dataset, free of outliers and irregularities, with underlying deterministic relationships between a handful of parameters. Smaller neural networks might be capable of capturing the underlying dependencies in the data, and regularization might turn out to be optional.
• The PSF and noise levels have a significant effect on the obtained images. Where possible, they should be used as inputs to the model.
• To understand how dramatic the effect of the noise is, we might want to establish an upper bound on model performance by training on the noiseless images first.
• Given the availability of noiseless images, models that take advantage of this additional information should be explored.

Model Building¶

In our model building we used multiple approaches:

Autoencoders¶

An autoencoder has two potential advantages:

• First, an autoencoder is learning from both the labels and the noiseless images and thereby might incorporate more information about the underlying relationships into the network.
• Secondly, an autoencoder can be used not only to predict the labels but also to denoise and reconstruct the images.

Furthermore, we take inspiration from Jun Yin's work with Conditional Autoencoders mentioned during the lectures.

• We experiment with several possibly related architectures.
• We begin by trying to replicate GalSim output through building a decoder first — clean images are a deterministic function of the labels.
• Then we invert the decoder architecture to build a basic encoder structure.

Neural architecture search¶

Using a small subset of the data we run AutoKeras, an AutoML tool, to quickly test vanilla CNNs, ResNets, and Xception networks with different complexities, regularization, and normalization parameters. The search is guided by Bayesian optimization with Gaussian Processes.

Grid search of hyperparameters¶

• We pick several key hyperparameters of the best model, expand their range, and evaluate the effect using a small portion of the data.
• This allows us to significantly reduce the model size, decrease the training time, while simultaneously maintaining or even improving predictive performance.

Denoising pipeline¶

Informed by a large gap between performance metrics for noiseless and noisy data, we test a two-stage pipeline described in Madireddy (2019) that uses a separate denoising network as the first step.

Autoencoder¶

The initial idea of using an autoencoder for the task occurred to us during the EDA stage. We recognized that the dataset available to us has both the parameters and the noiseless image. An autoencoder can make use of both to train a network and predict them at the same time. An autoencoder can thus solve two problems simultaneously and may provide more robustness to the model. We were further inspired by Jun Yin's comments about using a Conditional Autoencoder with parameters representing the latent space.

Decoder¶

We know that the simulated noiseless images generated by GalSim are a deterministic function of five parameters plus the PSF. As the first step, we looked for an architecture of the Decoder which can take 5 parameters as input and output an image of the galaxy, i.e. learn the GalSim deterministic function via a deep CNN architecture. Our reasoning behind this being that if we can figure out a network to generate correct images from the labels, inverting the same architecture may let us do the reverse, i.e. generate the parameters from a noiseless image. If we have such an encoder, we can then look at suitably modifying the same to generate the parameters from the noisy image as well. During this exercise we explored the following aspects:

• Normalizing the parameters (inputs)

  • Since the parameters are on different scales, particularly the Flux, we normalize all to $[0,1]$.

• Normalizing the noiseless images (output)

  • After experimenting with individual image level normalization and training set level normalization, we chose the latter, as it produces better results.

• We modify the ConvDecoder layer used in our DCGAN assignment

  • Increase the number of Conv2D layers to six.
  • Change the activation function to LeakyReLU, except for the last layer where we use sigmoid because of the $[0,1]$ normalization of the images.
  • We also investigate different striding and in the end choose a stride of $(2,2)$ .
  • Using Dropout with Conv2D reduces the effectiveness of the ConvDecoder and hence we set the dropout rate to $0$.

• We also experiment with the number of Dense layers and ultimately use two Dense layers

  • We use linear activation for the first layer and LeakyRelu for the second layer.
  • The reason for the first linear activation is to maximize the information captured from the input layer.
  • We use Dropout layers with a dropout rate of $0.50$ after both Dense layers.

• We use BinaryCrossentropy as our loss function

  • The last layer of our network uses sigmoid activation function to align with the $[0,1]$ normalization of the images.
  • We empirically found binary cross-entropy to work better compared to MSE on our dataset. Our hypothesis is that might be either due to the distribution of pixel values that are close to either 0 or 1, or due to a less pronounced vanishing gradient problem.

Our final decoder architecture¶

Encoder¶

As the first step towards building the encoder, we reverse the decoder architecture above as follows:

• Normalize the noiseless images and parameters similar to the decoder above
• The ConvDecoder layer is replaced by ConvEncoder from our DCGAN assignment.
  • Conv2D layers are replaced with Conv2DTranspose layers with a stride of $(2,2)$
  • All of the layers use LeakyReLU activation
  • The dropout rate is set to $0$
• The output Dense layer uses sigmoid activation and the Dense layer immediately preceding it uses linear activation
  • We use Dropout layers with a dropout rate of $0.50$ after the intermediate Dense layers.
• Again use BinaryCrossentropy as our loss function together with sigmoid activation of the final layer

The above architecture produces good results on noiseless images, so we proceed to train the same on the noisy dataset. After multiple experiments with different datasets, additional inputs, batch level normalization, additional image preprocessing (as recommended by Jun Yin during the office hours), we make the following modifications to our architecture:

• We do away with the normalization of the input images and resort to a BatchNormalization layer in the encoder
• We learn that Sigma & PSF can be independently computed and hence are available as inputs
• Image input is now a 7-channel input layer with these channels:
  1. Original Noisy Image: $I$
  2. Noisy Image divided by Sigma ($\sigma$), the level of Gaussian noise of the image: $\frac{I}{\sigma}$
  3. Squared Noisy Image divided by $1.41$ times Sigma squared: $\frac{I^{2}}{1.41\sigma^{2}}$
  4. Noisy Image divided by the Variance of pixel values of the image: $\frac{I}{\sigma_p^{2}}$
  5. Noisy Image divided by the product of the Squared Variance and Sigma: $\frac{I}{\sigma_p^{4}.\sigma}$
  6. PSF Image: $P$
  7. PSF Squared Image: $P^{2}$
• Apart from the image input layer, we add a secondary input layer as follows:
  1. Sigma ($\sigma$), the Gaussian noise level of the image
  2. The PSF value for the full-width-half-max of the Moffat profile ($\text{psf}_r$)
  3. Standard Deviation ($\sigma_p$) of pixel values of the image
  4. Mean ($\mu_p$) of pixel values of the image

Our final encoder architecture¶

Model Parameters¶

Autoencoder: Training¶

Training Process¶

Our training process is divided into $3 \times 3$ parts: $3$ datasets times $3$ steps. For each dataset, we follow these steps:

1. Train the Encoder
2. Train the Decoder
3. Train the Autoencodder to fine-tune the weights

We find that $3$-step training achieves better results than training the Autoencoder in one go.

Training Hyper-parameters¶

Training Time¶

Validation Results¶

Dataset 1: Fixed Gaussian Noise Level and Fixed PSF¶

Dataset 2: Varying Gaussian Noise Level and Fixed PSF¶

Dataset 3: Varying Gaussian Noise Level and Varying PSF¶

Autoencoder: Source Code¶

The notebook detailing the training process and the code can be found in a supplemental notebook named autoencoder_training.ipynb, submitted as part of this project.

Neural Architecture Search¶

We perform the initial neural architecture search with the help of AutoKeras, an AutoML tool from the Keras team. On a small subset of the data:

• We test different type of architectures: regular CNNs, ResNets, and Xception models
• Vary the specifics of the architectures: depth and the width of the network, the convolutional kernel size, striding, type of pooling, different activation functions
• We test whether data normalization and dropout regularization are necessary
• As well as briefly check the effect of three optimizers: Adam, Adagrad, and SGD

Our hypothesis is that due to the simulated and homogeneous nature of the dataset, the performance on a small subset will be indicative of the performance on the full dataset.

The output of the search demonstrates that:

• ResNets perform significantly worse on our data compared to regular CNNs and Xception models
• Normalization and specifically standardization of the images to zero mean and standard deviation of one is beneficial
• Dropout is generally unnecessary, as it negatively impacts predictive performance and increases the training time
• Adam converges faster than the alternative optimizers and produces the lowest mean-squared error

The resulting best model chosen by AutoKeras has Xception architecture with 18 depthwise separable convolutional layers, 8 residual connections, 16 convolutional filters in the entry flow, 384 filters in the separable convolutions, a kernel size of 5, selu activation and global average pooling.

We implement this architecture in code and instantiate the model:

Grid Search of Hyperparameters¶

Like any tool, AutoKeras is not perfect. We use it as a starting point to select one possible architecture that works reasonably well on our data. The resulting model is relatively large, containing around 3.5 mln parameters and trains slowly. In an effort to optimize the model further, we identify two key hyperparameters that determine model complexity the most and run a grid search on them. We choose to keep the network quite deep, as selected by the neural architecture search. This setting does not influence the total number of parameters of Xception architecture significantly while permitting it to capture highly non-linear relationships.

We then interactively analyze the results of the grid search by slicing the data along each axis. Please trust the notebook to display the visualization tool. The help section of the tool lists the various ways of using it:

The performance of many Xception models is quite similar. Still, based on the visualization of the grid search, sep_num_filters=64 produce low values of the loss with a very small variance. We combine it with a low number of filters in the initial convolutional layers conv2d_num_filters=16 to obtain our best model. This ends up significantly reducing model complexity from 3.5 mln parameters to just over 109K, while at the same time improving performance and substantially cutting down the training time:

Below we show an abbreviated version of the model's architecture with just two residual blocks for the sake of readability. The central section of two separable convolution blocks and the following addition is repeated six more times in the final model:

Model Training: Xception CNN¶

We proceed to train the model on the full dataset with 200,000 observations. We use 90% of the data for training and 10% for validation. We train with early stopping monitoring validation loss, which is MSE for the five target labels. To establish an upper bound on model performance, we first train the model on noiseless images. We then progress to noisy images from the three datasets of increasing complexity. For the final dataset with varying PSF, we additionally supply the model with PSF images as input.

Training on a dataset of 200,000 images with early stopping takes around three hours on JupyterHub, i.e. on Nvidia Tesla K80 GPU. After training the models, we print out the resulting RMSE for each of the target labels in their original units:

Model performance in terms of RMSE drops as we add noise and increase data complexity. By far the largest decrease in predictive ability happens when we transition from noiseless images to noisy ones. RMSE is around 20 times worse for noisy images compared to noiseless ones, even when the level of Gaussian noise is fixed. Informed by this finding, as well as by Madireddy (2019) we set to explore whether state-of-the-art denoising might be helpful to reduce the gap between the quality of parameter estimation on noisy and noiseless data.

Denoising Pipeline¶

In their paper Modular Deep Learning Analysis of Galaxy-Scale Strong Lensing Images Madireddy et al. explore multi-stage denoising and deblurring pipeline that precedes the classification and regression networks. In their experiments, a modern image reconstruction neural network applied to noisy images significantly improved model performance on the final task.

We follow the same approach by implementing EDSR (Enhanced Deep Residual Networks) and RDN (Residual Dense Networks), two state-of-the-art algorithms for image denoising and restoration. We train these models on a separate dataset of simulated galaxies, then apply them to our noisy inputs and proceed with an Xception CNN model. At full capacity, both denoising networks do reduce the RMSE on the validation set, although not by much. The two-stage pipeline approach raises the overall complexity of the method, as well as drastically increases the total number of parameters. A possible direction of future research would be to try to combine the denoising network and the CNN into a single model.

The implementations of denoising algorithms are available in our GitHub repo for future reference.

Evaluation¶

We will first analyze our models visually and then proceed to quantitative evaluation. One way we can visualize the predictions of the Xception model is by generating an image of the galaxy with GalSim using the estimated parameters:

Even though the noiseless and the reconstructed images look very similar, there are in fact significant differences between the two, which matter. For instance, our model is substantially underestimating pixel intensity in the very center of the galaxy. In other words, the model is underestimating the Sérsic index. The remaining parameters are also somewhat off, but to a lesser degree:

Our autoencoder model produces its own reconstructed images, which we can analyze in a similar fashion:

From a performance perspective, there are cases when the model does a reasonably good job of reconstructing the images. There are also situations where the differences are very apparent. Overall, we find that most of the error is concentrated in the center of the galaxy, as seen in the plots of residuals. Furthermore, we discover that the distortion is largely due to errors in parameter estimation, i.e. the encoder, as demonstrated by the plots on the log-scale.

We can also visualize the correspondence between the true and the estimated parameters across the whole validation set:

Our model seems to have no difficulties estimating galaxy flux from noisy images. The deviation from a 45-degree line is minimal in the first plot. In addition to that, images with higher flux also have a higher signal-to-noise ratio. Still, the model is estimating almost the whole range of flux values equally well.

The correspondence between the true and the predicted values for the other parameters is much weaker. Especially problematic is the Sérsic index, which demonstrates the highest level of dispersion around the diagonal line. For all of these parameters there is a clear relationship between SNR and the precision of the estimates: the higher the signal-to-noise ratio, the better our model can estimate the parameter and vice versa.

To better understand the patterns in the errors, we will plot those against the true values:

In the case of flux, there seems to be a slight tendency to underestimate the value of the parameter at the higher end of the spectrum. For the other parameters, we see a common trend: the parameters are overestimated at the lower values and underestimated at the higher values. This relationship is mainly visible for the images with low signal-to-noise ratios. It is much less pronounced for images with high SNR.

Overall, the distributions of errors vaguely resemble bell-like and mostly symmetric shapes, which are centered around zero. That, however, does not exclude the presence of systematic errors in particular regions of each parameter, as we've seen previously:

Benchmark against the Cramér–Rao bound (CRB)¶

So far our model evaluation was performed on the validation set. To get an additional insight into the performance of the models, we will use two separate test sets, generated with a fixed set of parameters:

The images within each dataset differ only by the realizations of the noise. We will be computing the bias and the standard deviations for our two main models:

The resulting performance metrics on the test datasets look as follows:

When we train a network, our objective is to minimize the loss on a validation dataset. However, this does not help us judge how close or far we are from the best results we can achieve. A neural network is nothing but a "function estimator", with multiple parameters, and our network is estimating 5 different parameters. We also know that the total variance of any unbiased estimator is bounded by the Cramér–Rao lower bound. Here we can see our results vis-a-vis this lower bound.

Firstly, our estimator is biased, and in theory, the variance can be made smaller at the cost of increasing the bias. At the same time, our biases are mostly less than the standard deviation, and the standard deviations of our model's errors are very close to the Cramér-Rao bounds.

Baseline Model¶

The purpose of the baseline model is to benchmark a new technique against a conventional one. In our case, the conventional approach is to carry out non-linear optimization for 5 parameters with the noisy galaxy image as the input.

We use scipy.optimize.minimize to perform optimization with L-BFGS-B algorithm. Nelder-mead & BFGS do not support bounds and lead to generation of "invalid" values of parameters during the minimization process when we use GalSim to produce model images.

The objective function we minimize takes the noisy image, iteratively generates new images using GalSim, and calculates the loss. We compute the loss as the (negative) log-likelihood of the Image given the current iteration of parameter values, i.e, the $\log$ of $P(Image | Parameters)$. We take the negative of this to get the minimization objective.

Overall, the baseline model does a reasonable estimation and comes close to the Cramér-Rao bound. However, an iteration of minimization / parameter estimation takes about 35 seconds on JupyterHub, whereas the neural networks can process a million images in under 15 minutes. We include the code for the baseline model as part of our submission.

Limitations¶

We need to consider the limitations of our model and the approach in general. Besides the demonstrated bias, which is significant for some of the parameters, the most important fact is that our models have been trained and evaluated on synthetic data:

• The limited range of values the models have been trained on might lead to unexpected or overconfident results outside of that range of values.
• Beyond the synthetic dataset, we currently do not have an understanding of the behavior of the model. Many questions remain open:
  • Can we train a model on synthetic data and expect at least some of the learning to be transferable to the images of real galaxies? What is the amount of covariate shift between the distributions of generated and real images?
  • If the models need to be retrained on real data, were would the labels for supervised learning be coming from?
• Lastly, a neural network is inherently a black box tool. Compared to traditional estimation methods, deep learning algorithms might hinder the interpretation of results. This may be especially important for natural sciences. The networks might also be failing silently, without us even noticing that.

Eventually, our models are proof of concept, and the applicability to real-life astrophysical tasks needs to be thoroughly investigated.

Future Work¶

We would like to conclude by outlining some of the directions for possible future research. This includes:

• Grad-CAM and Saliency maps: examining these may provide insights for improving the estimation of parameters like Sérsic Index and Sérsic Radius.
• Attention layers: focusing a network's attention on particular regions on the image might prove useful for estimating galaxy parameters.
• Additional training: can help improve the performance of some of our models.
• Ensuring robustness: by incorporating uncertainty into our point estimates using Bayesian or frequentist methods we can make sure that the models can safely be used on out-of-distribution data.
• Real galaxy images: training and testing on real data will demonstrate whether our proof-of-concept approach is applicable to actual astrophysical tasks.

References¶

1. A. W. Graham, S. P. Driver, 2005, A concise reference to (projected) Sersic $R^{1/n}$ quantities, including Concentration, Profile Slopes, Petrosian indices, and Kron Magnitudes
2. D. Tuccillo, M. Huertas-Company, E. Decencière, S. Velasco-Forero, H. Domínguez Sánchez, P. Dimauro, 2017, Deep learning for galaxy surface brightness profile fitting
3. S. Madireddy, N. Li, N. Ramachandra, P. Balaprakash, S. Habib, 2019, Modular Deep Learning Analysis of Galaxy-Scale Strong Lensing Images
4. François Chollet, 2017, Xception: Deep Learning with Depthwise Separable Convolutions
5. Peter Erwin, 2014, Imfit: A Fast, Flexible New Program for Astronomical Image Fitting
6. I. Breda, P. Papaderos, J. M. Gomes, S. Amarantidis, 2019, A new fitting concept for the robust determination of Sérsic model parameters
7. B. Lim, S. Son, H. Kim, S. Nah, K. M. Lee, 2017, Enhanced Deep Residual Networks for Single Image Super-Resolution
8. Y. Zhang, Y. Tian, Y. Kong, B. Zhong, Y. Fu, 2020, Residual Dense Network for Image Restoration