Selection functions of strong lens finding neural networks
Convolution Neural Networks trained for the task of lens finding with similar
architecture and training data as is commonly found in the literature are
biased classifiers. An understanding of the selection function of lens finding
neural networks will be key to fully realising the potential of the large
samples of strong gravitational lens systems that will be found in upcoming
wide-field surveys. We use three training datasets, representative of those
used to train galaxy-galaxy and galaxy-quasar lens finding neural networks. The
networks preferentially select systems with larger Einstein radii and larger
sources with more concentrated source-light distributions. Increasing the
detection significance threshold to 12$\sigma$ from 8$\sigma$ results in 50 per
cent of the selected strong lens systems having Einstein radii
$\theta_\mathrm{E}$ $\ge$ 1.04 arcsec from $\theta_\mathrm{E}$ $\ge$ 0.879
arcsec, source radii $R_S$ $\ge$ 0.194 arcsec from $R_S$ $\ge$ 0.178 arcsec and
source S\'ersic indices $n_{\mathrm{Sc}}^{\mathrm{S}}$ $\ge$ 2.62 from
$n_{\mathrm{Sc}}^{\mathrm{S}}$ $\ge$ 2.55. The model trained to find lensed
quasars shows a stronger preference for higher lens ellipticities than those
trained to find lensed galaxies. The selection function is independent of the
slope of the power-law of the mass profiles, hence measurements of this
quantity will be unaffected. The lens finder selection function reinforces that
of the lensing cross-section, and thus we expect our findings to be a general
result for all galaxy-galaxy and galaxy-quasar lens finding neural networks.
A. Herle, C. M. O'Riordan, S. Vegetti
# Selection functions of strong lens finding neural networks

###### Abstract

Convolution Neural Networks trained for the task of lens finding with similar architecture and training data as is commonly found in the literature are biased classifiers. An understanding of the selection function of lens finding neural networks will be key to fully realising the potential of the large samples of strong gravitational lens systems that will be found in upcoming wide-field surveys. We use three training datasets, representative of those used to train galaxy-galaxy and galaxy-quasar lens finding neural networks. The networks preferentially select systems with larger Einstein radii and larger sources with more concentrated source-light distributions. Increasing the detection significance threshold to \(12\sigma\) from \(8\sigma\) results in 50 per cent of the selected strong lens systems having Einstein radii \(\theta_{\rm E}\geq 1.04\) arcsec from \(\theta_{\rm E}\geq 0.879\) arcsec, source radii \(R_{S}\geq 0.194\) arcsec from \(R_{S}\geq 0.178\) arcsec and source Sersic indices \(n^{\rm S}_{\rm Sc}\geq 2.62\) from \(n^{\rm S}_{\rm Sc}\geq 2.55\). The model trained to find lensed quasars shows a stronger preference for higher lens ellipticities than those trained to find lensed galaxies. The selection function is independent of the slope of the power-law of the mass profiles, hence measurements of this quantity will be unaffected. The lens finder selection function reinforces that of the lensing cross-section, and thus we expect our findings to be a general result for all galaxy-galaxy and galaxy-quasar lens finding neural networks.

keywords: gravitational lensing: strong; methods: data analysis, statistical; techniques: image processing

## 1 Introduction

Strong gravitational lensing allows one to address a broad range of cosmological and astrophysical questions, from the nature of dark matter (e.g Vegetti et al., 2018; Ritondale et al., 2019; Gilman et al., 2020; Hsueh et al., 2020) and galaxy evolution (e.g. Sonnenfeld et al., 2019; Mukherjee et al., 2021; Rizzo et al., 2021) to measuring the Hubble constant and other cosmological parameters (e.g. Birrer & Treu, 2021; Collett & Auger, 2014). The field is set to profit enormously from ongoing and upcoming large-sky surveys (with e.g. Euclid, the Square Kilometre Array and the Vera Rubin Observatory), as the number of known strong gravitational lens systems is expected to increase by many orders of magnitude (Collett, 2015; McKean et al., 2015). An understanding of the selection function of these surveys is required for a correct interpretation of the scientific results from strong gravitational lensing analyses of the larger sample.

There are two sides to the problem of selection bias in the newly discovered strong lens systems: the effect of the lensing cross-section, and that of the lens finding method used to identify the lenses. Sonnenfeld et al. (2023) focused on the former and showed that strong lenses are a biased subset of the true population with respect to the lens and source parameters. The effect of the lens finding methodology, however, is still not well understood.

The task of finding strong gravitational lens systems amongst the large number of objects imaged in future surveys is a challenging one, especially as strong lenses are, by their nature, very rare. Hence, automation in lens finding is expected to play an important role in the identification of these systems.

Neural networks are a class of machine learning models that are often applied in astronomy to search for rare objects. Lens finder neural networks are those trained to sift through large volumes of image data specifically to find strong gravitational lens candidates. Over the years, Convolution Neural Networks (CNNs) have emerged as the state of the art for lens finding. These models are essentially classifiers, trained to distinguish between _lens_ and _non-lens_ gravitational systems. For example, Lanusse et al. (2018) developed CMU DeepLens, which is a CNN trained to find strong gravitational lens systems on LSST data. The ResNet architecture implemented in their work became the standard for this task, and was adapted for the DESI Legacy survey by Huang et al. (2020) and Huang et al. (2021). CNNs have also been used to find strong lenses in the Kilo-Degree Survey (Petrillo et al., 2017, 2019, 2019), HST/ACS data of the COSMOS field (Pourrahmani et al., 2018), PANSTARS (Canameras et al., 2020), HSC (Canameras et al., 2021, 2023) and LOFAR (Rezaei et al., 2022).

In this paper, we present the first systematic study of the strong lens finder selection function. We develop lens finder neural networks in a similar fashion as commonly done in the lens finding literature and then characterise the selection effects that they introduce into the recovered sample of strong gravitational lenses. In particular, we are interested in identifying the characteristics of a strong gravitational lens system that drive the classification. Finally, we discuss the implications of these selection effects for different scientific applications of strong gravitational lensing.

The paper is organised as follows. In Section 2 we discuss the
mathematical formalism used. In Section 3 we describe the image simulation process for the three datasets used in this work. In Section 4 we detail the machine learning training and interpretability framework. In Section 5 we present our results on the selection function. In Section 6 and 7 we discuss our results and summarise our conclusions.

## 2 Selection bias formalism

From Sonnenfeld et al. (2023), the probability \(\mathrm{P_{SL}}\) of selecting a sample of strong lens systems for a given selection criterion S can be expressed as:

\[\mathrm{P_{SL}}(\Psi_{I},\Psi_{S}|S)\propto\mathrm{P_{I}}(\Psi_{I})\mathrm{P_{ s}}(\Psi_{S})\mathrm{P_{sl}}(\Psi_{I},\Psi_{S}|S)\;. \tag{1}\]

Here, \(\Psi_{I}\) and \(\Psi_{s}\) are the set of parameters describing the lens galaxy mass and the background source light distributions, respectively. \(\mathrm{P_{I}}\) and \(\mathrm{P_{s}}\) are the corresponding probability density distributions. \(\mathrm{P_{sel}}\) encapsulates the probability that a specific combination of lens and source produce a strong lens system and that this is found in the survey. It can be further separated into two components:

\[\mathrm{P_{sel}}(\Psi_{I},\Psi_{S}|S)=\mathrm{P_{det}}(\Psi_{I},\Psi_{S}) \mathrm{P_{find}}(\Psi_{I},\Psi_{S}|S). \tag{2}\]

where \(\mathrm{P_{det}}\) is the probability that multiple images of the source form as clearly distinct features in the survey image, and \(\mathrm{P_{find}}\) is the probability that this image is correctly classified as a strong lens system. The latter depends on the identification procedure adopted. By focusing on the first term, Sonnenfeld (2022) and Sonnenfeld et al. (2023) have quantified how strongly the properties (\(\Psi_{I}\) and \(\Psi_{s}\)) of the lenses and background sources are biased with respect to the general population of galaxies. The goal of this paper is to quantify \(\mathrm{P_{find}}\) for the case in which a CNN is used to identify lens systems in a given survey.

## 3 Data

The Euclid mission is expected to find about \(10^{5}\) strong gravitational lens systems (Collett, 2015). In this work, we focus on neural networks trained only with single-band data because most of the gravitational lens systems that will be found in the Euclid Wide Survey will be observed using the Visual imager (VIS) instrument.

VIS is a broadband optical instrument with a resolution of 0.16 arcsec, which is about three times better than that of the three infrared bands of the Near Infrared Spectrometer and Photometer (NISP) instrument (O'Riordan et al., 2023; Euclid Collaboration et al., 2022). Since the Einstein radius of the strong gravitational lenses expected to be found by Euclid peaks around 0.5 arcsec (Collett, 2015), the arcs and lens light are expected to be blended together within the NISP instrument. Hence, these features will be much better resolved by VIS, at the cost of losing colour information. The absence of multi-band data can be crucial in informing the decision rationale that CNNs arrive at during training, and may have a significant effect on the type of strong gravitational lens systems that will be identified. Moreover, Petrillo et al. (2019) found that the best performing network in terms of purity and completeness was the one trained on single-band data. Similarly, Lanusse et al. (2018) showed that CMU DeepLens is a successful tool for finding strong gravitational lens systems, and that even without colour information, it is able to learn enough about the lensed-arc morphology to solve the classification problem.

We make three datasets, \(\mathcal{D_{A}}\), \(\mathcal{D_{B}}\) and \(\mathcal{D_{C}}\), with \(10^{6}\) images for training and \(2\times 10^{5}\) for testing each. The samples are split evenly between two distinct classes, _lens_ and _non-lens_, which contain all and none of the lens systems, respectively.

Fig. 1 and 2 show examples of the two classes from \(\mathcal{D_{B}}\) and \(\mathcal{D_{C}}\), and Table 1 summarises the simulation parameter distributions for all datasets.

We refer the reader to the following sections for more details on the properties of each dataset. Briefly, the datasets \(\mathcal{D_{A}}\) and \(\mathcal{D_{B}}\) have extended sources, while the sources in \(\mathcal{D_{C}}\) are unresolved. This is done to understand how the selection function, \(\mathrm{P_{find}}\), differs between galaxy-galaxy and galaxy-quasar lens systems. The lens galaxies in \(\mathcal{D_{A}}\) have a simple analytical model for the lens-light distribution, while those in \(\mathcal{D_{B}}\) and \(\mathcal{D_{C}}\) are characterised by a higher level of complexity. This is done in order to quantify how the lens light model complexity affects the neural network selection function.

We add complexity to the lens rather than the source light because distinguishing between these components becomes non-trivial when the lens contains features resembling lensed source emission (e.g. arcs). The lensed emission is typically sufficiently distinct from analytic Sersic components because of its distorted morphology, such that introducing structure in the source will only make the task of lens finding easier. Moreover, using analytical source models allows us to more easily quantify how their parameters affect the network selection function for different choices of the lens light distribution.

### Source light

We use Sersic profiles (Sersic, 1963; Ciotti and Bertin, 1999) for the source-light model in the case of \(\mathcal{D_{A}}\) and \(\mathcal{D_{B}}\). A point-like source model is used for \(\mathcal{D_{C}}\). We refer the reader to Table 1 for more detailed information.

### Lens light

The lens light distribution in \(\mathcal{D_{A}}\) is created using a Sersic profile with an index and effective radius sampled uniformly from between 1.0 and 4.0 and between 0.05 and 0.5 arcsec, respectively. \(\mathcal{D_{B}}\) and \(\mathcal{D_{C}}\) use images generated with a Generative Adverserial Network (GAN, see Holzschuh et al., 2022, for more details). The GAN was trained to generate images that imitate those created from the SKIRT code on the IllustrisTNG simulation (Springel et al., 2018; Rodriguez-Gomez et al., 2019). The resulting datasets consist, therefore, of realistic early- (i.e. Sersic-like) and late-type galaxy images, several of which contain complex features like star-forming clumps, spiral arms and satellite galaxies (see Fig. 1 for examples of lens systems taken from \(\mathcal{D_{B}}\)). Structures of this kind can easily be confused for lensed-source emission when dealing with single-band data, making the problem of lens-finding more challenging. In this respect, \(\mathcal{D_{A}}\) represents the simplest formulation of the problem for the network to solve.

The lens magnitudes in \(\mathcal{D_{A}}\) and \(\mathcal{D_{B}}\) are sampled from a uniform distribution, and the source magnitudes are chosen such that they are dimmer than the lens. In \(\mathcal{D_{C}}\), we instead use fluxes defined relative to the source. In practice, the lens flux is set such that it is 5 to 100 times brighter (sampled uniformly in this range) than the source light. We chose these values based on the distribution of magnifications of the lens mass models (see Section 3.4 for more details) in this dataset, which peaks at \(\mu\approx 3\), thus ensuring that the number of systems where the lens light is completely absent in the image are negligible. These values are also chosen such as to keep the flux per pixel area the same, as the source light is concentrated into a point.
Figure 1: A representative sample from \(\mathcal{D}_{\mathcal{B}}\), ranked in increasing order of SNR from the top left to the bottom right in each panel. Note the range of complexities of the lens light models in the sample. In some cases, distinguishing between lens light and lensed source emission is trivial, but in cases where the lens light model contains complex features like spiral arms, this is not the case. Examples of lenses are shown on the left, and examples of non-lenses are shown on the right. The simulation parameters are sampled as shown in Table 1.

Figure 2: A representative sample from \(\mathcal{D}_{\mathcal{C}}\), ranked in increasing order of SNR from the top left to the bottom right in each panel. The images are simulated as described in Section 4, with the inclusion of contaminant galaxies and a variable PSF FWHM. Examples of lenses are shown on the left, and examples of non-lenses are shown on the right. The simulation parameters are sampled as shown in Table 1.
### Light contaminants

We do not include any light contaminants in \(\mathcal{D}_{\mathcal{R}}\) and \(\mathcal{D}_{\mathcal{B}}\). On the other hand, \(\mathcal{D}_{\mathcal{C}}\) does contain light contribution from nearby field galaxies. These are sampled from the GAN dataset and are placed randomly on the image, while ensuring that they lie outside the Einstein radius of the system and at least a set minimum distance apart from each other. These contaminants are 1 to 100 times brighter than the source light and their number in each image is sampled uniformly between 0 and 4. We sample their redshift between 0.2 and 5.0 from a probability density distribution based on the comoving volume at these redshifts. Examples of gravitational lens systems from \(\mathcal{D}_{\mathcal{C}}\) are shown in Fig. 2.

### Lens mass

For all three datasets, an elliptical power law is assumed for the lens mass model with the slope, \(\gamma\), ranging uniformly between 1.8 and 2.2 (Koopmans et al., 2006). In the case where the lens light is a GAN image (i.e. \(\mathcal{D}_{\mathcal{B}}\) and \(\mathcal{D}_{\mathcal{C}}\)), we use the image moments of the latter to align the position angle and axis-ratio of the mass profile with those of the lens light, thus ensuring that the light and mass distributions are consistent with each other. Moreover, the Einstein radius is scaled proportionally to the total flux of the selected GAN image. Thus, \(\mathcal{D}_{\mathcal{B}}\) and \(\mathcal{D}_{\mathcal{C}}\) have the same simulation parameters as \(\mathcal{D}_{\mathcal{R}}\), except for the Einstein radius and lens axis-ratio, for which the distributions are shown in the left and centre panels of Fig. 3, respectively.

External shear with no preferred direction and with strength sampled from a uniform distribution is applied to all datasets. Additionally, for \(\mathcal{D}_{\mathcal{C}}\), we further ensure that the numbers of doubly- and quadruply-imaged systems are equal within the lens class. This is done by first calculating the area inside the inner caustic for the specific axis-ratio chosen. We then find the radius of a circle that would be twice this area and sample the source position uniformly from within this circle.

### Noise

In order to add noise to the data, we follow the definition of signal-to-noise ratio (SNR) by O'Riordan et al. (2019):

\[S_{\mathrm{T}}=\frac{\Sigma_{i}^{N}m_{i}\mathrm{d}_{i}}{\sigma_{d}\sqrt{ \Sigma_{i}^{N}m_{i}}}\,, \tag{3}\]

where \(\sigma_{d}\) is the standard deviation of the sky noise, \(\mathrm{d}_{i}\) is the pixel value at index \(i\), and \(m_{i}\) is the masking variable defined such that \(m_{i}=1\) for pixels corresponding to lensed source-light and 0 otherwise. Pixel values that lie within twice the Einstein radius of the lens subtracted image of the system are considered to be source-light for \(\mathcal{D}_{\mathcal{R}}\) and \(\mathcal{D}_{\mathcal{B}}\), and in the case of \(\mathcal{D}_{\mathcal{C}}\), this is considered to be all pixels with a value greater than 1.5 times the standard deviation of the lens subtracted image.

The SNR of each image in \(\mathcal{D}_{\mathcal{R}}\) and \(\mathcal{D}_{\mathcal{B}}\) is set by the specific combination of source and lens magnitudes, and the resulting distribution is shown in the right panel of Fig. 3. In the case of \(\mathcal{D}_{\mathcal{C}}\), we first sample a value for the SNR from a uniform distribution (see Table 1), then add uncorrelated Gaussian noise to the images with a standard deviation calculated from Eqn. (3).

### Point-Spread Function

We use a circular Gaussian function for the Point-Spread Function (PSF). The Full Width at Half Maximum (FWHM) has a value of 0.16 arcsec for \(\mathcal{D}_{\mathcal{R}}\) and \(\mathcal{D}_{\mathcal{B}}\). For \(\mathcal{D}_{\mathcal{C}}\), we employ a variable FWHM, which is sampled uniformly between 0.16 and 0.3 arcsec. A FWHM

\begin{table}
\begin{tabular}{l c c c} \hline Parameter & \(\mathcal{D}_{\mathcal{R}}\) & \(\mathcal{D}_{\mathcal{B}}\) & \(\mathcal{D}_{\mathcal{C}}\) \\ \hline Source Light & Sérsic & Sérsic & Point-like \\ Source Radius, \(R_{\mathrm{S}}\) (arcsec) & \(\mathcal{U}(0.05,0.3)\) & \(\mathcal{U}(0.05,0.3)\) & - \\ Source Sérsic Index, \(n_{\mathrm{Sc}}^{\mathrm{S}}\) & \(\mathcal{U}(1,4)\) & \(\mathcal{U}(1,4)\) & - \\ Source Axis-Ratio, \(q_{\mathrm{S}}\) & \(\mathcal{U}(0.5,1.0)\) & \(\mathcal{U}(0.5,1.0)\) & - \\ Source Redshift & 2.0 & 2.0 & 2.0 \\ Source Apparent Magnitude, \(M_{\mathrm{VIS}}^{\mathrm{S}}\) & \(\mathcal{U}(M_{\mathrm{VIS}}^{\mathrm{S}},24.0)\) & \(\mathcal{U}(M_{\mathrm{VIS}}^{\mathrm{S}},24.0)\) & - \\ Source Flux & - & - & 1.0 \\ Lens light & Sérsic & GAN & GAN \\ Lens mass axis-ratio, \(q\) & \(\mathcal{U}(0.5,1.0)\) & - & - \\ Lens power-law slope, \(\gamma\) & \(\mathcal{U}(1.8,2.2)\) & \(\mathcal{U}(1.8,2.2)\) & \(\mathcal{U}(1.8,2.2)\) \\ Einstein Radius, \(\theta_{\mathrm{L}}\) (arcsec) & \(\mathcal{U}(0.5,2.0)\) & - & - \\ Lens redshift & 0.8 & 0.8 & 0.8 \\ Lens Sérsic Index, \(n_{\mathrm{Sc}}^{\mathrm{L}}\) & \(\mathcal{U}(1,4)\) & - & - \\ Lens Apparent Magnitude, \(M_{\mathrm{VIS}}^{\mathrm{L}}\) & \(\mathcal{U}(18.0,22.0)\) & \(\mathcal{U}(18.0,22.0)\) & - \\ Lens Flux, \(\mathrm{F}_{\mathrm{L}}\) & - & & \(\mathcal{U}(5.0,100.0)\) \\ Shear Strength & \(\mathcal{U}(0.0,0.1)\) & \(\mathcal{U}(0.0,0.1)\) & \(\mathcal{U}(0.0,0.1)\) \\ Contaminants & No & No & Yes \\ Contaminant Flux, \(F_{\mathrm{C}}\) & - & - & \(\mathcal{U}(1.0,100.0)\) \\ Pixel Size (arcsec) & 0.1 & 0.1 & 0.1 \\ Field-of-View (arcsec) & 10 & 10 & 10 \\ PSF\({}_{\mathrm{Iwhm}}\) (arcsec) & 0.16 & 0.16 & \(\mathcal{U}(0.16,0.3)\) \\ \hline \end{tabular}
\end{table}
Table 1: Parameter distributions used to create all datasets. From top to bottom: the parameters describing the source light distribution, the mass and light properties of the lens, distribution for the light contaminants and the observational set up. The superscripts S and L denote source and lens properties respectively.
of 0.16 arcsec corresponds to that of the VIS instrument aboard the Euclid telescope, and 0.3 arcsec represents one that is about 2 times worse than this.

## 4 Methods

### Training phase

During the training phase, the datasets were split into batches of \(10^{3}\) images each. ResNet18 CNNs were trained for 800 epochs each on \(\mathcal{D_{A}}\) and \(\mathcal{D_{B}}\). For \(\mathcal{D_{C}}\), we also included dynamic learning rate scheduling, wherein the learning rate is decreased by a factor \(10^{-0.5}\) if the test loss remains stagnant for 20 epochs. These values were chosen after extensive hyper-parameter tuning. The convergence criteria is set to be three consecutive drops in learning rate without a change in the test loss. With this setup, the network is trained for 260 epochs.

For all datasets, several data augmentation techniques were used during training, namely: random cropping to \(80\times 80\) pixels, random vertical and horizontal flips, rotation and erasing. The learning rate was set to 0.001 and the Adam optimizer was used to minimize the binary cross entropy loss. The images were normalised such that each pixel value is between 0 and 1 before being passed to the neural network.

Table 2 shows the differences in the final testing accuracy and loss for each of the three networks, as well as the True Positive Rate (TPR), False Positive Rate (FPR) and Area Under the Receiver Operator Characteristic (AUROC) curve for each network with corresponding dataset version. All three networks achieve very high classification accuracy, and have converged. This is also reflected in the AUROC values: the closer the AUROC is to 1, the closer the network is to a perfect classifier; thus the AUROC curve is a measure of the quality of the classifier for the particular dataset it is trained on.

The small decrease in network performance between \(\mathcal{D_{A}}\) and \(\mathcal{D_{B}}\) is a consequence of the fact that the network is exposed to a larger fraction of lens-light models and negative examples (galaxies in the _non-lens_ class) that have complex structure which can be confused for lensed arcs. We note that the higher test accuracy of the network trained on \(\mathcal{D_{A}}\) does not imply a superior performance, rather, it indicates the relative simplicity of the lens light model.

### Detection significance of a _lens_ image

We use the inverse error function to convert from probabilities to detection significance:

\[C=\sqrt{2}\ \mathrm{erf}^{-1}(p)\, \tag{4}\]

where \(C\) is the network lens detection significance and \(p\) is the probability that an image belongs to the _lens_ class. More details can be found in the Appendix A1. In a lens finding campaign, a detection significance threshold is chosen such that all classifications made by the network above this \(\sigma\)-cut constitute a strong lens system. The distributions of the outputs for the three different networks on their respective datasets, after converting to significance, are shown in the right panel of Fig. 4. Note that the distributions are roughly Gaussian and peak at different values for each dataset.

### Kullback-Leibler divergence

We would like to quantify the difference between the ground truth population of lenses and sources (\(\Psi_{I}\), \(\Psi_{S}\)) and the population recovered by the lens finder. To do this, we employ the Kullback-Leibler (KL) divergence, which is defined as:

\[D_{KL}(p||q)=-\int\ p(x)\ln\left[\frac{q(x)}{p(x)}\right]dx. \tag{5}\]

Here, \(p(x)\) is the probability distribution function of the true parent distribution (all lenses in the testing dataset) for the parameter \(x\), and \(q(x)\) is that of the selected sample at a particular detection

Figure 3: From left to right: distribution of lens Einstein radius, axis-ratio and SNR of the images for the datasets \(\mathcal{D_{B}}\) (orange histograms) and \(\mathcal{D_{C}}\) (green histograms). The bins are in log-space for the panel on the right.
significance threshold. We estimate these distribution functions at a specific \(\sigma\)-cut using Kernel Density Estimation (KDE) with a Gaussian kernel. We use Scott's rule (Scott, 2015) to estimate a reasonable bandwidth for the KDE.

The purpose of using the KL-divergence as a metric is two-fold: (i) from an interpretability point of view, the parameters that show an increase in the KL-divergence at high detection significance thresholds are clearly important for the networks to classify a system as a _lens_ or a _non-lens_, lending insight into what the networks have learnt in order to solve the lens finding problem. (ii) For physical parameters that are measured in a strong lensing survey, the KL-divergence indicates by how much the parameters inferred will deviate from the parent distribution for a given detection threshold.

It must be noted here that we have neglected the correlations between the different parameters when calculating the KL-divergence, which is akin to marginalising over all parameters except the one being considered in the calculation. This means that the variation of the KL-divergence for a specific parameter as the detection significance threshold becomes stricter indicates by how much this parameter will differ from the true sample, if we were only interested in measuring this specific parameter in a given survey. A more thorough analysis, which is beyond the scope of this work, would require the covariance of the parameters to be taken into account.

We also note that as the detection significance threshold increases, the number of lens images drops exponentially, as shown in the left panel of Fig. 4. Thus, there is an increase in the KL-divergence due to sampling noise, which needs to be accounted for in order to disentangle it from the increase in the KL-divergence that is a result of the selected sample differing from the truth. To this end, we sample the true distribution with different numbers of total samples (N) corresponding to each \(\sigma\)-cut, and then calculate the KL-divergence of this sample relative to the the entire true distribution of \(10^{5}\) systems. This is repeated several thousands of times to get enough realisations to account for sampling variance.
model and, therefore, the network does not need to learn complex concepts to distinguish between lensed and un-lensed emission. The more complex lens-light distributions in \(\mathcal{D}_{B}\) cause the network to not only be more sensitive to the lens axis-ratio and the source radius (we see a steep rise in KL-divergence already at \(>7\sigma\) for \(R_{\rm S}\)), but also to the Sersic index of the source. In particular, it is more likely to select sources with larger radii (\(R_{\rm S}\)), Sersic indices (\(n_{\rm SC}^{\rm s}\)) and axis-ratios (\(q_{\rm S}\)). Specifically, 50 per cent of the selected sample have \(R_{\rm S}\geq 0.194\) arcsec, \(n_{\rm SC}^{\rm S}\geq 2.62\) and \(q_{\rm S}\geq 0.758\) at a 12-\(\sigma\) cut as opposed to the ground truth values of \(R_{\rm S}\geq 0.175\) arcsec, \(n_{\rm SC}^{\rm s}\geq 2.51\) and \(q_{\rm S}\geq 0.751\). All these properties lead to larger and more concentrated source light distributions, making the lensed arcs more easily distinguishable from the lens light.

Additionally, all three networks tend to identify as lens systems those with a slightly more elliptical lens mass distribution at very strict detection significance thresholds, (i.e. \(>12\sigma\) for the networks trained on \(\mathcal{D}_{\mathcal{A}}\) and \(\mathcal{D}_{B}\)). An increase in the detection significance cut from 8 to 12\(\sigma\) causes 50 per cent of the selected lens systems to have \(q\leq 0.736\) and \(q\leq 0.713\) from \(q\leq 0.750\) and \(q\leq 0.723\), for \(\mathcal{D}_{\mathcal{A}}\) and \(\mathcal{D}_{B}\) respectively (with associated true values of \(q\leq 0.750\) and \(q\leq 0.723\)). The KL-divergence plot for \(\mathcal{D}_{C}\) (Fig. 10) shows more sensitivity (steep rise in KL-divergence at C \(>6\sigma\)) to the lens axis-ratio compared to the extended source lens finders. At a 7-\(\sigma\) cut, half of the lens systems have lens axis-ratios \(q\leq 0.670\) as opposed to the ground truth values of \(q\leq 0.723\). Moreover, this network is more efficient at selecting images with lower PSFwhm values. At lower PSFwhm, it becomes easier to distinguish between the lensed point-source and the lens/contaminant galaxies. However, the inclusion of stars as contaminants will make this more difficult, and the preference for lower PSFwhm may be negligible with stars as contaminants.

For this dataset, we also calculate the KL-divergence as a function of \(\sigma\)-cut for the flux of the lens \(F_{\rm L}\) and the contaminants, \(F_{\rm C}\), which are proxies for how bright these components are with respect to the source. We find that the KL-divergence is negligible at all values of \(\sigma\)-cut. This points to the fact that the network has learnt to distinguish between the lens/contaminant light and source light very effectively.

Interestingly, the increase in KL-divergence for the lens power-law slope for all three networks is consistent with the sampling noise, leading us to conclude that there is no preference for mass models of a particular slope.

## 6 Discussion

Our main findings indicate that the sample of strong gravitational lens systems identified in wide-sky surveys with CNNs are biased with respect to both lens and source properties. This could have important implications for the many scientific applications of strong gravitational lensing. In this section, we discuss the implications of these biases.

The CNN selection function, coupled with the lensing cross-section, will lead to samples of deflectors which are biased towards higher masses. This will effectively limit our prospects of studying the properties of lens galaxies with strong lensing to only the most massive objects. Essentially, the blurring by the PSF produces a lower limit on the Einstein radii of lens systems that can be found in a survey. This effect will be increased by the selection function of the neural networks.

Interestingly, we find no selection bias in favour of specific values of the slope, \(\gamma\), of the lens mass density profiles. Assuming that this result is confirmed with other CNNs, we can conclude that the large samples of new lens systems will allow one to measure this parameter and its evolution with redshift. This will provide a valuable probe of galaxy evolution and feedback models (e.g. Koopmans et al., 2006; Mukherjee et al., 2021).

Strong gravitational lensing acts as a cosmic telescope, allowing the study of high-redshift galaxies at higher physical resolution and SNR than usually possible otherwise. We have shown that networks
trained on images with complex and more realistic lens light models prefer larger source radii and Sersic indices. Sources with these properties are bigger and more concentrated, and thus produce arcs which are easier to distinguish from the lens light. A similar selection bias also occurs from the strong lensing cross-section (Serjeant, 2012; Hezaveh et al., 2012), and the neural networks will reinforce this effect. Hence, attempts to interpret the properties of lensed galaxies in the context of galaxy evolution models (as for example in Oldham et al., 2017, 2017; Stacey et al., 2021) from the large sample of discovered lenses will need to additionally account for the effect of the neural network selection function.

Time-delay cosmography can be used as an additional probe of the Hubble constant (Refsdal, 1964; Birrer et al., 2020; Birrer & Treu, 2021; Shajib et al., 2023). However, this requires the breaking of the mass-sheet degeneracy (Falco et al., 1985) by incorporating independent mass measurements from stellar dynamics. Velocity dispersion measurement uncertainties limit the overall precision on \(H_{0}\approx 10\) per cent (Kochanek, 2020; Schneider & Sluse, 2013). These uncertainties are too large to address the Hubble tension in a meaningful way. To overcome this issue, Birrer & Treu (2021) proposed a joint analysis using mass models and kinematic properties of galaxy-galaxy lenses as a prior on the mass profiles for the analysis of the

Figure 5: Distributions of the parent and selected sample for detection significance thresholds of 3, 7 and 10\(\sigma\) for the network trained on \(\mathcal{D}\): \(\theta_{\rm E}\) is the Einstein radius, \(q\) is the axis-ratio of the lens mass, \(\gamma\) is the power-law slope, \(M_{\rm VIS}^{\rm L}\) and \(M_{\rm VIS}^{\rm S}\) are the apparent magnitude of the lens and source respectively, \(R_{\rm S}\) is the source radius, \(n_{\rm Sc}^{\rm S}\) is the Sersic index of the source, \(q_{\rm S}\) is the axis-ratio of the source and SNR is the signal-to-noise ratio.
time-delayed galaxy-quasar lenses. This method inherently assumes that two types of strong gravitational lens systems are drawn from the same deflector parent population (see also Sonnenfeld, 2021). However, Sonnenfeld et al. (2023) pointed out that in the case of quadruply-imaged quasars, the lens galaxies tend to have larger ellipticities and halo masses for a given stellar mass. Similarly, we find that the neural network trained to find lensed quasars has a preference for higher ellipticity lens profiles as these have wider caustics and thus a higher chance of producing four image systems. More importantly, the networks trained to find lensed extended sources show a much weaker preference for higher ellipticities. Hence, the approach proposed by Birrer & Treu (2021, see also Gomer et al. 2022; Birrer et al. 2020) is more likely to introduce additional systematic errors as opposed to constraints. These can be accounted for with an understanding of the selection function of the lens finder used, which can be obtained with an analysis of the type done in this work.

### Limitations of our work

Our image simulations already account for several characteristics of real data that is often ignored in the lens finding literature, like complexity in the lens light models, field galaxies and a variable

Figure 6: Distributions of the parent and selected sample for detection significance thresholds of \(3,7\) and \(10\sigma\) for the network trained on \(\mathcal{D}_{\mathcal{B}}\). \(\theta_{\rm E}\) is the Einstein radius, \(q\) is the axis-ratio of the lens mass, \(\gamma\) is the power-law slope, \(M_{\rm VIS}^{\rm L}\) and \(M_{\rm VIS}^{\rm S}\) are the apparent magnitude of the lens and source respectively, \(R_{\rm S}\) is the source radius, \(n_{\rm Sc}^{\rm S}\) is the Sersic index of the source, \(q_{\rm S}\) is the axis-ratio of the source and SNR is the signal-to-noise ratio.
PSF FWHM. However, the training data used in this work could be further improved by the inclusion of complexities that: (i) make the data more diverse (e.g an elliptical PSF model with varying ellipticities, stars and non-lensed quasars as contaminants, complex source light models, gamma-rays, CCD artefacts) and (ii) make the task of separating lensed emission from source light more difficult (for example ring galaxies). We have seen that complexities like the former require the network to become adept at ignoring contaminants, as with the network trained on \(\mathcal{D}_{\mathcal{C}}\). For the latter, we expect a more complex selection function when trained on more realistic data (as with the neural network trained on \(\mathcal{D}_{\mathcal{B}}\)).

We have confined this study to the case of single-band data, thus it is unclear how training with multi-band data will influence the selection function of the neural networks. It is possible that in this scenario, colour related selection biases may be introduced in addition to the ones listed in this paper. Moreover, we have considered only the case of the ResNet18 architecture, as it is often the best performing network for lens finding (Canameras et al., 2023). How the selection function changes for different CNN architectures and for different ML models is an interesting question which is beyond the scope of this work.

Note that our datasets are split evenly between the two classes. A real survey will have far more non-lens than lens images. This is another example of the class imbalance problem in the machine learning literature. A real sample would have at most 1 lens system for every 100 objects, and a neural network trained on a dataset with a class ratio of 1:99 could achieve a 99 per cent accuracy by simply classifying every object as a non-lens. In order to alleviate this issue, the machine learning community has developed many techniques that centre around the theme of artificially making the class ratio closer to 1:1, which can be achieved by over-sampling the minority class or under-sampling the majority class during training. Since we use simulated data, we can circumvent this issue by creating a dataset with an equal number of lens and non-lens systems. Using a realistic class ratio could alter the notions that the networks learn during training. This may further exacerbate the selection effects that we outline in this paper, but we do not study the effect of class imbalance here.

Figure 7: Distributions of the parent and selected sample for detection significance thresholds of 3, 4, 5 and 6\(\sigma\) for the network trained on \(\mathcal{D}_{\mathcal{C}}\). \(\theta_{\text{E}}\) is the Einstein radius, \(q\) is the axis-ratio of the lens mass, \(\gamma\) is the power-law slope, \(F_{\text{L}}\) and \(F_{\text{C}}\) are the lens and contaminant flux respectively and SNR is the signal-to-noise ratio.
Figure 8: The increase in the Kullback-Leibler Divergence (\(D_{KL}\)) as the detection significance threshold is increased for the network trained on \(\mathcal{D}_{\mathcal{H}}\). The blue lines show the \(D_{KL}\) that is due to sampling noise, with the increasingly shaded regions depicting 1, 2 and 3\(\sigma\) contours.

Figure 9: The increase in the Kullback-Leibler Divergence (\(D_{KL}\)) as the detection significance threshold is increased for the network trained on \(\mathcal{D}_{\mathcal{B}}\). The blue lines show the \(D_{KL}\) that is due to sampling noise, with the increasingly shaded regions depicting 1, 2 and 3\(\sigma\) contours.
## 7 Conclusions

In this paper, we quantified the selection function of the machine learning models that will likely be used to find the strong lens systems in future surveys. We focused our efforts on the ResNet18 architecture trained on three different datasets, and found that the classification task of lens finding leads to a selection bias on the parameters of the identified sample of strong lenses.

This results from the fact that lens finding neural networks are more efficient at finding lenses with certain properties, making them more likely to be above any detection significance threshold chosen for a given survey. Moreover, samples of lenses found by neural networks in the first pass are used as training data for the next iteration of lens finding networks. This might have the effect of further exacerbating the selection bias in each iteration.

We have shown that neural networks are most sensitive to the Einstein radius of the system as they preferentially select strong lens systems with larger values of \(\theta_{\rm E}\). In addition, the networks are biased towards bigger sources with more concentrated light distributions. Galaxy-quasar lens finding neural networks also show a stronger preference for more clinical lens mass distributions than those trained to find galaxy-galaxy lens systems. We also find that the networks show no preference for any values of the lens power-law slope.

Lens finding neural networks reinforce the biases introduced by the lensing cross-section. Our results clearly show that an analysis of the selection effects of lens finding neural networks is a key additional step that needs to be incorporated into any systematic attempt to find strong gravitational lenses in upcoming surveys.

## Acknowledgements

AH thanks the Max Planck Computing and Data Facility (MPCDF) for computational resources and support. AH also thanks Daniel Grun, Matteo Guardiani and Philipp Frank for useful insights and discussions. SV thanks the Max Planck Society for support through a Max Planck Lise Meitner Group, and acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (LEDA: grant agreement No 758853).

## Data Availability

The data used in this paper are available from the corresponding author on request.
