Strong Lens Discoveries in DESI Legacy Imaging Surveys DR10 with Two Deep Learning Architectures

TL;DR

This study employs deep learning models to identify strong gravitational lens candidates across nearly the entire southern extragalactic sky in DESI Legacy Imaging Surveys DR10, resulting in the discovery of 811 new candidates and a total of 3868 across multiple papers.

Contribution

First to use ensemble deep learning models for comprehensive lens candidate search in DESI Legacy Surveys DR10 covering ~14,000 deg².

Findings

Discovered 811 new strong lens candidates.

Combined with previous papers, identified a total of 3868 candidates.

Applied ensemble neural networks with visual inspection for candidate ranking.

Abstract

We have conducted a search for strong gravitational lensing systems in the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys Data Release 10 (DR10). This paper is the fourth in a series of searches (following Huang et al. 2020; Huang et al. 2021; Storfer et al. 2024, Paper I, II, & III respectively). This is the first catalog of lens candidates covering nearly the entirety of the extragalactic sky south of declination $\delta\approx +32$ deg, all of it observed by the DECam, covering $\sim$14,000 $deg^2$. We impose a $z$-band magnitude cut of < 20 in AB magnitude. We deploy a Residual Neural Network and EfficientNet as an ensemble trained on a compilation of known lensing systems and high-grade candidates as well as nonlenses in the same footprint. The predictions from these two base models are aggregated using a meta-learner. After applying our ensemble to the survey data, we exclude known candidates and systems, and use our own visual inspection portal to rank images in the top 0.01 percentile of all neural network recommendations. We have found 811 new lens candidates. These include 484 new candidates in the Legacy Surveys DR9 footprint, all parts of which have been searched for strong lenses at least once before, either by our group or others. Combining the discoveries from this work with those from Paper I (335), II (1210), and III (1512), we have discovered a total of 3868 new candidates in the DESI Legacy Surveys.

Tables

Table 1. Table 1: Object Counts in DECaLS DR10 z < z\leavevmode\nobreak\ <\leavevmode\nobreak\ 20 mag

SER		DEV		REX	EXP	Total by Region
$21,482,748$		$7,774,926$		$9,514,266$	$4,774,571$	$43,546,511$

Table 2. Table 2: Lenses Counts in the Training Sample

SER		DEV		REX	EXP
960		266		112	34

Table 3. Table 3: Strong Lens Candidates by Type

Grade	A		B		C	Total by Type
SER	65		71		357	493
DEV	9		12		100	121
REX	11		18		133	162
EXP	5		3		27	35
Total by Grade	90		104		617	811

Table 4. Table 4: Strong Lens Candidates

Grade	A		B		C	Total by Type
Within DR9	21		51		412	484
Outside DR9	69		53		205	327
Total by Grade	90		104		617	811

Table 5. Table 5: Purity for Our Neural Network Ensemble

The Tractor Type	SER		DEV		REX	EXP
Purity by Type	5		18		18	80

Table 6. Table 6: Model size, training time, and AUC comparison.

	EfficientNet	ResNet
Trainable Parameters	20,542,883	194,433
Training Time [mins]	61	56
Highest AUC	0.9987	0.9984

Equations

$${\mathcal{L}}_{CE}=-\sum _{i=1}^N{y}_i\log {\hat{y}}_i+(1-y_i)\log (1-{\hat{y}}_i)$$

$$\mathrm{TPR}=\frac{\mathrm{TP}}{\mathrm{P}}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}$$

$$\mathrm{FPR}=\frac{\mathrm{FP}}{\mathrm{N}}=\frac{\mathrm{FP}}{\mathrm{FP}+\mathrm{TN}}$$

Full text

galaxies: high-redshift – gravitational lensing: strong

1 Introduction

Strong gravitational lensing systems are a powerful tool for cosmology. They have been used to study dark matter in galaxies and clusters (e.g., Kochanek, 1991; Broadhurst et al., 2000; Koopmans & Treu, 2002; Bolton et al., 2006; Koopmans et al., 2006; Massey et al., 2010; Grillo et al., 2015; Shu et al., 2017), and are uniquely suited to probe substructure in cluster and galaxy scale lenses, as well as line-of-sight low-mass halos and test the predictions of the cold dark matter (CDM) model beyond the local universe (e.g., Vegetti & Koopmans, 2009; Huang et al., 2009; Vegetti et al., 2010; Meneghetti et al., 2020, 2023; Çağan Şengül et al., 2021; Wagner-Carena et al., 2022; Fagin et al., 2024). Searches and spectroscopic follow-up of lensed quasars (e.g., Dawes et al., 2023; Sheu et al., 2024; He et al., 2023) have enabled time-delay cosmography, offering an independent route to measure $H_{0}$. Multiply lensed supernovae (SNe) are also ideal for measuring time delays and $H_{0}$ because of their well-characterized light curves, and in the case of Type Ia, with the added benefit of standardizable luminosity (e.g., Refsdal, 1964; Treu, 2010; Oguri & Marshall, 2010), provided microlensing can be accurately characterized (Yahalomi et al., 2017). Furthermore, lens models can be constructed after the SNe have faded (Ding et al., 2021). In the last decade, strongly lensed supernovae, both core-collapse (Rodney et al., 2016; Chen et al., 2022) and Type Ia (e.g., Pierel et al., 2023, 2024; Pascale et al., 2024) have been discovered. Sheu et al. (2023) conducted a retrospective search for strongly lensed supernovae in the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys (Dey et al., 2019), and found seven promising candidates. Time-delay $H_{0}$ measurements from multiply imaged supernovae (e.g., Pierel & Rodney, 2019; Suyu et al., 2020; Huber et al., 2021; Pascale et al., 2025), combined with measurements from distance ladders (e.g., Freedman et al., 2020; Riess et al., 2021) and lensed quasars (e.g., Wong et al., 2019; Millon et al., 2020; Birrer et al., 2020), can be an important test of the tension between $H_{0}$ measured locally and the value inferred from the Cosmic Microwave Background (CMB; Planck Collaboration et al., 2020). In addition, magnified (but not multiply-lensed) SNe Ia were identified (Patel et al., 2014; Nordin et al., 2014; Rubin et al., 2018) and used to test the lens models. Finally strong lensing systems can be used to constrain the properties of dark energy (e.g., Jullo et al., 2010; Li et al., 2024), especially using systems with multiple source planes (e.g., Collett & Auger, 2014; Linder, 2016; Sharma & Linder, 2022; Caminha et al., 2022).

The introduction of neural networks to identify gravitational lens candidates in imaging surveys has been transformative and continues to show success in identifying lensing candidates (e.g., Jacobs et al., 2017; Metcalf et al., 2018; Jacobs et al., 2019a, b; Cañameras et al., 2020, 2021; Rojas et al., 2022; Shu et al., 2022; Zaborowski et al., 2023; González et al., 2025). In our recent work, we discovered over 3000 new strong lenses (Huang et al., 2020, 2021; Storfer et al., 2024, hereafter referred to as Papers I, II, and III, respectively) in the DESI Legacy Imaging Surveys (Dey et al., 2019) by using residual neural networks trained on observed images.

In this paper, we present results from our fourth search for strong lenses in the Legacy Surveys, now using DR10. We provide an overview of the observations in §2. In §3, we describe the construction of the training sample, our neural network ensemble, and our optimized visual inspection method. New lens candidates are presented in §4. We discuss our discoveries in §5 and conclude in §6.

2 Data

2.1 The Tractor

First, we briefly describe The Tractor (Lang et al., 2016), a forward modeling algorithm that performs probabilistic astronomical source detection and typing, and constructs the source catalogs for the Legacy Surveys. Source extraction is done on pixel-level data, taking as input the individual images from multiple exposures in multiple bands, with different seeing in each. The Tractor treats the fitting process as a $\chi^{2}$ minimization problem. A detected source is retained if the initial penalized $\chi^{2}$ is improved by 25111For more details, see https://www.legacysurvey.org/dr10/description/. The Tractor models detected sources as the better of either a point source (“PSF”) or round exponential (“REX”) galaxy. A detected source can be further classified as de Vaucouleurs (DEV; Sérsic index $n=4$) or exponential (EXP; $n=1$) profile over REX/PSF if such a fit improves the $\chi^{2}$ by 9. The classification becomes a Sérsic profile (SER) with an improvement in the $\chi^{2}$ by another 9 over DEV/EXP.

2.2 Observations

The Legacy Surveys are composed of three surveys: the Dark Energy Camera Legacy Survey (DECaLS), the Beijing-Arizona Sky Survey (BASS), and the Mayall $z$-band Legacy Survey (MzLS). DECaLS is observed by the Dark Energy Camera (DECam; Flaugher et al., 2015) on the 4-m Blanco telescope, which covers $\sim 9000$ deg2 of the sky in the range of $-18^{\circ}\lesssim\delta\lesssim+32^{\circ}$. As in DR9, DR10 also contains DECam data reprocessed from the Dark Energy Survey (DES; Dark Energy Survey Collaboration et al., 2016) for $\delta\lesssim-18^{\circ}$, which provides a footprint of \sim 5000\leavevmode\nobreak\deg2. The newly observed area in DR10 consists of other DECam surveys: the DELVE (Drlica-Wagner et al., 2021) and DeROSITAS (Zenteno et al., 2025) programs. The DECam surveys will hereafter be referred to in their entirety as DECaLS.222Note that the footprint of DECaLS has expanded from DR7 to DR10, as we show later in this work. BASS/MzLS are observed in the $g$ and $r$ bands by the 90Prime camera (Williams et al., 2004) on the Bok 2.3-m telescope and in the $z$-band by the Mosaic3 camera (Dey et al., 2019) on the 4-meter Mayall telescope. Together, BASS/MzLS cover the same $\sim 5000$ deg2 of the northern subregion of the Legacy Surveys. Since DR10 includes new imaging data only from DECam, and not in this northern sub-region, this work performs a lens search only in DECaLS.

The Legacy Surveys is imaged with a median $5\sigma$ PSF depths of 23.95 ($g$), 23.54 ($r$), and 22.50 ($z$) AB mag. The median FWHMs for the delivered images are: $1.29^{\prime\prime}$ ($g$), 1.18*′′* ($r$), and 1.11*′′* ($z$) for DECaLS (Dey et al., 2019).

Figure 2.2 displays the $z$-band depth map of objects classified as SER with $z<20.0$ mag in the Legacy Survey DR10. The figure also outlines the coverage of the DECaLS footprint in DR9, illustrating the overlapping regions as well as the additional areas newly observed in DR10. SER is the most common galaxy type in this magnitude regime. Table 1 shows the total counts for each galaxy type (SER, DEV, REX, and EXP) with $z$ $<\leavevmode\nobreak\ 20.0$ mag for the DECaLS region.

3 The Training Sample, Neural Network Ensemble, and Deployment

3.1 Training Sample

As in Papers I, II and III, we continue to use observed images of both lenses and nonlenses for training. Our training sample consists of 1372 lenses and high quality candidates, 869 of which come from Paper I, II, and III, with the rest selected from the following searches arranged by publication date: Zitrin et al. (2012); More et al. (2012); Anguita et al. (2012); Brownstein et al. (2012); Gavazzi et al. (2014); More et al. (2016); Carrasco et al. (2017); Diehl et al. (2017); Wong et al. (2018); Johnson et al. (2018); Sonnenfeld et al. (2018, 2019); Jacobs et al. (2019b); Petrillo et al. (2019); Jaelani et al. (2020); Cañameras et al. (2020); Sonnenfeld et al. (2020); Chan et al. (2020); Stein et al. (2022); Talbot et al. (2021); Rojas et al. (2021); Cañameras et al. (2021); Li et al. (2021); Jaelani et al. (2021); Wong et al. (2022); O’Donnell et al. (2022); Shu et al. (2022).

To select the lenses for the training sample, we inspect the $\sim 12,000$ known lenses and candidates reported in literature that overlap with the DECam ($\delta<+32^{\circ}$). We apply a uniform set of selection criteria (similar to those in Paper I, II, and III) based solely on the Legacy Surveys images for these systems, regardless of where they were originally discovered. The locations of the selected systems are shown in Figure 3.1.

We apply cuts to the training sample based on $z$-band magnitude ($z<20.0$ mag), location ($\delta$ $<+32^{\circ}$), and depth (at least three exposures333Papers I-III refer to “exposures” as “passes”. In this work, we use “exposures” instead for consistency with Legacy Survey documentation. in the $g$, $r$, and $z$ bands). The $z$-band magnitude cut was initially tested in Paper I, where we discovered that including objects fainter than the 20th magnitude yields diminishing returns for identifying lens candidates.

Paper III excluded EXP objects from the training sample since lens galaxies rarely fall under that type. However, the model in Paper III was able to perform out-of-sample classifications and correctly identified EXP candidates during deployment that meet the selection criteria for the training sample in this work. Hence, we now include all types (SER, DEV, REX, and EXP) in the training sample. Table 2 breaks down the counts per The Tractor type present in the training sample.

For the selection of nonlenses, for each type (SER, DEV, REX) in DECaLS, we bin by the number of exposures in $z$-band. This is to prevent potential bias by the neural net based on the number of exposures (see Paper II). We then select nonlenses randomly, but keep the proportionality to the lenses the same in each bin ($\sim$ 100:1). This ratio has been steadily increased from Paper I through Paper III (where this ratio was $\sim$ 33:1 in Paper III) to better reflect the expected nonlens-to-lens ratio during deployment, which is estimated to be $\sim$ 10,000:1. We cross-reference this nonlens selection with the $\sim 12,000$ known lenses and candidates found in the literature to remove potential contaminants (i.e., cutouts that contain a lens or a lens candidate). This is made possible implementing spherimatch444https://github.com/technic960183/spherimatch (Hsu, 2025), a package intended for cross-referencing catalogs in spherical coordinates. The package has a time complexity of $\mathcal{O}(N\log N)$, making it well-suited for handling large catalogs, as required in this work. In an effort to remove possible additional contamination, we deploy the trained ResNet from Paper III to this selection and inspect the images that receive a probability greater than 0.4 (the same threshold used in Paper III). After visually inspecting the recommendations from this model, we remove 9 potential contaminants from this selection. The contaminants we remove are comparable to or better than Grade D in Paper III (also, see § 3.3.1). This results in a total of 134,182 nonlenses in the training sample. This is 10$\times$ larger than training sample in our first lens search (Paper I).

The final dataset of lenses and nonlenses is randomized and divided into 70% for training and 30% for validation, while maintaining the $\sim$ 100:1 nonlens-to-lens ratio in each split. For both the training and validation catalogs, we extract image cutouts of 101 × 101 pixels, corresponding to an angular size of $\sim$ 26*′′* × 26*′′*. The exact same training and validation sets are used to train our two models described below in § 3.2.

3.2 The Neural Network Models

In this section, we describe the two neural network models used in our lens search, a residual network (ResNet; § 3.2.1) and a so-called “EfficientNet” (§ 3.2.2).

3.2.1 ResNet

The “shielded” residual network (ResNet) model, as described in Paper II, is an improvement upon the architecture presented in Lanusse et al. (2018, L18). The “shielding” layers (Szegedy et al., 2014) perform 1$\times$1 convolutions that reduce the dimensionality (i.e., the number of layers) of the output from each ResNet block.

We use the cross-entropy loss function,

[TABLE]

where $y_{i}$ is the label for the $i$th image (1 for lens and 0 for nonlens), $\hat{y}_{i}\in[0,1]$ is the model predicted probability, and $N$ is the number of images in one training step (the same as batch size, given below). During training, the ResNet attempts to minimize the loss while the Area Under the Curve (AUC) metric is used throughout training to identify the best-performing model. Figure 3.2.1 shows the loss minimization process alongside the Receiver Operating Characteristic (ROC) curve for the model that achieved the highest validation AUC during training.

The ResNet model contains a total of 194,433 trainable parameters and was trained with a batch size of 2048 ($\times$4 more than the EfficientNet, § 3.2.2), an initial learning rate (LR) of  5$\times 10^{-4}$, and a learning rate decay schedule that reduced the LR by half at the 80th epoch. The ResNet is trained on 4 NVIDIA A100 GPUs on one GPU node on the Perlmutter supercomputer at the National Energy Research Scientific Computing Center (NERSC)555A U.S. Department of Energy Office of Science user facility, https://www.nersc.gov/. . Using TensorFlow’s MirroredStrategy, training on 4 GPUs is straightforward. The multi-GPU implementation results in a speed gain of $3.7\times$ (remarkably close to the naive expectation of $4\times$): the training time per epoch is reduced from 78 to 21 seconds.

In addition to the cross-entropy loss, we further assess the performance of our trained ResNet using the ROC, which plots the True Positive Rate (TPR) against the False Positive Rate (FPR), where a “Positive” (P) indicates a lens and a “Negative” (N) a nonlens. The definitions of TPR and FPR are given as:

[TABLE]

and

[TABLE]

where True Positives (TP) are lenses correctly identified, False Positives (FP) are nonlenses incorrectly identified as lenses, True Negatives (TN) are nonlenses correctly rejected, and False Negatives (FN) are lenses incorrectly rejected.

Random classifications produce a diagonal line on the ROC curve with an AUC of 0.5, while a perfect classifier achieves an AUC of 1. During training, the validation AUC is calculated at each epoch, and the model with the highest validation AUC is selected for deployment.

3.2.2 EfficientNet Model

For the first time with respect to Papers I through III, we integrate the EfficientNetV2, referred to in this work as the “EfficientNet,” which improved the EfficientNetV1 presented in Tan & Le (2019). EfficientNetV1 comes from the search for a principled method for scaling Convolutional Neural Network (CNN) parameters (width of layers, depth of layers, resolution scaling) in order to achieve high accuracy, even with a small model that can be run on mobile devices. The primary building block for EfficientNetV1 is the mobile inverted bottleneck MBConv from Sandler et al. (2018), with a squeeze and excitation layer. The MBConv block’s input and output are thin “bottle-necked” layers as opposed to the standard residual block, which uses expanded representations in the input. In an MBConv block, the input is first expanded with a 1$\times$1 convolution which increases the number of channels. Lightweight depthwise convolution is performed on each channel of the expanded input separately, followed by a pointwise (1$\times$1) convolution, which builds features via linear combinations of channels. The output is sent into a squeeze-and-excitation block, an attention mechanism which gives weights to each feature map. The squeeze and excitation block’s output is put through another pointwise (1$\times$1) convolution, this time to reduce the number of channels back down (This is similar to the “shielding” layers we added to the ResNet model; § 3.2.1). Compared to a standard convolution, this process requires less memory and fewer computations. However, because depthwise convolutions suffer from poor parallelization, they are slower to train using modern GPUs. To improve speed, EfficientNetV2 introduces Fused-MBConv blocks, which replace the expansion and depthwise convolution with a regular 3$\times$3 convolution. While this is more expensive in terms of parameters and computations, because of parallelization, using FusedMBConv in early layers results in a $\sim$40% increase in training speed when using a GPU, and only a $\sim$4% increase in parameters compared to EfficientNetV1.

Unlike our ResNet model, the EfficientNet was pre-trained. In this work, we fine-tune the EfficientNet model using the same cross-entropy loss function as for the ResNet (Equation 1). The training results are shown in Figure 3.2.2.

The EfficientNet model contains a total of 20,542,883 trainable parameters and was trained with a batch size of 512, an initial learning rate (LR) of 3.88 $\times 10^{-4}$, and a learning rate decay schedule that reduced the LR by half at the 130th epoch. Training was also conducted on 4 GPUs at NERSC (the same as in § 3.2.1). With a per GPU batch size of 512, the training time per epoch improves from 72 to 23 seconds, or, an approximately $3\times$ speed gain, slightly less than the naive expectation of 4$\times$. This results in a 3 hour reduction of training time, from 4.7 hours to train a 160 epoch model down to 1.3 hours.

Figure 3.2.2 shows the AUC and loss progression for the training and validation sets, and the ROC curve for our best-performing model, with an AUC of 0.9987, which is comparable to the ResNet’s best AUC.

The divergent behavior of the loss and AUC curves in Figure 3.2.2 is worth noting. After epoch 50, where the model achieves the highest validation AUC, the training and validation loss continue to improve (i.e., decrease), unlike the validation AUC, which deteriorates. The validation AUC fluctuates near 0.995 before starting to decline around epoch 80. Figure 3.2.2 addresses this discrepancy by comparing the probability distributions of an EfficientNet at 50 epochs vs. at 160 epochs, showing that the 160-epoch model puts 283 more images in its smallest prediction probability bin (0, 0.04) than the 50-epoch model. Due to the significant class imbalance in the training sample (negative examples outnumber positive examples by 100:1), the model is able to achieve an overall lower loss by pushing these images into the lowest probability bin. As indicated by the lower panel of Figure 3.2.3, most of the additional 283 images came from images assigned a probability in the adjacent (low) probability bins. The vast majority of images that were not in the top or bottom bin (predicted probability between 0.04 and 0.96) were already correctly identified as TNs by the 50-epoch model (300 of 404 images were TNs). Thus, pushing them into the lowest bin improves (lowers) the loss but does not improve the AUC. On the other hand, 104 of these images are TPs, so pushing them to the lowest bin should negatively affect our performance. The lowest bin of the 50-epoch model is 99.98% TNs, but the same bin of the 160-epoch model is only 99.94% TNs. This reduction in purity will increase the loss. However, given the much larger number of TNs that got pushed into the lowest bin, the overall effect of redrawing the boundary by the 160-epoch model will improve (lower) the loss, but hurt (lower) the AUC. Similarly, some lenses that were assigned to the second highest probability bin by the 50-epoch model are now assigned to the highest probability bin by the 160-epoch model. This will improve (lower) the loss but will not affect the AUC. This logic is why the AUC is the definitive metric to select the best-performing model for deployment.

3.2.3 The Meta-learner

The ResNet and EfficientNet models, described in §  3.2.1 and § 3.2.2 respectively, generate independent predictions for each image during deployment. To combine these predictions, we employ feature weighted stacking (Coscrato et al., 2020), which consists of a meta-learner that aggregates the probabilities produced by the two base models. This meta-learner is implemented as a simple one-layer neural network with 300 nodes, designed to take the probabilities of the base models as its input and output a single probability.

The meta-learner is trained using the probability outputs of the base models on the training dataset. During validation, the meta-learner is evaluated by feeding it the predictions of the base models on the validation set. The primary objective of the meta-learner is to find an optimal combination of the probabilities of the two base models.

Given that the two base models are trained on the same dataset, their predictions are inherently correlated. This correlation constrains the meta-learner’s ability to substantially outperform simple averaging of the base models’ predictions, as illustrated by the AUCs in Figure 3.2.3. However, the meta-learner still optimizes the weighting of the base model contributions, which provides a more systematic approach than manual averaging. Figure 3.2.3 shows that the meta-learner achieves a higher overall AUC compared to each individual base model.

To further enhance the use of two or more models, strategies such as training the base models on different subsets of the data could be explored, as doing so may produce more diverse predictions by the base models. For example, Shu et al. (2022) trained models using lenses at different redshifts, which resulted in the models finding different candidates.

3.3 Visual Inspection

3.3.1 Visual Inspection Scheme

We choose a probability threshold of $\geq 0.9867$ from the meta-learner, resulting in 5680 recommendations from approximately 43 million images (the top 0.01 percentile). Out of 5680 pre-selected images, we remove 1102 recommendations present in the $\sim$12,000 known candidates and systems from literature. Our chosen probability threshold and elimination of known systems from visual inspection ensures that only the most confident predictions are considered and only new discoveries are reported in this paper. This final set of 4578 recommendations is then visually inspected to identify small blue or red galaxy/galaxies next to (typically) a red galaxy/galaxies that exhibit the following traits:

• •

are typically 1 - 5*′′* away from the red galaxy/galaxies

• •

have low surface brightness

• •

curve toward the red galaxy/galaxies

• •

have counter/multiple images with similar colors (especially in Einstein-cross like configuration)

• •

are elongated (including semi- or nearly full rings)

We assign grades based on the number of traits and the level of clarity of each trait exhibited by a recommendation. The criteria used by human graders in this work, built upon the ones used in Papers I through III, are the following:

• •

Grade A: We have a high level of confidence in these candidates. Many of them have one or more prominent arcs, usually blue. The rest have one or more clear arclets, sometimes arranged in multiple-image configurations with similar colors, also typically blue. However, there are clear cases with red arcs as well.

• •

Grade B: Giant arcs tend to be fainter than those for Grade A. Likewise, the putative arclets tend to be smaller and/or fainter, or isolated (without counter images).

• •

Grade C: They generally have features that are even fainter and/or smaller than what is typical for Grade B candidates, but are nevertheless suggestive of lensed arclets. Counter images are often not present or hardly discernible. In a number of cases, the angular scales of the candidate systems are comparable to or only slightly larger than the seeing. Therefore, to attain a higher level of certainty for some of these candidates, higher spatial resolution, deeper data, or spectroscopic observations would be required.

Co-authors JI, BK, CS, and XH performed generous initial inspection of the recommended images using the grading scheme above. CS and XH provided final inspection grades, and we report the average and difference of the results of these two graders.

3.3.2 Visual Inspection Portal

To optimize the grading process, we have developed a visual inspection portal. Given a list of recommendations, the portal displays each co-added image along with its individual grz channels with “Previous”/“Next” buttons that enable scrolling through individual recommendations. The cutout image for each recommendation is displayed with text-boxes to enter a grade and an optional comment, which are written into a CSV file. Each image is also linked to the Legacy Surveys’ sky viewer666https://www.legacysurvey.org/viewer for zoom and other options, e.g., to provide larger scale context for a possible group/cluster scale lens. Figure 3.3.2 shows the user interface of the portal.

With respect to Papers I–III, this is the first time that individual grz bands are used in addition to the co-added images to visually inspect recommendations. Typical lenses are early-type galaxies with redder spectral energy distributions (SEDs), often peaking at longer wavelengths, which makes them more prominent in the z band. Conversely, lensed sources, which are often star-forming galaxies or quasars, tend to have bluer SEDs, making their arcs or multiple images more prominent in the g band. The r band provides an intermediate perspective, capturing both the redder lenses and bluer sources with moderate contrast. By examining the separate bands, it is possible to detect subtle features that may otherwise be overlooked in a combined RGB image, as it has been noted by other lens searches such as Storfer et al. (2025).

4 Results

In this section, we first present the lens candidates identified from the visual inspection in this work (§ 4.1). We then provide a summary of all lens candidates discovered by our group in the Legacy Surveys DR7, 8, 9, and 10 (§ 4.2).

4.1 Lens Candidates in DR10

In this section, we present the strong lens candidates discovered in Legacy Surveys DR10 from this work. We deployed our trained models on $\sim$ 43 million cutouts with $\geq 3$ exposures in $g$, $r$ $z$-bands, centered on all galaxy types (non-PSF) with $z$-band mag $<20.0$ mag in the DECaLS footprint.

As described in Section 3.3.1, the top 5680 recommendations with highest probabilities are selected, and then 1102 known candidates and systems are removed, which come from: Stein et al. (2022); O’Donnell et al. (2022); Jaelani et al. (2020); Cañameras et al. (2020); Petrillo et al. (2019); Huang et al. (2020, 2021); Wong et al. (2022); Sonnenfeld et al. (2018); Wong et al. (2018); Li et al. (2021); Sonnenfeld et al. (2019); Rojas et al. (2021); Jacobs et al. (2019a); Talbot et al. (2021); Storfer et al. (2024); More et al. (2016, 2012); Sonnenfeld et al. (2020); Diehl et al. (2017); Gavazzi et al. (2014); Savary et al. (2022); Wong et al. (2022); Cañameras et al. (2021); Carrasco et al. (2017); Chan et al. (2020)

Resulting from visual inspection, we report 811 new lens candidates. These candidates contain 90 Grade As, 104 Bs, and 617 Cs. Figure 4.1 shows three examples for each of the four categories of new lens candidates. The locations in the sky of all new candidates are shown in Figure 4.1 and Table 3 breaks them down by The Tractor type.

Our search demonstrates that we can find a large number of new strong lens candidates in a footprint that has been mined for lenses before in certain sub-regions, repeatedly. Remarkably, $\sim 60\%$ of our 811 new candidates lie within the DR9 footprint, which is an area that has been the focus of multiple previous lens searches. Table 4 breaks down the distribution of these candidates between the previously searched DR9 region and the newly added area in DR10. Figure 4.1 shows a subset of the highest-quality (A-grade) candidates found within DR9. The abundance of new high-quality candidates from DR9 in this work may be owing to our use of separate bands for visual inspection (see §3.3.2), which previous searches within DR9 lack.

The purity of systems by type is shown in Table 5. These values are computed by scaling the grade distribution of new candidates to known ones, under the assumption that both new and rediscovered candidates follow the same grade distribution. This is due to the fact that only recommendations that did not match with known candidates and systems were inspected. The SER, DEV, and REX types exhibit improved purity compared to Paper III, where, respectively, one in every 23, 22, and 27 recommendations was a candidate. The only exception to this trend is galaxies typed as EXP by The Tractor, for which Paper III reported a purity of 39, whereas this work finds a purity of 80.

4.2 Lens Candidates in DR7, 8, 9, and 10

From our four searches in the Legacy Surveys, we have found a grand total of 4,869 new candidates (Paper I in DR7 among DEV and COMP in DECaLS only, Paper II in DR8, Paper III in DR9, and this work in DR10). The entire catalog of these lenses can be found on our project website777https://sites.google.com/usfca.edu/neuralens/publications/lens-candidates-inchausti-et-al-2025. The positions of all candidates on the sky are shown in Figure 4.2.

5 Discussion

This search in DR10 employs two neural network architectures: a custom ResNet (§3.2.1) and an EfficientNet (§3.2.2). Table 6 compares the sizes, training times, and AUCs of these two models. Notably, both trained in approximately the same amount of time ($\sim 1$ hour), despite the EfficientNet having $\sim$100 times more trainable parameters than ResNet. This stark difference in size, combined with the similar training durations, highlights EfficientNet’s scalability. Conversely, our custom ResNet is impressive in that it achieves virtually identical AUC performance as EfficientNet despite having far fewer trainable parameters.

Figure 5 compares the output probabilities assigned by the ResNet and the EfficientNet for new candidates graded A, B, and C. Although both models achieve similar AUCs, they exhibit slight differences in their probability distributions. The EfficientNet assigns a broader range of probabilities, extending down to $\sim 0.86$, whereas the ResNet probabilities only drop to $\sim 0.92$. This suggests that the ResNet is more confident in its predictions, while the EfficientNet is more conservative. Despite these differences, the marginal histograms show that both models concentrate candidates in the upper probability range, with all grades peaking at the highest probability bin.

While other lens searches, such as Shu et al. (2022), have implemented two different classifiers simultaneously, we advance this technique by using a meta-learner to aggregate the probabilities of our base models (see §3.2.3). This method of combining probabilities results in an overall AUC higher than that of the individual models, as shown in Figure 3.2.3. While in this case, the performance of the meta-learner is the same as the simple average of the predictions from the two models, this nevertheless provides a systematic approach to optimally combine the results from multiple models in future lens searches.

The use of individual grz bands for visual inspection (see §3.3.2) has proven to be highly effective. As demonstrated in Figure 3.3.2, different bands can help distinguish the lens from the source. This capability is particularly valuable in cases where the Einstein radius is comparable to the seeing.

6 Conclusions

We have carried out a search for strong gravitational lensing systems in the DESI Legacy Surveys Data Release 10 (DR10). This is the fourth paper in a series of lens searches in the DESI Legacy Surveys, following Paper I (Huang et al., 2020, DR7 DECaLS, DEV and COMP only), Paper II (Huang et al., 2021, DR8, DEV, COMP, and REX), and Paper III (Storfer et al., 2024, DR9, SER, DEV, REX, and EXP). We use a neural network ensemble consisting of a custom ResNet (Lanusse et al., 2018; Paper II) and an EfficientNet (Tan & Le, 2019), trained on observed lenses and nonlenses. We apply our trained models to $\sim$ 43 million non-PSF (i.e., SER, DEV, REX, and EXP, extracted and typed by The Tractor) cutout images with at least three exposures in each of the $grz$ bands and a $z$-band magnitude cut of $<20.0$ mag for the galaxy at the center of each image.

We hold a high standard in grading these candidate systems. We report 811 new lens candidates with a grade breakdown of 90 As, 104 Bs, and 617 Cs. Combining all four searches, we have found a grand total of new 3868 strong lens candidates. Grade D systems are not included in the counts, but they can be found on our project website (URL provided in §4.2).

In this search we have made improvements in four major areas: machine learning methods, training sample, neural network training hardware, and visual inspection:

• •

We have implemented a neural network ensemble composed of models with two distinct architectures whose predictions are optimally combined using a meta-learner, as opposed to a single architecture (Papers I-III).

• •

We continue to use observed images for training as in Papers I–III, but we have increased the size of our training sample by increasing the nonlens-to-lens ratio from 33:1 to 100:1, respectively, from Paper III to this work. The total training sample size thus has increased from 66,545 to 135,554 between Paper III and this work.

• •

We have trained both models on four NVIDIA A100 GPUs on a single GPU node, whereas for Papers I-III, the training was done on one GPU. Using 4 GPUs results in approximately 3.7 and $3\times$ speed gain for the ResNet and EfficientNet, respectively. Our training and hyperparameter optimization for two models with a training sample that is twice as large as for Paper III takes advantage of the speed afforded by the parallelized training.

• •

We use both the combined RGB and individual grz band images during visual inspection. This makes it easier to identify systems that have small Einstein radius and/or faint lensed arc(s), but sufficient color difference between the lens and lensed arcs. This work also uses a specialized portal that streamlines the process by enabling the inspection and annotation of candidates efficiently.

With the improvements and innovations outlined above, we find an increase in purity in this work where 1 in every 5 SER, 18 DEV, 18 REX, and 80 EXP neural network recommendations is a lens candidate.

To our knowledge, the additional footprint covered by DR10 compared to DR9 (see Figure 2.2) has not been searched for lenses. In this sub-region, we have found 327 lens candidates. However, in all other parts of the DECaLS footprint, multiple lens searches have been performed, by our group or others. And yet, we have found 484 new candidates in these parts. Therefore, systematic comparisons of lens searches to investigate the demographics (the “selection function”) of lens candidates (especially those with A and B grades) found and missed by various searches covering the same footprint would be beneficial to astrophysics and machine learning. Work in this direction has begun (e.g., Knabel et al., 2020; More et al., 2024) and the search results from this work make it clear that more such studies are needed.

This work, together with other searches, has clearly demonstrated that machine learning approaches are highly effective in discovering large numbers of high-quality strong lensing candidates. This will certainly continue to be the case for future surveys such as the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST), Euclid, and the Nancy Grace Roman Space Telescope. The discovery of strong lenses will therefore continue to accelerate and the potential for using strong lensing to tackle major astrophysical and cosmological questions has never been greater.

7 Acknowledgement

This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231 and the Computational HEP program in The Department of Energy’s Science Office of High Energy Physics provided resources through the “Cosmology Data Repository” project (Grant #KA2401022) X.H. acknowledges the University of San Francisco Faculty Development Fund. A.D.’s research is supported by National Science Foundation’s National Optical-Infrared Astronomy Research Laboratory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation.

This paper is based on observations at Cerro Tololo Inter-American Observatory, National Optical Astronomy Observatory (NOAO Prop. ID: 2014B-0404; co-PIs: D. J. Schlegel and A. Dey), which is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at The Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovacão, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Enérgeticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, the National Optical Astronomy Observatory, the University of Nottingham, The Ohio State University, the OzDES Membership Consortium the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.