Three-dimensional Reconstruction of Weak-lensing Mass Maps with a Sparsity Prior. I. Cluster Detection

In order to reconstruct high-resolution, high-S/N mass maps, we incorporate prior information on the density contrast field into the reconstruction by modeling the density field as a sum of basis atoms in a "dictionary":

$$\begin{eqnarray}&&\delta ={\boldsymbol{\Phi }}\cdot x,\end{eqnarray} \tag{ 2 }$$

where Φ is the matrix operator transforming from the projection coefficient vector x to the density contrast field δ. The column vectors of Φ are the basis "atoms" of the model dictionary and denoted as ϕ_s .

There have been a few studies that adopt different dictionaries for 3D weak-lensing map reconstructions: (i) Simon et al. (2009) perform a reconstruction in Fourier space, which is equivalent to representing the density contrast field with sinusoidal functions. However, sinusoidal functions are not localized in configuration space, and the density contrast field is not sparse in Fourier space. Therefore, sparsity priors cannot be directly applied to this model dictionary for high-resolution mass-map reconstructions. (ii) Leonard et al. (2014) model the density contrast field with starlets (Starck et al. 2015). However, the starlet dictionary does not account for the angular scale difference at different lens redshifts, and is not specifically designed to model the clumpy mass distribution in the universe. It is worth exploring other dictionaries for our purpose of weak-lensing mass reconstruction.

In the standard cosmological model, dark matter is concentrated in roughly spherical "halos", which have the NFW density profile (Navarro et al. 1997). Motivated by this fact, we generate a model dictionary using NFW atoms with N typical scale radii r_s (s = 1,...,N) in comoving coordinates. An atom has the NFW surface density profile on the transverse plane (Takada & Jain 2003) and the Dirac δ function in the line-of-sight direction. Following Leonard et al. (2014), we neglect the size of halos along the line of sight since the resolution scale of the reconstruction is much larger than the halos.

The multiscale NFW atom is defined as

$$\begin{eqnarray}&&{\phi }_{s}({\vec{r}}_{\theta },z)=\displaystyle \frac{f}{2\pi {r}_{s}^{2}}\,F(| {\vec{r}}_{\theta }| /{r}_{s})\,{\delta }_{D}(z),\end{eqnarray} \tag{ 3 }$$

where ${\vec{r}}_{\theta }$ is the projection of the comoving position on the transverse plane,

$$\begin{eqnarray}F(x)=\left\{\begin{array}{ll}-\tfrac{\sqrt{{c}^{2}-{x}^{2}}}{(1-{x}^{2})(1+c)}+\tfrac{\mathrm{arccosh}\left(\tfrac{{x}^{2}+c}{x(1+c)}\right)}{{\left(1-{x}^{2}\right)}^{3/2}} & (x\lt 1),\\ \tfrac{\sqrt{{c}^{2}-1}}{3(1+c)}\left(1+\tfrac{1}{c+1}\right) & (x=1),\\ -\tfrac{\sqrt{{c}^{2}-{x}^{2}}}{(1-{x}^{2})(1+c)}+\tfrac{\arccos \ \left(\tfrac{{x}^{2}+c}{x(1+c)}\right)}{{\left({x}^{2}-1\right)}^{3/2}} & (1\lt x\leqslant c),\\ 0 & (x\gt c),\end{array}\right.\end{eqnarray} \tag{ 4 }$$

$f=1/[\mathrm{ln}(1+c)-c/(1+c)]$, and c is the halo concentration. For simplicity, we fix c = 4 for the NFW atoms in our dictionary. In the present paper, we adopt a hard truncation on the NFW profile at a radius equals cr_s . Studying the influence of different truncation forms (Oguri & Hamana 2011) on the mass-map reconstruction is left to our future work.

Compared to other model dictionaries, our NFW dictionary is motivated by a physical consideration on the clumpy mass distribution in the universe. Furthermore, the multiscale NFW atoms are set up in comoving coordinates, with an account of the scale difference in the angular coordinates for halos at different lens redshifts. The corresponding NFW atom in the angular separation coordinates is

$$\begin{eqnarray}&&{\phi }_{s}(\vec{\theta },z)=\displaystyle \frac{f{\chi }^{2}(z)}{2\pi {r}_{s}^{2}}\,F(| \vec{\theta }| \chi (z)/{r}_{s})\,{\delta }_{D}(z),\end{eqnarray} \tag{ 5 }$$

where χ(z) is the comoving distance to redshift z.

For the NFW dictionary, the transform from the projection coefficient vector to the density contrast field of Equation (2) is written as

$$\begin{eqnarray}&&\delta (\vec{r})=\sum _{s=1}^{N}\int {d}^{3}r^{\prime} {\phi }_{s}(\vec{r}-\vec{r^{\prime} }){x}_{s}(\vec{r^{\prime} }).\end{eqnarray} \tag{ 6 }$$

To simplify the notation, we compress the projection coefficient vectors on the basis atoms with multiple scales into a single column vector:

$$\begin{eqnarray}&&x=\left(\begin{array}{c}{x}_{0}\\ {x}_{1}\\ ...\\ {x}_{N}\end{array}\right),\end{eqnarray} \tag{ 7 }$$

and write the dictionary transform operator into a row vector:

$$\begin{eqnarray}&&{\boldsymbol{\Phi }}=\left(\int {d}^{3}r\,{\phi }_{0}(\vec{r})\,\int {d}^{3}r\,{\phi }_{1}(\vec{r})\,...\,\int {d}^{3}r\,{\phi }_{N}(\vec{r})\right).\end{eqnarray} \tag{ 8 }$$

For an additional test and comparison, we also construct a dictionary with the point-mass atoms which are represented by the 3D Dirac function as

$$\begin{eqnarray}&&{\phi }_{\mathrm{PM}}(\vec{\theta },z)={\delta }_{D}({\theta }_{1})\,{\delta }_{D}({\theta }_{2})\,{\delta }_{D}(z).\end{eqnarray} \tag{ 9 }$$

The 2D profiles of the NFW atoms and the point-mass atom on the transverse plane are shown in Figure 1. The corresponding 1D profiles are plotted in Figure 2. The profiles are smoothed with a Gaussian kernel and pixelized on linearly spaced grids. Details on the smoothing and pixelizing operations are described in Section 2.3.2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We define the forward transform operator A = P · Q · Φ where P represents systematic effects and Q represents the physical lensing effect. With Equations (1) and (2), we write the transform from the coefficient vector x to the observed lensing-shear field as

$$\begin{eqnarray}&&\gamma ={\boldsymbol{A}}\cdot x+\epsilon .\end{eqnarray} \tag{ 10 }$$

The lensing-transform operator can be expressed as

$$\begin{eqnarray}&&{\boldsymbol{Q}}={\int }_{0}^{{z}_{s}}{{dz}}_{l}\,K({z}_{l},{z}_{s})\int {d}^{2}{\theta }^{{\rm{{\prime} }}}\,\overrightarrow{{\rm{D}}}(\overrightarrow{\theta }-\overrightarrow{{\theta }^{{\rm{{\prime} }}}}).\end{eqnarray} \tag{ 11 }$$

K(z_l , z_s ) is the lensing kernel between lens redshift z_l and source redshift z_s (Bartelmann & Schneider 2001), which is given by

$$\begin{eqnarray}K({z}_{l},{z}_{s})=\left\{\begin{array}{ll}\tfrac{3{H}_{0}{{\rm{\Omega }}}_{m}}{2c}\tfrac{{\chi }_{l}{\chi }_{{sl}}(1+{z}_{l})}{{\chi }_{s}E\left({z}_{l}\right)} & ({z}_{s}\gt {z}_{l}),\\ 0 & ({z}_{s}\leqslant {z}_{l}),\end{array}\right.\end{eqnarray} \tag{ 12 }$$

where E(z) is the Hubble parameter as a function of redshift, in units of H₀.

$$\begin{eqnarray}&&D(\vec{\theta })=-\displaystyle \frac{1}{\pi }{\left({\theta }_{1}-i{\theta }_{2}\right)}^{-2}\end{eqnarray} \tag{ 13 }$$

is the Kaiser–Squires kernel (Kaiser & Squires 1993), which decays proportionally to the inverse square of the distance. Here θ_1,2 are the two components of the angular position vector $\vec{\theta }$.

The top panel of Figure 3 shows the lensing kernels as a function of source redshift for lenses at different lens redshifts. The bottom panel of Figure 3 shows the correlation between the lensing kernels. Each lensing kernel is highly nonlocal, and the lensing kernels of different lens redshifts are strongly correlated, as can be seen in the bottom panel. These inherent properties of the lensing kernels result in strong correlation between the column vectors that constitute the forward transform matrix A . As will be discussed in Section 2.4, the strong correlation makes it challenging to reconstruct mass maps with high resolution in the line-of-sight direction.

Zoom In Zoom Out Reset image size

Shear measurement deviates from the true, physical shear owing to a variety of systematic effects in real observations. In this section, we discuss the influence of several major systematics on the lensing-shear measurement, and describe the corresponding transform operator by decomposing into three parts of R (photometric redshift uncertainties), W (smoothing), and M (masking).

2.3.1. Photometric Redshift Uncertainty

Photometric redshifts (photo-z) are estimated with a limited number of broad optical and near-infrared bands in the current generation weak-lensing surveys. For example, nine bands are used for the KIDS+VIKING survey (Hildebrandt et al. 2020), and five bands for the DES and HSC Survey. Unlike the high-precision spectroscopic redshift (spect-z) estimation, the photo-z estimation suffers from considerable statistical uncertainties.

Figure 4 shows the histogram of the best-fit estimates of the Machine Learning for Photo-z (MLZ; Carrasco Kind & Brunner 2013) algorithm ⁸ for galaxies in tract ⁹ 9347 from the photo-z catalog (Tanaka et al. 2018) of the first-year HSC data release (Aihara et al. 2018). The galaxies are divided into 10 bins according to the photo-z best-fit estimates. Figure 5 shows the average probability density function (PDF) for galaxies in each redshift bin.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the presence of photo-z uncertainties, a source galaxy with a best-fit photo-z estimate z_s has a posterior probability P(z∣z_s ) of being actually located at redshift z. This means that there is a possibility P(z∣z_s ) that the galaxy image is actually distorted by the shear at redshift z. Therefore, the photo-z uncertainty smears the lensing kernel statistically, which we model by a smearing operator

$$\begin{eqnarray}&&{\boldsymbol{R}}=\int {{dz}}_{s}P(z| {z}_{s}).\end{eqnarray} \tag{ 14 }$$

Figure 3 shows the lensing kernels and their correlations for source redshifts with photo-z uncertainties demonstrated in Figure 5. Compared to the lensing kernels of spect-z that converge to zero at source redshifts lower than the lens redshift, the lensing kernels of photo-z do not converge to zero at the low source galaxy redshifts. This is simply because the galaxies with photo-z estimated lower than the lens redshifts may actually be located at higher redshifts. In addition, the photo-z uncertainty only slightly increases the correlations between lensing kernels at different lens plane.

2.3.2. Smoothing

The source galaxies are not uniformly distributed in the sky, with substantial fluctuations in the number density. We first smooth the lensing-shear measurements, and pixelize the smoothed shear field onto a regular grid. After these procedures, the fast Fourier transform can be directly conducted on the transverse plane in each redshift bin to compute the convolution operation in Equation (1). Another benefit of smoothing is that it reduces bias arising from the aliasing effect in the pixelization process since the smoothing kernel reduces the amplitude of the shear signal at high frequency.

We convolve the lensing shear measured from galaxy images with a 3D smoothing kernel $W(\vec{\theta },z)$. The expectation of the lensing-shear field after smoothing is

$$\begin{eqnarray}&&{\gamma }_{\mathrm{sm}}(\vec{\theta },z)=\displaystyle \frac{{\sum }_{i}W(\vec{\theta }-{\vec{\theta }}_{i},z-{z}_{i}){\gamma }_{L}({\vec{\theta }}_{i},{z}_{i})}{{\sum }_{i}W(\vec{\theta }-{\vec{\theta }}_{i},z-{z}_{i})},\end{eqnarray} \tag{ 15 }$$

where z_i and θ_i refer to the photometric redshift and transverse position of the ith galaxy in the galaxy sample, respectively. γ_L refers to the shear field before the smoothing.

The smoothing kernel $W(\vec{\theta },z)$ can be decomposed into a transverse component ${W}_{t}(\vec{\theta })$ and a line-of-sight component W_l (z) as

$$\begin{eqnarray}&&W(\vec{\theta },z)={W}_{t}(\vec{\theta })\,{W}_{l}(z).\end{eqnarray} \tag{ 16 }$$

We use an isotropic 2D Gaussian kernel and a 1D top-hat kernel to smooth the shear field in the transverse plane and in the line-of-sight direction, respectively. These components are given by

$$\begin{eqnarray}\begin{array}{rcl}{W}_{t}(\vec{\theta }) & = & \frac{1}{2\pi {\beta }^{2}}\exp \left(-\frac{| \vec{\theta }| }{2{\beta }^{2}}\right),\\ {W}_{l}(z) & = & \left\{\begin{array}{ll}1/{\rm{\Delta }}z & (| z| \lt {\rm{\Delta }}z/2),\\ 0 & \mathrm{otherwise},\end{array}\right.\end{array}\end{eqnarray} \tag{ 17 }$$

where we set β = 1.5' in this paper. By definition, the smoothing kernel is normalized as

$$\begin{eqnarray}&&\int {d}^{3}r\,W(\vec{r})=1.\end{eqnarray} \tag{ 18 }$$

Assuming that the galaxy number distribution varies slowly at the smoothing scale, the smoothed shear can be written into

$$\begin{eqnarray}&&{\gamma }_{\mathrm{sm}}={\boldsymbol{W}}\cdot {\gamma }_{L},\end{eqnarray} \tag{ 19 }$$

where W is the smoothing operator defined as

$$\begin{eqnarray}&&{\boldsymbol{W}}=\int {d}^{3}r^{\prime} \,W(\vec{r}-\vec{r^{\prime} }).\end{eqnarray} \tag{ 20 }$$

We note that another widely used scheme is to average the shear measurements in each pixel. Such a scenario is equivalent to resampling the shear field smoothed with a 3D top-hat kernel with the same scale as the pixels.

After smoothing, we pixelize the shear field on N_x × N_y × N_z grids, where N_x and N_y are the number of pixels of the two orthogonal axes of the transverse plane, and N_z is the number of pixels in the line-of-sight direction. We denote by γ_α the smoothed shear measurements recorded on the pixel with index α, where 1 ≤ α ≤ N_x × N_y × N_z . The grids on the transverse planes are equally spaced with a pixel size of 1'. In the line-of-sight direction, we set binning with an equal galaxy number as shown in Figure 4.

Similarly, we pixelize the projection coefficient vector x onto an N_x × N_y × N_l grid. In this paper, the projection coefficient vector is pixelized in equal spacing ranging from redshift 0.01 to redshift 0.85. Here, we use N_l to denote the number of the lens planes and x_β to denote the projection field element with index 1 ≤ β ≤ N_x × N_y × N_l × N, where the last N is the number of NFW dictionary frames representing different physical scale radii. The corresponding pixelized elements of the forward transform matrix A is denoted as A_{α β} .

In order to see the effect of smoothing clearly, we plot the histograms of the statistical shear measurement errors due to shape noise and sky variance for the galaxies in tract 9347 of the first-year HSC shear catalog (Mandelbaum et al. 2018). In Figure 6, we compare the error on an individual galaxy basis and that of individual pixels after the smoothing and pixelization procedures described in this section. The standard deviation of the individual pixel errors is much smaller than for the galaxies owing to the smoothing. We also note that even though the shear measurement error for the galaxies does not follow a Gaussian distribution, the error after smoothing and pixelization is well described by a Gaussian distribution, which is simply explained by the central limit theorem.

Zoom In Zoom Out Reset image size

2.3.3. Masking

In real observations, shear measurements can be performed in a finite region of the sky, and the boundaries often have complicated geometries. Moreover, many isolated subregions near the bright stars are masked out.

We define the masking-window function as

$$\begin{eqnarray}M(\vec{r})=\left\{\begin{array}{ll}1 & {n}_{\mathrm{sm}}\geqslant 1,\\ 0 & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 21 }$$

where n_sm is the smoothed galaxy number density. We define the masking operator as

$$\begin{eqnarray}&&{\boldsymbol{M}}=\int {d}^{3}r^{\prime} \,M(\vec{r^{\prime} }){\delta }_{D}(\vec{r}-\vec{r^{\prime} }),\end{eqnarray} \tag{ 22 }$$

where ${\delta }_{D}(\vec{r})$ is 3D Dirac delta function.

Taking into account the systematics discussed previously in Sections 2.3.1–2.3.3, the systematic operator is

$$\begin{eqnarray}&&{\boldsymbol{P}}={\boldsymbol{M}}\cdot {\boldsymbol{W}}\cdot {\boldsymbol{R}}.\end{eqnarray} \tag{ 23 }$$

The projection coefficients can be estimated by optimizing a penalized loss function. An estimator is generally defined as

$$\begin{eqnarray}&&\hat{x}=\mathop{\arg \,\min }\limits_{x}\left\{\displaystyle \frac{1}{2}{}_{{\rm{\Sigma }}}\parallel (\gamma -{\boldsymbol{A}}x){\parallel }_{2}^{2}+\lambda C(x)\right\},\end{eqnarray} \tag{ 24 }$$

where ${}_{{\rm{\Sigma }}}\parallel (\gamma -{\boldsymbol{A}}x){\parallel }_{2}^{2}$ is the l² norm of residuals weighted by the inverse of the covariance matrix Σ of the shear measurement error, which measures the difference between the prediction and the data. C(x) is the penalty term measuring the deviation of the coefficient estimate x from the prior assumption. The estimation with the "penalty" term prefers parameters that are able to describe the observation with a specified prior information. The penalty parameter λ adjusts the relative weight between the data and the prior assumption in the optimization process.

There have been a few studies that adopt different penalties for 3D weak-lensing map reconstructions: (i) Simon et al. (2009) propose to use the Wiener filter, which is also known as l² ridge penalty ($C=\parallel x{\parallel }_{2}^{2}$), to find a penalized solution in Fourier space. Oguri et al. (2018) apply the method of Simon et al. (2009) to the first-year data of the HSC Survey (Aihara et al. 2018). It is found that the density maps reconstructed by the method suffer from significant line-of-sight smearing with a standard deviation of σ_z = 0.2 − 0.3. (ii) Leonard et al. (2014) use the GLIMPSE algorithm, which adopts a derivative version of the l¹ LASSO penalty ($C=\parallel x{\parallel }_{1}^{1}$) to find a sparse solution in the starlet dictionary (Starck et al. 2015). GLIMPSE reduces the smearing by adopting a "greedy" coordinate descent algorithm that forces the structure to grow only on the most related lens redshift plane.

In this paper, we use another derivative version of the LASSO penalty—the adaptive LASSO penalty. We first normalize the column vectors of the forward transform matrix in Section 2.5. We then introduce the loss function with the adaptive LASSO penalty in Section 2.5.1. Finally, we find the minimum of the loss function in Section 2.5.2 with the FISTA algorithm (Beck & Teboulle 2009).

The l² norm of the ith column vector of A weighted by the inverse of the noise covariance matrix is defined as

$$\begin{eqnarray}&&{{ \mathcal N }}_{i}={\left({{\rm{\Sigma }}}^{-1}\right)}_{\alpha \beta }{A}_{i\alpha }{A}_{i\beta }.\end{eqnarray} \tag{ 25 }$$

The column vectors have different weighted l² norms. Since the gradient descent algorithm, which will be used to solve Equation (24) in Section 2.5.2, takes each basis vector equally, we normalize the column vectors before performing the density map reconstruction to boost the convergence speed of the gradient descent iteration. The normalized forward transform matrix and projection parameters are given by

$$\begin{eqnarray}\begin{array}{rcl}{A}_{\alpha \beta }^{{\prime} } & = & {A}_{\alpha \beta }/{{ \mathcal N }}_{\alpha }^{\tfrac{1}{2}},\\ {x}_{\beta }^{{\prime} } & = & {x}_{\beta }{{ \mathcal N }}_{\beta }^{\tfrac{1}{2}}.\end{array}\end{eqnarray} \tag{ 26 }$$

2.5.1. Adaptive LASSO

The LASSO algorithm uses the l¹ norm of the projection coefficient vector to regularize the modeling. The estimator is defined as

$$\begin{eqnarray}&&{\hat{x^{\prime} }}^{\mathrm{LASSO}}=\mathop{\arg \,\min }\limits_{x}\left\{\displaystyle \frac{1}{2}{}_{{\rm{\Sigma }}}\parallel (\gamma -{\boldsymbol{A}}^{\prime} \cdot x^{\prime} ){\parallel }_{2}^{2}+\lambda \parallel x^{\prime} {\parallel }_{1}^{1}\right\},\end{eqnarray} \tag{ 27 }$$

where $\parallel \cdot {\parallel }_{1}^{1}$ and $\parallel \cdot {\parallel }_{2}^{2}$ refer to the l¹ norm and l² norm, respectively, and λ is the penalty parameter for the LASSO estimation.

The LASSO algorithm searches and selects the parameters that are relevant to the measurements, and simultaneously estimates the values of the selected parameters. It has been shown by Zou (2006) that when the column vectors of the forward transform matrix A ' are highly correlated, the algorithm cannot select the relevant atoms from the dictionary consistently. In addition, the estimated parameters are often biased owing to the shrinkage in the LASSO regression. We note that for the density map reconstruction problem here, the column vectors are highly correlated even in the absence of photo-z uncertainties since the lensing kernels for lenses at different redshifts overlap significantly, i.e., they are highly correlated as shown in Figure 3. Therefore, the LASSO algorithm cannot precisely determine the consistent mass distribution in redshift, and the reconstructed map suffers from smearing in the line-of-sight direction even in the absence of noise.

Figure 7 shows an example reconstruction results for a single halo with mass M₂₀₀ = 10¹⁵ h⁻¹ M_⊙ at redshift 0.35. ¹⁰ The shear measurement and photo-z uncertainties are not included in this simulation. We find significant smearing of the mass distribution with the LASSO algorithm (left panel of Figure 7).

Zoom In Zoom Out Reset image size

To overcome the problem, Zou (2006) proposes the adaptive LASSO algorithm, which uses adaptive weights to penalize different projection coefficients in the l¹ penalty. The adaptive LASSO algorithm performs a two-step process. In the first step, the standard (nonadaptive) LASSO is used to estimate the parameters. Let us denote the preliminary estimation as ${\hat{x^{\prime} }}_{\mathrm{LASSO}}$. In the second step, the preliminary estimate is used to calculate the nonnegative weight vector for penalization as

$$\begin{eqnarray}&&\hat{w}=\displaystyle \frac{1}{| {\hat{x^{\prime} }}_{\mathrm{LASSO}}{| }^{\tau }},\end{eqnarray} \tag{ 28 }$$

where we set the hyperparameter τ to 2 . The adaptive LASSO estimator is then given by

$$\begin{eqnarray}&&\hat{x^{\prime} }=\mathop{\arg \,\min }\limits_{x^{\prime} }\left\{\displaystyle \frac{1}{2}{}_{{\rm{\Sigma }}}\parallel (\gamma -{\boldsymbol{A}}^{\prime} \cdot x^{\prime} ){\parallel }_{2}^{2}+{\lambda }_{\mathrm{ada}}\parallel \hat{w}\circ x^{\prime} {\parallel }_{1}^{1}\right\},\end{eqnarray} \tag{ 29 }$$

where "◦" refers to the element-wise product. λ_ada is the penalty parameter for the adaptive LASSO, which does not need to be the same as the penalty parameter for the preliminary LASSO estimation λ. The adaptive weights enhance the shrinkage in the soft thresholding for the coefficients with smaller amplitudes, whereas the weights suppress the shrinkage for the coefficients with larger amplitudes.

To simplify the equations in the following, we rewrite the loss function with the Einstein notation as

$$\begin{eqnarray}\begin{array}{rcl}L(x^{\prime} ) & = & \displaystyle \frac{1}{2}{\left({{\rm{\Sigma }}}^{-1}\right)}_{\alpha \beta }\left({\gamma }_{\alpha }^{* }-{A}_{\alpha i}^{{\prime} * }{x}_{i}^{{\prime} }\right)({\gamma }_{\beta }-{A}_{\beta j}^{{\prime} }{x}_{j}^{{\prime} })\\ & & +\,{\lambda }_{\mathrm{ada}}\hat{{w}_{\beta }}| {x}_{\beta }^{{\prime} }| ,\end{array}\end{eqnarray} \tag{ 30 }$$

and the quadruple term in the loss function as

$$\begin{eqnarray}&&G(x^{\prime} )=\displaystyle \frac{1}{2}{{\rm{\Sigma }}}_{\alpha \beta }^{-1}({\gamma }_{\alpha }^{* }-{A}_{\alpha i}^{{\prime} * }{x}_{i}^{{\prime} })({\gamma }_{\beta }-{A}_{\beta j}^{{\prime} }{x}_{j}^{{\prime} }).\end{eqnarray} \tag{ 31 }$$

2.5.2. FISTA

Beck & Teboulle (2009) propose the Fast Iterative Soft Thresholding Algorithm (FISTA) to solve the LASSO problem. Since the loss functions of the LASSO and the adaptive LASSO differ only in their penalty terms, FISTA is also applicable to solve the adaptive LASSO problem in a straightforward manner. In this paper, we apply FISTA to solve both the preliminary LASSO estimation and the adaptive LASSO estimation.

We first explain the LASSO preliminary estimation. The coefficients are initialized as ${x}_{i}^{\left(1\right)}=0$. According to FISTA, we iteratively update the projection coefficient vector x. Taking the nth iteration as an example, a temporary update is first calculated as

$$\begin{eqnarray}&&{x}_{i}^{{\prime} \left(n+1\right)}={\mathrm{ST}}_{\lambda }\left({x}_{i}^{{\prime} \left(n\right)}-\mu {\partial }_{i}G({x}^{{\prime} \left(n\right)})\right),\end{eqnarray} \tag{ 32 }$$

where ST is the soft-thresholding function defined as

$$\begin{eqnarray}&&{\mathrm{ST}}_{\lambda }\left(x^{\prime} \right)=\mathrm{sign}(x^{\prime} )\max \left\{| x^{\prime} | -\lambda ,0\right\}.\end{eqnarray} \tag{ 33 }$$

The soft thresholding is a part of the LASSO algorithm, which selects the modes with amplitude greater than λ, and shrinks the selected estimation by λ. The coefficient μ is the step size of the gradient descent iteration, and ${\partial }_{i}G({x}^{{\prime} \left(n\right)})$ is the ith element of the gradient vector of G at point ${x}^{{\prime} \left(n\right)}$:

$$\begin{eqnarray}&&{\partial }_{i}G({x}^{{\prime} \left(n\right)})={{\rm{\Sigma }}}_{\alpha \beta }^{-1}{\mathfrak{R}}\left\{{A}_{\alpha i}^{{\prime} * }({\gamma }_{\beta }-{A}_{\beta j}^{{\prime} }{x}_{j}^{{\prime} })\right\},\end{eqnarray} \tag{ 34 }$$

where ${\mathfrak{R}}\left\{\bullet \right\}$ returns the real part of the input vector. FISTA requires an additional update using a weighted average between ${x}^{{\prime} \left(n+1\right)}$ and ${x}^{{\prime} \left(n\right)}$ as

$$\begin{eqnarray}\begin{array}{rcl}{t}^{\left(n+1\right)} & = & \displaystyle \frac{1+\sqrt{1+4{\left({t}^{\left(n\right)}\right)}^{2}}}{2},\\ {x}^{{\prime} \left(n+1\right)} & \leftarrow & {x}^{{\prime} \left(x+1\right)}+\displaystyle \frac{{t}^{\left(n\right)}-1}{{t}^{\left(n+1\right)}}({x}^{{\prime} \left(n+1\right)}-{x}^{{\prime} \left(n\right)}),\end{array}\end{eqnarray} \tag{ 35 }$$

where the relative weight is initialized as ${t}^{\left(1\right)}=1$.

FISTA converges as long as the gradient descent step size μ satisfies

$$\begin{eqnarray}&&0\lt \mu \lt \displaystyle \frac{1}{\parallel {{\boldsymbol{A}}}^{\dagger }\cdot {{\boldsymbol{\Sigma }}}^{-1}\cdot {\boldsymbol{A}}\parallel },\end{eqnarray} \tag{ 36 }$$

where the operation ∥•∥ returns the spectrum norm of the input matrix. The spectral norm is estimated by simulating a large number of random vectors with l² norms equal to 1 with different realizations. The matrix operator A ^† · Σ⁻¹ · A is subsequently applied to each random vector to yield a corresponding transformed vector. The spectral norm is determined as the maximum l² norm of the transformed vectors.

As summarized below, we first initialize the projection coefficients as zero, and use FISTA to find the global minimum of the LASSO loss function. Thus obtained global minimum is the preliminary one. We then use the preliminary LASSO estimation to weight the coefficients and calculate the adaptive LASSO loss function. Finally, we set the preliminary LASSO estimation as the "warm" start of the adaptive LASSO estimation, and FISTA again to find the global minimum of the adaptive LASSO loss function, which is our final solution.

Algorithm 1. Our Algorithm

Input γ: pixelized complex 3D shear field

Output δ: 3D array of density contrast

1: Normalize column vectors of A

2: Estimate step size μ and ${\boldsymbol{\Sigma }}$

3: Initialization:

4: ${x}^{{\prime} \left(1\right)}=0$

5: $\hat{w}=1$

6: ${t}^{\left(1\right)}=1$,i = 1, j = 1

7: while $j\leqslant 2$ do

8: while $i\leqslant {N}_{\mathrm{iter}}$ do

9: # soft thresholding

10: ${x}_{i}^{{\prime} \left(n+1\right)}={\mathrm{ST}}_{\hat{w}\lambda }\left({x}_{i}^{{\prime} \left(n\right)}-\mu {\partial }_{i}G({x}^{{\prime} \left(n\right)})\right)$

11: # FISTA algorithm

12: ${t}^{\left(n+1\right)}=\tfrac{1+\sqrt{1\,+\,4{\left({t}^{\left(n\right)}\right)}^{2}}}{2}$

13: ${x}^{{\prime} \left(n+1\right)}\leftarrow {x}^{{\prime} \left(x+1\right)}+\tfrac{{t}^{\left(n\right)}-1}{{t}^{\left(n+1\right)}}({x}^{{\prime} \left(n+1\right)}-{x}^{{\prime} \left(n\right)})$

14: $i=i+1$

15: end while

16: Reinitialization:

17: $\hat{w}=| {\hat{x^{\prime} }}_{\mathrm{LASSO}}{| }^{-2}$, $\lambda \leftarrow {\lambda }_{\mathrm{ada}}$

18: ${\hat{x^{\prime} }}^{\left(1\right)}={x}^{{\prime} \left({N}_{\mathrm{iter}}\right)}$

19: ${t}^{\left(1\right)}=1$, i = 1

20: $j=j+1$

21: end while

22: $\delta ={\boldsymbol{\Phi }}{{ \mathcal N }}^{-\tfrac{1}{2}}{x}^{{\prime} \left({N}_{\mathrm{iter}}\right)}$

Download table as:  ASCIITypeset image

We use halos with a variety of masses and redshifts in the range 10¹⁴ h⁻¹ M_⊙ < M < 10¹⁵ h⁻¹ M_⊙, and 0.05 < z < 0.85, respectively. We divide the parameter space into eight redshift bins and eight mass bins with equal separation. We randomly shift the input halo redshift and halo mass from the bin center by a small amount in order to avoid repeatedly sampling at the exact same halo mass and redshift.

The concentration parameter c_h of a NFW halo is determined as a function of the halo mass (M₂₀₀) and redshift (z_h ) according to Ragagnin et al. (2019):

$$\begin{eqnarray}&&{c}_{{\rm{h}}}=6.02\times {\left(\displaystyle \frac{{M}_{200}}{{10}^{13}{M}_{\odot }}\right)}^{-0.12}{\left(\displaystyle \frac{1.47}{1.+{z}_{{\rm{h}}}}\right)}^{0.16}.\end{eqnarray} \tag{ 37 }$$

The weak-lensing-shear fields of the NFW halos are calculated according to Takada & Jain (2003). The shear distortions are applied to 100 realizations of galaxy catalogs with the HSC-like shear measurement error and photo-z uncertainty.

The mock galaxy catalogs are generated using the HSC S16A shear catalog (Mandelbaum et al. 2018). We use the galaxies in a 1 square degree field at the center of tract 9347 (Aihara et al. 2018). The average galaxy number density in this region is $22.94\,{\mathrm{arcmin}}^{-2}$. The positions of galaxies are randomized to distribute homogeneously in the 1 square degree stamp. We randomly assign redshift for each galaxy following the MLZ photo-z PDF (Tanaka et al. 2018).

We simulate the HSC-like shear measurement error due to shape noise and sky variance with different realizations by randomly rotating the galaxies in the shear catalog. The shear measurement errors on an individual galaxy level and an individual pixel level are demonstrated in Figure 6. The standard deviation map of the noise is demonstrated in Figure 8.

Zoom In Zoom Out Reset image size

In this subsection, we test the performance of our algorithm by adopting models where the matter density field is represented by multiscale NFW atoms. The dictionary is constructed with three frames of different NFW scale radii in comoving coordinates: 0.12 h⁻¹ Mpc, 0.24 h⁻¹ Mpc, and 0.36 h⁻¹ Mpc. The truncation radii are set to four times the scale radii for the atoms in the dictionary, i.e., we assume c_h = 4. Note that each frame of our dictionary fixes the scale radius in comoving coordinates and thus the NFW atoms have different angular sizes when placed at different redshift.

We test the algorithm by varying the penalty parameter for the LASSO estimation with λ = 3.5, 4.0, and 5.0. The corresponding penalty parameters for the final adaptive LASSO estimations are set to λ_ad = λ^τ+1. Here, both the LASSO estimation and our adaptive LASSO estimation select the pixels with S/Ns greater than λ in each gradient descent iteration, and the local density is estimated for the selected pixels with a shrinkage of the estimation amplitude. The LASSO estimation shrinks the density amplitudes by λ for every selected pixel, whereas the adaptive LASSO estimation suppresses the shrinkage if the preliminary estimation for the pixel is greater than λ, by down-weighting their penalties, and otherwise it enhances the shrinkage.

We would like to reconstruct with the resolution limit set by the Gaussian-smoothing kernel with a standard deviation of 15 and the 1' pixel scale as described in Section 2.3.2. In addition, we smooth the reconstructed density field with the same Gaussian kernel in each lens redshift plane.

Figure 9 shows the 3D density maps reconstructed with different penalty parameters for a halo with M₂₀₀ = 10^15.02 h⁻¹ M_⊙ at redshift 0.164. The corresponding pixel-value histograms are shown in Figure 9. Our adaptive LASSO algorithm assigns zero to a fraction of the reconstructed pixels that are mostly noises, while retaining strong signals. The right panels demonstrate a case where almost all of the pixels with pure noise are set to zero when a large penalty parameter is adopted. It is important to note that the reconstructed density maps are not compromised by line-of-sight smearing, i.e., the reconstructed lens is localized in redshift.

Zoom In Zoom Out Reset image size

Following Lanusse et al. (2016), we normalize the detected peaks in the lth (l = 1...20) lens redshift plane to account for the peak amplitude difference arising from the difference in the norm of the lensing kernels in different redshift bins:

$$\begin{eqnarray}&&{\delta }_{\mathrm{peak}}^{n}(\vec{\theta },{z}_{l})={\delta }_{\mathrm{peak}}(\vec{\theta },{z}_{l})/{{ \mathcal R }}_{l}^{\tfrac{1}{2}},\end{eqnarray} \tag{ 38 }$$

where the normalization matrix is defined as

$$\begin{eqnarray}&&{{ \mathcal R }}_{l}=\sum _{s}{K}^{2}({z}_{l},{z}_{s}),\end{eqnarray} \tag{ 39 }$$

where K(z_l , z_s ) is the lensing kernel. In Figure 10, we show the histograms of the normalized peaks with different penalty parameters. There, we stack the histograms from 100 realizations of all the halos sampled in the redshift–mass plane. In addition, we generate 1000 realizations of pure noise catalogs and perform the reconstruction using the noise catalogs in order to examine the noise properties. The histograms of the normalized peaks detected from the pure noise catalogs are shown in Figure 10 along with the best-fit Gaussian functions of the noise-peak histograms.

Zoom In Zoom Out Reset image size

We find that the number counts including both true and false peaks decrease as the penalty parameter λ increases. In addition, the standard deviation of noise peaks decreases as λ increases. As a result, with λ = 5.0, we find a clearer excess in the positive peak counts compared with the noise-peak histograms, especially at the high density contrast. This is expected because a higher penalty parameter prefers a sparser solution; more peaks originating from the noise are removed than those from real clusters at the high density contrast.

We demonstrate the position offsets of the detected halos compared to the input positions in Figure 11. The left panel shows the 2D number histograms stacked over all the simulations as a function of the angular distance and the redshift offsets. We see clear clustering of the peaks close to the position of the input halo on the stacked number histogram. For each simulation, we identify positive peaks closest to the input position (in pixels). If a closest peak is located inside the region denoted with the dashed box in the left panel of Figure 11, we regard it as a "true" peak detection. Other identified peaks, which include both positive and negative peaks, are judged to be false detection.

Zoom In Zoom Out Reset image size

The right panel of Figure 11 shows the estimated redshift of the "true" detections averaged over the simulations (with different noise realizations) for each halo as a function of the input redshift. As shown, the estimated redshifts are slightly lower than the ground truth by Δz ∼ 0.03 for halos in the low-redshift range (z ≤ 0.4). For halos at 0.4 < z ≤ 0.85, the relative redshift bias is below 0.5%. The standard deviation of the redshift estimation is 0.092.

In order to reduce the number of false detections, we select peaks as galaxy clusters if the peak values are greater than an empirically determined threshold. The threshold is set in units of the standard deviation of the noise peaks. We use different detection thresholds (1.5 and 3.0σ) to detect galaxy clusters from the mass maps reconstructed with different penalty parameters. Figures 12 and 13 show the detection rates for halos in the redshift–mass plane for different penalty parameters (λ = 3.5 and λ = 5.0). The corresponding number of false detections per square degree as a function of detection threshold are also demonstrated in the figures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As shown, the false peak density is successfully reduced for relatively large detection threshold, but the detection rate of the halo also decreases. After a few experiments, we have decided to set the detection threshold to 1.5σ and set the penalty parameter λ to 5.0 since the combination suppresses the false detection to 0.022 while keeping a high halo detection rate. In summary, our method is able to detect halos with minimal mass of 10^14.0 h⁻¹ M_⊙, 10^14.7 h⁻¹ M_⊙, and 10^15.0 h⁻¹ M_⊙ for the low (z < 0.3), median (0.3 ≤ z < 0.6), and high (0.6 ≤ z < 0.85) redshift ranges, respectively.

Using the detection rate measured from our simulations, we are able to predict the number density of detected clusters by assuming the halo mass function of Tinker et al. (2008). We use HMF (Murray et al. 2013), an open-source package, ¹¹ to calculate the halo mass function. The predicted halo detection number density for the λ = 5 and 1.5σ detection threshold that was set up is shown in Figure 14. The resulting cluster number density is 0.49 deg⁻², which is much higher than the false detection rate of 0.022 deg⁻². This cluster number density corresponds to 78.4 clusters for the first-year HSC shear catalog (Mandelbaum et al. 2018) with a survey area of ∼ 160 deg². The expected number detected is slightly higher than the number of 2D cluster detections (63 detected clusters) for the first-year HSC shear catalog (Miyazaki et al. 2018). Furthermore, our 3D detection method provides an accurate redshift estimation for individual clusters. In contrast, the redshift information is not provided from the 2D mass-map reconstruction.

Zoom In Zoom Out Reset image size

We perform an additional test by substituting the default NFW dictionary with point mass. This test may indicate some certain limitations of our method when applied to a case with very compact (point) objects, which, however, is unlikely to happen in actual lensing observations.

We set the penalty parameter for the preliminary LASSO to λ = 3.5 and 5.0. Pramanik & Zhang (2020) propose to incorporate group information into different adaptive LASSO penalty weights by setting the weights for projection coefficients in a group to the average of the adaptive weights in the same group. We assume that the neighboring pixels in the same redshift plane belong to the same structural group (e.g., galaxy cluster and void), and smooth the amplitude of the preliminary LASSO estimation in each lens redshift plane with a top-hat filter of comoving diameter r_c = 0.25 h⁻¹ Mpc. Let us denote the amplitude of the preliminary LASSO estimation with the smoothing as $| {\hat{x}}^{\mathrm{LASSO}}{| }_{\mathrm{sm}}$, and adopt the penalty weights given by $\hat{w}=1/| {\hat{x}}^{\mathrm{LASSO}}{| }_{\mathrm{sm}}^{\tau }$.

Figure 15 shows the reconstruction result with the point-mass dictionary. Interestingly, several "discrete" masses are assigned to at different redshift bins in the neighboring region of the true halo center. In contrast, as we have seen in Figure 9, the NFW dictionary manages to recover a consistent mass distribution. The problem of the point-mass dictionary originates from the fact that the profile of the point-mass atom in the transverse plane is much more compact than the profile of the input halo, especially when placed at low redshifts.

Zoom In Zoom Out Reset image size