High Performance Computing for gravitational lens modeling: single vs double precision on GPUs and CPUs

TL;DR

This paper demonstrates that using high performance computing techniques, including mixed precision algorithms, significantly accelerates gravitational lens modeling on CPUs and GPUs, reducing energy consumption and enabling cost-effective scientific analysis.

Contribution

The authors develop a mixed precision algorithm for lens modeling that maintains accuracy while greatly improving performance on CPUs and GPUs, including consumer-grade GPUs.

Findings

Speedup factors of up to 170 on CPU clusters.

Energy consumption reduced by up to 98%.

Consumer GPUs achieve near high-end GPU performance with mixed precision.

Abstract

Strong gravitational lensing is a powerful probe of cosmology and the dark matter distribution. Efficient lensing software is already a necessity to fully use its potential and the performance demands will only increase with the upcoming generation of telescopes. In this paper, we study the possible impact of High Performance Computing techniques on a performance-critical part of the widely used lens modeling software LENSTOOL. We implement the algorithm once as a highly optimized CPU version and once with graphics card acceleration for a simple parametric lens model. In addition, we study the impact of finite machine precision on the lensing algorithm. While double precision is the default choice for scientific applications, we find that single precision can be sufficiently accurate for our purposes and lead to a big speedup. Therefore we develop and present a mixed precision algorithm which only uses double precision when necessary. We measure the performance of the different implementations and find that the use of High Performance Computing Techniques dramatically improves the code performance both on CPUs and GPUs. Compared to the current LENSTOOL implementation on 12 CPU cores, we obtain speedup factors of up to 170. We achieve this optimal performance by using our mixed precision algorithm on a high-end GPU which is common in modern supercomputers. We also show that these techniques reduce the energy consumption by up to 98%. Furthermore, we demonstrate that a highly competitive speedup can be reached with consumer GPUs. While they are an order of magnitude cheaper than the high-end graphics cards, they are rarely used for scientific computations due to their low double precision performance. Our mixed precision algorithm unlocks their full potential. The consumer GPU delivers a speedup which is only a factor of four lower than the best speedup achieved by a high-end GPU.

Figures

Tables

Table 1. Table 1: Theoretical maximum computing performance for our used CPU and GPU models. These values can only serve as a rough indicator of expected performance, as the real application performance will depend on many parameters such as the used algorithm and its implementation. We list the base frequency for the CPU while we use the boost frequency for the GPUs, as the CPU typically reaches the boost frequency only on a few cores and and not on all cores simultaneously. We compute the CPU maximum computing performance using the following formula: Two operations per cycle × \times frequency × \times AVX vectorization × \times number of cores [ Besl , 2013 ] . Note that the AVX factor for SP is two times larger than for DP. We use the same formula, but without the AVX factor, for GPUs. Graphics cards have a different number of cores for SP and DP computations and thus a different maximum performance. The number of GPU cores listed in the table is for SP computations. Due to the lower number of DP cores, the Nvidia GTX 1080 Ti’s GP102 GPU has thirty-two times less performance in DP computations than in SP [e.g., Harris , 2016 ] , while the P100’s GP100 GPU and the V100’s GV100 GPU are two times slower [ Nvidia Corporation , 2016 , 2017b ] . The number of cores and the frequencies are taken from Intel Corporation [ 2014 ] and Nvidia Corporation [ 2017a , b , 2016 ] .

	Intel Xeon E5-2680 v3	Nvidia GTX 1080 Ti	Nvidia P100	Nvidia V100
	12 cores, 2.50 GHz	3584 cores, 1582 MHz	3584 cores, 1480 MHz	5120 cores, 1530 MHz
Double precision	240 GFLOPS	354 GFLOPS	5304 GFLOPS	7833 GFLOPS
Single precision	480 GFLOPS	11340 GFLOPS	10609 GFLOPS	15667 GFLOPS

Table 2. Table 2: Benchmark results of the gradient computation implementations for a HFF-like lens. We show the double and mixed precision run times for the HPC-optimized CPU version using 12 cores and the three GPU models. The mixed precision algorithm for the GPUs uses a hybrid CPU/GPU implementation. Column three presents the measured run time advantage of mixed precision over double precision. Note that the mixed precision algorithm requires a substantial amount of additional computations compared to the double precision algorithm, as it must check for which pixels single precision is accurate enough and for which ones the gradient must be recomputed in double precision. Despite this overhead, the mixed precision implementation is the fastest for all four hardware devices. The fourth column shows the speedup of the mixed precision implementation with respect to the best currently achievable speed of LENSTOOL .

	Run time	Run time	Run time reduction	Speedup factor
	Double Precision	Mixed Precision	Double $\to$ Mixed Precision	compared to
	(seconds)	(seconds)	(%)	LENSTOOL
CPU	8.1	7.40	9	9
GTX	6.8	1.58	77	42
P100	0.71	0.58	18	115
V100	0.46	0.39	15	171

Table 3. Table 3: Energy comparison of the different gradient computation implementations for one run of the benchmark. We estimate the energy to solution by multiplying the Thermal Design Power (TDP) of the different hardware devices with the respective benchmark run times. In the case of the mixed precision GPU implementations, which use a hybrid CPU/GPU approach, we add the TDPs of the GPU and the CPU. The TDP values are taken from Intel Corporation [ 2014 ] and Nvidia Corporation [ 2016 , 2017a , 2017b ] . The last column shows the percentage of energy saved by using the double or mixed precision algorithm instead of the current LENSTOOL software.

	Hardware	Energy to solution	Energy to solution	Energy saved	Energy saved compared
	TDP	Double Precision	Mixed Precision	Double $\to$ Mixed Precision	to LENSTOOL
	(Watt)	(Joule)	(Joule)	(%)	(%)
LENSTOOL	120	8016	-	-	-
CPU	120	972	888	9	88 / 89
GTX	250	1700	585	66	79 / 93
P100	300	213	244	-15	97 / 97
V100	300	138	164	-19	98 / 98

Equations

$$\mathbf{\beta }=\mathbf{\theta }-\mathbf{\alpha }(\mathbf{\theta }),$$

$$\mathbf{\alpha }=\nabla \psi ,$$

$$\psi (\mathbf{\theta })=\frac{1}{\pi }{\int }_{{\mathbb{R}}^2}{\text{d}}^2{\theta }^{\prime }\kappa ({\mathbf{\theta }}^{\prime })\ln |\mathbf{\theta }-{\mathbf{\theta }}^{\prime }|,$$

$$\kappa (\mathbf{\theta })=\frac{\Sigma (\mathbf{\theta })}{{\Sigma }_{\text{crit}}},$$

$${\Sigma }_{\text{crit}}=\frac{c^2}{4\pi G}\frac{D_{\text{OS}}}{D_{\text{OL}}{D}_{\text{LS}}},$$

$$\Sigma (\mathbf{\theta })=\int \text{d}z\rho (\mathbf{\theta },z),$$

$$\Sigma (\mathbf{\theta })=\frac{{\sigma }_v^2}{2G{D}_{\text{OL}}|\mathbf{\theta }|},$$

$$\kappa (\mathbf{\theta })=\frac{{\theta }_E}{2|\mathbf{\theta }|},$$

$$\displaystyle\theta_{\text{E}}=4\pi\Big{(}\frac{\sigma_{v}}{c}\Big{)}^{2}\frac{D_{\text{LS}}}{D_{\text{OS}}},$$

$$|\mathbf{\alpha }|={\theta }_{\text{E}}.$$

$${\chi }^2=\sum _i\sum _j\frac{(x_{\text{obs},ij}-x_{ij}{)}^2}{{\sigma }_{ij}^2}$$

$$\Psi =\frac{D_{\text{OS}}}{D_{\text{LS}}}\psi$$

$${\beta }_1={\theta }_1-\frac{D_{\text{LS}}}{D_{\text{OS}}}(\nabla \Psi (\mathbf{\theta }){)}_1,$$

$${\beta }_2={\theta }_2-\frac{D_{\text{LS}}}{D_{\text{OS}}}(\nabla \Psi (\mathbf{\theta }){)}_2.$$

$$\psi (\mathbf{\theta })=\sum {\psi }_{\text{cluster}}(\mathbf{\theta })+\sum {\psi }_{\text{galaxy}}(\mathbf{\theta }),$$

$$\nabla \Psi (\mathbf{\theta })=\sum \nabla {\Psi }_{\text{cluster}}(\mathbf{\theta })+\sum \nabla {\Psi }_{\text{galaxy}}(\mathbf{\theta }).$$

$$R=\sqrt{{\theta }_{1,ϵ}^2+{\theta }_{2,ϵ}^2},$$

$$\displaystyle\phi=\arctan\Big{(}\frac{\theta_{2,\epsilon}}{\theta_{1,\epsilon}}\Big{)},$$

$${\theta }_{1,ϵ}=\sqrt{a_{1,ϵ}}{\theta }_1,$$

$${\theta }_{2,ϵ}=\sqrt{a_{2,ϵ}}{\theta }_2,$$

$$a_{1,ϵ}=1-ϵ,$$

$$a_{2,ϵ}=1+ϵ.$$

$${\psi }_{ϵ}(\mathbf{\theta })=\psi (R,\varphi ),$$

$${\alpha }_{1,ϵ}(\mathbf{\theta })=|\mathbf{\alpha }(R,\varphi )|\sqrt{a_{1,ϵ}}\cos (\varphi ),$$

$${\alpha }_{2,ϵ}(\mathbf{\theta })=|\mathbf{\alpha }(R,\varphi )|\sqrt{a_{2,ϵ}}\sin (\varphi ).$$

$$\displaystyle\big{(}\nabla\Psi_{\epsilon}\big{)}_{1}=(1-\epsilon)~{}b_{0}~{}\frac{\theta_{1}}{R},$$

$$\displaystyle\big{(}\nabla\Psi_{\epsilon}\big{)}_{2}=(1+\epsilon)~{}b_{0}~{}\frac{\theta_{2}}{R},$$

$$b_{0}=4\pi\Big{(}\frac{\sigma_{v}}{c}\Big{)}^{2}.$$

$$x=\sigma \times {\bar{x}}_2\times 2^e,$$

$${\delta }_1({\theta }_1,{\theta }_2)$$

$${\delta }_2({\theta }_1,{\theta }_2)$$

$$|0.03\text{arcsec}-{\delta }_i({\theta }_1,{\theta }_2)|<0.0092\text{arcsec}$$

$$|\Delta x_{\text{source},1}({\theta }_1,{\theta }_2)|$$

$$=|\Delta x-{\delta }_1({\theta }_1,{\theta }_2)|,$$

$${\delta }_3({\theta }_1,{\theta }_2)$$

$${\delta }_4({\theta }_1,{\theta }_2)$$

$$|{\delta }_i({\theta }_1,{\theta }_2)|>0.0104\text{arcsec}$$

Full text

††journal: Astronomy and Computing

1 Introduction

The $\Lambda$CDM cosmology standard model describes our universe with great precision, but it also introduces two unknown quantities, Dark Energy and Dark Matter. They dominate the energy density of the universe [e.g., Planck Collaboration et al., 2016], but their physical nature has so far remained elusive. Consequently their study is one of the prime targets of cosmological research.

Strong gravitational lensing is a unique tool for cosmology, as it is sensitive to the total (baryonic and dark) matter density and thus it probes Dark Matter directly [see e.g. Kneib and Natarajan, 2011, Schneider et al., 2006, for reviews]. Its application has led to constraints on cosmological parameters [e.g., Jullo et al., 2010, Bonvin et al., 2017] and the Dark Matter self-interaction cross-section [e.g., Randall et al., 2008, Bradač et al., 2008]. In addition, the magnification effect of a strong gravitational lens can be used to study the early universe and to constrain its reionization [e.g., Atek et al., 2015, Ishigaki et al., 2015]. Strong lensing requires deep, high-resolution data and indeed the field has prospered thanks to programs like the Cluster Lensing And Supernova survey with Hubble [CLASH, Postman et al., 2012] and the Hubble Frontier Fields [HFF, Lotz et al., 2017].

Future missions like Euclid, the James Webb Space Telescope (JWST), the Large Synoptic Survey Telescope (LSST), and the Wide Field Infrared Survey Telescope (WFIRST) will provide a large amount of excellent data sets for lensing. These will enable the lensing community to further push the boundaries of cosmological knowledge. This, however, will only be feasible if we are able to efficiently harvest the wealth of information available in the data. This will be a challenge, e.g. due to the amount of data available or the high quality of the data, which permits the creation of lens models with a high level of detail and precision, but also requires more computing time. Gravitational lensing software and pipelines will have to be ready to process these data sets in a reasonable amount of time.

Therefore we are currently redesigning the strong lensing software LENSTOOL111Open source software publicly available at https://projets.lam.fr/projects/lenstool/wiki [Jullo et al., 2007, Kneib et al., 1996]. LENSTOOL has been successfully used to model many strong lensing galaxy clusters with high precision [see e.g. Jauzac et al., 2014, 2015, Limousin et al., 2016, for recent lens models] and has been serving the lensing community for more than two decades. In a recent comparison of strong lensing modeling software it has performed very well [Meneghetti et al., 2016]. However, the HFF data sets provided the greatest number of lensing constraints so far and this posed a computing challenge for LENSTOOL. It took several weeks to compute a single HFF lens model and several different lens models from different priors are required to find the best fitting model.

The new version is designed to meet this computation challenge by using High Performance Computing (HPC) methods. The LENSTOOL algorithms are very well suited for massive parallelism and we employ this technique to accelerate the computations. While we focus on lensing by galaxy clusters, a recent publication by Tessore et al. [2016] has shown that massive parallelism holds also great promise for the modeling of galaxy lenses. In this paper, we discuss the central lensing algorithm of LENSTOOL and in particular the performance-critical computation of deflection potential gradients. We have implemented the gradient computation algorithm using two different hardware types in order to be able to compare performance. The first version is a highly optimized and parallelized CPU code and the second version uses Graphics Processing Unit (GPU) acceleration.

During the development phase, we have asked ourselves the question: Can we do even better by using single precision instead of the commonly used double precision? The computing power of both CPUs and GPUs is higher for single precision [see e.g. Eijkhout et al., 2016, Besl, 2013], so we can expect a significant performance improvement. The downside is that this might lead to an error in our results if single precision is not precise enough for our computations. Therefore we use error propagation to compute the impact of single precision on the results of the central lensing algorithm. In addition, we measure and compare the single and double precision performance of both CPU and GPU implementations.

The paper is organized as follows: Section 2 gives a concise introduction to strong gravitational lensing and the LENSTOOL algorithms. It also presents the CPU and GPU implementations. Section 3 introduces the single and double precision floating-point representations and investigates if single precision is precise enough for our computations. We present and compare the performance measurements of the single and double precision CPU and GPU implementations in section 4. We discuss our results in section 5 and conclude in section 6.

2 Accelerating lensing with massive parallelism

2.1 Strong gravitational lensing

Galaxies and galaxy clusters are so dense that they locally deform space-time. As a result, they can act as a lens for background objects, which are magnified and distorted or even multiply imaged. Lensing also changes the locations at which we observe the lensed images on the sky so that they are typically not coincident with the locations at which we would observe the background sources in the absence of lensing. In practice, we can only observe the lensed images of a background source, but not the background source itself. However, the position of the background source on the sky can be calculated with the lens equation [see e.g. the reviews Kneib and Natarajan, 2011, Bartelmann and Schneider, 2001, for a derivation],

[TABLE]

where the two dimensional vectors $\boldsymbol{\beta}$, $\boldsymbol{\theta}$, and $\boldsymbol{\alpha}$ describe respectively the location of the source in the source plane, the location of the lensed image in the image plane, and the scaled deflection angle. Note that these quantities are angles. In the case of multiple images, the lens equation has more than one solution $\boldsymbol{\theta}$ for a fixed value of $\boldsymbol{\beta}$ [e.g., Bartelmann and Schneider, 2001]. The lens equation is derived under the assumption that we have only one lens, that the gravitational field is weak enough so that the field equations of General Relativity can be linearized, that we can use the Born approximation, and that the physical extent of the lens is small compared to the angular diameter distances between observer and lens, $D_{\text{OL}}$, and lens and source, $D_{\text{LS}}$.

The background objects are typically extended sources like galaxies. The shape of the lensed images will differ from the shape of the source, since the light coming from the object at coordinate $\boldsymbol{\beta}^{\prime}$ will be lensed slightly differently than the light coming from the object at coordinate $\boldsymbol{\beta}^{\prime\prime}$ [e.g., Bartelmann and Schneider, 2001]. Therefore we can use the lens equation to compute $\boldsymbol{\theta}$ for each coordinate $\boldsymbol{\beta}$ of the object and thus the shape of the lensed image due to distortion and magnification.

The scaled deflection angle $\boldsymbol{\alpha}$ is the gradient of the deflection potential $\psi$,

[TABLE]

and $\psi$ depends on the dimensionless projected surface mass density $\kappa$,

[TABLE]

where $\Sigma(\boldsymbol{\theta})$ is the projected surface mass density,

[TABLE]

and we defined the critical projected surface mass density $\Sigma_{\text{crit}}$. Here $\rho$ is the mass density, $c$ is the speed of light, $G$ is the gravitational constant, and $D_{\text{OS}}$ is the angular diameter distance between observer and source. We can see from these equations that $\boldsymbol{\alpha}$ and thus the strength of the lensing effect depend on the projected surface mass density. Therefore lensing probes the total surface mass density of the lens, including baryonic and Dark Matter components.

The value of $\kappa$ is a good indicator to distinguish the so-called “weak” and “strong” lensing regimes. In the case of weak lensing, the lensed image appears slightly magnified and distorted, and in the case of strong lensing, the image is strongly magnified and distorted and multiple images appear. A mass distribution which has $\kappa\geq 1$ somewhere produces multiple images for some source positions $\boldsymbol{\beta}$ [e.g., Bartelmann and Schneider, 2001]. In the case of cluster lensing, the strong lensing area and thus the multiple images are typically located in the central regions of the cluster, where the projected surface mass density is large enough [e.g., Kneib and Natarajan, 2011].

We will illustrate gravitational lensing with an example. We will look at a simple lens model, the Singular Isothermal Sphere (SIS), which we will use for the remainder of this paper as it has a relatively simple mathematical expression and is thus very instructive. The projected surface mass density is

[TABLE]

where $\sigma_{v}$ is the line-of-sight velocity dispersion of the “particles” (e.g. galaxies in a galaxy cluster), which are assumed to be in virial equilibrium [e.g., Bartelmann and Schneider, 2001]. Thus we have

[TABLE]

where we defined the Einstein deflection angle $\theta_{\text{E}}$. Using equations 2 and 3, we find that the magnitude of the scaled deflection angle is constant,

[TABLE]

We see that the lens equation has infinitely many solutions for $\boldsymbol{\beta}=\mathbf{0}$, namely each point on the circle with radius $\theta_{\text{E}}$. Therefore a background source at this location will be strongly lensed into a perfect Einstein ring.

2.2 Strong lensing algorithm

2.2.1 Overview

LENSTOOL models strong lensing galaxy clusters by using parametric models of the large-scale cluster halos and the galaxy-scale halos. In a typical merging cluster, we have two large-scale halos and hundreds of galaxy halos. Depending on the chosen parametric model, we have several free parameters such as $x$ and $y$ position, velocity dispersion, etc. for each halo. It is possible to constrain the range of the free parameters or to reduce their number, e.g. by assuming a scaling relation like the Faber-Jackson relation [Faber and Jackson, 1976] for galaxy-scale halos [Natarajan et al., 1998]. Nevertheless, the best lens model will still be hidden in a massive, high dimensional parameter space. LENSTOOL uses BayeSys3222Publicly available at http://www.inference.org.uk/bayesys/, a Bayesian Markov Chain Monte Carlo (MCMC) software package, to sample this parameter space [see Jullo et al., 2007, for a detailed description]. For each parameter combination probed by the MCMC, LENSTOOL computes the goodness of fit of the corresponding lens model given the observational data. It does this by modeling the lens with the given set of parameters and, using this model, lensing the observed multiple images into the source plane and subsequently back into the image plane, see figure 1. If the probed lens model is close to the true matter distribution, the re-lensed multiple image positions will be close to the observed multiple image positions and the goodness of fit parameter

[TABLE]

will be small [Jullo et al., 2007]. We denote the observed position of multiple image $j$ of multiple image system $i$ with $x_{\text{obs},ij}$, the re-lensed position with $x_{ij}$, and the error budget of the position with $\sigma_{ij}$. Since the parameter space probed by the MCMC is massive, it typically takes several weeks of computation time to find the best model for lenses with HFF-like data.

There are two ways to speed up the computation. The first is to speed up the MCMC, e.g. by parallelizing it. The second way is to speed up the $\chi^{2}$ computation. In this paper, we will focus on accelerating a crucial part of it, the gradient computation. Since we will have to take a very precise look at the algorithm when we compute the impact of single and double precision on its result, we will now present a detailed description.

2.2.2 Gradient computation in the $\chi^{2}$ algorithm

Before we present the $\chi^{2}$ algorithm, we reformulate the lens equation 1 by introducing

[TABLE]

and making the gradient dependence explicit:

[TABLE]

As a result, we only have to compute the constant $D_{\text{LS}}/D_{\text{OS}}$ once instead of for every image pixel. Note that the deflection potential at position $\boldsymbol{\theta}$ is a superposition of all cluster-scale and galaxy-scale deflection potentials $\psi_{\text{cluster}}$ and $\psi_{\text{galaxy}}$,

[TABLE]

[see e.g. Jullo et al., 2007] and as a result we have

[TABLE]

We see that the lens equation is computationally cheap to evaluate once the total gradient $\nabla\Psi$ is known. The computation of $\nabla\Psi$, however, involves potentially complicated gradient calculations for hundreds of potentials and as we will see in the next paragraph, it has to be computed for every pixel in our image. The Hubble Space Telescope Advanced Camera for Surveys (HST ACS) produces images with $4096\times 4096$ pixels at a pixel scale of $\approx 0.05~{}\text{arcsec/pixel}$ [Avila et al., 2017], which we can typically upsample to 0.03 arcsec/pixel [Lotz et al., 2017], so that HFF images have a total of $\approx 6730\times 6730~{}\text{pixels}\approx 45~{}\text{million~{}pixels}$. This shows that the computation of $\nabla\Psi$ is computationally expensive and an excellent target for speedup with HPC parallelism methods.

The $\chi^{2}$ computation is now performed as follows. We compute $\nabla\Psi$ for each pixel of the image plane. Then we loop over each multiple image $j$ in each multiple image system $i$. For each multiple image, we use equation 13 to compute the source coordinates, $\beta_{ij,1}$ and $\beta_{ij,2}$. Subsequently, we determine the barycenter of the sources of a given multiple image system $i$. If we are close to the true lens model, all multiple images will be mapped onto approximately the same source location, but in general the locations of the predicted sources can differ substantially, which makes it necessary to use the barycenter. In the next step, we re-lens the barycenter back into the image plane to obtain the locations of the multiple images predicted by the lens model. However, the lens equation cannot easily be inverted, so we have to find the locations in a different way. First, we divide each pixel in the image plane into two triangles, see figure 1. We do this because lensing always maps triangles onto triangles, but not squares onto squares. Second, we lens each triangle into the source plane by using equation 13 and we check if the barycenter is inside this triangle in the source plane. If it is, a predicted multiple image location in the image plane is found. Once we have found the locations of all predicted multiple images for all multiple image systems, we compute the $\chi^{2}$ according to equation 11.

The gradient calculations will naturally differ for different chosen parametric models. As an example, we present the gradient computation for a generalized form of the SIS, the pseudo-elliptical SIS (henceforth called SIE), in algorithm 1. It is necessary to generalize the parametric model, as we want to use this algorithm to model any SIS lens configuration by simply choosing the appropriate number of lenses and parameter values. We expand our treatment of the SIS in subsection 2.1 by following the procedure in Golse and Kneib [2002]. We introduce the pseudo-ellipticity of the deflection potential, $\epsilon$, and the coordinate system

[TABLE]

with

[TABLE]

Note that we call $\epsilon$ a pseudo-ellipticity, because the resulting elliptical shapes will only correspond to ellipses with classical ellipticty $\epsilon^{\prime}=1-b/a$, where $a$ and $b$ are the semi-major and semi-minor axes of the ellipse, for small values of $\epsilon$ [Golse and Kneib, 2002]. Therefore we assume in the following $\epsilon\ll 1$. The advantage of using a pseudo-elliptical parametric model is that it leads to relatively simple analytic expressions of the derived lensing quantities [Golse and Kneib, 2002]. Now we can simply calculate the values of the pseudo-elliptical deflection potential $\psi_{\epsilon}$ at location $\boldsymbol{\theta}$ by using the relation [Golse and Kneib, 2002]

[TABLE]

and analogous for $\Psi_{\epsilon}$. The resulting pseudo-elliptical shape is stretched along the $\theta_{1}$-axis, so that we have $\Phi=0$, where $\Phi$ is the counter-clockwise angle between the semi-major-axis and the $\theta_{1}$-axis. Algorithm 1 extends this approach to potentials with $\Phi\neq 0$ by using rotations. We obtain the following equations for the scaled deflection angle [Golse and Kneib, 2002],

[TABLE]

We can now combine the equations 2, 9, 10, 12, and 20 to derive the gradient expressions for the SIE,

[TABLE]

where we introduced the parameter

[TABLE]

The presented equations for the SIE always reduce to the previously presented equations for the spherical SIS for $\epsilon=0$.

2.3 CPU and GPU implementations

We implement a performance-optimized CPU version of the gradient computation in C++333C++ is a programming language standardized by the International Organization for Standardization, public website: https://isocpp.org/ by using the following techniques. First, we structure our data in the Structures of Arrays (SoA) format instead of the Arrays of Structures (AoS) format, see figure 2. To illustrate the difference, we take a look at the internal representation of five SIS potentials using SoA and AoS. In the AoS format, they are stored as an array comprised of five different data structures. Each data structure corresponds to a SIS potential and it contains the respective data of the SIS like $\theta_{\text{center},1}$, $\theta_{\text{center},2}$, and $\sigma_{v}$. In the SoA format, the potentials are stored in one data structure which consists of different arrays. Each array corresponds to a SIS parameter like $\theta_{\text{center},1}$ and array element number one of the $\theta_{\text{center},1}$-array would correspond to the $\theta_{\text{center},1}$ location of SIS number one, element number two to the $\theta_{\text{center},1}$ location of SIS number two etc. As a result, the SIS parameters occupy contiguous parts of the memory, which is usually beneficial for vectorized computations [e.g., Eijkhout et al., 2016, Besl, 2013]. Second, we use Advanced Vector Extensions (AVX) technology available on the latest CPU generations to harvest their built-in vectorization potential. For example, Besl [2013] obtained a significant speed-up by using SoA and AVX. CPU cores with AVX technology use registers with a width of 256 bits to process 8 single precision or 4 double precision values simultaneously, see figure 3 [Besl, 2013]. Note that the same operation has to be performed for each of the simultaneously processed data values. Therefore AVX is a Single Instruction Multiple Data (SIMD) parallelism technique [Eijkhout et al., 2016]. Third, we parallelize the computation using Open Multi-Processing (OpenMP)444OpenMP is an application programming interface managed by the non-profit OpenMP Architecture Review Board, public website: http://www.openmp.org on the outermost loop of algorithm 1. Each core of the multi-core CPU will now work on computing the total gradient for its assigned pixel and thus we compute the gradients for several pixels in parallel.

We implement the GPU version of the algorithm with CUDA555CUDA is a parallel computing platform and programming model for general computing on GPUs managed by Nvidia Corporation, public website: https://developer.nvidia.com/cuda-zone. First, we structure our data again in the SoA format. Second, we use the massively parallel architecture of GPUs to parallelize the gradient computation. Modern GPUs have many Streaming Multiprocessors (SM), which in turn consist of many Streaming Processors (SP), so the total amount of processor cores is computed by multiplying the two [e.g., Eijkhout et al., 2016]. The number of cores available depends on the GPU model, for example the Nvidia Tesla P100 (henceforth called P100) possesses 3584 cores for single precision computations [Nvidia Corporation, 2016]. In addition, GPUs are designed to be extremely efficient at switching between threads, where all threads in a single block of threads execute the same instruction [Eijkhout et al., 2016]. Therefore we can effectively use many more threads than we have GPU cores. Different blocks of threads can be processed independently. This GPU parallelism is called Single Instruction Multiple Thread (SIMT) [Eijkhout et al., 2016]. We use GPU threads to parallelize the outermost loop of algorithm 1. Each GPU thread computes the total gradient for its assigned pixel. Therefore we can compute the gradients for thousands of pixels simultaneously.

3 Finite machine precision errors in strong lensing

3.1 Single and double precision

Modern computers usually store real numbers in the IEEE 754 single precision floating-point representation (henceforth called SP) or the IEEE 754 double precision floating-point representation (henceforth called DP) [Institute of Electrical and Electronics Engineers, 2008, see e.g. Goldberg [1991] for an overview of floating-point arithmetic]. A real number $x\in\mathbb{R}$ in decimal representation is thus stored in a binary format,

[TABLE]

where the integer $e$ is the exponent, the sign $\sigma$ equals $+1$ or $-1$, and $\bar{x}_{2}$ is a binary number satisfying $(1)_{2}\leq\bar{x}_{2}<(10)_{2}$ [Institute of Electrical and Electronics Engineers, 2008]. Note that the binary number $\bar{x}_{2}$ consists of several integer digits $d\in\{0,1\}$, i.e. $\bar{x}_{2}=d_{0}.d_{1}d_{2}\dots d_{p-1}$. In the remainder of this paper we will denote the binary format by using the subscript $2$, so $(1)_{2}$ and $(10)_{2}$ correspond to the numbers $1$ and $2$ in decimal representation. For example, the number 2.25 would correspond to $\sigma=+1$, $\bar{x}_{2}=(1.001)_{2}$, and $e=(1)_{2}$666In practice, the leading bit of $\bar{x}_{2}$ would be implicit and $e$ would be stored as a biased exponent, but we can ignore such intricacies here to simplify the presentation. The number of digits in $\bar{x}_{2}$ is called the precision $p$ of the representation. According to IEEE 754, SP has a precision of $p=24$ digits and an exponent $-126\leq e\leq 127$, while DP has $p=53$ and $-1022\leq e\leq 1023$. SP values are stored using 4 bytes (= 32 bits) and DP values using 8 bytes (= 64 bits) [Institute of Electrical and Electronics Engineers, 2008]. As a result, DP can store a number $x$ with higher accuracy than SP, but this comes at the price of increased memory consumption and usually also reduced computing performance [e.g., Besl, 2013, Eijkhout et al., 2016]

Both DP and SP have only a limited amount of memory available and thus their accuracy is limited. We define the machine epsilon $\epsilon$ as the difference between 1 and the next larger number that can be stored using the given representation [Eijkhout et al., 2016]. For SP and DP we thus have respectively $\epsilon=2^{-23}\approx 1.2\times 10^{-7}$ and $\epsilon=2^{-52}\approx 2.2\times 10^{-16}$. These errors are so small that they might seem unimportant at first, but they will be magnified by the different computing operations performed in the course of an algorithm, so that they can become very large and relevant once the final result is obtained.

To illustrate this point, we now look at a hypothetical calculator777This illustration is inspired by an example in the lecture notes of Catalin Trenchea, available online at http://www.math.pitt.edu/~trenchea/math1070/MATH1070_2_Error_and_Computer_Arithmetic.pdf. For simplicity, it does not use SP or DP, but a decimal number representation with 6 digits precision and no exponent. We compute a relatively simple function, $f(x)=x\times(\sqrt{x+1}-\sqrt{x})$. For $x=50,000$, the result from the hypothetical calculator is 100, while the true result is 111.8, so we have a relative error of more than 10%. To understand this behavior, we take a look at the different steps which the calculator has to perform. It computes $\sqrt{50,001}$ and rounds the result to 6 digits (result: 223.609) and then it repeats these steps for $\sqrt{50,000}$ (result: 223.607). Therefore we have two rounding errors, but they are very small. However, now the calculator subtracts two almost equal numbers to obtain $000.002$, so only the last number of the result is a significant digit. We have lost a lot of accuracy which we cannot recover. The subsequent multiplication does not increase the error, but it propagates it into the final result. This example illustrates that even with the high precision available in modern computers, the result of a sufficiently long and complex algorithm can be significantly affected by the chosen number representation.

DP permits a much higher accuracy than SP and therefore it is tempting to simply use it for all computations. However, this accuracy comes at the price of computing performance. As shown in table 1, this is particularly true for GPUs. While the theoretical maximum computing performance of a modern CPU decreases by a factor of two, the peak performance of a consumer GPU like the Nvidia GeForce GTX 1080 Ti (henceforth called GTX) drops by two orders of magnitude. This is a significant problem for GPU-accelerated scientific software, where SP is often not accurate enough. To ameliorate this issue, graphics card manufacturers introduced new hardware specifically designed to improve the DP performance. The P100 and its recently released successor, the Nvidia Tesla V100 (henceforth called V100) achieve half of their SP performance when using DP. However, these special purpose GPUs are much more expensive than regular consumer GPUs like the GTX, typically by an order of magnitude. Table 1 shows that the SP performance of a high-end consumer GPU is comparable to the SP power of the special purpose GPUs. Thus, if it is possible to use SP instead of DP in our lensing algorithm, we would not only significantly increase the code performance on both CPUs and GPUs, but we might also be able to achieve a close to optimal performance with relatively cheap hardware.

3.2 Computing finite precision errors for strong lensing

We will now show that using SP in algorithm 1 is accurate enough for a large fraction of the image pixels. We restrict ourselves again to the SIE model. It is possible to generalize these results to other parametric models, but the fraction of the image for which SP is accurate enough will vary and has to be computed for each model independently.

The lens equation 13 maps the triangular pixels in the image plane onto triangular pixels in the source plane. We assume a HFF pixel size of 0.03 arcsec and we maximize the lensing effect by using $D_{\text{LS}}/D_{\text{OS}}=1$. As a result, the lens equation is now a simple subtraction of $\nabla\Psi_{\epsilon}(\theta_{1},\theta_{2})$. We now look at an observed multiple image in the image plane. Note that our ability to locate the multiple image is observationally constrained by the size of the image pixels, so there is an observational error budget on the image location of half a pixel. In addition, the algorithm lenses both the triangular pixel and the image into the source plane. It is possible that their respective errors due to machine precision have the same magnitude but the opposite sign, and therefore the error budget shrinks by another factor of two. As a result, the value of $\nabla\Psi_{\epsilon}$ can be considered accurate enough if the error $E$ is smaller than a quarter of a pixel. Thus our upper limit for the gradient error is $E_{i}\leq 7.5\times 10^{-3}~{}\text{arcsec}$, where $i=1,2$.

However, this error budget does not yet account for the magnification effect of strong lensing. Background sources and distance scales appear magnified when they are strongly lensed and consequently distance scales in the image plane like pixel sizes will be be de-magnified when they are mapped into the source plane. The resulting error budget for $\nabla\Psi_{\epsilon}$ becomes thus $E_{i}\leq 7.5\times 10^{-3}~{}\text{arcsec}/M_{i}$, where $M_{i}$ is the magnification along the $\theta_{i}$-axis.

In A we derive an upper bound for the error of $\nabla\Psi_{\epsilon}$ due to finite machine precision. We assume that the lens is a strong lensing cluster modeled with two cluster-scale SIE halos. The SP upper error bound along the $\theta_{i}$-axis is $\Delta(\nabla\Psi_{\epsilon})\leq 2.3\times 10^{-3}~{}\text{arcsec}$ if we use the following approach. As discussed in A, we compute the gradients with SP except in pixel grids of $400\times 400~{}\text{pixels}$ around cluster-scale halos and $20\times 20~{}\text{pixels}$ around galaxy-scale halos, where we use DP. This corresponds to approximately 1% of all image pixels. As a result, SP is accurate enough for each of the remaining 99% of the image pixels if the respective magnification along both $\theta_{i}$-axes is $M_{i}\leq 3.26$. In strong lensing, we typically measure the magnification of the area of a multiple image and not the magnification along an axis. The measured values are typically single digits [see e.g. Jauzac et al., 2015, for magnification values of a HFF cluster]. While these values cannot easily be converted to axis-magnifications due to the the typically arc-like shape of strongly magnified images, they strongly suggest that SP will be accurate enough for a large fraction of the image. However, strong lensing clusters have critical lines where the magnification diverges. In the case of the SIS, this critical line is the Einstein ring. While the magnification does not become infinite in practice [see e.g. Bartelmann and Schneider, 2001, for a detailed discussion], it can become very large and thus SP will no longer be accurate enough. Consequently, we can use SP for a large fraction of the image, but we also need to implement a mechanism which ensures that we compute the gradients with DP whenever SP is not enough due to high magnification.

3.3 Fixing the missing accuracy close to critical lines

We add the missing accuracy close to critical lines as follows. First, we compute $\nabla\Psi_{\epsilon}$ for each pixel in the image using the approach presented in the previous subsection. Second, we compute for each pixel

[TABLE]

which is computationally cheap because we have already computed the gradient values for all pixels. For the HST ACS, we have a pixel height and width $\Delta x=0.03~{}\text{arcsec}$. Note that $\delta_{i}$ corresponds to the change of the pixel length along the $\theta_{i}$-axis due to lensing into the source plane. Third, we recompute $\nabla\Psi_{\epsilon}$ in DP for all pixels where

[TABLE]

for $i=1$ or $i=2$, which implies that $M_{i}>3.26$. We derive this condition by computing the pixel length in the source plane $\Delta x_{\text{source}}$ along the $\beta_{1}$-axis,

[TABLE]

where we used the lens equation. An analogous relation holds for the $\beta_{2}$-axis. Note that taking the absolute value of $\Delta x_{\text{source}}$ is necessary because lensing can change the image parity [see e.g. the review Kneib and Natarajan, 2011]. For the assumed HFF pixel scale, the condition that $M_{i}>3.26$ translates into $\Delta x_{\text{source},i}<0.0092~{}\text{arcsec}$. The lensing effect along the $\theta_{1}$- and $\theta_{2}$-axis is shown in figure 4 and the criterion in equation 25 is illustrated in the top part of figure 5.

We can assume that each source is lensed along a chosen $\theta_{i}$-axis, as this can be achieved by a simple change of the image plane coordinate system. However, the shape of the image plane pixels is not invariant under such a transformation, as figure 4 illustrates. Therefore we need two additional criteria. We define

[TABLE]

and we recompute $\nabla\Psi_{\epsilon}$ in DP if

[TABLE]

for $i=3$ or $i=4$. This case is illustrated in the bottom part of figure 5.

In summary, we compute $\nabla\Psi_{\epsilon}$ in SP everywhere except in small patches centered on the origin of each lens as described in A and for the pixels where the criteria defined in equations 25 and 28 hold. This is illustrated in figure 6.

4 Performance measurements

We implement both the GPU and the CPU version of the gradient computation twice, once in SP and once in DP. The respective versions are identical up to the change in precision. In addition, we implement the mixed precision algorithm for both types of hardware. In the first step, this algorithm computes the result for each pixel in SP. In the second step, it checks which results are not accurate enough and recomputes these with DP. For this purpose, the algorithm uses the criterions developed in the previous section. The GPU implementation of the mixed precision algorithm uses asynchronous computations and load balancing for the second step, i.e. the computation of the most expensive DP calculations can be dispatched asynchronously between the CPU and the GPU. We adjust the load balancing for the different GPU models. As table 1 shows, this is particularly interesting for the GTX: The use of a hybrid CPU/GPU approach alleviates the very low DP performance of this card and drastically reduces the impact of the DP computations on the overall run time.

In the next step, we want to measure the performance gain of using HPC methods in strong lensing. For this purpose, we measure the time which the different software implementations require to compute the gradient of a HFF-like cluster lens for each pixel of a Hubble image. We have repeated this measurement several times and find that the benchmark results are stable, i.e. they do not significantly vary in different runs with the same setup. We also compute the gradients with the current LENSTOOL software, which serves as a reference. We assume a $\Lambda$CDM cosmology with $H_{0}=70~{}\text{km/(s~{}Mpc)}$, $\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$. We use an image with $6730\times 6730~{}\text{pixels}$ and a pixel scale of 0.03 arcsec/pixel to simulate images from the HST ACS. The galaxy cluster consists of two cluster-scale and 700 galaxy-scale halos like in the HFF cluster Abell 2744 [e.g., Jauzac et al., 2015]. We model the lens using SIE halos. The lens redshift is 0.3 and all sources are at the same redshift $z_{\text{source}}=2.0$. The velocity dispersion determined by Jauzac et al. [2015] for one cluster-scale halo of Abell 2744 is approximately 1200 km/s and for a galaxy halo it is roughly 150 km/s. However, these values correspond to the velocity dispersion parameter of the parametric lens model chosen in Jauzac et al. [2015], which is not identical to measured line-of-sight velocity dispersions of galaxies. The exact conversion must be computed numerically, but for our purposes a rough agreement is enough, so we can use a conversion factor of 0.85 [Elíasdóttir et al., 2007]. This leads to $\sigma_{v}\approx 1000~{}\text{km/s}$ and we choose this value and a pseudo-ellipticity of the potential $\epsilon=0.05$ for the first cluster scale halo. For the second large-scale halo, we use $\sigma_{v}=700~{}\text{km/s}$ and $\epsilon=0.04$. We model the galaxy halos by allowing $\sigma_{v}$ and $\epsilon$ to randomly vary between 10 and $15~{}\text{km/s}$ and 0 and 0.15, respectively. Note that the galaxy halo velocity dispersions are approximately a factor ten smaller than the ones used in Jauzac et al. [2015], because the magnitude of the scaled deflection angle for a SIE does not decrease with distance from the lens center, as it does for more realistic parametric lens models. Thus we need to decrease the velocity dispersion to limit the lensing effect of individual galaxies at large separations from the galaxy.

Figure 7 presents the gradient values in $\theta_{1}$ direction for the cluster lens. Figure 8 shows that the error resulting from our mixed precision algorithm is within the allowed limit for each pixel. Figure 9 illustrates the hypothetical error of a pure SP algorithm close to a cluster-scale halo. We see that the area in which we re-compute $\nabla\Psi_{\epsilon}$ shown in figure 6 covers nicely the area in which the SP error is largest.

Table 2 and figure 10 present the benchmark performance of the different gradient computation implementations. The speedup of the HPC-optimized codes with respect to the current LENSTOOL software is considerable. The indicated LENSTOOL performance is obtained by using all 12 CPU cores. It is thus the best currently achievable speed, as LENSTOOL cannot be run on multiple computer nodes and its performance is thus limited by the number of CPUs available on a single node. In addition, the mixed precision implementations are consistently faster than the DP ones, which validates our approach. The HPC CPU version reduces the run time of the benchmark by one order of magnitude and the GPU implementations by up to two while keeping the error within the allowed bounds. The consumer-grade GTX card displays the largest performance gain with respect to the DP computation. Figure 11 demonstrates that our HPC-optimized CPU software scales almost perfectly with the number of cores available on a single node using multi-threading. Figure 12 compares the benchmark performance of the Nvidia GTX with the high-end GPU based on the same Pascal GPU architecture (Nvidia P100). The P100 is an order of magnitude more expensive than the GTX. The P100 is considerably faster when only DP is used, which is expected as it was designed for this purpose. However, as soon as the mixed precision algorithm is used, the GTX reduces its run time dramatically and the performance comes close to the P100.

5 Discussion

The performance measurements in the last section demonstrate clearly the value of HPC methods for strong lensing software. They lead to a speedup of one to two orders of magnitude, depending on the chosen hardware. In addition, our measurements show that GPUs are perfectly suited for the massively parallel lensing calculations. As expected, the high-end GPUs have a big performance advantage in DP computations, but our mixed precision algorithm and the hybrid CPU/GPU approach bring the consumer GPU’s performance very close to its more expensive siblings. Note that the use of mixed precision also benefits the benchmark performance of the the high-end GPUs and of the CPU, but not on the same scale. Mixed precision thus leads to a performance benefit regardless of used hardware while also delivering accurate results.

Furthermore, table 3 demonstrates that the use of HPC methods dramatically reduces the energy consumption. We estimate the required energy to solution of the respective gradient computation implementations by multiplying the Thermal Design Power (TDP) of the used hardware with the time to solution. Note that we use only the TDP of the CPU and the GPUs for the energy to solution computations and we neglect the power consumption of other components which is typically considerably lower [see e.g. Cumming et al., 2014, for a detailed energy-efficiency study of a computing cluster]. Energy savings of up to 98% are possible compared to the current LENSTOOL software. The HPC techniques are thus friendly to the environment and lower the electricity bill of the computing cluster. While the use of the mixed precision algorithm further reduces the energy consumption in the case of the CPU and the GTX, its use increases the required energy for the P100 and the V100. This is due to the hybrid CPU/GPU approach in the mixed precision implementation for the GPUs. In the case of the GTX, the decrease in run time can handily offset the additional power consumption of the CPU, while this is not the case for the high-end GPUs.

It is possible to generalize this approach to other commonly used parametric lens models like the Navarro-Frenk-White (NFW) profile [Navarro et al., 1997, 1996] or the dual Pseudo Isothermal Elliptical mass distribution (dPIE) [e.g., Elíasdóttir et al., 2007, Kassiola and Kovner, 1993]. However, this will lead to much more complicated and longer gradient computation algorithms than the one for the SIE studied in this paper. As a result, it is possible that the resulting error due to finite machine precision will become bigger, too. While the exact fraction of the image for which SP is accurate enough must be explicitly calculated and studied for the respective model, it is possible that it will be considerably lower than for the SIE. Therefore the performance gain from using mixed precision might shrink accordingly. A separate study will be necessary to determine whether the expected performance gain is worth the effort.

Finally, it is necessary to discuss legal aspects of the GPU drivers provided by the graphics chip manufacturer, Nvidia. The drivers are required to use the hardware and their use is subject to certain terms and conditions, which Nvidia has recently updated. These new terms might be interpreted to prohibit the use of consumer-grade graphics cards like the GTX in computing clusters, which is the typical deployment in the scientific community. Thus the customers would be effectively forced to buy the much more expensive high-end GPUs, even though the cheaper graphics cards might be fully sufficient for the intended application. The authors of this paper are strongly concerned about this development, in particular given the limited financial budgets of academic research worldwide. Therefore they have contacted Nvidia and were informed in writing that Nvidia has no intentions to prohibit the use of the cheaper consumer-grade cards for the non-commercial purposes of researchers. The authors urge Nvidia to formalize this permission for scientific use by including it in the terms and conditions or preferably to remove this restriction altogether, thus allowing everyone to use the GPU which best suits their respective needs.

6 Conclusion

In this paper, we demonstrate the value of High Performance Computing techniques for strong lensing software. We study a performance-critical part of the widely used LENSTOOL lens modeling software, namely the deflection potential gradient computations of the $\chi^{2}$ calculation algorithm. We present and discuss an optimized CPU version with Advanced Vector Extensions and OpenMP and a GPU implementation in CUDA for a SIE lens model.

In addition, we calculate the impact of finite machine precision on the strong lensing algorithm. We demonstrate for the SIE model that single precision is accurate enough for a large part of the image. We develop a mixed precision algorithm which allows us to use single precision for performance while computing critical parts of the image in double precision.

Finally, we measure the computing performance for a galaxy cluster lens similar to the Hubble Frontier Fields. We find that our HPC techniques accelerate the computation by an order of magnitude on CPUs and by up to two orders of magnitude on GPUs. In addition, they reduce the energy consumption by up to 98%. The mixed precision approach delivers the best performance for every type of hardware while providing accurate results. It also permits to harness the full potential of a consumer-grade GPU, which can achieve a competitive benchmark performance for a small fraction of the monetary cost of a high-end GPU.

7 Acknowledgments

MR thanks Yves Revaz for fruitful discussions of GPU-accelerated computing. GF gratefully acknowledges support from the EPFL Faculté des Sciences de Base. This work was supported by EPFL through the use of the facilities of its Scientific IT and Application Support Center. The authors gratefully acknowledge the use of facilities of the Swiss National Supercomputing Centre (CSCS) and they thank Colin McMurtrie and Hussein Harake for their continued support. This research made use of matplotlib [Hunter, 2007], Inkscape, Astropy [Astropy Collaboration et al., 2013], TeX Live, Wolfram Alpha, and NASA’s Astrophysics Data System.