Data simulation¶

Cosmology¶

This work follows in the general footsteps of data analysis of weak lensing surveys like HSC Hikage et al. (2019) and KiDS Köhlinger et al. (2017)Asgari, Marika et al. (2021), as we aim to replicate their methodologies. Throughout the work, we assume a fiducial cosmology of fixed parameter values as shown in Tab. %s. In particular, all data maps will be generated from a power spectrum following this cosmology. We will refer to this power spectrum as the fiducial power spectrum $C(l)$. A leading modelling choice comes with the redshift distribution $n(z)$. We model it as a Smail-type distribution Smail et al. (1995)Kirk et al. (2012),

The choice of parameters has been made to emulate bin 5 of the KiDS1000 survey Hildebrandt, H. et al. (2021), with parameters $\alpha=3.5$, $\beta=4.5$, $z_0=1$.

Map making pipeline¶

When dealing with maps of finite size $(L, L)$ and pixel resolution $(N,N)$, Fourier space is going to have boundaries, just as real space does. These limits are given by

Therefore, a map of size $(10^{\circ},10^{\circ})$ and a grid $64\times64$ will have limits $(l_{\text{min}}, l_{\text{max}}) = (36, 2304)\text{rad}^{-1}$. Now that we have the physical range for the power spectrum, we can generate data.

To generate a GRF, we employ the algorithm mentioned in Eq. (13). However, in making a lognormal field, the matter is a bit more complicated. Any GRF on which we apply the lognormal transformation $\mathcal{L}_\kappa$, from Eq. eq:L_k, becomes a lognormal field. The reason we do not just lognormal transform any field is given by the fact that they would not recover the fiducial power spectrum $C(l)$. Our goal is therefore to find a transfer GRF $\tilde{\kappa}^G$ which, when transformed lognormally, gives rise to a lognormal field $\kappa^L$ that recovers $C(l)$. To do so, we follow the sequence of transformations illustrated in the top branch of Fig. 6.

Noise and mask¶

As all data is being simulated, we only take into account the so-called shape noise, which is due to the intrinsic distribution of ellipticities and angle formed with respect to us. GPs will treat each pixel of the map as a random variable Gaussianly distributed with standard deviation given by Croft et al. (2017),

We use values of $\sigma_e=0.26$, $n_g= 4, 10, 30, 100 \text{ galaxies}/\text{arcmin}^2$, and a pixel area given by the pixel resolution squared, $A_{px}=(L/N)^2$. For our highest-resolution run, we have $N=64$ and use $n_g = 10 \text{ galaxies}/\text{arcmin}^2$, which results in a noise standard deviation of $\sigma_\text{noise} \sim 0.0088$. The mask being used as seen in Fig. 7, covers approximately $10\%$ of the patch Grewal et al. (2024). We keep the noise and mask random seeds fixed throughout the work.

Kernel¶

The kernel of a Gaussian process is given by a function of the form $k(x,y)$. It takes in two points, $x$ and $y$, and returns the value of their correlation. In our case, specifically, $x$ and $y$ will be two points in a grid of shape $(N,N)$, and the kernel function will be the convergence angular autocorrelation function. As all code used for this paper is written in JAX Bradbury et al. (2018) we opt for the use of the library tinygp Foreman-Mackey et al. (2024) for all Gaussian process computations. tinygp allows for the use of a custom kernel with a custom evaluate method, which takes two points on the grid and returns the correlation value. We follow with two Python pseudocodes of our kernel implementations.

class kernel_Hankel:
    def __init__(self, cl, N, L):
        self.w   = Hankel(cl, N, L)
        self.r   = L / N
        
    def evaluate(self, x, y):
        theta    = self.r * sqrt(sum(x - y))
        return self.w(theta)

class kernel_FFT:
    def __init__(self, cl2D, N, L):
        M        = sqrt(cl2D.size)
        self.w2D = abs(ifft2(cl2D)) * M**2 / L**2
        self.r   = M / N
        
    def evaluate(self, x, y):
        d0       = self.r * abs(x[0] - y[0])
        d1       = self.r * abs(x[1] - y[1])
        return self.w2D[int(d0)][int(d1)]

Hankel transform¶

In order to build a Gaussian process kernel, we need to find the correlation function that best describes the data. The way we do this is by computing the angular power spectrum and transforming it in the corresponding angular correlation function. Let’s explore the previously discussed flat-sky relation between angular power spectrum and correlation function given by Eq. (8). This particular integral of a Bessel function is also known as a zeroth order Hankel transformation,

We will explore two methods for computing this integral: the integration method Sec. ?? and the FFTlog method Sec. ??.

Integration¶

Integration is the most straight forward way to evaluate the integral, but it requires to implement an algorithm for the approximation of the Bessel function $J_0$. The advantages of this method:

• it is easy to integrate over the correct $l$-range, from $l_{\text{min}}$ and $l_{\text{max}}$, a freedom that we do not have with the FFTlog.

The disadvantages:

• integration is computationally slow, especially when dealing with highly oscillatory behaviour introduced by the Bessel function, which requires fine sampling.

FFTlog¶

The FFTlog method Simonović et al. (2018) is a fast implementation of the Hankel transformation. In fact it simplifies Eq. (23) by assuming a power decomposition of $C(l)$. Such a decomposition is achievable by taking the fast Fourier transformation (FFT) in $\log k$, hence the name FFTlog. The power spectrum becomes,

Substituting in Eq. (23),

Take $x=l\theta$ and $s_\alpha-1=\upsilon + i\eta_\alpha+1$,

Lastly we recognise $\int_0^\infty dx x^{s-1}J_0(x) = \frac{2^{s-1}}{\pi}\sin(\pi s/2) [\Gamma(s/2)]^2$, a Mellin transform. Using this tabulated result we conclude that the correlation is given by the sum,

For the Mellin transform identity to hold analytically, the integral bounds have to go from 0 to $\infty$. Although computationally we don’t need to consider such a wide range, we still have to broaden the integration limits to something larger than our $l$-range; so as to avoid ringing effects. Experimentally we have found that extending the $l$-range between $l_{\text{min}}/4$ and $l_{\text{max}}$, is enough to compensate for such effects. The advantages of this method:

• the FFTlog method is ultimately much faster than the integration method.

The disadvantages:

• the widening of the $l$-range needed to avoid ringing effects inevitably adds power to the correlation function resulting in a higher variance and overall amplitude.

Fourier transform¶

The second transformation that we can use to go from power spectrum to correlation function, is the Fourier transformation discussed in Eq. (7). It is a more fundamental relation than the Hankel transformation as it does not assume that $l$ is radially symmetric. In our case, low resolution, coupled with a square map, means that the radial assumption might not apply. Note that, the lowest resolution grid we use is $32\times32$, with $L/N \sim 19$ arcmin; which is one order of magnitude worse than today’s weak lensing data, around $\sim 3$ arcmin according to Zhou et al. (2023). In the following Sec. ?? we will find experimentally that the Hankel methods described above do not work well with our GP setup. For this reason, we explore this method of conversion between angular power spectrum and correlation function,

As computers can only deal with discrete functions, it is important to note that we will be performing a discrete Fourier transformation (DFT). We will use the widely known FFT algorithm to compute DFTs. In particular we use JAX’s implementation of the 2D FFT algorithm, jax.numpy.fft.ifft2.

Let us find the relation between continuous Fourier transformation and DFT. The definitions of continuous and discrete Fourier transformations in 1D, are respectively:

First of all, a DFT is dimensionless. Secondly, it is discrete and bounded. We can therefore rewrite Eq. (27) using the substitution $k(p)=2\pi \frac{p}{L}$ to discretise k-space,

Applying the definition of DFT as seen in Eq. eq:DFT, it follows that in 1D

Which means that the correlation function is given by a backwards normalised inverse DFT to be scaled by a factor $(N/L)^d$, where $d$ is the dimension of the considered space. In our case, $d=2$.

In order to perform a FFT in 2D, we will need a two dimensional extension of the angular power spectrum. We make use of its radial symmetry with respect to $\bm{l}=(l_x,l_y)$ and create a 2D grid of shape $(M,M)$ as shown in Fig. 8.

Now that we have a 2D power spectrum we can take the inverse two dimensional fast Fourier transformation to obtain a 2D correlation function. In practice, we will test two grids, which we name full-range FFT method and half-range FFT method.

Full-range FFT¶

The full-range FFT is defined by a grid of shape $(2N,2N)$. Since the grid is radial, it will be centered. This implies that if we want to keep information of all $N$ modes,

the grid will have to be at least of shape $M=2N$. The advantages of this method:

• it keeps information on the full range of modes.

The disadvantages:

• it introduces rounding errors, as the field has shape $(N,N)$ and not all of the possible distance combinations of such a grid are covered bya grid of shape $(2N,2N)$.

Half-range FFT¶

The half-range FFT is defined by a grid of shape $(N,N)$. The advantages of this method:

• has the perk of having no shape mismatch between field and correlation function.

The disadvantages:

• it loses half of the $l$-range, missing information on small scales.