im3shape: A maximum-likelihood galaxy shear measurement code for cosmic gravitational lensing
We present and describe im3shape, a new publicly available galaxy shape measurement code for weak gravitational lensing shear. im3shape performs a maximum likelihood fit of a bulge-plus-disc galaxy model to noisy images, incorporating an applied point spread function. We detail challenges faced and choices made in its design and implementation, and then discuss various limitations that affect this and other maximum likelihood methods. We assess the bias arising from fitting an incorrect galaxy model using simple noise-free images and find that it should not be a concern for current cosmic shear surveys. We test im3shape on the GREAT08 Challenge image simulations, and meet the requirements for upcoming cosmic shear surveys in the case that the simulations are encompassed by the fitted model, using a simple correction for image noise bias. For the fiducial branch of GREAT08 we obtain a negligible additive shear bias and sub-two percent level multiplicative bias, which is suitable for analysis of current surveys. We fall short of the sub-percent level requirement for upcoming surveys, which we attribute to a combination of noise bias and the mis-match between our galaxy model and the model used in the GREAT08 simulations. We meet the requirements for current surveys across all branches of GREAT08, except those with small or high noise galaxies, which we would cut from our analysis. Using the GREAT08 metric we we obtain a score of Q=717 for the usable branches, relative to the goal of Q=1000 for future experiments. The code is freely available from https://bitbucket.org/joezuntz/im3shape.
keywords:methods: statistical, methods: data analysis, techniques: image processing, cosmology: observations, gravitational lensing: weak, dark energy
Weak gravitational lensing has the most potential for constraining the cosmological parameters that describe the dark universe (Albrecht et al. 2006; Peacock & Schneider 2006). As dark energy slows the growth of dark matter structures, these in turn affect the appearance of distant galaxies, distorting them slightly as gravitational fields deflect their light. This distortion is known as cosmic shear. Its measurement from individual galaxy images is crucial to the field of weak gravitational lensing, as the amount of distortion is directly related to the amount of intervening dark matter. However, the stretch induced by cosmic shear is typically on the order of a few percent, which renders accurate shape measurement an intricate but also intriguing problem.
Several upcoming and future observational surveys plan to capitalise on the potential of cosmic shear for discovering the nature of dark energy. These include imminent ground-based projects (the Kilo-Degree Survey, Hyper Suprime-Cam111http://www.naoj.org/Projects/HSC/HSCProject.html and the Dark Energy Survey: DES222http://www.darkenergysurvey.org), ground-based telescopes under construction (the Large Synoptic Survey Telescope: LSST333http://www.lsst.org), and future space telescopes (Euclid444http://sci.esa.int/euclid and WFIRST555http://wfirst.gsfc.nasa.gov).
It is perhaps surprising that a major challenge for these experiments is simply measuring the shapes of the galaxy images. The difficulties arise because the galaxies involved are: small, less than the Point Spread Function (PSF) size in many cases; noisy, with a signal-to-noise ratio (S/N) of ; often irregular; and convolved with a PSF which itself induces ellipticity in the galaxy that is typically greater than the gravitational shear we wish to measure. To fulfill the potential of cosmic shear these problems must be overcome, allowing precise and accurate measurements of the underlying gravitational shear.
1.1 Existing methods
The best known and most commonly-used method in weak lensing shape estimation remains that described by Kaiser et al. (1995), Luppino & Kaiser (1997) and Hoekstra et al. (1998), known as KSB or KSB+ (see also Schneider, 2006 for an overview). KSB uses weighted quadrupole moments on galaxy images and PSF images to provide an estimator of gravitational shear. Many implementations of this approach, and derived methods adapted for space-based data (e.g. Rhodes et al., 2000, Schrabback et al., 2010), have been used to estimate gravitational shear over more than a decade. An important discovery of the Shear TEsting Programme (STEP) competitions was the sensitive dependence of the accuracy of KSB upon the details of its implementation (see Section 1.2, below; Heymans et al., 2006; Massey et al., 2007).
Recent modifications to KSB have sought to improve accuracy by adding higher-order terms to series approximations in the method (e.g. Okura & Futamase, 2009; Viola et al., 2011; Melchior et al., 2011; Okura & Futamase, 2012). More fundamentally different shape measurement approaches have also been proposed. A relatively early example was the ‘commutator’ method of Kaiser (2000); although this method has not found broad adoption in practice, the work contained a number of insights that were to inform the later development of the field.
Gaussianization, or Re-Gaussianization, approaches (e.g., REGLENS: see Mandelbaum et al., 2012, Section 4.3; Mandelbaum et al., 2005; Hirata & Seljak, 2003) first convolve galaxy images with an additional kernel designed to make the final PSF as Gaussian, and isotropic, as possible. Such images then closely match many of the implicit assumptions in weighted quadrupole moment approaches such as KSB, so moments can then be used to extract shear estimates with reduced bias. Gaussianization methods remain competitive in terms of systematic bias (Mandelbaum et al., 2005; Mandelbaum et al., 2012) and have been used successfully in a number of studies in the Sloan Digital Sky Survey (e.g. Huff et al., 2011; Reyes et al., 2010; Mandelbaum et al., 2006; Hirata et al., 2004).
A family of methods known commonly as shapelets construct a model of PSF-convolved galaxy images as a sum of orthonormal Gauss-Lagurre polynomial functions (e.g. Bernstein & Jarvis, 2002; Refregier & Bacon, 2003; Massey & Refregier, 2005; Kuijken, 2006; Massey et al., 2007; Nakajima & Bernstein, 2007). Effectively the ellipticity of derived shapelet models can then be used to provide an estimator of gravitational shear. Melchior et al. (2010) raised some conceptual concerns with the approach, and while similar (perhaps more drastic) conceptual problems with KSB did not affect its popularity there are fewer examples of shapelet shear estimation from real data in the literature (although see Velander et al., 2011).
The number of proposed shape measurement methods is steadily growing. Interest in the problem has undoubtedly been stimulated by impending survey data of unprecedented statistical power, and by the blind analysis contests we describe in Section 1.2, below. Examples of innovative recent proposals the competitions described in the next section, and include Fourier Domain Null Testing (Bernstein, 2010) and shear estimation via lookup table (e.g., ‘MegaLUT’; see the review in the appendices of Kitching et al., 2012).
The final family of shape measurement methods are commonly known as model-fitting methods, and involve the comparison of simple parametric galaxy models (convolved by the PSF) to imaging data. The galaxy model ellipticity, among other parameters inferred in the model fit, make up the shear estimator.
The lensfit (Miller et al., 2007; Kitching et al., 2008; Miller et al., 2012) method seeks to perform Bayesian inference on galaxy model parameters and performed well in the GRavitational lEnsing Accuracy Teasting (GREAT08) Challenge: see Section 1.2, below; (Bridle et al., 2010). The more recent DeepZot method that also used a model-fitting algorithm proved the eventual winner of GREAT10 (Kitching et al., 2012).
The im2shape public code (Bridle et al., 2002)666http://www.sarahbridle.net/im2shape/ modelled both the galaxy and the PSF as a sum of Gaussians, inspired by Kuijken (1999). It performed MCMC sampling and took the mean of the of the ellipticity values of the samples as the shear estimate. It was used by Bardeau et al. (2005); Bardeau et al. (2007) to measure galaxy cluster masses.
1.2 Shear measurement contests
Researchers have come together to assess these many shear measurement methods in different regimes. The Shear TEsting Programme (STEP) was a joint blind analysis of simulated data by groups from within the shear measurement community (Heymans et al., 2006; Massey et al., 2007). It showed that the shear measurement problem is far from trivial, but demonstrated that the cosmological constraints set at that time were not limited by shear measurement accuracy. The first two GRavitational lEnsing Accuracy Testing (GREAT) Challenges posed the problem to computer scientists and provided sufficient simulations to assess whether methods were suitable for future surveys (Bridle et al., 2008; Kitching, 2010). They showed that there are some regimes where existing methods work well enough, but others where further progress is necessary (Bridle et al., 2010; Kitching et al., 2012). The third challenge in the series, GREAT3,777http://great3challenge.info is currently in progress towards a challenge opening and simulation data release in 2013. In addition to a large simulation dataset, GREAT3 is leading the collaborative development of an open source extragalactic image simulation software toolkit,888https://github.com/GalSim-developers/GalSim further advancing the STEP and GREAT programme of providing the community with a common reference for comparison and improvement of galaxy shape measurement methods.
This paper presents a successor code to the im2shape code described above, called im3shape. It uses (by default) a sum of two Sersic profiles for the galaxy, similar to lensfit (Miller et al., 2007; Kitching et al., 2008; Miller et al., 2012), and a Moffat profile for the PSF. It calculates the maximum likelihood parameter values and corrects for noise bias (see, e.g., Refregier et al., 2012; Melchior & Viola, 2012; Hirata et al., 2004; Bernstein & Jarvis, 2002) using the calibration scheme described in Kacprzak et al. (2012).
In Section 2 we discuss the requirements for shape measurement for current and future data sets. Then in Section 3 we present the im3shape code that we are releasing, the details of its operation and the choice of parameters it uses. Section 4 is concerned with noise-free simulations, to verify the code in ideal circumstances and to explore what biases the use of an incorrect model can introduce to the results. We present the results of tests of the code on the GREAT08 simulation set in Section 5, and conclude in Section 6. We give more details of the code implementation in Appendix A, of the termination criteria in Appendix B and overview the analysis process in Appendix C. Appendix D gives the equations and details used for analytic calculation of likelihood gradients with respect to the parameters.
Following Heymans et al. (2006), the accuracy of shear measurement methods are generally quantified in the literature in terms of the multiplicative error and additive error on the true shear , such that
where is the shear estimate and the subscript refers to the two shear components. It is assumed there is no cross contamination between, for example, and .
The requirements on and for cosmic shear are computed in Amara & Réfrégier (2008) as a function of the survey area, depth and galaxy number density (see also Kitching et al., 2009). Amara & Réfrégier (2008) consider a two-point statistical analysis of the shear field, and the requirements are set so that the systematic error equals the statistical error, thus providing an upper limit on the allowed bias.
The GREAT08 (Bridle et al., 2008) challenge quantifies the accuracy of shear measurement methods using a single number called the quality factor, . The equation relating to and is given in Voigt & Bridle (2010). A related but different definition is used in GREAT10 (Kitching et al., 2012) but we stick to the GREAT08 definition here because we use the GREAT08 simulations.
In Table 1 we show the requirements on , and for three sets of survey parameters, chosen to represent current, upcoming and future cosmic shear surveys. Examples of surveys with requirements comparable to the ‘current’ set of requirements are the Canada-France-Hawaii Telescope Lensing Survey (CFHTLenS: Heymans et al., 2012), and likely the first results from KIDS, DES and HSC (de Jong et al., 2012; The Dark Energy Survey Collaboration, 2005; Takada, 2010). The upcoming survey requirements correspond roughly to the full DES and HSC surveys. The future survey requirements correspond to LSST and Euclid (e.g. Chang et al., 2012; Laureijs et al., 2011).
3 Im3shape methodology & code
The code we have developed and present in the remainder of this paper is called im3shape. It belongs to the family of model-fitting methods: it forward-fits a galaxy model to each data image it is supplied with and reports the parameters of the best fitting model, including the ellipticity components. For finding the best fit we use a maximum likelihood approach, i.e. we search for those model parameters that minimize , the summed square difference between our generated model and the actual pixel data, weighted according to a user-supplied weight image. For details on the maxmimum likelihood estimation see Appendix A and D.
Throughout we assume that we have selected a square portion of the image, with stamp_size pixels on a side.
3.1 Code layout
im3shape is a modular C code, with a significant amount of Python glue code that enables the setting up of new models and their options automatically. Key processes like optimization and galaxy model generation routines are encapsulated in complex structures with simple interfaces, rather than complicated functions. This architecture makes it particularly easy to design new galaxy or star models to be fit to data. Since the process of fitting to a large number of postage stamps is inherently trivial to parallelize, the core of the code is all single-threaded.
3.2 Model choice
Using the flexible systems of im3shape we have experimented with fitting a variety of models, mainly to the images in the GREAT08 challenge.
We chose by default a two-component Sérsic, with a de Vaucouleurs bulge and an exponential disc, with the same centroid and ellipticities. In addition we simplify the fitting by fixing the ratio of scale radii, at unity.999This is similar but not identical to the model adopted in the GREAT08 simulations. im3shape should therefore work reasonably well on those images, and we note that this model was originally chosen to be a reasonable description of real galaxies.
Within that general model we can choose different parameterizations for some components - these are discussed in detail in Appendix A.
Our complete set of parameters is shown in Table 2. The parameters marked as fixed were kept at default values during our standard analysis but can be varied easily if desired or to study model bias.
|Disc half-light radius|
|Bulge peak flux|
|Disc peak flux|
|Bulge-disc size ratio||✓|
|Disc Sérsic index||✓|
|Bulge Sérsic index||✓|
3.3 Model fitting
Generating the model image described in the previous Section is not trivial: we need to ensure that the model is at high enough resolution to capture sharp central galaxy features. This is discussed in Section 4.1.
Assuming images of the galaxy model can be accurately generated (see Section 4.1), maximizing the computed likelihood (i.e. minimizing the weighted, squared difference between the non-linear model and the data) is technically and computationally challenging. Reaching the current level of speed and accuracy in maximization within im3shape required considerable exploration of maximization algorithms and their settings. We therefore give a brief overview of our findings in Appendix B, and an overview of the whole process of image generation and fitting in Appendix C.
With tuned parameter settings the maximizer, which uses a Levenberg-Marquardt algorithm, works well. Figure 1 shows example likelihoods from the first galaxy of set 3 in GREAT08. It can be seen that the point is indeed a good minimum (and further searches can confirm it to be a global one). The asymmetry in some of the parameters is part of the cause of the noise bias phenomenon (Kacprzak et al., 2012).
4 Noise-free simulations
Generating the model image described in the previous Section is not trivial; for example, we need to ensure that the model is at high enough resolution to capture sharp central galaxy features without aliasing. We also need to understand how best-fitting models may be biased in cases where they cannot represent the data perfectly. In this Section we investigate the behaviour and performance of im3shape on noise-free simulations to find optimal parameter settings for model generation, and analyse in a simple way what potential biases can arise when the model is insufficiently flexible to describe the data.
4.1 Upsampling biases
The im3shape software uses the Discrete Fourier Transform (DFT) to render images of convolved galaxy profiles: this can be done with speed thanks to the Fast Fourier Transform (FFT) algorithm, freely available in highly efficient implementations (Frigo & Johnson, 2005). However, representing the convolution of telescope PSFs with analytic, continuous, non-bandlimited functions such as the Sérsic profile requires careful numerical approximation, as can (if the PSF is also not bandlimited) the subsequent rendering of convolved profiles into pixellated images (see, e.g., Marks, 2009). Decisions must be made about how and where to make approximations in the representation of these profiles, balancing precision against computational cost.
Upsampling, which is the more accurate rendering of model images using higher resolution (upsampled) intermediate images for internal operations such as convolution (see, e.g., Rowe et al., 2012), is a key tool to eliminate defects in modelling due to aliasing. We therefore investigate the requirements on the upsampling parameters used by the image generation code in im3shape. These three resolution parameters are:
, called upsampling in parameter files - the convolution between the galaxy model and the PSF model is performed on a grid of pixels in size, i.e. each image pixel is divided into sub-pixels across. The flux within each sub-pixel is computed assuming the intensity is constant within it and equal to the intensity at the centre of the sub-pixel, for all pixels except the central pixels (see (2)).
, called n_pixels_to_upsample in parameter files - within a central galaxy region of image sub-pixels the flux is integrated on a even finer grid within each sub-pixel (instead of assuming a constant intensity across each sub-pixel).
, called in parameter files - determines the number of integration intervals within a sub-pixel in the central image region; therefore the overall resolution in the central galaxy region is times the data image resolution.
The required values of the these resolution parameters depend on the survey (see section 2). In figure 2 we plot the multiplicative bias as a function of the parameter for different values of and with . (The additive bias is negligible in comparison to the multiplicative bias because we use a circular PSF for this investigation.) Results are shown for a pure bulge and a pure disc with and 1.14, which are similar to the fiducial and small size branches used in GREAT08 (see Bridle et al., 2010, including for size definition). The bias is computed using a ring-test (Nakajima & Bernstein, 2007) with 4 angles, degrees.
The multiplicative bias is obtained by shearing each galaxy in the ring by an amount along the -axis and along the line . The PSF is a circular Moffat profile with FWMH=2.85 pixels and shape parameter . The postage stamp size (stamp_size in parameter files) is 25 pixels. We find that the requirements on the resolution parameters do not change if we increase the stamp size to 45 pixels. The bulge-to-disk scale size ratio, , is fixed to 1 during the fitting.
We find that for an upcoming survey (medium green shaded area in figure 2) we require 5, , and 7. For the results are unstable.
4.2 Model bias
Having established that the resolution parameters described above are high enough to remove any bias from aliasing effects, we now turn to the issue of model or external bias arising from fitting an incorrect model to the data. If a forward model approach could be perfectly implemented and the posterior distribution accurately characterized for each galaxy (rather than being reduced to a single point-estimate like the maximum likelihood) and combined appropriately, then this would be the only remaining source of bias for simple images like those in GREAT08. In practice, there may be other problems such as neighboring galaxies and image artefacts.
Very realistic galaxy models are almost certainly out of reach of forward methods for the forseeable future - modelling spiral arms, for example, is hard enough for a single high-resolution image, let alone noisy ones. We defer testing on extremely realistic models (for example using the SHERA simulation code of Mandelbaum et al., 2012, or Galsim), and instead apply a necessary but not sufficient test - fitting a simple but incorrect model using im3shape.
Again, we perform these tests without noise, to separate the issue of noise bias on maximum-likelihood estimates from external bias; of course, there may be an interaction between the two. As we only touch upon these issues here, we refer the interested reader to Voigt & Bridle (2010) and Melchior et al. (2010) for a more extensive disucssion.
Except where otherwise noted, our galaxy simulations in this section used , equal component scale radii, and half the flux in each of the bulge and disc components . The PSF is a Moffat with , , and FWHM.
4.2.1 Sérsic Indices
We first test a simple difference between the input (true) model and the fitted one - a different index of the Sérsic profile. Our input (true) model is generated with a single component, a Sérsic profile with index , however, we fit our standard bulge plus disc model.
Figure 3 demonstrates this bias: it shows the pre- and post-PSF convolved simulated galaxy images, and the residuals between the simulated galaxy image and best-fit galaxy model, for two different values of the true Sérsic index. The images are noise-free and the correct PSF model has been used in the fits. The error on the true galaxy ellipticity is approximately for , and for .
The residual image shows different behaviour along the and directions, which is what causes the ellipticity to be mis-estimated. As discussed in Voigt & Bridle (2010), the galaxy size is worst estimated along the direction where the galaxy is smallest relative to the PSF. Therefore the ellipticity is incorrect.
Our model used in im3shape is correct when or , since at these points the true model is a subset of the im3shape model. We generate simulations with our fiducial sets of galaxy parameters, except for the Sérsic index, which we vary. We apply a ring test to the simulations to assess the multiplicative and additive biases; the former is shown in figure 4, a figure for the latter is omitted as it is negligible in the noise-free case. As expected, the bias becomes zero when the true and fitted model coincide at and .
Figure 4 suggests that this kind of model bias has the potential to be a problem for upcoming surveys but is only marginally relevant for current generation of experiments. Allowing for additional flexibility of the model could resolve this potential issue. Further tests are needed on more realistic galaxy images to assess the true impact on upcoming surveys.
4.2.2 Different bulge and disc ellipticity
In figure 5 the multiplicative bias arising from a slightly more realistic model difference is shown - a bulge and disc model with components that are not co-elliptical. This is of interest because it is a difference between the im3shape model and the GREAT08 simulations. The input image bulge ellipticity101010We define ellipticity as with from Table 2 and Appendix A. is fixed at 0.3 and the disc ellipticity varied so that . The difference between the bulge and disc ellipticity, , is varied between -0.25 and 0.25. Typically, bulges are rounder than discs, i.e . The model fit to the data is our standard two-component co-elliptical profile.
Once again a ring test is applied to assess the bias, and the noise is zero. The other galaxy parameters are set to the fiducial values described above, except that the different curves show different relative amplitudes of the bulge and disc components. The bias is zero when either is small or the galaxy is completely bulge- or disc-dominated, since in that case the secondary component does not affect the fit at all. The effect is greatest when the components have comparable total fluxes, and this bias too is small for the current generation of experiments, if bulges are rounder than discs.
4.2.3 Incorrect Radius Ratio
As depicted in Table 2, the ratio between the half-light radii of the bulge and the disc is fixed by default in the code. This is partly to reduce the number of parameters (since fewer parameters generally means easier optimization), but also because this particular parameter causes problems with convergence in some regions of its range.
This is clearly a potential source of model bias - the ratio of these radii varies in real galaxies (when bulges and discs can be clearly identified). Furthermore, in GREAT08 this value is around . We again apply a ring test to noise-free simulations to check how severe this problem is. If we alter the ratio without changing the disc radius, then the bulge radius changes accordingly. To separate the effect of fitting small galaxies from this model bias we keep the ratio of the convolved image to PSF size, i.e. , fixed at three different values corresponding to GREAT08 sizes.
Figure 6 shows the results of this test. If real galaxies have smaller bulge scale radii than discs then the biases are mostly smaller than the requirements for upcoming surveys. If they are larger then the bias would be acceptable for current surveys but problematic for upcoming surveys, with a sign which depends on the galaxy size. For the GREAT08 fiducial branch (, true ) we expect a slight over-estimate of the shear. This will combine with the slight underestimate expected from figure 5 for bulges which are rounder than discs, and may partially cancel out to bring results down to be compatible with the requirements for upcoming surveys. We recommend a program of testing on realistic images for calibration of upcoming surveys.
We tested im3shape on the Gravitational LEnsing Accuracy Testing 2008 (GREAT08) Challenge (Bridle et al., 2008, 2010) which has a sufficient number of galaxies to test methods at the accuracy required for future surveys (unlike STEP1 and STEP2). Furthermore, it features a constant shear across the images which greatly facilitates the computation of success metrics (unlike GREAT10) and therfore performance evaluation. The GREAT3 Challenge is still under construction (Mandelbaum, Rowe et al. in prep).
The GREAT08 Challenge includes two blind competitions: one with low noise, which has a signal-to-noise ratio111111Here, we adopted the GREAT08 definition of signal-to-noise ratio defined by (A11) in Bridle et al. (2010), Appendix A, page 14. of 200 (‘LowNoise_Blind’) and one with more realistic signal-to-noise ratio around 20 (‘RealNoise_Blind’). In this section we show the results of im3shape on these simulations.
5.1 Low Noise Blind
The LowNoise_Blind challenge consists of 15 image files, each of which has 10,000 galaxy postage stamps of 3939 pixels. All the galaxies are a mixture of bulge and disk components which are co-centered but not co-aligned. The misalignment angle is drawn from a Gaussian with a standard deviation of 20 degrees. We therefore expect a model bias to apply in these cases as described in section 4.2. The signal-to-noise ratio of 200 is large enough that noise bias described in Kacprzak et al. (2012) is expected to be negligible. The PSF and galaxy size parameters are described in detail in Bridle et al. (2010). In summary, there are 5 image files of each of 3 different galaxy sizes. This allows us to show how the success metrics depend on galaxy size.
The GREAT08 images are simulated at essentially infinite resolution by using photon shooting for the PSF and galaxy profiles (Bridle et al., 2010). im3shape makes its models at high resolution (see section 4.1). For this test we ran at two different values of this resolution, with an upsampling for a coarse run and for an improved, higher resolution run that took about twice as long to complete. We are interested to see both results on the small LowNoise_Blind dataset and would like to use the coarser resolution for the much larger RealNoise_Blind dataset.
In figure 7 we show the im3shape results on the LowNoise_Blind GREAT08 data. At high resolution the resulting -factor scores reach level for all branches, a level at which the statistical noise in the challenge dominates - scores above this level are consistent with zero error. We reach this score in spite of a model bias: the the GREAT08 simulated galaxies have misaligned isophotes, whereas we fit co-elliptical ones, and the ratio of bulge to disc size ratio is different.
At the coarser resolution we are within the requirements for future surveys for all but the multiplicative bias on small galaxies, which is nonetheless acceptable for the current generation. We therefore use the coarser resolution in the remainder of this paper.
We also show the results submitted by others at the time of GREAT08. Significant developments have been made in other codes in the meantime, but we can conclude that at that time there were no methods as effective as im3shape in this high signal-to-noise regime. The closest competitor is the stacking method by AL, which has not yet been developed into a method that can be used on realistic data. We should however note that the model in the im3shape code is very well matched to the GREAT08 simulations, in that they both use exponential disc plus de Vaucouleurs bulge galaxies. An investigation into this topic is beyond the scope of this work, and will be addressed by future challenges including GREAT3.
In running this analysis we determined which parameters in the code can be tweaked to minimize runtime while retaining high accuracy. The most important such parameters are those controlling the details of the LevMar optimization and convergence behavior. In subsequent work we have found that tuning these parameters is quite specific to the noise levels and other parameters of the galaxies in question, and we have not been able to draw any general conclusions about them.
5.2 Real Noise Blind
The GREAT08 RealNoise_Blind Challenge is considerably larger, with 2700 sets of galaxy stamps. The 2700 sets are divided into 9 branches which span a range of simple galaxy models and PSF ellipticities and a range of realistic galaxy sizes and noise levels. The fiducial branch has a Signal-to-Noise a factor ten lower than in LowNoise_Blind but the same galaxy size as the central branch of LowNoise_Blind and the same PSF and galaxy type (co-centred but not co-elliptical bulge plus disc). Two branches explore the same two extrema in galaxy size as LowNoise_Blind; two branches vary the Signal-to-Noise by a factor either side of the fiducial value; two branches vary the PSF ellipticity by a factor of two and in direction; and two branches explore variants on the galaxy model: bulge or disc and bulge plus disc with an offset in the center. For more details on the challenge, see Bridle et al. (2008).
We tuned the maximization parameters described in Appendix B.1 to improve the mean likelihood of all the galaxies in a single set; on the basis that if a set of minimisation parameters did not find the best likelihood then it must not have converged fully. We note that the mean ellipticity obtained for this single set was quite sensitive to the maximization parameters, and therefore we recommend users of the im3shape code should perform a similar tuning on new galaxy samples.
5.3 Noise Bias Calibration
At the noise levels in the RealNoise_Blind Challenge we do expect that we need to correct for the noise bias effect described in Refregier et al. (2012), Kacprzak et al. (2012) and references therein. We describe our approach to noise bias calibration in this subsection.
In Kacprzak et al. (2012) we ran im3shape on noisy simulations with a range of input parameters to find the behaviour of the noise bias with galaxy properties. We then attempted to remove the bias by matching the measured galaxy properties to the simulation input parameters to predict the expected bias, which we subtracted off. We ran this method on branch 3 of our GREAT08 analysis, which has the same galaxy model as im3shape, and found , with and which is nearly satisfactory for current experiments but well outside our requirements for upcoming experiments.
This straightforward correction fails to do better because the observed properties are themselves noisy (and thus noise-biased), and so our bias estimate is itself biased. In practice, we have found it ineffective to fully calibrate galaxy ellipticity measurements one-by-one using their individual properties, and turn instead to a population-based approach in which a complete class of galaxies is calibrated collectively.
The approach that we take is discussed in Appendix E. Briefly, we use the deep data from the low-noise part of the challenge to provide distributions of the galaxy properties, and then perform simulations to compute the mean multiplicative and additive biases for those properties. We then apply these mean biases to the shear estimates from each set.
To predict biases for a given parameter set we simulate galaxies at grid points with scale radii and bulge flux fractions , and with ellipticity from to . At each grid point we fit a cubic polynomial to the bias as a function of , and store the coefficients. To get the bias for intermediate parameters we use a Gaussian radial basis function interpolator. Since the simulations are costly we simulate only at , and extrapolate to other values using the result from Refregier et al. (2012) and Kacprzak et al. (2012) that the bias scales as the inverse of the square of the signal-to-noise.
There are several factors about the structure of GREAT08 that are unrealistic in this context. Firstly, not every branch in the main sample has a corresponding low-noise version, whereas in real surveys it should be possible to match populations with more care in most cases. This would tend to make noise calibration harder on GREAT08. Secondly, the sizes of galaxies in GREAT08 are not drawn from a population - they have a single true value. This would tend to make calibration more effective than is realistic. We also note that the information needed to fully perform the calibration in this way was not public at the time of the challenge.
Where possible we match branches using the corresponding low-noise branch; where none exists we use the nearest approximation. See Appendix E for more details.
Q-factor scores and and values for the pre- and post-noise bias calibration results are shown in figure 8. Even before the calibration the bias on the galaxies was small enough that we achieved a quality factor , with and small enough for upcoming surveys. At lower signal to noise, as we expect, the pre-calibration Q factor drops to below 100.
After the noise bias calibration procedure described above the Q factors and the values reach the levels required for upcoming experiments for most of the branches of the challenge. The results are stable with PSF type. The m values are somewhat more variable but, depending on the exact branch, can reach the required levels. In particular the value is good in the b/d branch, where where we expect no model bias.
Our most encouraging result is the stability across the galaxy type branches. The only branch for which we use the correct model is the first galaxy type, where the simulation model is single-component bulge or disc. For all other branches we expect some model somewhat like the ones discussed in section 4.2, but as in that section the effect is not critical for current surveys.
The branches for small and noisy galaxies are clearly more problematic, and we would be forced to remove such galaxies from any current analysis. We believe the extrapolation fails badly in the high-noise case because of a slowdown in the magnitude of the bias with ellipticity, which does not follow the relation that we used. Very noisy galaxies have much broader scatter in measured ellipticity, and when these values start to push up to the edge of the space, at , this acts as an bias which partially counteracts the usual one. Missing this effect can lead to extrapolated values over-correcting the bias.
We have presented a public galaxy shape fitting code, im3shape, for weak lensing shear measurement. Our code efficiently finds the maximum likelihood fit of our model to any given galaxy, and does so extremely well: we find little sign of any bias due to poor minimization for simple simulated galaxies, following careful tuning of the model image resolution settings and other parameters. The code is modular and easily extensible, and we would welcome community contribution to any part of it121212http://bitbucket.org/joezuntz/im3shape. The technical problems in generating models and fitting them to galaxies are extensive, with particular problems of resolution, parameterization, and non-linear numerical optimization. We discuss our solutions to these problems in the appendices to this paper.
We find that we need to generate model images at five times the data resolution to reach our desired accuracy for upcoming surveys. We also find that to accurately capture small scale features in de Vaucouleurs cores we must sample a small central region even more finely - at times the data resolution.
Our code fits a simple bulge plus disc galaxy model to galaxy images, and we have made an exploration of model bias, which can arise when true galaxies fail to be as neatly analytic in form as our model ones. We find that reasonable deviations between the simple true and assumed models cause only minor bias, not significant for current data and early results from the upcoming generation of surveys. The bias will, however, attain greater significance as upcoming surveys reach their full statistical potential. Our tests were rather limited in that we considered only simple differences from our canonical model - deviations in Sersic index, and differences between the size and ellipticity of the two model components. We look forward to performing more realistic tests with the GalSim tool set in future work (Mandelbaum, Rowe et al. in prep).
We have demonstrated the utility of our method and code by testing on the GRaviational lEnsing Accuracy Testing 2008 (GREAT08) Challenge in both low and realistic noise regimes. We find extremely good results for the low noise simulations, with multiplicative and additive shear measurement biases well within the requirements for upcoming (Stage III) surveys, for the medium sized and large galaxies. For practical reasons we chose to compromise on the resolution we used in the fits and thus our results for the smallest galaxies are only sufficiently good for current surveys.
In Kacprzak et al. (2012) we calculated the degradation of shear measurements due to noise on images and discussed correction schemes. We do not expect this to be significant at the signal-to-noise of the LowNoise simulations (). We expect the main difficulty to arise from resolution and the differences between our fitted model and the GREAT08 galaxy models. Both our galaxy models and the GREAT08 low noise simulations assume galaxies are made of a de Vaucouleurs bulge plus an exponential disc. The difference is that we assume these two components are co-centric and co-elliptical with the same size ratio, whereas in GREAT08 low noise the two components have different ellipticities, are not co-aligned and have a bulge scale radius which is approximately 50% larger than the bulge scale radius. We find that there is little resulting problem arising from model bias alone, as expected from our simpler investigations.
We ran im3shape on the 27 million galaxies in the realistic noise part of the GREAT08 Challenge, which we chose over GREAT10 because it is simpler and the variable shear fields used in the latter were not relevant to this method.
At this signal-to-noise () we do expect significant noise bias due to pixel noise in the images (Kacprzak et al., 2012), leading to a multiplicative shear measurement error of a few per cent. We found that the correction scheme tested in (Kacprzak et al., 2012) does not work well on GREAT08 due to the bias-on-bias problem described in that work. We therefore deviated from the methods discussed there, which corrected each galaxy individually, in favour of a global correction using the knowledge of galaxy properties measured in the corresponding low noise sample. The leading order correction to shear power spectrum measurements is expected to simply depend on the mean multiplicative bias of a population, which can be corrected globally. Since the noise bias depends on galaxy properties like size, and galaxies with similar properties can be clustered, then it is possible this scheme may not be satisfactory for high precision corrections.
We scale the noise bias correction from the fiducial value according to signal-to-noise, as expected in the analytic calculations of Refregier et al. (2012). This works well for the lower noise images but less well for those at higher noise levels as also noted in (Kacprzak et al., 2012). As stated in that work we are not satisfied with the performance of our optimiser in this highest noise regime and would cut such images from cosmic shear analyses at this point, or calibrate using a more extensive suite of population matched simulations. A similar concern applies to the smallest galaxies.
For all other branches of GREAT08 we attain additive and multiplicative biases which are within the requirements for current surveys, and have quality factor metrics which are as good or better than other non-stacking methods at the time of the GREAT08 Challenge.
Further work appears to be required before im3shape can be applied to the final years of upcoming Stage III cosmic shear surveys. We speculate that our residual shear measurement bias is due to the combined effects of using the incorrect galaxy model in the fit, and the noise on the images. Ideally an improved shear measurement method could be found that did not suffer from these types of bias. Alternatively the calibration scheme needs to be tailored more accuractly to the analysis in question, for example by using a suite of low noise images with known shears.
To be fully applicable to real data there are a number of extensions required to the code, for example (i) the ability to handle multiple exposures of the same patch of sky (ii) the capacity to deal with or flag neighboring objects contaminating the shear measurement.
Furthermore, we have not investigated here the effect of postage stamp size or the effect of a weight map cutting out parts of the galaxy image. These developments are ongoing but testing them is beyond the scope of the present work.
We thank Steve Gull, Jean-Paul Kneib, Eduardo Cypriano, Phil Marshall, John Bridle, Patrick Hudelot, Sebastien Bardeau, David Hogg, Alexandre Refregier, Adam Amara, Filipe Abdalla, Niall Maccrann, Emmanuel Bertin, Jim Bosch, Mike Jarvis, Fabrizio Sidoli, Gary Bernstein, Catherine Heymans and Tom Kitching for helpful discussions, and Dugan Witherick for patient and extensive computational help.
JZ, TK, SB, BR, LV and MH acknowledge support from European Research Council in the form of a Starting Grant with number 240672.
- Amara & Réfrégier (2008) Amara A., Réfrégier A., 2008, MNRAS, 391, 228
- Bardeau et al. (2005) Bardeau S., Kneib J.-P., Czoske O., Soucail G., Smail I., Ebeling H., Smith G. P., 2005, A&A, 434, 433
- Bardeau et al. (2007) Bardeau S., Soucail G., Kneib J.-P., Czoske O., Ebeling H., Hudelot P., Smail I., Smith G. P., 2007, A&A, 470, 449
- Bernstein (2010) Bernstein G. M., 2010, MNRAS, 406, 2793
- Bernstein & Jarvis (2002) Bernstein G. M., Jarvis M., 2002, Astronomical Journal, 123, 583
- Bridle et al. (2010) Bridle S. et al., 2010, MNRAS, 405, 2044
- Bridle et al. (2008) Bridle S., et al., 2008, ArXiv 0802.1214
- Bridle et al. (2010) Bridle S., et al., 2010, MNRAS, 405, 2044
- Bridle et al. (2008) Bridle S. et al., 2008, ArXiv 0802.1214
- Bridle et al. (2002) Bridle S. L., Kneib J.-P., Bardeau S., Gull S. F., 2002, in Natarajan P., ed., The Shapes of Galaxies and their Dark Halos. pp 38–46
- Chang et al. (2012) Chang C. et al., 2012, ArXiv e-prints: astro-ph/1206.1378
- de Jong et al. (2012) de Jong J. T. A., Verdoes Kleijn G. A., Kuijken K. H., Valentijn E. A., 2012, Experimental Astronomy, p. 34
- Frigo & Johnson (2005) Frigo M., Johnson S. G., 2005, Proceedings of the IEEE, 93, 216
- Heymans et al. (2006) Heymans C., et al., 2006, MNRAS, 368, 1323
- Heymans et al. (2012) Heymans C. et al., 2012, MNRAS, 427, 146
- Hirata & Seljak (2003) Hirata C., Seljak U., 2003, MNRAS, 343, 459
- Hirata et al. (2004) Hirata C. M. et al., 2004, MNRAS, 353, 529
- Hoekstra et al. (1998) Hoekstra H., Franx M., Kuijken K., Squires G., 1998, Ap. J., 504, 636
- Huff et al. (2011) Huff E. M., Eifler T., Hirata C. M., Mandelbaum R., Schlegel D., Seljak U., 2011, ArXiv astro-ph/1112.3143
- Kacprzak et al. (2012) Kacprzak T., Zuntz J., Rowe B., Bridle S., Refregier A., Amara A., Voigt L., Hirsch M., 2012, ArXiv 1203.5049
- Kaiser (2000) Kaiser N., 2000, Ap. J., 537, 555
- Kaiser et al. (1995) Kaiser N., Squires G., Broadhurst T., 1995, Ap. J., 449, 460
- Kitching et al. (2009) Kitching T. D., Amara A., Abdalla F. B., Joachimi B., Refregier A., 2009, MNRAS, 399, 2107
- Kitching et al. (2012) Kitching T. D., et al., 2012, MNRAS, 423, 3163
- Kitching et al. (2008) Kitching T. D., Miller L., Heymans C. E., van Waerbeke L., Heavens A. F., 2008, MNRAS, 390, 149
- Kitching (2010) Kitching T. e. a., 2010, ArXiv 1009.0779
- Kuijken (1999) Kuijken K., 1999, A&A, 352, 355
- Kuijken (2006) Kuijken K., 2006, A&A, 456, 827
- Laureijs et al. (2011) Laureijs R. et al., 2011, ArXiv e-prints: astro-ph/1110.3193
- Luppino & Kaiser (1997) Luppino G. A., Kaiser N., 1997, Ap. J., 475, 20
- Mandelbaum et al. (2006) Mandelbaum R., Hirata C. M., Broderick T., Seljak U., Brinkmann J., 2006, MNRAS, 370, 1008
- Mandelbaum et al. (2012) Mandelbaum R., Hirata C. M., Leauthaud A., Massey R. J., Rhodes J., 2012, MNRAS, 420, 1518
- Mandelbaum et al. (2005) Mandelbaum R. et al., 2005, MNRAS, 361, 1287
- Marks (2009) Marks R. J., 2009, Handbook of Fourier Analysis & Its Applications, 1st edn. Oxford University Press
- Massey et al. (2007) Massey R., et al., 2007, MNRAS, 376, 13
- Massey & Refregier (2005) Massey R., Refregier A., 2005, MNRAS, 363, 197
- Massey et al. (2007) Massey R., Rowe B., Refregier A., Bacon D. J., Bergé J., 2007, MNRAS, 380, 229
- Melchior et al. (2010) Melchior P., Böhnert A., Lombardi M., Bartelmann M., 2010, A&A, 510, A75+
- Melchior & Viola (2012) Melchior P., Viola M., 2012, MNRAS, 424, 2757
- Melchior et al. (2011) Melchior P., Viola M., Schäfer B. M., Bartelmann M., 2011, MNRAS, 412, 1552
- Miller et al. (2012) Miller L. et al., 2012, ArXiv e-prints: astro-ph/1210.8201
- Miller et al. (2007) Miller L., Kitching T. D., Heymans C., Heavens A. F., van Waerbeke L., 2007, MNRAS, 382, 315
- Nakajima & Bernstein (2007) Nakajima R., Bernstein G., 2007, Astronomical Journal, 133, 1763
- Okura & Futamase (2009) Okura Y., Futamase T., 2009, Ap. J., 699, 143
- Okura & Futamase (2012) Okura Y., Futamase T., 2012, Ap. J., 748, 112
- Press et al. (1992) Press W. H., Teukolsky S. A., Vettering W. T., Flannery B. P., 1992. Cambridge: Cambridge Univ. Press
- Refregier & Bacon (2003) Refregier A., Bacon D., 2003, MNRAS, 338, 48
- Refregier et al. (2012) Refregier A., Kacprzak T., Amara A., Bridle S., Rowe B., 2012, ArXiv 1203.5050
- Reyes et al. (2010) Reyes R., Mandelbaum R., Seljak U., Baldauf T., Gunn J. E., Lombriser L., Smith R. E., 2010, Nature, 464, 256
- Rhodes et al. (2000) Rhodes J., Refregier A., Groth E. J., 2000, Ap. J., 536, 79
- Rowe et al. (2012) Rowe B., Bacon D., Massey R., Heymans C., Haeussler B., Taylor A., Rhodes J., Mellier Y., 2012, ArXiv e-prints: astro-ph/1211.0966
- Schneider (2006) Schneider P., 2006, in Meylan G., Jetzer P., North P., Schneider P., Kochanek C. S., Wambsganss J., eds, Saas-Fee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro. pp 269–451
- Schrabback et al. (2010) Schrabback T., et al., 2010, A&A, 516, A63
- Takada (2010) Takada M., 2010, in Kawai N., Nagataki S., eds, American Institute of Physics Conference Series Vol. 1279, American Institute of Physics Conference Series. pp 120–127
- The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration 2005, ArXiv Astrophysics e-prints
- Velander et al. (2011) Velander M., Kuijken K., Schrabback T., 2011, MNRAS, 412, 2665
- Viola et al. (2011) Viola M., Melchior P., Bartelmann M., 2011, MNRAS, 410, 2156
- Voigt & Bridle (2010) Voigt L. M., Bridle S. L., 2010, MNRAS, 404, 458
Appendix A Modelling and parameterization
a.1 Choice of ellipticity parameters
The performance of all the optimizers we have tried is heavily dependent on the particular choice of parameterization. In particular, the choice of how to parameterize the shear components is particularly important since they are our target. We define ellipticity to match the effect of shear on images: . The two obvious parameterizations based on this, and , were both found to be problematic.
Using the two separate ellipticities and as parameters creates a non-linear boundary at , the edge of the unit circle.
While most optimization methods cleanly handle linear boundaries (for example simple limits on parameters), non-linear ones are less well served. We tried simply assigning an extremely low likelihood to any proposed sample beyond this boundary, but fast optimization methods use likelihood differences to determine step sizes, so they typically behave wildly if large negative values are used to signal edges. Since the images of model galaxies diverge as simply making the models of these galaxies is also prone to numerical problems.
Using the polar coordinates as parameters works well at high ellipticity where we use a simple parameter limit on , but fails at lower ellipticities. As the likelihood becomes a plateau in since the coordinates are degenerate. Optimizers spent an extremely long time wandering around the circle.
We found a reliable parameter set which is a variant of the first parameterisation, . We use as the parameters varied by the optimizer, but map them to a disc wherever they stray beyond it. We mirror the parameter space radially about some just below unity (we found 0.95 worked well), so that as the input increases above the maximum, the value given to the minimiser first decreases and eventually reaches zero, and turns negative so that the signs of the output and are flipped. Although this mirror structure causes repetition in the likelihood surface and multiple modes we have not found this to be an issue in practice.
a.1.1 Other parameters
There are various radius parameters we could choose for each of the components. The choice is particularly important for the sharp bulge profiles, as their full-width half-maxima are extremely small and so they are not suitable parameters.
A better choice that works for both components is the radius which encloses half the total flux of the complete profile. Rather than using two radius parameters, we the disc radius and the ratio between the radius of the bulge and disc. We find (see section 4.2) that no large model bias is generated by fixing the latter parameter provided it is approximately correct. We therefore use a fixed value for any given data set. For our GREAT08 runs this value was arbitrarily set to unity, to avoid trying to match too closely to the GREAT08 truth values.
A likelihood plateau occurs for small objects like stars, where the shears are irrelevant. In these regimes the fitted radii of the objects tended to fail to converge as they varied around zero. We solved this problem by setting a small transition radius sub-pixels and a replacing the radius with when and otherwise. The minimum value radius actually used is then sub-pixels.
We have experimented with marginalization over the amplitude components of the models. Since the model is linear in these parameters their likelihood through a given slice in the other parameters is Gaussian and we can analytically marginalize over them. We have found this to be slightly problematic in general - the marginalization is not particularly numerically stable and since we have two components we need to perform the PSF convolution twice if we sum them after marginalization. These two effects meant it was in fact slower to use this trick, provided that we have a reasonable initial estimate of the ellipticity, particularly as our gradient-based minimizer finds the linear parameters very simple to optimize.
If we move the modelling into Fourier space we could further marginalize over the centroid parameters and of the model - combining these dimensionality reductions could well yield a significant speed up. The downside of this is decreased flexibility in the modelling. We defer further discussion of this issue.
a.1.3 Point spread function (PSF)
im3shape has the general facility to use an arbitrary image for PSF, but by default we read a catalog of either Moffat or Airy function parameters and use that to build a PSF image. It is built at the same upsampled resolution as the image itself, and transformed and stored in Fourier space for faster convolution with trial images.
When these images are made the convolution is done using the FFTW package Frigo & Johnson (2005).
a.1.4 Fourier Sampling Kernel
Our model for the pixels in an image is really a constant value across the surface of the pixel. But as far as a Discrete Fourier transform is concerned a 2D array represents a grid of function samples from the image. To account for this difference our Fourier-space model of the image needs to include convolution with a square top-hat function the size of a pixel - this is true even when simulating at the same resolution as the data. Since we are also convolving with a PSF we can fold together these convolutions into a single Fourier-space product. The Fourier transform of a 2D top-hat is given by:
where the image dimensions are and the upsampling factor is .
The model image in im3shapeis made at a resolution several times higher than the data resolution, and the convolution with the PSF also happens at this high resolution. Only then is the model downsampled to compare to the data. This procedure ensures that small and sharp features in the galaxy can be captured.
We find, though, that the central regions of galaxy bulges with the standard de Vaucouleurs profile are so sharp that an even higher resolution is needed to accurately model them. Rather than generate the whole image at this extreme resolution we only use it for the central pixels where it is relevant.
There are therefore three different resolution parameters to consider: the upsampling (the ratio of the model to image pixel size in the outer image region), the size of the central higher resolution region, and the resolution increase in this central region (compared to the outer region). The values that we used for these parameters are discussed in more detail in section 4.1.
Appendix B Numerical optimization
The optimization problem of finding the best model parameters can be tackled with a choice of algorithms; we have found that the best performing is the Levenberg-Marquadt (LevMar) method. This method uses more information than methods (such as Powell’s method or Nelder-Meald) which operate only on the total likelihood as a function of the input parameters: it uses the complete mismatch map, (i.e. the per pixel) for each parameter set considered.
The LevMar method uses gradients (analytic or numerical) to optimize; it is therefore not possible to signal a failure of a particular parameter set (which can occur when parameters are near their limits) using NaNs or other flag values. While we do try to put in place a parameterization that prevents this from happening (see below), we also signal the failure by setting the for the failing model to 2 per pixel, a value large enough to be rejected but hopefully not hugely problematic.
b.1 Termination criteria
The LevMar method depends on a number of manually set thresholds that determine when the numerical optimization should terminate. In particular, it stops and reports the parameters of the best model fit when at least one of the following criteria is met:
One of the partial derivatives of the objective function () with respect to the model parameters falls below a threshold , i.e. , where denotes one of the model parameters. Note that this termination criterion can be sensitive to the changes of image pixel scaling and changes in the signal-to-noise ratio (SNR).
The relative change of the norm of the estimated parameter vector is smaller than , i.e. , where is the current iteration index. This termination criterion can be sensitive to individual parameter scaling.
The value of falls below ; again, this termination criterion is sensitive to image scaling and the S/N ratio.
The relative change of falls below , i.e. , where is the objective value at the current iteration and is the total number of pixels. This termination criterion was implemented by the authors and is not included in the LevMar implementation we are using. It can be sensitive to the step size that LevMar is taking and thus will be affected by the levmar_init_mu parameter - the initial damping of the normal equations solved in the Levenberg - Marquadt algorithm.
The maximum number of iterations is reached (the minimization has failed).
Appendix C Analysis Overview
Before analyzing a collection of galaxy images we:
Read a collection of options specifying the model to be used, parameter ranges, and various options from a parameter file
Build a model object encapsulating the model we will fit to the galaxies, including its likelihood and prior functions
Before analyzing each individual image we then:
Generate the point-spread image based on PSF parameters; compute its FT.
Generate the sinc function for the specified upsampling (the FT of a top-hat).
Multiply these two Fourier space quantities together to generate the Fourier kernel.
Generate a weight map based on the image noise levels.
Compute a starting estimate for the galaxy properties from defaults and weighted quadrupole moments.
We then pass the model likelihood function to a minimizer object which runs the chosen algorithm. At each evaluation of the posterior, for our standard model, we:
Check if the parameters are outside the desired ranges and give up early if not.
Build high-resolution (up-sampled) real-space models of each component of the galaxy (e.g. bulge and disc).
Replace the central few pixels of each model with summed up pixels from even higher resolution, since this is the region where the image changes quickly.
Sum and convolve the components with the precomputed kernel.
Down-sample the image to the data resolution.
Compute the likelihood from the weight map and model image.
Save the model image for use in the optimizer and residual analysis.
Appendix D Details on the maximum likelihood estimation
The maximization problem outlined in section B is based on the following model:
where denotes the observed galaxy image, , the convolution of the PSF and a square top-hat pixel window function (see section A.1.4), the (higher resolution) galaxy model with parameter vector , and the downsampling operator, which averages into lower resolution pixels. As detailed in section 3 we assume a two-component Sérsic galaxy model with the parameters (see Table 2 for a detailed parameter description)
where for the disc component and for the bulges, and denotes the sheared squared distance
with the centered coordinate vector and the covariance matrix
For additive Gaussian noise, i.e. maximizing the likelihood is equivalent to minimizing
with respect to the model parameters , where denotes the residual image between the observed galaxy image and the downsampled, PSF convolved galaxy model image , i.e.
For simplicity we assume to be uncorrelated white noise and so consider only diagonal terms of the covariance matrix, i.e. ; we neglect any correlation in noise of adjacent pixels.
For minimizing (7), im3shape uses the free software package levmar131313Available at http://www.ics.forth.gr/~lourakis/levmar/, a C/C++ implementation of the Levenberg-Marquardt (LM) algorithm, a standard technique for solving non-linear least-squares problems (Press et al., 1992). The LM method is based on a linear approximation of in the neighbourhood of , i.e.
where is a small perturbation of the parameter vector and denotes the Jacobian matrix of :
The Jacobian (10) can be numerically approximated, but to speed convergence of the minimization procedure it is generally recommended to compute it analytically. For our model (3) the Jacobian reads with respect to the model parameters
Note that for clarity and the ease of exposition we omitted the image upsampling within a central image region as described in section 4.1 in the derivation above, i.e. we considered the special case where is set to one. The generalization to is straightforward since it only involves an additional summation over all sub-pixels to yield the flux for each pixel within the central image region.
Appendix E Noise bias calibration using deep data
Noise bias is the difference between the expectation of the maximum-likelihood result found by a model fitter like im3shape, which is a biased estimator, and the true underlying value, coming from the non-linear mapping between parameters and image. Its typical value is a few percent, and we need to remove it to reach our target accuracy levels.
The size of the bias depends sensitively on the true parameters of the galaxy, and if these were known we could remove noise bias completely. However, we have access only to noisy estimates of these parameters and therefore our estimate of the noise bias is itself noisy, and biased. We can think of this as a bias-on-bias problem, and we find that failing to account for it means we significantly miss our targets.
We can get around this problem if there is a subset of our observational data that is deeper than the bulk of our sample, and therefore of greater signal-to-noise, to a degree sufficient that it has negligible noise bias. This is often the case in real surveys that seek to detect high redshift supernovae, and it is approximated in the GREAT08 challenge with LowNoise_Blind data set (although the galaxy types in this set do not match those in the main sample exactly). In this case we do not calibrate each galaxy individually, but instead find a mean bias for the population and apply it en masse.
Biases and will afflict each of our galaxies, depending on the true galaxy properties . We can calculate mean values of this bias and , and apply these evenly to all the galaxies, such that the mean galaxy ellipticity (which is our goal) is correct.
The mean of the bias across the population is given by:
We find the distribution using fits to the deep data, and using simulations - for each point in a grid in the parameter space we simulate many galaxies and determine the value . We can then interpolate between these values to do the integral, and finally apply the mean to all the galaxy estimates.
A similar process is used for the bias, except that varies directly with the PSF, so we use this information - we calculate the assuming a fiducial PSF with ellipticity aligned with the direction. The applied value to apply to each galaxy is then:
Appendix F Results tables
Tables 3-6 give the GREAT08 , , , and values for im3shape with (I3) and without (I3-U) noise bias calibration applied. Also shown are the scores for the other entrants at the time of the challenge, as described in Bridle et al. (2010).