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”:

where $\mathbf{\Phi}$ is the matrix operator transforming from the projection coefficient vector $x$ to the density contrast field $\delta$. The column vectors of $\mathbf{\Phi}$ are the basis “atoms” of the model dictionary and denoted as $\phi_{s}$.

There have been a few studies that adopt different dictionaries for $3$D 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 $\delta$ 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

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

$f=1/[\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

where $\chi(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

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

and write the dictionary transform operator into a row vector:

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

The $2$D profiles of the NFW atoms and the point-mass atom on the transverse plane are shown in Figure 1. The corresponding $1$D 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.

We define the forward transform operator $\mathbf{A}=\mathbf{P}\cdot\mathbf{Q}\cdot\mathbf{\Phi}$ where $\mathbf{P}$ represents systematic effects and $\mathbf{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

The lensing transform operator can be expressed as

$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

where $E(z)$ is the Hubble parameter as a function of redshift, in units of $H_{0}$.

is the Kaiser-Squares kernel (Kaiser & Squires, 1993), which decays proportional to the inverse-square of the distance. Here $\theta_{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 $\mathbf{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.

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 $\mathbf{R}$ (photometric redshift uncertainties), $\mathbf{W}$ (smoothing), and $\mathbf{M}$ (masking).

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

where ${}_{\Sigma}\mathinner{\!\left\lVert(\gamma-\mathbf{A}x)\right\rVert}_{2}^{2}$ is the $l^{2}$ norm of residuals weighted by the inverse of the covariance matrix $\Sigma$ 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 $\lambda$ 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 $3$D weak-lensing map reconstructions: (i) Simon et al. (2009) propose to use the Wiener filter, which is also known as $l^{2}$ ridge penalty ($C=\mathinner{\!\left\lVert x\right\rVert}^{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 Hyper Suprime-Cam Survey (Aihara et al., 2018). It is found that the density maps reconstructed by the method suffer from significant line of sight smearing with standard deviation of $\sigma_{z}=0.2-0.3$. (ii) Leonard et al. (2014) use GLIMPSE algorithm, which adopts a derivative version of the $l^{1}$ LASSO penalty ($C=\mathinner{\!\left\lVert x\right\rVert}^{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^{2}$ norm of the $i$th column vector of $\mathbf{A}$ weighted by the inverse of the noise covariance matrix is defined as

The column vectors have different weighted $l^{2}$ norms. Since the gradient descent algorithm, which will be used to solve Equation (24) in Section 2.5.2, takes each column 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

We use halos with a variety of masses and redshifts in the range $10^{14}\text{h}^{-1}M_{\odot}<M<10^{15}\text{h}^{-1}M_{\odot}$, 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_{\rm h}$ of a NFW halo is determined as a function of the halo mass ($M_{200}$) and redshift ($z_{\rm h}$) according to Ragagnin et al. (2019):

The weak-lensing shear fields of the NFW halos are calculated according to Takada & Jain (2003). The shear distortions are applied to one hundred 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~{}\rm{arcmin}^{-2}$. The positions of galaxies are randomized to distribute homogeneously in the one-square degree stamp. We randomly assign redshift for each galaxy following the MLZ photo-$z$ probability distribution function (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 error on individual galaxy level and individual pixel level are demonstrated in Figure 6. The standard deviation map of the noise is demonstrated in Figure 8.

In this section, 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 coordinate: $0.12~{}\text{h}^{-1}$ Mpc, $0.24~{}\text{h}^{-1}$ Mpc, and $0.36~{}\text{h}^{-1}$ Mpc. The truncation radii are set to four times the scale radii for the atoms in the dictionary, i.e., we assume $c_{\rm 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 with varying the penalty parameter for the LASSO estimation with $\lambda=3.5$, $4.0$, and $5.0$. The corresponding penalty parameters for the final adaptive LASSO estimations are set to $\lambda_{\rm{ad}}=\lambda^{\tau+1}$. Here, both the LASSO estimation and our adaptive LASSO estimation select the pixels with signal-to-noise ratios greater than $\lambda$ 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 $\lambda$ for every selected pixels, whereas the adaptive LASSO estimation suppresses the shrinkage if the preliminary estimation for the pixel is greater than $\lambda$, by down-weighting their penalties, and otherwise it enhances the shrinkage.

We’d like to reconstruct with the resolution limit set by the Gaussian smoothing kernel with a standard deviation of $1\farcm 5$ and the $1\arcmin$ pixel scale as described in Sections 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 $3$D density maps reconstructed with different penalty parameters for a halo with $M_{200}=10^{15.02}~{}\text{h}^{-1}M_{\odot}$ 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 parameteris 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.

Following Leonard et al. (2014), we normalize the detected peaks in the $l$th ($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:

where the normalization matrix is defined as

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.

We find that the number counts including both true and false peaks decrease as the penalty parameter $\lambda$ increases. In addition, the standard deviation of noise peaks decreases as $\lambda$ increases. As a result, with $\lambda=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 $2$D 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 the pixel unit). 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.

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 $\Delta z\sim 0.03$ for halos in the low-redshift range ($z\leq 0.4$). For halos at $0.4<z\leq 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\sigma$ and $3.0\sigma$) 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 ($\lambda=3.5$ and $\lambda=5.0$). The corresponding number of false detections per square degree as a function of detection threshold are also demonstrated in the figures.

As shown, the false peak density is successfully reduced for relatively large detection threshold, but the detection rate of halo also decreases. After a few experiments, we have decided to set the detection threshold to $1.5\sigma$ and set the penalty parameter $\lambda$ 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}\text{h}^{-1}M_{\odot}$, $10^{14.7}\text{h}^{-1}M_{\odot}$, $10^{15.0}\text{h}^{-1}M_{\odot}$ for the low ($z<0.3$), median ($0.3\leq z<0.6$) and high ($0.6\leq 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⁴⁴4https://github.com/halomod/hmf, to calculate the halo mass function. The predicted halo detection number density for the setup $\lambda=5$ and $1.5\sigma$ detection threshold is shown in Figure 14. The resulting cluster number density is $0.49$ $\deg^{-2}$, which is much higher than the false detection rate of 0.022 deg^-2. 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 $\sim 160\deg^{2}$. The expected number of detection is slightly higher than the number of $2$D cluster detections ($63$ detected clusters) for the first year HSC shear catalog (Miyazaki et al., 2018b). Furthermore, our $3$D detection method provides an accurate redshift estimation for individual clusters. In contrast, the redshift information is not provided from the $2$D mass map reconstruction.

We perform an additional test by substituting the default NFW dictionary with point-mass. This test may indicate some certain limitation 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 $\lambda=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_{\rm c}=0.25~{}\text{h}^{-1}$ Mpc. Let us denote the amplitude of the preliminary LASSO estimation with the smoothing as $\mathinner{\!\left\lvert\hat{x}^{\rm{LASSO}}\right\rvert}_{\rm{sm}}$, and adopt the penalty weights given by $\hat{w}=1/\mathinner{\!\left\lvert\hat{x}^{\rm{LASSO}}\right\rvert}_{\rm{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.