{'fraction_non_alphanumeric': 0.03457849570923776, 'fraction_numerical': 0.019687026754164564, 'mean_word_length': 4.684648493543759, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pixelated Reconstruction of Gravitational Lenses using Recurrent Inference Machines

Modeling strong gravitational lenses in order to quantify the distortions in the images of background sources and to reconstruct the mass density in the foreground lenses has traditionally been a difficult computational challenge. As the quality of gravitational lens images increases, the task of fully exploiting the information they contain becomes computationally and algorithmically more difficult. In this work, we use a neural network based on the Recurrent Inference Machine (RIM) to simultaneously reconstruct an undistorted image of the background source and the lens mass density distribution as pixelated maps. The method we present iteratively reconstructs the model parameters (the source and density map pixels) by learning the process of optimization of their likelihood given the data using the physical model (a ray-tracing simulation), regularized by a prior implicitly learned by the neural network through its training data. When compared to more traditional parametric models, the proposed method is significantly more expressive and can reconstruct complex mass distributions, which we demonstrate by using realistic lensing galaxies taken from the cosmological hydrodynamic simulation IllustrisTNG.


Introduction
Strong gravitational lensing is a natural phenomenon through which multiple distorted images of luminous background objects, i.e. early-type star-forming galaxies, are formed by massive foreground objects along the line of sight (e.g., Vieira et al., 2013;Marrone et al., 2018;Rizzo et al., ICML 2022 Workshop on Machine Learning for Astrophysics, Baltimore, Maryland, USA, 2022. Copyright 2022 by the author(s). 2020; Sun et al., 2021). These distortions are tracers of the distribution of mass in foreground objects, independent of the electromagnetic behaviour of these overdensities. As such, this phenomenon offers a powerful probe of the distribution of dark matter and its properties outside of the Milky Way (e.g., Dalal & Kochanek, 2002;Treu & Koopmans, 2004;Hezaveh et al., 2016;Gilman et al., 2020;.
Lens modeling is the process of inferring the parameters describing both the mass distribution in the foreground lens and the light emitted by the background source. This has traditionally been a time-and resource-consuming procedure.
A common practice to model the mass of lensing galaxies is to assume that the density profiles follow simple parametric forms, e.g., a power law ρ ∝ r −γ . These profiles generally provide a good fit to low-resolution data and are easy to work with due to their small number of parameters (e.g., Koopmans et al., 2006;Barnabè et al., 2009;Auger et al., 2010). However, as high-resolution and high signal-to-noise ratio (SNR) images become available, lens analysis with simple models requires introducing additional parameters representing the complexities in the lensing galaxies and their immediate environments (e.g., Sluse et al., 2017;Wong et al., 2017;Birrer et al., 2019;Rusu et al., 2020;Li et al., 2021). This approach becomes intractable as the quality of images increases. For example, no simple parametric model of the Hubble Space Telescope (HST) Wide Field Camera 3 (WFC3) images of the Cosmic Horseshoe (J1148+1930) -initially discovered by Belokurov et al. (2007) -has been able to model the fine features of the extended arc (e.g., Bellagamba et al., 2016;Cheng et al., 2019;Schuldt et al., 2019).
In this work, we develop a method for pixelated strong gravitational lensing mass and source reconstruction, allowing it to reconstruct complex distributions. Our method is based on the Recurrent Inference Machine (RIM, Putzky & Welling, 2017), which proposes to learn an iterative inference algorithm, moving away from hand-chosen inference algorithms and hand-crafted priors. In this framework, the prior is implicit in the dataset used to train the neural network. We also present a new architecture based on the original RIM to allow the inference of pixelated maps for this highly non-linear and under-constrained problem.

Data
The background source brightness distributions are taken from the Hubble Space Telescope (HST) COSMOS field Scoville et al., 2007), acquired in the F814W filter. A dataset of mag limited (F814W < 23.5) deblended galaxy postage stamps  was compiled as part of the GREAT3 challenge (Mandelbaum et al., 2014). The data is publicly available (Mandelbaum et al., 2012), and the preprocessing is done through the open-source software GALSIM . The final set has 13 321 galaxy images cropped to 128 2 pixels and with a flux greater than 50 photons cm −2 s −1 . We split this set into a training set (90%) and a test set (10%) before data augmentation and denoising with an autoencoder (Vincent et al., 2008).
The projected surface density maps (convergence) of lensing galaxies were made using the redshift z = 0 snapshot of the IllustrisTNG-100 simulation (Nelson et al., 2019) in order to produce physically realistic realizations of dark matter and baryonic matter halos. We selected 1604 profiles, split into a training set (90%) and a test set (10%), with the criteria that they have a total dark matter mass of at least 9 × 10 11 M . We then collected all dark matter, gas, stars, and black hole particles from the profiles. We then compute smoothed projected surface density maps with an adaptive Gaussian kernel following the prescriptions from Aubert et al. (2007) and Rau et al. (2013). The final training set is composed of 3 different projections (xy, xz and yz) of each profiles, rendered on a pixelated grid with a resolution of 0.55 kpc/h and 128 2 pixels. Several data augmentation rounds, including rescaling each convergence maps randomly to produce a set with an Einstein radius uniformly distributed θ E ∼ U(0.5 , 2.5 ), were used to increase the number of profiles to 50 000 for training a VAE and the RIM.

Data Augmentation with VAEs
When working with limited data, data augmentation is crucial to ensure that the trained model is robust against perturbations -like rotations of images -which are not implicitly included as symmetries in the architecture of the model. We trained a variational auto-encoder (VAE, Kingma & Welling, 2013) for data augmentation of the source maps and another VAE for the convergence maps.
Direct optimization of the ELBO loss for VAEs can prove difficult because the reconstruction term could be relatively weak compared to the Kullback Leibler (KL) divergence term (Kingma & Welling, 2019). To alleviate this issue, we follow the work of Bowman et al. (2015) and Kaae Sønderby et al. (2016) in setting a warm-up schedule for  The fine-tuning procedure is able to recover this subhalo because of its strong signal in the lensed image and reduces the residuals to noise level.
the KL term in the ELBO loss, starting from β = 0.1 up to β max . Following the work of Lanusse et al. (2021), we also introduce an 2 penalty between the input and output of the bottleneck fully-connected layers to encourage an identity map between them. This regularisation term is slowly removed during training.

Raytracing
Simulations of lensed images are produced using a raytracing code, which maps the brightness distribution of background sources to the observed coordinates. The foreground pixel coordinates θ i and the source pixel coordinates β i are related by the lens equation where α is the deflection angle. The deflection angle is calculated from the projected surface density field κ (also commonly refered to as the convergence) by the integral The intensity of a pixel in a simulated lensed image is obtained by bilinear interpolation of the source brightness distribution at the coordinate β i . The integral in equation (2) is computed in near-linear time using Fast Fourier Transforms (FFT).
A blurring operator -i.e., a convolution by a point spread function -is then applied to the lensed image to replicate the response of an imaging system. This operator is implemented as a GPU-accelerated matrix operation since the y Observation Source Convergencê blurring kernels used in this paper have a significant proportion of their energy distribution encircled inside a small pixel radius.

Recurrent Inference Machine
The RIM (Putzky & Welling, 2017) is a form of learned gradient-based inference intended to solve inverse problems of the form where y is a vector of noisy lensed images, f is a function encoding the physical model, x is a vector of parameters of interest, and N is a vector of additive noise. This framework has been applied in the context of linear inverse problems, where the function f can be represented in a matrix form, in particular in the cases of under-constrained problems for which the prior on the parameters x, p(x), is either intractable or hard to compute Lønning et al., 2019). The use of the RIM to solve non-linear inverse problems was first investigated in (Modi et al., 2021). In our case, the inverse problem mapping function F is non-linear w.r.t. the convergence parameters which is a consequence of the non-linearity of equation (1).
The governing equation for the RIM is a recurrent relation that takes the general form where p(y | x) is an isotropic gaussian likelihood function, characterized by the noise standard deviation σ. By minimizing the weighted mean squared loss the RIM learns to optimize the parameters x given a likelihood function. Unlike previous work (Andrychowicz et al., 2016;Putzky & Welling, 2017;Morningstar et al., 2018;Lønning et al., 2019), the data vector y -or observation -is fed to the neural network in order to learn the initialization of the parameters,x (0) = g ϕ (0, y, 0), as well as their optimization. We found in practice that this significantly improves the performance of the model for our problem and it avoids situations where the model would get stuck in local minima at test time due to poor initialization.
The RIM used in this work is designed based on a U-net architecture (Ronneberger et al., 2015). The most important aspect of our implementation is the use of Gated Recurrent Units (Cho et al., 2014) placed in each skip connection which guides the reconstruction independently at different levels of resolution. The gradient of the likelihood is computed using automatic differentiation. Following Modi et al. (2021), we preprocess the gradients using the Adam algorithm (Kingma & Welling, 2013).

Fine-tuning
Once the RIM is trained, we can treat the RIM optimization procedure as a baseline estimator of the parameters x given a noisy observation y. We now concern ourselves with a strategy to improve this estimator. This is important for observations with high SNR, for which the estimator must be extremely accurate to model all the fine features present in the arcs. The fine-tuning objective is to minimize directly the likelihood over each time steps of the RIM: log p(y |x (t) ) + log p(ϕ) . (6) This objective function makes no use of labels, meaning that only an observation y is required to fine-tune the RIM. This allows us to use this objective at test time, at which point the RIM is trained to reconstruct this specific observation.
We use elastic weight consolidation (EWC, Kirkpatrick et al., 2016) as prior over the model parameters. To compute the Fisher matrix in EWC, it is necessary to sample from a distribution of lensing systems that are conditioned on the observation. This is accomplished by sampling the latent space of both the source VAE and the convergence VAE near the latent code of the baseline prediction of the RIM. More details regarding this procedure are given in appendix E. maps with traditional maximum-likelihood methods using analytical profiles would require significant user input and considerable computational resources due to large parameter degeneracies. The mock observations as well as the reconstructed lensed images are also shown, alongside the residuals of the reconstructions and the χ 2 statistic. The first 3 reconstructions have statistically significant residuals, with no pixels exceeding the 5σ threshold. The last reconstruction, which is arguably the most complex to perform, has a few pixels reaching 5σ.

Results and Discussion
In addition to a visual inspection of the reconstructed sources and convergences, we compute the coherence spectrum to quantitatively assess the quality the reconstructions γ(k) = P 12 (k) where P ij (k) is the cross power spectrum of images i and j at the wavenumber k. Figure 4 shows the mean value and the 68% inclusion interval of those spectra for the convergence and source maps in a test set of 3000 examples. The fine-tuning procedure, shown in red, is able to improve significantly the coherence of the baseline background source, shown in black, at all scales. The coherence spectrum of the convergence remains unchanged by the fine-tuning procedure. Still, we note that many examples in the dataset showcase significant improvement which we illustrate in Figure 1. The results obtained here demonstrate the effectiveness of machine learning methods for inferring pixelated maps of the distribution of mass in lensing galaxies. Since this is a heavily under-constrained problem, stringent priors are needed to avoid overfitting the data, a task that has traditionally been difficult to accomplish (e.g., Saha & Williams, 1997). The model proposed here can implicitly learn these priors from a set of training data.
The flexible and expressive form of the reconstructions means that, in principle, any lensing system (e.g., a single simple galaxy, or a group of complex galaxies) could be analyzed by this model, without any need for pre-determining the model parameterization. This is of high value given the diversity of observed lensing systems, and their relevance for constraining astrophysical and cosmological parameters.
The source code, as well as the various scripts and parameters used to produce the model and results is available as open-source software under the package Censai 1 . The model parameters, as well as the convergence maps and the background sources used to train these models, the test set examples and the reconstructions results are available as open-source datasets hosted by Zenodo 2 . This research made use of Tensorflow (Abadi et al., 2015), Tensorflow-Probability (Dillon et al., 2017), Numpy , Scipy , Matplotlib (Hunter, 2007), scikit-image (Van der Walt et al., 2014), IPython (Pérez & Granger, 2007), Pandas (Wes McKinney, 2010; pandas development team, 2020), Scikit-learn (Pedregosa et al., 2011), Astropy (Astropy Collaboration et al., 2013; and GalSim .

C. VAE Architecture and optimisation
For the following architectures, we employ the notion of level to mean layers in the encoder and the decoder with the same spatial resolution. In each level, we place a block of convolutional layers before downsampling (encoder) or after upsampling (decoder). These operations are done with strided convolutions like in the U-net architecture of the RIM.

D. RIM architecture and optimisation
The notion of link function Ψ : Ξ → X , introduced by Putzky & Welling (2017), is an invertible transformation between the network prediction space ξ ∈ Ξ and the forward modelling space x ∈ X . This is a different notion from preprocessing, discussed in section 2.1, because this transformation is applied inside the recurrent relation 4 as opposed to before training. In the case where the forward model has some restricted support or it is found that some transformation helps the training, then the link function chosen must be implemented as part of the network architecture as shown in the unrolled computational graph in Figure 7. Also, the loss L ϕ must be computed in the Ξ space in order to avoid gradient vanishing problems when Ψ is a non-linear mapping, which happens if the non-linear link function is applied in an operation recorded for backpropagation through time (BPTT).
For the convergence, we use an exponential link function with base 10:κ = Ψ(ξ) = 10 ξ . This Ψ encodes the non-negativity of the convergence. Furthermore, it is a power transformation that leaves the linked pixel values ξ i normally distributed, thus improving the learning through the non-linearities in the neural network. The pixel weights w i in the loss function (5) are chosen to encode the fact that the pixel with critical mass density (κ i > 1) have a stronger effect on the lensing configuration than other pixels. We find in practice that the weights encode this knowledge in the loss function and improved both the empirical risk and the goodness of fit of the baseline model on early test runs.
For the source, we found that we do not need a link function -its performance is generally better compared to other link function we tried like sigmoid and power transforms -and we found that the pixel weights can be taken to be uniform, i.e. w i = 1 M .
In the first optimisation stage, we trained 24 different architectures from a small set of valid hyperparameters previously identified for approximately 4 days (wall time using a single Nvidia A100 gpu). Following this first stage, 4 architectures were deemed efficient enough to be trained for an additional 6 days. y g ϕ

Learned initializationξ
(1) Figure 7. Unrolled computational graph of the RIM. Operations along solid arrows are being recorded for BPTT, while operations along dashed arrows are not. The blue arrows are only used for optimisation during training. During fine-tuning or testing, the loss is computed only as an oracle metric to validate that our methods can recover the ground truth.

E. Fine-Tuning
We follow the work of Kirkpatrick et al. (2016) to define a prior distribution over ϕ that address the issue of catastrophic forgetting (McCloskey & Cohen, 1989;Ratcliff, 1990): where diag(I(ϕ D )) is the diagonal of the Fisher information matrix encoding the amount of information that some set of gravitational lensing systems from the training set similar to the observed test task carries about the baseline RIM weights ϕ D -the parameters that minimize the empirical risk over the training dataset D. The Lagrange multiplier λ is tuning our estimated uncertainty about the neural network weights for the particular task at hand.
We can understand the need for a conditional sampling distribution by looking at the posterior of the RIM parameters. Suppose we are given a training set D and a test task T which are conditionally independent given ϕ and have uniform priors, then the posterior of the RIM parameters ϕ can be rewritten using the Bayes rule as The sampling distribution in this expression appears as the conditional p(T | D), which can also be viewed as the set of examples from the training set similar to the test task by rewriting it as p(T | D) ∝ p(D | T ). The EWC term is then derived by a Laplace approximation of the prior p(ϕ | D) around ϕ D , which we also take to be proportional to the training loss, the likelihood of each time steps and an 2 loss Each reconstruction is performed by fine-tuning the baseline model on a test task composed of an observation vector, a PSF and a noise amplitude. In practice, fine-tuning the test set of 3 000 examples can be accomplished in parallel so as to be done in at most a few days by spreading the computation on ∼ 10 Nvidia A100 GPUs (or 10 hours on ∼ 100 GPUs). Each reconstruction uses at most 2000 steps, which turns out to be approximately 20 minutes (wall-time) per reconstruction. Early stopping is applied when the χ 2 reaches noise level (χ 2 = ν). Figure 8. Examples similar to the test task shown in Figure 3. They are sampled from the latent space of the source VAE and convergence VAE near the RIM baseline latent codeẑ (T ) and used to estimate diag(I(ϕ D )).