A Differentiable Perturbation-based Weak Lensing Shear Estimator

Upcoming imaging surveys will use weak gravitational lensing to study the large-scale structure of the Universe, demanding sub-percent accuracy for precise cosmic shear measurements. We present a new differentiable implementation of our perturbation-based shear estimator (FPFS), using JAX, which is publicly available as part of a new suite of analytic shear algorithms called AnaCal. This code can analytically calibrate the shear response of any nonlinear observable constructed with the FPFS shapelets and detection modes utilizing auto-differentiation (AD), generalizing the formalism to include a family of shear estimators with corrections for detection and selection biases. Using the AD capability of JAX, it calculates the full Hessian matrix of the non-linear observables, which improves the previously presented second-order noise bias correction in the shear estimation. As an illustration of the power of the new AnaCal framework, we optimize the effective galaxy number density in the space of the generalized shear estimators using an LSST-like galaxy image simulation for the ten-year LSST. For the generic shear estimator, the magnitude of the multiplicative bias $|m|$ is below $3\times 10^{-3}$ (99.7% confidence interval), and the effective galaxy number density is improved by 5%. We also discuss some planned future additions to the AnaCal software suite to extend its applicability beyond the FPFS measurements.


INTRODUCTION
Weak gravitational lensing refers to the weak but coherent distortions in the shapes of distant background galaxies caused by the foreground inhomogeneous mass distortion.It is a unique and essential cosmological tool, providing crucial insight into the Universe's most mysterious constituents: dark matter and dark energy.The fundamental principle behind this tool lies in Einstein's general theory of relativity, which posits that the presence of mass can distort the fabric of space-time, causing light to follow curved trajectories (see Bartelmann & Schneider 2001;Refregier 2003a;Massey et al. 2010;Kilbinger 2015 for reviews on weak lensing).
The weak lensing distortion, which is termed "shear", is a measure of the anisotropic stretching of the image of a galaxy due to the intervening mass distribution.Accurate measurement of the weak lensing shear is vital for understanding the structure and evolution of the universe.Due to its subtle nature, this task is inherently complex and is compounded by observational challenges, including the point-spread function (PSF) from telescope optics and atmospheric effects (Paulin-Henriksson et al. 2008;Liaudat et al. 2023); model bias from unrealistic assumptions on galaxy morphology (e.g., Bernstein 2010); bias from image noise due to the non-linearity in shear estimators (e.g., Refregier et al. 2012); bias due to anisotropies in our selected galaxy sample induced by the sample selection (Kaiser   ⋆ xiangchl@andrew.cmu.eduet al. 1995) and detection (Sheldon et al. 2020), and bias arises from blending and deblending of galaxies (Li et al. 2018).We refer interested readers to Mandelbaum (2018) for a comprehensive review of systematics in shear estimation.Despite these difficulties, METADETECTION (Huff & Mandelbaum 2017;Sheldon & Huff 2017;Sheldon et al. 2020), an empirical shear estimator, has been developed to measure shear from blended galaxy images with subpercent accuracy without relying on any calibration from external image simulations.In addition to METADETECTION, a Bayesian method BFD (Bernstein et al. 2016) reaches subpercent accuracy for isolated galaxies, but there is a residual percent-level multiplicative bias for blended galaxies.
In addition to these shear estimation methods, there has been a long tradition of developing purely analytic 1 shear estimators (Zhang 2008;Zhang et al. 2015;Li et al. 2018Li et al. , 2022b;;Li & Mandelbaum 2023).FPFS is a perturbation-based analytic shear estimator without any assumption regarding galaxy morphology (Li et al. 2018).Specifically, Li et al. (2018) constructs an ellipticity estimator with a set of linear observables that is obtained by convolving the image with a set of kernels (e.g., shapelet functions: Refregier 2003b; Massey & Refregier 2005) after PSF deconvolution (Zhang 2008).Li et al. (2018) derive the linear shear responses of the linear 1 By analytic here, we mean that all the formulae used to measure and calibrate the shear are derived directly from first principles and do not involve any empirical steps.
observables using the shear responses of these kernels to obtain the shear response of the ellipticity.Li et al. (2022b) corrects for secondorder noise bias using the second-order derivatives (Hessian matrix) of the ellipticity with respect to the linear observables and the covariance matrix of the measurement error on the linear observables induced by image noise.Finally, Li & Mandelbaum (2023) interpreted the detection and selection processes as a weighting of image pixels after convolution with a set of kernels, and corrected for the detection bias in the FPFS shear estimator.Li & Mandelbaum (2023) showed that the shear estimation bias in FPFS is less than 0.5% of the shear signal in the presence of blending.The FPFS shear estimator is more than a hundred times faster than the widely used METADETECTION (Sheldon et al. 2020) due to its analytic features2 .
However, the formalism in that work neglected a few relatively small elements in the Hessian matrix to simplify the analytic derivation of the noise bias correction.This leaves a residual multiplicative bias of about -0.4 percent relative to the true shear for faint galaxies.In addition, it is too complicated to analytically extend the ellipticity and detection / selection weights to a general form and derive the shear response and noise bias for them.These limitations motivate our work in this paper.
JAX (Bradbury et al. 2018) is a powerful computational framework developed by Google that aims to offer a new perspective on numerical computing.It extends the ability of the Python programming language to automatic differentiation of Python and NumPy functions.To be more specific, we can use a JAX-based framework to compute derivatives of the functions with respect to the input variables.In this paper, we put the FPFS shear estimator into the JAX framework, and use the auto-differentiation (AD) feature to improve the noise bias correction in the FPFS shear estimator.In addition, we generalize the FPFS shear estimator to a more generic form and improve its precision using the AD capabilities of JAX.In order to provide a simple interface through one repository that enables users to seamlessly update as new versions of analytic shear estimators become available, we introduce a novel analytical framework for shear calibration, termed AnaCal3 .This system is dedicated to improving weak lensing shear calibration by creating easy-to-use tools that accurately turn measured galaxy shapes into shear distortions.For better manageability, modularity, and scalability, we have organized our codebase into two separate repositories.The first repository handles the analysis of galaxy images to generate shape catalogs and is accessible at FPFS GitHub Repository4 .The second repository is dedicated to calculating shear based on these shape catalogs, and can be found at the impt GitHub Repository5 .However, users of this method can simply interact with the AnaCal repository.
In this paper, we first introduce the FPFS shear estimator in section 2. After that, we present two applications of our shear estimator in section 3: first, improving the noise bias correction and second, using a more generic shear estimator to improve precision.Finally, we summarize these results and the future outlook in section 4.

FPFS shear estimator
In this subsection, we briefly review the FPFS shear estimator developed in Li et al. (2018Li et al. ( , 2022b)); Li & Mandelbaum (2023).A weak lensing shear distortion caused by a foreground inhomogeneous mass distribution changes the light profiles of distant background galaxies.Here is the Jacobian matrix of the mapping from the lensed sky to the true sky.The parameters (g 1 , g 2 ) are employed to quantify shear distortion within the framework of the celestial coordinate system.Specifically, g 1 accounts for the elongation of the image in the horizontal direction, while g 2 is responsible for the stretching of the image along an axis inclined at a 45-degree angle to the horizontal.In this paper, we set the lensing convergence (Bartelmann & Schneider 2001) to zero to simplify the notation.
The FPFS ellipticity, a spin-2 observable, which negates under a 90-degrees rotation (see appendix B of Li & Mandelbaum 2023 for more details), is used to measure weak lensing shear.It is typically described as a complex number, g = g 1 + ig 2 .For this work, it is defined with the FPFS shapelet modes (denoted as M nm ): The weighting parameter C, introduced by Li et al. (2018), is a constant for a galaxy sample.Different choices for C assign different relative weights to galaxies with different brightnesses.Moreover, it ensures that the noise bias correction remains finite and manageable (Li et al. 2022b).M nm are polar shapelet modes measured from the observed noisy image after PSF deconvolution in Fourier space (Li et al. 2018): where f p k is the observed (PSF-convolved, noisy) image in Fourier space and p k is the PSF image in Fourier space, and χnm is the Fourier transform of the polar shapelet basis function with radial quantum number "n" and angular quantum number "m" (also known as spin number).The polar shapelet basis functions (Massey & Refregier 2005;Bernstein & Jarvis 2002) are defined as where are the Laguerre polynomials and σ h is the smoothing scale of shapelets and the detection kernel.n can be any non-negative integer, and m is an integer between −n and n in steps of two.σ h is set to 0. ′′ 52 in this paper, which optimize the effective galaxy number density of the algorithm for simulated images with seeing size of 0. ′′ 8 .
In addition to shapelets, we use eight detection modes (v i , where i = 0 . . .7) to detect galaxies, and the detection modes are computed for each pixel as where the coordinate center is set to the pixel center.The detection kernels ψ * i (k) are defined in Fourier space: where x i are vectors with length equal to the image pixel side length and orientations pointing towards eight directions (i = 0 . . .7) separated by π/4 .The detection modes measure the difference in brightness between this pixel and its eight adjacent pixels.In practice, we convolve images with the detection kernels to compute detection modes for each pixel, and use a smooth step function to select peaks with all of these detection modes above a threshold as galaxy candidates.In our current detection algorithm, galaxies with centers proximate to pixel edges or corners may be identified twice or four times, respectively.However, each of these detections receives a lesser weight, approximately 0.5 or 0.25 depending on the number of detections.Furthermore, the centroid of each detection is located at the peak pixel's center.Given that the shear response of centroid refinement cannot be analytically derived with our current methods, it prevents us from refining the estimated galaxy centroid.We intend to delve deeper into these aspects and optimize them in future work on the detection algorithm.
More specifically, we use a galaxy weight w to select galaxies with flux and size above corresponding thresholds (Li et al. 2022b), and detect galaxies from images (Li & Mandelbaum 2023).The galaxy weight is required to be second-order differentiable for the noise bias correction (Li & Mandelbaum 2023).The weight can be constructed as a product of either truncated sine functions (see equation ( 46 where T sel 0 is used to select bright galaxies with high signal-to-noise ratio (SNR), T sel 2 is used to select well-resolved large galaxies, and T det is used to detect galaxies.Effectively, we assign a weight to each pixel, and the weight reflects the possibility of the peak being the center of a galaxy.
Our detection algorithm is a simplified version of the one used by SExtractor (Bertin & Arnouts 1996).It operates by identifying peak pixels that exceed the sky background threshold.However, unlike the SExtractor peak detection algorithm, our method does not employ an adaptive process to refine galaxy centers and galaxy boundaries.This is an "ambiguous detection", where detected peaks can have larger offsets from true galaxy center compared to detections with centroid refinement.Despite this limitation, our algorithm excels in providing an accurate analytic response to shear.This response is crucial for correcting biases in shear estimation, making our approach both effective and reliable in this respect.
As shown in equation ( 9) of Li & Mandelbaum (2023), assuming that the histogram of image noise is symmetric with respect to zero after subtracting the background level, the expectation value of the noisy sheared observable (with a tilde) is related to the noiseless sheared one (without accent tag) as and, following equation ( 6) of Li & Mandelbaum (2023), the expectation value of the sheared noiseless observable (without accent tag) is related to the intrinsic unsheared noiseless one (with a bar) as q i;ρ ≡ ∂q i /∂g ρ is the shear response of the linear observable q i .K q i q j ≡ ⟨δq i δq j ⟩ n is the covariance matrix of the linear observables that are used by the FPFS shear estimator, which is demonstrated in Fig. 1.Here ⟨ • ⟩ is the average over noiseless or noisy galaxies; whereas ⟨ • ⟩ n is the average over noise realizations.In Equations ( 8) and ( 9), q i are the array of shapelets and detection modes.We note that, in principle, q i can be any linear observables as long as we derive the shear responses and noise covariance matrix for these linear observables.We refer interested readers to Li & Mandelbaum (2023) for a more detailed discussion.
It is worth noting that, when deriving the algorithm, we operate under the assumption that PSFs are oversampled.This implies that galaxies, as extended objects, are also oversampled.According to the Nyquist-Shannon sampling theorem, a continuous signal can be accurately reconstructed from its sampled, discrete counterpart without any information loss, and the signal amplitude beyond the Nyquist wave number of our sampling remains sufficiently low in Fourier space.Consequently, after sampling the continuous signal, there is no aliasing bias on the Discrete Fourier Transform of the sampled signal below the Nyquist wave number.This assumption is generally valid and practical for ground-based astronomical observations.In practice, we compute the DFT directly from images using a sufficiently large postage stamp size so that the wave-number resolution is adequately high to ensure that the integrals in equations ( 3) and ( 5) can be effectively approximated by a Riemann sum.Furthermore, it's important to note that image noise, which arises due to the quantum nature of light, can actually be undersampled even in ground-based images.This is because image noise is not smeared by the PSF, leading to its susceptibility to the aliasing effect in observed images.For accurate noise bias correction, it becomes essential to estimate the two-point correlation function of the aliased noise.This estimation is derived from blank pixels that do not contain any detections.Under the condition that PSFs are undersampled for single exposure images of space missions, both the observed stellar images and galaxy images experience aliasing bias.This bias predominantly affects large |k| near the Nyquist wave number in Fourier space.It is important to note that this aliasing effect cannot be straightforwardly eliminated through PSF deconvolution.We refer readers to

Shear Response and Noise Bias Functionals
Following Li & Mandelbaum (2023), we correct for noise bias in the weighted ellipticity using equation ( 8), and the noise bias functional is denoted as S , where The noise bias functional returns the noise bias correction function for an input function of weighted ellipticity.Additionally, the functionals for the two diagonal elements of the shear response matrix of the weighted ellipticity are given by where µ = 1, 2 (Li & Mandelbaum 2023).The two off-diagonal terms in the response are spin-4 and therefore have zero mean.Moreover, following Li & Mandelbaum (2023), we can use equation ( 8) to correct for noise bias in the shear response measurement, and the correction term is S (R µ ) .Finally, the shear estimator is given by ĝµ The numerator and the denominator in equation ( 12) are evaluated and averaged over the galaxies in the galaxy sample.Similar to Li et al. (2022b); Li & Mandelbaum (2023), we use we µ = wẽ µ −S wẽ µ and R µ = R µ wẽ µ −S R µ wẽ µ to denote the corresponding observables after noise bias correction.
In this paper, we put the shear response functional (R) and the noise bias functional (S ) into the JAX ecosystem, in order to utilise the AD in JAX to derive the shear response and noise bias correction functions for any definition of ellipticity, selection and detection.

APPLICATIONS
In this section, we present two applications of the AD code to the FPFS method, resulting in an improved shear estimator within the AnaCal framework.In section 3.1, we improve the noise bias correction by including all the terms in the second-order noise bias correction.In section 3.2, we extend the FPFS ellipticity to a more generic form and find the optimal shear estimator from among the more general family of estimators considered.
We quantify the bias in the shear estimation using multiplicative (m 1,2 ) and additive (c 1,2 ) biases (Huterer et al. 2006).Specifically, the estimated shear, g 1,2 , is related to the true shear, g 1,2 , as The LSST-like image simulation used in this paper is introduced in Sheldon et al. ( 2023), and we give a brief review in appendix A. We show a simulated g-band and a griz-band coadded image with the ten-year LSST noise level in Fig. 2. In this work, we use the coadded LSST ten-year images in griz bands to perform detection and shear estimation6 .Our image simulations are divided into 5000 subfields corresponding to 590 square degrees in total.We use all of them for section 3.1 and a subset 100 of them for section 3.2.
Although the simulation package (Sheldon et al. 2023) can be used to generate images with stars, bad pixels, bright-star masks and image coaddition, we focus on testing our updated shear estimation algorithm with the basic galaxy-only image simulations in this paper.We will test our algorithm with more realistic setups (including these additional systematics) in future work.

Second Order Noise Bias Correction
In Li & Mandelbaum (2023), we analytically derived the shear responses R µ (we), the noise bias corrections S (we), and S R µ (we) for the FPFS weighted ellipticity.However, when deriving the noise bias correction, we set many elements in the covariance matrix of measurement error (as shown in the bottom panel of Fig. 1) to zero to simplify our derivation.We found a multiplicative bias in the estimated shear ranging from −0.3% to −0.5% with 1-σ error of ∼0.1%, and the magnitude of the bias increases as the image noise level increases, which suggests that the bias originates from insufficient noise bias correction.
In this section, we compare the multiplicative biases between the analytic version of our FPFS shear estimator and the JAX version on the LSST-like image simulations.We note that the difference between these two codes is that the JAX version (publicly available within the AnaCal framework) includes the complete secondorder noise bias correction while the analytic version neglects several terms in the covariance matrix of the measurement error.2) with C = 7.6 .Red lines are the full noise bias correction with JAX, and blue lines are the analytic noise bias correction (labeled as ANA).The errorbars show the 1σ and 3σ uncertainties, and the shaded region is the LSST ten-year requirement on control of residual multiplicative biases in shear (The LSST Dark Energy Science Collaboration et al. 2018).Use of the full noise bias correction from JAX makes the shear estimator noisier, but its multiplicative bias is consistently below the LSST requirement (99.7% confidence level).In this paper, we use SNR> 12 as our default selection, as is indicated with the vertical dashed line.
Fig. 3 shows the multiplicative and additive biases as functions of the SNR cut at the lower end.The resolution cut is set to 0.05 and the lower end following Li & Mandelbaum (2023).We find that there is a negative multiplicative bias for the analytic correction, which is consistent with the findings in Li & Mandelbaum (2023).In contrast, for the JAX correction, the multiplicative bias is consistent with zero, and its magnitude is less than 0.3 × 10 −3 , the LSST ten-year requirement.

Extension of FPFS
To improve the precision of shear estimation with an optimal weighting scheme, we extend the FPFS ellipticity in equation ( 2) to a more generic form: where σ 0 and σ 2 are the 1σ measurement errors on M 00 and M 00 + M 20 due to image noise.α and β are power-law exponents.If we set c 0 = C/σ 0 , c 2 = 25, α = 1 and β = 0 7 , the extended ellipticity defined in equation ( 14) reduces to the original FPFS ellipticity defined in equation ( 2).We note that this extension of the FPFS ellipticity gives us more freedom to adjust the relative weights between galaxies with different properties.We can use the JAX framework to minimize the statistical error (including intrinsic shape noise and shape measurement error) in the extended hyperparameter space (α, β, c 0 , c 2 ).
To study the dependency of the statistical error on the hyperparameters and find the minima of the statistical error function, we use 100 subfields of simulated noisy blended galaxy images (see Fig. 2 for an example).Each simulated image has different realizations of galaxy positions, galaxy population and image noise realizations, and it covers 0.12 square degrees.We measure shear from each image following the process in Appendix A and use these 100 shear measurements to derive the statistical error.After normalizing the statistical error according to the area of each simulated image, we obtain the per component statistical error on shear estimation for a region of one square arcmin, which is denoted as σ g .Then the effec-7 Note, once β = 0, c 2 does not matter in the ellipticity definition.Figure 5.The effective galaxy number density as a function of c 0 and c 2 for a fixed SNR lower limit of 12.We refer readers to equation ( 14) for the definitions of c 0 and c 2 .The optimal configuration is indicated with a black point.We fix the other two hyperparameters to the optimal ones: α = 0.83, β = 0.18 .We refer readers to equation ( 14) for the definitions of α and β .The optimal configuration is indicated with a black point.We fix the other two hyperparameters to the optimal ones: c 0 = 6, c 2 = 25 .
tive number density is estimated as where 0.26 is the per component root-mean-square (RMS) of intrinsic shape noise for the reGauss shear estimation method (Hirata & Seljak 2003) that is widely used in the weak lensing community.Although different methods have different shape noise RMS, we normalize the effective number density with the reGauss intrinsic shape noise RMS (0.26) so that we can compare our estimation of effective number density with other methods that adopt the same normalization and the predictions using reGauss (e.g., The LSST Dark Energy Science Collaboration et al. 2018).In summary, we first estimate the statistical uncertainty in shear estimation for a one square arcmin region.Then we derive the equivalent galaxy number for the shear if it were to be estimated with reGauss.Note that we do not estimate the effective number density by counting the galaxy number (with weights) after detecting galaxies from images (Chang et al. 2013;The LSST Dark Energy Science Collaboration et al. 2018); that process assumes that the statistical error in shear estimation for each galaxy is uncorrelated.Finally, we use the Nelder-Mead minimizer (Nelder & Mead 1965) implemented in scipy (Virtanen et al. 2020) to find the minimum of the statistical error (σ g ), or equivalently the maximum of the effective number density (n eff ) -which occurs at parameter values c 0 = 6, c 2 = 25, α = 0.83, β = 0.18 .At this point, the effective galaxy number density is 24.8 arcmin −2 .
During the optimization, we fix the SNR cut to 12, and the resolution cut to 0.05 to select galaxies that are sufficiently bright and well-resolved for shear inference.In Fig. 4, we show the multiplicative and additive biases for the estimator with the optimal hyperparameters.We find that the magnitude of the multiplicative bias is consistently below 0.3%, the LSST science requirement (The LSST Dark Energy Science Collaboration et al. 2018).The magnitude of the additive bias is consistently below 5 × 10 −5 .
In Figs. 5 and 6, we show the n eff on 2D slices through the highdimensional hyperparameter space, corresponding to the (c 0 , c 2 ) and (α, β) planes.The maximum of the effective number density is ∼5% larger than that for the original FPFS shear estimator.We find that the effective number density is much more sensitive to the choice of power-law indices (α, β) than to the constants (c 0 , c 2 ) since the power law indices have a more significant influence on the relative weights between galaxies with different brightnesses.Additionally, the effective number density varies rapidly in the β = 2 α direction in the (α, β) plane, and it varies slowly in the orthogonal direction.Moreover, the JAX framework facilitates the flexibility of the shear estimator by allowing both the weight and ellipticity to be defined as any functions of the linear modes.This capability expedites the exploration and optimization of effective number density within the space of possible functions while still ensuring the method is unbiased.
Additionally, we show the effective number density as a function of SNR cut in Fig. 7.We find that for an SNR lower limit exceeding 15, the effective number density significantly decreases.However, the plot shows that there is relatively little value in decreasing the SNR lower limit well below 15, since galaxies at this faint end are dominated by the measurement error caused by image noise.This finding explains our choice of a default lower limit of SNR> 12.However, we note that the absolute value of the effective number density may not be realistic since we have not precisely calibrated the number density in the simulations to any real dataset.We primarily focus, therefore, on relative changes in the number density.Figs.5-7 are used to demonstrate the possibility of optimizing the effective number density from the image level, but not to predict the effective number density of the LSST final year survey.
To make sure that our shear estimator is stable, we show the multiplicative bias on the 1D slice lines through this four-dimensional parameter space in appendix B. We find that the multiplicative bias is not significantly dependent on the hyperparameters, and it is less than 3 × 10 −3 within a wide range of values in the hyperparameter space.

Benchmarking
We benchmark the code on a Intel(R) Xeon(R) Gold 6240R CPU @ 2.40GHz using a simulated region covering one square degree with ten-year LSST depth to compress images to catalogs, and report the results as follows: (i) The FPFS code (https://github.com/mr-superonion/FPFS/) takes 0.18 CPU hour to detect and process the images and generate shear catalogs; (ii) The impt code (https://github.com/mr-superonion/imPT/) takes 0.045 CPU hour to derive shear from the catalogs.
(iii) The analytic code in Li & Mandelbaum (2023) takes 0.004 CPU hours to derive shear from the catalogs.
In summary, with the JAX-based version of FPFS that is available through AnaCal, it would take less than four thousand CPU hours to infer shear from the LSST ten-year coadds across the entire survey area.

SUMMARY AND OUTLOOK
In this paper, we show the results of implementing the FPFS shear estimator within the JAX ecosystem to utilise the automatic differentiation (AD) functionality in JAX to improve the second-order noise bias correction.Furthermore, we extend the FPFS shear estimator to a generic form and use the AD to derive shear response and noise bias correction for that estimator.
We test the algorithm with the LSST-like image simulations (Sheldon et al. 2023) using the coadded image in the griz-bands, and we report that the effective number density is improved by ∼5% compared to the original FPFS shear estimator.We find that the magnitude of the multiplicative bias is consistently less than 0.3% (LSST ten-year requirement on the control of multiplicative bias) within the 3σ uncertainties, illustrating the promise of this shear estimator for LSST.
In the future, as a further step towards a practical use of the application in real survey data, we will apply the shear estimator to image simulations with stars, bad pixels, and bright star masks as in Sheldon et al. (2023).This will enable confirmation that its performance carries over even in the presence of these additional challenges, along with opportunities to further develop the method to face any issues that this exercise reveals.
In order to provide a simple interface through one repository that enables users to seamlessly update as new versions (even with substantially different formalism) become available, we introduce a novel analytical framework for shear calibration, termed AnaCal.This framework is devised to bridge various analytic shear estimators that have been developed (e.g., FPFS) or are anticipated to be created in the future.We intend to develop a suite of analytical shear estimators capable of inferring shear with subpercent accuracy, all while maintaining minimal computational time.Below, we outline several possible avenues for future shear estimation method development.
• Enhancing Shear Estimation through Higher-Order Moments: Typically, the shear estimation process employs the lowest-order (specifically, second-order) spin-2 moments.Nonetheless, as suggested in Refregier & Bacon (2003), the utilization of fourth-order moments holds potential for advancing shear inference.We plan to incorporate this fourth-order shear estimator into the FPFS framework, aiming to amplify the SNR in shear inference by leveraging both second and fourth-order moments while providing additional opportunities for PSF systematics mitigation.
• Adaptive Kernel Shear Estimator: The existing FPFS shear estimator utilizes a fixed-kernel approach; more specifically, the kernel (namely, the Gaussian kernel in shapelets) does not dynamically adjust according to the size and shape of each observed galaxy.In future work, we aim to augment the precision of shear estimations by developing an adaptive kernel approach (e.g., Hirata & Seljak 2003), thereby improving the SNR of the estimated shear.that R 1 + = R 1 − which is true for our simulation since the input galaxy sample is the same for the images with positively and negatively distorted galaxies.The galaxies in each orthogonal galaxy pair and galaxies with different applied shears are selected and weighted independently.After that, we apply our shear estimator to the selected sample to test its performance.In addition, we divide our full galaxy sample into 5000 subfields (corresponding to about 590 square deg in total), and the errors on the means of multiplicative bias and additive bias are estimated using the standard deviation of these 5000 measurements.

APPENDIX B: STABILITY OF THE ESTIMATOR
We show the stability of our shear estimator with different hyperparameters (α, β, c 0 , c 2 ) in Fig. B1.Specifically, we check the multiplicative bias on the 1D slices along these four hyperparameters separately on a wide range covering the optimal parameter.The results show that the accuracy is not significantly dependent on the hyperparameters, and the multiplicative bias is below 3 × 10 −3 at the 3σ level.

Figure 1 .
Figure 1.The covariance matrix (K v i v j in equation (8)) of the measurement errors on linear observables (e.g., seven shapelet modes and eight detection modes with corresponding sixteen shear responses) due to the image noise.The shapelet modes from left to right are M 00 , M 20 , real and imaginary components of M 22 , M 40 , and real and imaginary components of M 42 .As shown, the elements in the lower-left block have the largest magnitude dominating the noise bias correction.The upper panel shows the full covariance matrix, and the lower panel shows the covariance matrix used for the analytic noise bias correction in Li & Mandelbaum (2023), which has several sub-blocks that were set to 0.
) of Li & Mandelbaum 2023) or truncated polynomial functions (Hazimeh et al. 2020) of shapelet modes M nm and eight detection modes

Figure 2 .
Figure 2. The g-band image (left), which has the highest SNR, and the coadded image of griz-bands (right) of the LSST-like image simulation (Sheldon et al. 2023) of the same region using the expected ten-year LSST noise level.The region is a random cut-out from the simulated exposure, and it covers 3.4 × 3.4 square arcmins.The images are coadded with inverse variance weight based on background photon noise.
Kannawadi et al. (2021) andYamamoto et al. (2023b) for testing shear estimation code on undersampled image simulations that are representative to single exposure images from space missions.However, we are planning to apply the algorithm to coadded images derived from single exposure images, and we refer readers toHirata et al. (2023) andYamamoto et al. (2023a) for deriving oversampled coadded images from single exposure images.However, since the formalism inLi & Mandelbaum (2023) work neglected a few relatively small elements in the Hessian matrix to simplify the analytic derivation of noise bias correction, the final results showed a residual bias of about -0.4 percent relative to the true shear for faint galaxies.Here we use JAX's AD feature to automatically derive the full noise bias correction, to extend the ellipticity and detection / selection weights to a general form, and to derive the shear response and noise bias for them.

Figure 3 .
Figure 3.The bias (upper panel) and additive bias (lower panel) of the FPFS shear estimator defined in equation (2) with C = 7.6 .Red lines are the full noise bias correction with JAX, and blue lines are the analytic noise bias correction (labeled as ANA).The errorbars show the 1σ and 3σ uncertainties, and the shaded region is the LSST ten-year requirement on control of residual multiplicative biases in shear (The LSST Dark Energy ScienceCollaboration et al. 2018).Use of the full noise bias correction from JAX makes the shear estimator noisier, but its multiplicative bias is consistently below the LSST requirement (99.7% confidence level).In this paper, we use SNR> 12 as our default selection, as is indicated with the vertical dashed line.

Figure 4 .
Figure 4.The multiplicative bias (upper panel) and additive bias (lower panel) of the generalized FPFS shear estimator defined in equation (14) with the optimal parameters: (c 0 = 6, c 2 = 25, α = 0.83, β = 0.18).The errorbars show the 1σ and 3σ uncertainties, and the shaded region is the LSST tenyear requirement on control of residual multiplicative biases in shear.The multiplicative bias for the generic FPFS shear estimator is consistently below the LSST requirement (99.7% confidence level).

Figure 6 .
Figure6.The effective galaxy number density as a function of α and β .We refer readers to equation (14) for the definitions of α and β .The optimal configuration is indicated with a black point.We fix the other two hyperparameters to the optimal ones: c 0 = 6, c 2 = 25 .

Figure 7 .
Figure7.The effective galaxy number density as a function of SNR cut.We set the hyperparameters to c 0 = 6, c 2 = 25, α = 0.83, β = 0.18 , corresponding to those for the optimal shear estimator identified in this work.The vertical dashed line indicates the default SNR cut used in this paper.

Figure B1 .
Figure B1.The multiplicative bias as a function of hyperparameters (α, β, c 0 , c 2 ).We show the function in 1D slices along these hyperparameters separately in four panels.The errorbars show the 1σ and 3σ uncertainties, and the shaded region is the LSST ten-year requirement on the systematic control (The LSST Dark Energy Science Collaboration et al. 2018).The vertical dashed lines indicate the optimal setups.