© ESO, 2015

1. Introduction

Weak gravitational lensing is one of the most powerful cosmological probes used to study the nature of dark energy and its evolution with redshift. An important step towards achieving this goal relies on the accurate measurements of the cosmic shear, i.e., the tiny change in ellipticity of the observed galaxies due to foreground cosmic structures. The measurement must be performed statistically, using millions of galaxies whose images are affected by the convolution of the atmospheric and instrumental point spread function (PSF) and by noise. Reliable image processing techniques are required to correct for the PSF smearing and to extract the faint cosmic shear signal.

So far, only a handful of truly different shape measurement algorithms are available. They split into four broad categories. The first category, known as forward fitting methods, fits analytical models of galaxies in order to measure their shape, decontaminated from the PSF (e.g., Gentile et al. 2012; Kacprzak et al. 2012; Miller et al. 2007). The second category decomposes the image of galaxies on the basis of functions of different shapes and reconstructs the original image (e.g., Refregier 2003; Kuijken 2006), again decontaminated from the PSF. The third class of methods is based on the calculation of the moments of the light distribution of galaxies and of their associated PSFs (Kaiser et al. 1995; Bernstein 2010). More recently, methods based on machine learning have also started to be implemented (e.g., Gruen et al. 2010; Tewes et al. 2012).

Because the above methods show different levels of sensitivity to noise and to various sources of systematic errors, it is worth exploring additional ways of performing PSF correction in weak lensing measurements. One way of doing this is to carry out full image deconvolution of the data prior to the shape measurement, which is the work undertaken in this article. Our new algorithm proposes to describe the image distortion by the PSF in a matrix equation using linear algebra formalism. In this formalism, the deconvolution problem can be solved by applying the Hopfield neural network iterative scheme. The galaxy ellipticities are then measured using the second-order moments of the two-dimensional auto-correlation function (ACF) of deconvolved images.

We start by introducing the image deconvolution problem and we discuss the method based on Hopfield neural networks to solve it in Sect. 2. In Sect. 3 we present our new adaptive neuron updating rule for the astronomical image reconstruction and the shape measurement algorithm for practical implementation. Sect. 4 deals with the calibration of the Tikhonov regularization parameter on simulated data in order to find the best value. We present our application of the method to the data of the GREAT10 Galaxy challenge (Kitching et al. 2011, 2012) in Sect. 5 and we give a conclusion in Sect. 6.

2. Image deconvolution problem

An observed image can be modeled as a function A of the true source image augmented by a noise term ${\mathit{u}}_{\mathrm{0}}\mathrm{=}\mathit{A}\mathrm{\left(}{\mathit{u}}_{\mathrm{s}}\mathrm{\right)}\mathrm{+}\mathit{n,}$(1)where u₀ and u_s denote the observed image and the source image, respectively, and n is additive noise.

In the case of telescope images, A is commonly represented by a convolution with a known PSF constrained from the images of stars, and a regular sampling by the detector. The problem of image restoration is to recover an unknown source image u_s from a given observed, blurred image u₀ and PSF h. Eq. (1) takes the form ${\mathit{u}}_{\mathrm{0}}\mathrm{=}\mathit{h}\mathrm{\otimes }{\mathit{u}}_{\mathrm{s}}\mathrm{+}\mathit{n,}$(2)where ⊗ denotes the convolution operation, and h is the convolution kernel representing the PSF of the telescope.

The deconvolution problem described above is a well-known ill-posed problem, which does not have a unique solution, as the convolution and the noise effectively destroy some information (Allison 1979). One possible approach for solving this problem is regularization as the convolution and the noise effectively destroy the telescope images. A widely used regularization method consists in minimizing the Tikhonov functional (Tikhonov & Arsenin 1977), $\mathit{J}\mathrm{\left(}\mathit{u}\mathrm{\right)}\mathrm{=}\mathrm{\parallel }{\mathit{u}}_{\mathrm{0}}\mathrm{-}\mathit{h}\mathrm{\otimes }\mathit{u}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{+}\mathit{\lambda }\hspace{0.17em}\mathrm{\parallel }\mathit{D}\mathrm{\left(}\mathit{u}\mathrm{\right)}{\mathrm{\parallel }}^{\mathrm{2}}\mathit{,}\hspace{0.17em}$(3)where ∥ . ∥ denotes the L₂ norm, and D is a differential operator. The regularization parameter λ> 0 is introduced to control the trade-off between the image bias and image variance (Haykin 2008). The image that minimizes the functional J(u) above is an estimation of the intrinsic image u_s (Eq. (1)).

2.1. Toeplitz matrix formulation

In linear algebra, the convolution operator can be represented as a matrix-vector multiplication by constructing a special Toeplitz matrix for a particular PSF kernel (Gray 2006). The Tikhonov functional J(u) is expressed in a vector form as: $\mathit{J}\mathrm{\left(}{u}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\parallel }{{u}}_{\mathrm{0}}\mathrm{-}\mathit{H}\hspace{0.17em}{u}{\mathrm{\parallel }}^{\mathrm{2}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\mathit{\lambda }\hspace{0.17em}\mathrm{\parallel }\mathit{D}{u}{\mathrm{\parallel }}^{\mathrm{2}}\mathit{,}$(4)where u₀ is the observed image and, the variable u denotes an approximation of the sours image u_s, both represented as column vectors, and H is a Toeplitz matrix encoding the convolution.

The convolution operation is now represented as a matrix multiplication, where known PSF image, h(x,y), is converted into the corresponding Toeplitz matrix H. A Toeplitz matrix is a sparse diagonal-constant matrix. If one considers the discrete convolution of a n × n image I with a n × n PSF kernel h, the convolved image may be presented as ${\mathit{I}}_{\mathrm{conv}}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\mathrm{=}\sum _{{\mathit{x}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{\mathit{n}\mathrm{-}\mathrm{1}}\sum _{{\mathit{y}}^{\mathrm{\prime }}\mathrm{=}\mathrm{0}}^{\mathit{n}\mathrm{-}\mathrm{1}}\mathit{I}\left(\mathit{x}\mathrm{-}{\mathit{x}}^{\mathrm{\prime }}\mathit{,y}\mathrm{-}{\mathit{y}}^{\mathrm{\prime }}\right)\hspace{0.17em}\mathit{h}\left(\mathit{x}\mathrm{-}{\mathit{x}}^{\mathrm{\prime }}\mathit{,y}\mathrm{-}{\mathit{y}}^{\mathrm{\prime }}\right)\mathit{,}$(5)where x and y are the discrete indices of rows and columns, respectively, and I_conv(x,y) is the result of the convolution operation for the pixel at position (x,y).

Fig. 1 Toeplitz matrix creation schema. A 5 × 5 image is convolved by a 3 × 3 PSF kernel. The resulting Toeplitz matrix is 25 pixels on a side.

Figure 1 illustrates how to generate the first row of a Toeplitz matrix H from a kernel h. Each row of the Toeplitz matrix is a vector composed by “flattened”, row by row, kernel matrix coefficients corresponding to a particular pixel position in the image. For a n × n image, the corresponding Toeplitz matrix has a n² × n² dimension.

2.2. Application of HNN to the image deconvolution problem

Fig. 2 Architecture of the Hopfield network of four neurons. The neurons are presented as light blue disks numbered 1, 2, 3, and 4, xi(k), (i = 1,2,3,4) present the corresponding neuron state at time k (iteration number). Wij, (i,j = 1,2,3,4) are the weights of interconnections between neurons. b(k) denotes the bias at current time k. The output of the sum function f(k) in iteration k is the input vector for the activation function. The activation function uses the thresholding rule to calculate the neuron updates to change their states. Z-1 denotes the time delay of the iterative circuit. Since the network energy E converges, the loop stops with the output vector yi(k), (i = 1,2,3,4).

The Hopfield neural network (HNN) model is a widely-used approach in many different optimization problems (Hopfield 1982; Rojas 1996). Neural networks are decision making tools operating on statistical data, inspired by the functioning of a biological nervous system. Each neural network model has a set of neurons (nodes) characterized by the basic elements:

• weights (strength of the connections between neurons),

• sum function (calculates the total signal incoming from all connected neurons),

• activation function (calculates the net output).

The principle of a HNN is shown in Fig. 2. The Hopfield model is a recurrent model, which consists of only one layer of neurons all connected with each other, so that each neuron is the input and output at the same time, creating a self-feedback loop with delayed connection from a neuron to itself. The weight matrix between neurons are symmetrical, e.g., w_ij = w_ji. One of the advantages of HNN is that no prior information is necessary, so the neurons can be initialized with any non-negative values.

The image reconstruction problem described above (Eq. (2)) can be solved using a HNN whose neurons correspond to the pixels of the reconstructed image, as is described in the following. Sun (2000) proposed the modified HNN model, where the activation function is a threshold function, determined by the General Updating Rule (GUR), which is theoretically proven to converge to the best solution of Eq. (2). According to the GUR, the updated value of the neuron i at iteration k is calculated as ${\mathit{v}}_{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ {\mathit{v}}_{\mathit{i}}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{-}\mathrm{\Delta }{\mathit{v}}_{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,} & \mathit{f}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{<}\mathrm{-}{\mathit{TH}}_{\mathit{i}}\\ {\mathit{v}}_{\mathit{i}}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{+}\mathrm{\Delta }{\mathit{v}}_{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,} & \mathit{f}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{>}\hspace{0.17em}{\mathit{TH}}_{\mathit{i}}\\ {\mathit{v}}_{\mathit{i}}\mathrm{(}\mathit{k}\mathrm{-}\mathrm{1}\mathrm{)}\mathit{,} & \mathrm{otherwise}\mathit{,}\end{array}$(6)where Δv_i(k) is the update value which equal the unit of the image brightness, and TH and f are the threshold and sum function output, respectively.

Hopfield nets have a scalar estimator associated with each iteration of the network referred to as the energy function, E of the network. At each iteration the neurons are updated according to the GUR, decreasing the norm of the energy function. The self-feedback loop of HNN stops at the stable state for the network – when there are no neurons to update according the activation function Eq. (6).

The energy function of HNN in a matrix form is defined as (Sun 2000) $\mathit{E}\mathrm{\left(}{u}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}{{u}}^{\mathit{T}}\mathrm{(}\mathit{W}\mathrm{+}\mathit{\lambda }\hspace{0.17em}{\mathit{G}}^{\mathrm{)}}\hspace{0.17em}{u}\mathrm{-}{\mathit{Q}}^{\mathit{T}}{u}\mathit{,}$(7)with $\mathit{W}\mathrm{=}{\mathit{H}}^{\mathit{T}}\mathit{H}\mathit{,}\mathit{G}\mathrm{=}{\mathit{D}}^{\mathit{T}}\mathit{D,}\mathit{Q}\mathrm{=}{\mathit{H}}^{\mathit{T}}{{u}}_{\mathrm{0}}\mathit{,}$(8)where W is the synaptic weight matrix related to the PSF, the Q term is associated with a bias vector of HNN and the combination (W + λG) is the interconnection matrix between neurons, represented by image pixels.

On the other hand, expanding Eq. (4), we have $\mathit{J}\mathrm{\left(}{u}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{{u}}^{\mathit{T}}\mathrm{(}{\mathit{H}}^{\mathit{T}}\mathit{H}\mathrm{+}\mathit{\lambda }\hspace{0.17em}{\mathit{D}}^{\mathit{T}}\mathit{D}\mathrm{)}{u}\mathrm{-}\mathrm{(}{\mathit{H}}^{\mathit{T}}{{u}}_{\mathrm{0}}{\mathrm{)}}^{\mathit{T}}{u}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{{u}}_{\mathrm{0}}^{\mathit{T}}{{u}}_{\mathrm{0}}\mathit{,}$(9)using (8) in Eq. (9), we obtain $\mathit{J}\mathrm{\left(}{u}\mathrm{\right)}\mathrm{=}\begin{array}{c}\underset{\mathrm{􏽼}0\mathrm{􏽻􏽺}0\mathrm{􏽽}}{\frac{\mathrm{1}}{\mathrm{2}}{{u}}^{\mathit{T}}\mathrm{(}\mathit{W}\mathrm{+}\mathit{\lambda }\hspace{0.17em}\mathit{G}\mathrm{)}{u}\mathrm{-}{\mathit{Q}}^{\mathit{T}}{u}}\\ \mathit{E}\mathrm{\left(}{u}\mathrm{\right)}\end{array}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{{u}}_{\mathrm{0}}^{\mathit{T}}{{u}}_{\mathrm{0}}\mathit{.}$(10)One can see that the Tikhonov regularization functional J(u) and the energy function of the HNN network E(u) differ only in the term $\frac{\mathrm{1}}{\mathrm{2}}{{u}}_{\mathrm{0}}^{\mathit{T}}{{u}}_{\mathrm{0}}$ which does not depend on u. Therefore, by minimizing the HNN energy we can find the best estimation of a solution of Eq. (1).

The network is run as a series of iterative steps where the values of selected neurons (pixels) are updated so as to minimize the energy of the whole network. Image pixels are updated according to the condition Eq. (6), where the threshold value is given by ${\mathit{TH}}_{\mathit{i}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\hspace{0.17em}\begin{array}{c}\mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\end{array}{}^{}{\mathit{W}}_{\mathit{ii}}\mathrm{+}\mathit{\lambda }\hspace{0.17em}{\mathit{G}}_{\mathit{ii}}\begin{array}{c}\mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\end{array}{}^{}\hspace{0.17em}\mathit{,}\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathit{..}\mathit{,}{\mathit{N}}^{\mathrm{2}}\mathit{.}$(11)At each iteration the pixels exceeding the threshold are updated, minimizing the energy function. The network stops running when its total energy E reaches a global minimum, that is, when no pixel needs to be updated the image is assumed to be fully deconvolved.

2.3. Modification of the Hopfield neural network: adaptive neuron updating rule

In the application of HNN to natural images the updating value Δv_i(k) is an integer constant, usually equal to 1 (for instance, for the images where the pixel value varies from 1 to 256). However, astronomical images have a great range of pixel intensities. Such an image can be divided into two regions, the background area that contains very little signal, and the regions of astronomical objects themselves, where the dynamical range is large. These special properties of astronomical images must be taken into account to determine the best updating value.

For the HNN deconvolution the selection of the most accurate pixel updating step is very important as:

• a small updating value requires a large number of iterations, making the minimization process very slow;

• a large updating value can result in an instability of the HNN energy function, preventing the neural network from converging.

We propose a new adaptive neuron updating rule, which takes into account the specific properties of astronomical images and that accelerates the HNN deconvolution process. At each iteration k, the neuron update for a given pixel Δv_i(k), is proportional to the gradient of the energy function, which allows us to reduce the difference between the observed image and the image being reconstructed convolved with the PSF kernel, $\mathrm{\Delta }{\mathit{v}}_{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathit{c}\begin{array}{c}\mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\end{array}{}^{}\mathrm{\nabla }\mathit{E}\mathrm{\left(}\mathit{v}\mathrm{\right(}\mathit{k}\mathrm{\left)}\mathrm{\right)}\begin{array}{c}\mathrm{}\\ \mathrm{}\\ \mathrm{}\\ \mathrm{}\end{array}{}^{}\mathrm{=}\mathit{c}\hspace{0.17em}\mathrm{\left[}\mathrm{\right(}\mathit{W}\mathrm{+}\mathit{\lambda }\hspace{0.17em}\mathit{G}\mathrm{)}\hspace{0.17em}\mathit{v}\mathrm{(}\mathit{k}\mathrm{)}\mathrm{-}\mathit{Q}\mathrm{]}\mathit{,}$(12)where c is a constant which scales the updating values, and therefore also influences the number of required iterations. The set of neuron updating values Δv(k) is an array, where each element is associated with certain pixel in the image.

3. Galaxy shape measurement

The galaxy shape measurement process consists of three main steps:

• 1.

  Galaxy image reconstruction by iterative application of the HNN image deconvolution algorithm (Sect. 3.1) to “postage stamps” containing a single galaxy image each.

• 2.

  Computation of the two-dimensional ACF for each deconvolved galaxy image (Sect. 3.2).

• 3.

  Measurement of the galaxy ellipticities by calculating the quadrupole moments of their ACF images (Sect. 3.3).

3.1. Image deconvolution scheme

The deconvolution algorithm presented in Sect. 2 is implemented in the following way:

We represent each galaxy postage stamp I of dimension n × n pixels, as a “flattened” vector v (Sect. 2). For each galaxy vector v we:

• Generate a vector y of length n² by stacking the rows of the observed image u₀;

• Generate the Toeplitz matrix H using the PSF kernel (Gray 2006) and the matrix D using a high pass filter d, designed such that the summation of d_ij is equal to zero. In our method we use the kernel d proposed by Sun (2000): $\begin{array}{ccc}\mathit{d}\mathrm{=}\frac{\mathrm{1}}{\mathrm{100}}\hspace{0.17em}\left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathrm{;}& & \end{array}$

• Calculate the matrices W, G and Q (Eq. (8));

• Set c = 0.001 ∗ max(u₀) and λ = 0.0001;

• Set k = 0andX(0) = Y;

• Compute the threshold vector according to Eq. (11);

• Set k = 1 and repeat the following steps while k<M:

  • 1.

    Calculate the negative gradient vector of the energy function (Eq. (12)):$\mathit{f}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathit{Q}\mathrm{-}\mathrm{(}\mathit{W}\mathrm{+}\mathit{\lambda }\hspace{0.17em}\mathit{G}\mathrm{)}\hspace{0.17em}\mathit{v}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{;}$(13)

  • 2.

    Compute the neuron updating vector Δv(k) according to Eq. (12);

  • 3.

    Update the vector v(k) according to the GUR (Eq. (6)). If HNN stopped at the stable state, e.g., there is no pixel to update, then break the iteration loop;

  • 4.

    Correct the vector v(k) as follows: ${\mathit{v}}_{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ \mathrm{0}\mathit{,} & {\mathit{v}}_{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{<}\mathrm{0}\\ {\mathit{v}}_{\mathit{i}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,} & \mathrm{otherwise}\mathrm{;}\end{array}$(14)

  • 5.

    Check that E(k + 1) <E(k), otherwise break the iteration loop;

• Reshape the flattened vector v(k) into a two-dimensional image I^⋆.

3.2. The autocorrelation function in a shape measurement

Fig. 3 Illustration of the HNN deconvolution, using galaxies from the GREAT10 image processing challenge. The three rows show three different galaxies. In each row, are shown from left to right: the original GREAT10 postage stamp for each galaxy, the associated PSF, the HNN deconvolution, the residual image (difference between data and deconvolved images, reconvolved by the PSF), the ACF of the deconvolved galaxy, and the ACF clipped at the radius that encloses 90% of the total light.

Once the deconvolved galaxies, I^⋆, have been obtained, the next step is to measure their ellipticities. Our goal with the HNN deconvolution method is to allow shape measurements free of any underlying model for the light distribution of galaxies. Therefore the ellipticity measurement has to be carried out in a non-parametric way, e.g., using the second-order moments of the light distribution of I^⋆.

In practice, such a measurement is sensitive to noise and requires the estimation of the centroid of the galaxies. A safe and convenient way to circumvent both limitations is to measure the galaxy ellipticities using the ACF of the light distribution and not the galaxy images directly. Indeed, it has been shown by Miralda-Escudé (1991) and proven later by van Waerbeke et al. (1997) that the ellipticity of the ACF is the same as the ellipticity of the image itself. In addition, by construction, the ACF of any function is exactly centered at the center of the frame. We therefore measure the ellipticity of the galaxy ACF rather than the ellipticity of the galaxy itself.

Figure 3 shows examples of galaxies borrowed from the GREAT10 image processing challenge, along with their HNN deconvolution and ACF functions. To minimize edge effects caused by the square nature of the galaxy stamps, the ACF is clipped at the radius that encloses 90% of the total light. Such a cutoff radius removes nearly all edge effects and improves the calculation of the ellipticity through the second-order moments method.

3.3. Estimation of galaxy ellipticities

We use the ACF of the deconvolved galaxy images, I^⋆, to measure their complex ellipticity, e = e₁ + i e₂. The ellipticity components (e₁,e₂) are derived from the unweighted quadrupole moments Q, calculated around the center of the light distribution (x_c,y_c) ${\mathit{e}}_{\mathrm{1}}\mathrm{=}\frac{{\mathit{Q}}_{\mathit{xx}}\mathrm{-}{\mathit{Q}}_{\mathit{yy}}}{{\mathit{Q}}_{\mathit{xx}}\mathrm{+}{\mathit{Q}}_{\mathit{yy}}}\mathit{,} {\mathit{e}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{2}\hspace{0.17em}{\mathit{Q}}_{\mathit{xy}}}{{\mathit{Q}}_{\mathit{xx}}\mathrm{+}{\mathit{Q}}_{\mathit{yy}}}\mathit{,}$(15)where $\begin{array}{c}\\ {\mathit{Q}}_{\mathit{xx}}\mathrm{=}\mathrm{\int }{\mathit{I}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\mathrm{d}\mathit{x}\hspace{0.17em}\hspace{0.17em}\mathrm{d}\mathit{y}\mathit{,}\\ {\mathit{Q}}_{\mathit{xy}}\mathrm{=}\mathrm{\int }{\mathit{I}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathrm{c}}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{c}}\mathrm{)}\hspace{0.17em}\hspace{0.17em}\mathrm{d}\mathit{x}\hspace{0.17em}\mathrm{d}\mathit{y}\mathit{,}\\ {\mathit{Q}}_{\mathit{yy}}\mathrm{=}\mathrm{\int }{\mathit{I}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{x,y}\mathrm{\right)}\hspace{0.17em}\hspace{0.17em}\mathrm{(}\mathit{y}\mathrm{-}{\mathit{y}}_{\mathrm{c}}{\mathrm{)}}^{\mathrm{2}}\hspace{0.17em}\hspace{0.17em}\mathrm{d}\mathit{x}\hspace{0.17em}\mathrm{d}\mathit{y}\mathit{.}\end{array}$(16)In the present case, we use the ACF of the light distribution, so (x_c,y_c) are known exactly and coincide with the center of the image.

4. Tikhonov regularization parameter

The important question in Tikhonov regularization is how to determine the best parameter λ> 0 in order to find a solution u^∗ close to the true noise-free solution u. The parameter λ regulates the balance between the accuracy decreasing with λ, and the smoothness of the restored image increasing with λ. In order to find the best value of λ we applied the image deconvolution method described above to the synthetic data. Hereafter we discuss the calibration of λ based on simulations.

Each synthetic dataset consists of 10 000 galaxy images of a Sérsic profile convolved with a simulated Moffat-profile PSF. The noise added to the blurred images has a Gaussian distribution with mean zero and the corresponding signal-to-noise ratio (S/N) ranges from 10 to 40 dB. The S/N of the convolved noisy image is estimated as $\mathit{S}\mathit{/}\mathit{N}\mathrm{=}\mathrm{10}{\log }_{\mathrm{10}}\mathrm{(}{\mathit{\sigma }}_{\mathrm{S}}^{\mathrm{2}}\mathit{/}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{)}$, where ${\mathit{\sigma }}_{\mathrm{S}}^{\mathrm{2}}$ is a variance from the noise-free convolved image and σ² is the variance of the noise added to the convolved image.

Fig. 4 Correct transition rate as a function of signal-to-noise ratio for the adaptive update approach (black dashed line) and contents update approach (blue solid line). λ = 10-3.5.

To evaluate the efficiency of our algorithm we use the metrics called Correct Transition Rate (CTR) and discriminative signal-to-noise ratio (dS/N) described below (Sun 2000; Ben-Arie & Rao 1993). To recover the individual galaxy image a certain number of iterations is required. In each iteration, the states of some neurons change according the GUR. If the state transition for a given neuron j decreases the difference between the state of the neuron and the corresponding pixel value of the original image | v^j(k) − u^j |, the transition is called correct transition. The CTR is defined as $\mathrm{CTR}\mathrm{=}\frac{{\mathit{K}}_{\mathrm{c}}}{{\mathit{K}}_{\mathrm{t}}}\mathit{,}$(17)where K_c is the number of the correct transitions, and K_t is the total number of the neural state transitions. The CTR varies from 0 to 1 and measures the convergence rate of the algorithm. The higher the CTR is, the better the efficiency of the algorithm.

Figure 4 shows the CTR factor for both the proposed adaptive neuron updating and standard constant neuron updating. One can see that, whatever the S/N, the adaptive neuron updating approach demonstrates a higher convergence rate compared to the constant updating.

Fig. 5 Correct transition rate as a function of the Tikhonov regularization parameter λ for S/N = 10 (red dotted line), S/N = 20 (blue dashed line), S/N = 40 (black solid line).

Fig. 6 Discriminative signal-to-noise ratio as a function of the Tikhonov regularization parameter λ for S/N = 10 (red dotted line), S/N = 20 (blue dashed line), S/N = 40 (black solid line).

Fig. 7 Accuracy on the shape measurement (as inverse RMSD) as a function of the Tikhonov regularization parameter λ for two different methods: the adaptive neuron update (black dashed line) and constant neuron update (blue solid line). Figures from top to bottom correspond to S/N equal 10, 20 and 40 dB, respectively.

The CTR values for λ varying from 10^-6 to 10⁶ is shown in Fig. 5 the CTR increases with S/N, though for λ> 10^-2.5 the difference on CTR is negligible. The highest CTR value corresponds to λ = 10^-3.5.

To estimate the performance of the image reconstruction, we use the dS/N criterion estimated as $\mathrm{d}\mathit{S}\mathit{/}\mathit{N}\hspace{0.17em}\mathrm{\equiv }\mathrm{10}\hspace{0.17em}{\log }_{\mathrm{10}}\left(\frac{\mathrm{\parallel }{\mathit{u}}_{\mathrm{0}}\mathrm{-}\mathit{u}{\mathrm{\parallel }}^{\mathrm{2}}}{\mathrm{\parallel }{\mathit{u}}^{\mathrm{\ast }}\mathrm{-}\mathit{u}{\mathrm{\parallel }}^{\mathrm{2}}}\right)\mathit{,}$(18)where u is the original simulated galaxy image, u^∗ is the final restored image and u₀ denotes the convolved noisy image. As shown in Fig. 6, dS/N is higher for higher S/N and decreases with λ.

The main parameter which controls the performance of the proposed approach is the accuracy on shape measurement. We measured the ellipticities of the reconstructed galaxies using the technique explained in Sect. 3. The accuracy on shape measurement is estimated using the inverse Root Mean Square Deviation (RMSD) $\mathrm{1}\mathit{/}\mathrm{RMS}\mathrm{D}\hspace{0.17em}\mathrm{\equiv }\sqrt{\frac{\mathrm{2}\mathit{N}}{{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathrm{(}{\mathit{e}}_{\mathrm{1}\mathit{,i}}^{}\mathrm{-}{\mathit{e}}_{\mathrm{1}\mathit{,i}}^{\mathrm{\ast }}\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{e}}_{\mathrm{2}\mathit{,i}}^{}\mathrm{-}{\mathit{e}}_{\mathrm{2}\mathit{,i}}^{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{2}}}}$(19)where ${\mathit{e}}_{\mathrm{1}\mathit{,i}}^{}$, ${\mathit{e}}_{\mathrm{2}\mathit{,i}}^{}$ are the intrinsic ellipticities and ${\mathit{e}}_{\mathrm{1}\mathit{,i}}^{\mathrm{\ast }}$, ${\mathit{e}}_{\mathrm{2}\mathit{,i}}^{\mathrm{\ast }}$ are the estimated ellipticities calculated from (Eq. (15)); i denotes the galaxy number in the set with N = 10 000 galaxies.

Figure 7 shows the accuracy on the shape measurement for both adaptive neuron update and constant neuron update. One can see that for λ ≤ 10^-4 the adaptive neuron update gives better accuracy regardless the noise level of the images. For λ> 10^-4 the accuracy is sensitive to the noise level of the images. While on S/N level of 20 and 40 the adaptive neuron update still shows better performance, on noisier images with S/N = 10 it works less well than with constant neuron updating.

Taking into account the CPU time and the quality of the reconstructed images (Figs. 4–7) the best value of λ = 10^-4 is suggested.

5. Application to the GREAT10 data

5.1. The GREAT10 galaxy challenge and its dataset

The image deconvolution scheme described in the paper was used to participate in the Gravitational LEnsing Accuracy Testing 2010 challenge (GREAT10), an image processing competition that ran from December 2010 to September 2011 (Kitching et al. 2011, 2012) and whose goals are (i) to test existing shear measurements methods and (ii) to promote the development of new shear measurement techniques. The GREAT10 challenge is the continuation of the GREAT08 challenge (Bridle et al. 2009, 2010) and STEP programs (Heymans et al. 2006; Massey et al. 2007; Bridle et al. 2010), with an increasing degree of complexity.

The data of the GREAT10 challenge consist of 24 datasets of 200 simulated galaxy fields. Each image contains 10 000 noisy, PSF-convolved galaxy images, arranged on a 100 × 100 grid. Each galaxy field is an image stamp of 48 × 48 pixels.

The main innovation in the 2010 edition of the challenge is the inclusion of variable components for the intrinsic ellipticities for the galaxies and of the shear fields. The PSF fields are also spatially varying and are provided to the challenge participants both as functional form and as FITS images.

With the simulated data, participants to the challenge must measure the ellipticities of the galaxies, as they were prior to the convolution by the PSF and the addition of noise. The results are submitted to the challenge organizers either under the form of a catalogue of ellipticities or as a shear power spectrum. A “quality factor”, Q is then computed by comparing the submissions of participants to the true characteristics of the simulations, known only to the organizers of the challenge.

Each GREAT10 image set is designed to evaluate the ability of competing methods to deal with galaxy or PSF fields with different properties (e.g., constant versus varying field, size, S/N). We summarize in Table 2 the main PSF and galaxy characteristics for each of the GREAT10 dataset.

The exact metric used for assessing various characteristics of the methods have been described in the GREAT10 Galaxy challenge paper (Kitching et al. 2012). It consists of:

• A “raw” quality factor Q that measures the difference, averaged over all sets, between the reconstructed and true shear power spectra. The value of the raw Q factor is used to rank the different submissions to the challenge and to show the leaderboard of the competition.

• A quality factor Q_dn is also obtained, after removal of biases caused by finite S/N or inherent shape measurement method-related noise.

• A split of the total bias on the galaxy ellipticities into an additive and a multiplicative biases, respectively denoted as c and m.

• Additional parameters and ℳ ≃ m² + 2m, intended to measure the additive and multiplicative biases calculated at power spectrum level.

Fig. 8 Q factor comparison for all 24 sets of GREAT10 data for different method variants (see Table 1). The last line give the mean Q factor for all 24 sets.

5.2. GREAT10 results

While a general analysis of the GREAT10 results for all participating teams is available in Kitching et al. (2012), we focus here on a more specific analysis of the results obtained using our HNN deconvolution.

In order to reduce the CPU time, we used smaller postage stamps of the galaxy images: 19 × 19 , 23 × 23 and 29 × 29 pixels. Our submissions to GREAT10 were done under the names NN19 (official challenge submission) and NN23 (post-challenge submission) which stands for “Neural Network”, followed by the stamp size of the galaxy images, in pixels. In both submissions, we used the PSF under its functional form, but we reconstructed the PSF on a 15 × 15 postage stamp to carry out the deconvolution.

Fig. 9 Shear power spectra for the 24 datasets of the GREAT10 challenge, for our submission “23/den/star” (see Table 1). The red lines show our estimate of the shear power spectra while the green lines represent the true shear power.

We also made submissions after the deadline of the competition and in some cases, we applied image denoising to the data prior to the deconvolution, following Nurbaeva et al. (2011). The submissions with denoising are labeled “den” in Table 1 and the ones without denoising are labeled “raw”. In these submissions we used postage stamps with 23 pixels on-a-side and we used both the functional (labeled “func”) form for the PSF and the pixelated FITS PSFs (labeled “star”).

Fig. 10 Multiplicative bias per set, ℳ/2 on the shear power spectrum, as defined in Kitching et al. (2012). The ℳ/2 values are given on the left panel for our 23/raw/star submission and the right panel shows the 23/den/star submission.

Fig. 11 Additive bias per set, $\sqrt{\mathrm{𝒜}}$ on the shear power spectrum, as defined in Kitching et al. (2012). The $\sqrt{\mathrm{𝒜}}$ values are given on the left panel for our 23/raw/star submission and the right panel shows the 23/den/star submission.

5.3. Analysis of the GREAT10 results

The quality factors, Q, obtained for our submissions are summarized in Table 1 as well as in Fig. 8. Figure 9 shows the power spectra for the “23/den/star” submission, i.e., for a deconvolution using 23 × 23 pixels stamps, denoising, and the pixelated version of the PSF. Finally, Figs. 10 and 11 show the multiplicative bias ℳ/2, and additive bias $\sqrt{\mathrm{𝒜}}$ computed for the shear power spectra.

It is interesting to analyze first the influence of our different submissions on the Q factors shown in Fig. 8. Before doing so, we should note that our first submission, NN19, was affected by a misunderstanding of the centering convention in the GREAT10 challenge. This resulted in a very low Q factor for this submission, which we discard in the following discussion. In the other submissions, we note that increasing the stamp size used for the deconvolution significantly improves the Q factor. A stamp size of 29 pixels, however, results in large cpu times, of the order of 4 s per galaxy. This would make the method very hard to use in practice on data from a real survey.

All of our mean Q factors are between 70 and 90, meaning that on average the use of the functional form for the PSF or the pixelated PSF has little influence, perhaps because the signal-to-noise in the pixelated PSFs is high enough to make the functional form unnecessary. We also note that denoising improves the Q factor by 5 to 10%.

One can see that the good Q factors are obtained for the sets 10 and 11, which have a realistic signal-to-noise distribution. The same is observed for sets 17 and 18 which contain smooth galaxies. Among these 4 sets, we see no difference in terms of Q when the PSF is kept fixed or variable.

Finally, maybe the most interesting observations are for Set 20 which contains turbulent PSFs and for Sets 23 and 24, which contain offset bulge/disk galaxies. Both characteristics are typical for realistic ground-based surveys which will contain complex PSFs and galaxies. Obtaining our best Q factors for these sets is therefore very encouraging. The purpose of direct image deconvolution is indeed to enable the use of any arbitrary galaxy shape, independent of a simplistic underlying model. This gain in galaxy complexity is made to the prize of a lower Q factor for sets with low signal-to-noise galaxies and for small galaxies. This is not unexpected as deconvolution in general tend to enhance noise and to perform less well as sampling is getting worse, i.e., exactly for small and faint galaxies.

Obtaining good Q factors is necessary but not sufficient to allow weak lensing measurements. We also need to ensure that the shape measurement method does not introduce significant bias on the shear power spectrum. These biases are summarized in Figs. 9–11. It is interesting to see that, while the Q factors for the low signal-to-noise galaxies and small galaxies are poor, the biases on the shear power spectrum remain comparable with those of sets with high signal-to-noise and large galaxies. In addition, the biases on the offset bulge/disk galaxies remain low. Finally, we see that denoising decreases the multiplicative bias and slightly increases the additive bias and this almost independently of the properties of the galaxies and PSFs in each set.

6. Conclusion

We have introduced a new method to spatially deconvolve images, and we have applied it for the first time to the field of weak gravitational lensing. While the deconvolution process tends to enhance the noise, it has the advantage to allow ellipticity measurements with no assumption of an underlying light profile. The noise enhancement itself is in part circumvented by measuring the galaxy ellipticities on the autocorrelation function of the deconvolved image. The autocorrelation function has the advantage of being centered for each stamp, resulting in higher accuracy of the shape measurement.

We have then confronted our algorithm to the simulated data proposed in the GREAT10 image processing challenge. We find that image deconvolution performs well on complex galaxy shapes although with somewhat lower performances at low signal-to-noise and for small galaxies (compared with the PSF size). These results are almost independent of the PSF properties, i.e., shape and spatial variations. We also find that image denoising improves the precision of the measurement (Q factor) without affecting much the accuracy, measured by the multiplicative and additive biases on the shear power spectrum.

The above behavior of our method is observed for the specific data of the GREAT10 challenge and would need more investigations either on more complex simulations or on real data. The method remains cpu-intensive and requires in practice about 1 s per galaxy to converge properly. The results obtained on GREAT10 are nevertheless encouraging, as the method performs best on more complex galaxy shapes. Application to the future data of the GREAT3 challenge (Mandelbaum et al. 2013), which will include HST-like galaxies, will be a crucial test for direct image deconvolution.

Acknowledgments

We are grateful to the GREAT10 Coordination Team for organizing this stimulating challenge, in particular Tom Kitching for his help and for sharing the shear analysis code. We thank Matthew Luke Nichols for helpful discussions. Malte Tewes acknowledges support by the DFG grant Hi 1495/2-1. This work is supported by the Swiss National Science Foundation (SNSF). The GREAT10 challenge itself was sponsored by a EU FP7 PASCAL 2 challenge grant.

References

All Tables

All Figures

	Fig. 1 Toeplitz matrix creation schema. A 5 × 5 image is convolved by a 3 × 3 PSF kernel. The resulting Toeplitz matrix is 25 pixels on a side.
In the text

Fig. 2 Architecture of the Hopfield network of four neurons. The neurons are presented as light blue disks numbered 1, 2, 3, and 4, xi(k), (i = 1,2,3,4) present the corresponding neuron state at time k (iteration number). Wij, (i,j = 1,2,3,4) are the weights of interconnections between neurons. b(k) denotes the bias at current time k. The output of the sum function f(k) in iteration k is the input vector for the activation function. The activation function uses the thresholding rule to calculate the neuron updates to change their states. Z-1 denotes the time delay of the iterative circuit. Since the network energy E converges, the loop stops with the output vector yi(k), (i = 1,2,3,4).

In the text

Fig. 3 Illustration of the HNN deconvolution, using galaxies from the GREAT10 image processing challenge. The three rows show three different galaxies. In each row, are shown from left to right: the original GREAT10 postage stamp for each galaxy, the associated PSF, the HNN deconvolution, the residual image (difference between data and deconvolved images, reconvolved by the PSF), the ACF of the deconvolved galaxy, and the ACF clipped at the radius that encloses 90% of the total light.

In the text

	Fig. 4 Correct transition rate as a function of signal-to-noise ratio for the adaptive update approach (black dashed line) and contents update approach (blue solid line). λ = 10-3.5.
In the text

	Fig. 5 Correct transition rate as a function of the Tikhonov regularization parameter λ for S/N = 10 (red dotted line), S/N = 20 (blue dashed line), S/N = 40 (black solid line).
In the text

	Fig. 6 Discriminative signal-to-noise ratio as a function of the Tikhonov regularization parameter λ for S/N = 10 (red dotted line), S/N = 20 (blue dashed line), S/N = 40 (black solid line).
In the text

	Fig. 7 Accuracy on the shape measurement (as inverse RMSD) as a function of the Tikhonov regularization parameter λ for two different methods: the adaptive neuron update (black dashed line) and constant neuron update (blue solid line). Figures from top to bottom correspond to S/N equal 10, 20 and 40 dB, respectively.
In the text

	Fig. 8 Q factor comparison for all 24 sets of GREAT10 data for different method variants (see Table 1). The last line give the mean Q factor for all 24 sets.
In the text

	Fig. 9 Shear power spectra for the 24 datasets of the GREAT10 challenge, for our submission “23/den/star” (see Table 1). The red lines show our estimate of the shear power spectra while the green lines represent the true shear power.
In the text

	Fig. 10 Multiplicative bias per set, ℳ/2 on the shear power spectrum, as defined in Kitching et al. (2012). The ℳ/2 values are given on the left panel for our 23/raw/star submission and the right panel shows the 23/den/star submission.
In the text

	Fig. 11 Additive bias per set, $\sqrt{\mathrm{𝒜}}$ on the shear power spectrum, as defined in Kitching et al. (2012). The $\sqrt{\mathrm{𝒜}}$ values are given on the left panel for our 23/raw/star submission and the right panel shows the 23/den/star submission.
In the text