Strong gravitational lensing (SGL) refers to the noticeable distortion of a distant source’s image caused by the gravitational field of massive objects located between the observer and the source (Narayan and Bartelmann, 1996). This phenomenon produces distinctive features, such as Einstein rings, arcs, or multiple images. Strong gravitational lensing systems serve as a valuable tool for cosmological research, enabling the study of the distribution of mass, particularly dark matter, within lensing objects (e.g. Narayan and White, 1988; Treu, 2010; Massey et al., 2010; Grillo et al., 2015; Meneghetti et al., 2020). Time-delay measurements between multiple lensed images provide a geometric method to constrain the Hubble constant (e.g. Refsdal, 1964; Grogin and Narayan, 1996; Suyu et al., 2017; Wong et al., 2020). Additionally, the magnification effect of SGL facilitates the observation of intrinsically faint populations at high redshifts, offering insights into early star formation and galaxy assembly (e.g. Kelly et al., 2018; Welch et al., 2022; Palencia et al., 2024). ¹

Cluster-scale strong lensing systems, where massive galaxy clusters ( M 500 ≳ 1 0 14 M ⊙ ) act as gravitational lenses, produce distinctive arcs or ring structures, offering numerous scientific opportunities (e.g. Narayan et al., 1984; Kneib et al., 1996; Kneib & Natarajan 2011; Jullo et al., 2007; Richard et al., 2010; Postman et al., 2012; Coe et al., 2019). These systems provide valuable insights into both astrophysical and cosmological aspects. Their sensitivity to gravitational effects enables the mapping the cluster mass distributions (Kneib et al., 1993; Jullo and Kneib, 2009; Coe et al., 2012; Jauzac et al., 2015), constraining the properties of dark matter or cosmology (Jullo et al., 2010; Caminha et al., 2016), and probing the mechanisms of cluster formation (Merten et al., 2015). Moreover, the magnification of high redshift background sources enables the study of galaxy formation during the earliest epochs, which extended to the dark age (Kneib et al., 2004; Coe et al., 2013). Furthermore, previous studies demonstrate that light deflection information can constrain cosmological parameters, offering a powerful geometric tool to probe Dark Energy (Jullo et al., 2010).

Galaxy cluster-scale strong lensing systems are, however, rare phenomena (Faure et al., 2008; Oguri, 2010). The advent of wide-field sky surveys, such as the Dark Energy Spectroscopic Instrument (DESI) Legacy Imaging Surveys (Dey et al., 2019), the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (LSST) (Ivezić et al., 2019), and Euclid (Euclid Collaboration et al., 2025), has increased the opportunity to discover these rare systems. With the vast amount of data generated, many strong lensing events likely remain undiscovered. However, identifying them is a challenge (Metcalf et al., 2019). Traditional detection pipelines, which rely on visual inspection (Khullar et al., 2021) or parametric morphology filters (e.g., arc length and curvature thresholds) (Lenzen et al., 2004; Horesh et al., 2005; Alard, 2006; Seidel and Bartelmann, 2007; Xu et al., 2016), scale poorly to such massive data volumes. Manual vetting becomes prohibitively time-intensive, whereas rule-based algorithms suffer from high false-positive rates due to contaminants such as spiral arms, tidal debris, and edge-on galaxies.

Recent advancements in machine learning, particularly in Convolutional Neural Networks (CNNs), have shown remarkable potential for detecting strong gravitational lenses in astronomical imaging data. Trained on large-scale simulations of lensing systems, CNNs are highly effective at identifying subtle morphological features, such as faint arcs or distorted counterimages, which often evade traditional detection algorithms (Jacobs et al., 2019). Early applications of CNN-based classifiers, such as those employed in the Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP) (Jaelani et al., 2024), the Kilo-Degree Survey (KiDS) (Petrillo et al., 2017; Li et al., 2021) and the Dark Energy Spectroscopic Instrument (DESI) Legacy Surveys (Huang et al., 2020; Huang et al., 2021; Storfer et al., 2024), have achieved recall rates of ≥ 90 % for galaxy-scale lenses. These methods reduce false positive rates by an order of magnitude compared to rule-based approaches, demonstrating their effectiveness in addressing the challenges of lens discovery in massive datasets.

In this study, we present a CNN-based methodology designed to identify cluster-scale strong gravitational lensing systems within the cluster catalog compiled by Zou et al. (2021) from the DESI Legacy Imaging Surveys. Our method employs an 18-layer ResNet architecture, chosen for its demonstrated effectiveness in extracting features from complex astronomical images, to systematically detect strong lensing signatures within galaxy clusters. To train the model, we generate realistic cluster-scale lensing images using a mock program that combines DESI legacy imaging of real galaxy clusters with synthetic lensed background sources, incorporating variations in mass distribution, redshift, and observational conditions. The initial CNN model, trained on mocked data, is applied to the DESI fields to identify strong lensing candidates. The confirmed true positives, obtained through visual inspection, are subsequently incorporated into the training set for iterative retraining of the model. our pipeline successfully identified 485 cluster-scale strong lensing candidates, where 247 of them are newly found compared to other work in DESI data.

This article is organized as follows: In Section 2, we introduce the data we used, namely the galaxy cluster catalog from Zou et al. (2021), and then describe our mock program used to produce mock lensing images to generate our training set and validation set. Section 3 provides the construction of our neural network model and the training process, then in Section 4 we present the strong lensing candidates we found. Finally in Section 5, we discuss our findings and draw a conclusion. Throughout this paper, all RGB color composite images are generated using the g , r , and z band images from the DESI Legacy Imaging Surveys, which are assigned to the blue, green, and red color channels, respectively. Our catalog of cluster lens candidates is publicly available online at: https://github.com/zkzdyzg/Galaxy-Cluster-Scale-Strong-lenses-from-the-DESI-Legacy-Imaging-Surveys.

In this section, we first describe the 540,432 galaxy clusters identified in the DESI Legacy Imaging Surveys, as reported by Zou et al. (2021), in Subsection 2.1. These clusters may serve as input for our CNN-based methodology to identify strong gravitational lenses. Next, we provide details on the mock framework developed to generate simulated lens images for training in Subsection 2.2, along with the preprocessing steps applied to the data in Subsection 2.3.

In this section, we present the methodology of our work. First, we describe the convolutional neural network (CNN) employed and the associated training configuration. Next, we outline the methods employed to evaluate network training. Finally, we describe the iterative process adopted to enhance the results.

Our lens detection system utilizes a modified 18-layer Residual Neural Network (ResNet-18) architecture (He et al., 2015), reimplemented in PyTorch 1.10 with specific adaptations for processing astronomical images. This adaptation follows a similar approach to that used in Li et al. (2021), which has demonstrated efficient performance in comparable applications. The baseline model is reconfigured to accept 3-channel 223 × 223 pixel inputs, which correspond to the dimensions of our pre-processed galaxy cluster cut-outs that combine g , r , and z -band observations through chromatic stacking. The architecture consists of five residual block groups that progressively reduce spatial dimensions from 112 × 112, following the initial convolution and max-pooling, to 7 × 7 feature maps. Strided convolutions with a stride of 2 are applied at strategic layers to preserve morphological features critical for lens identification. Each residual block contains two 3 × 3 convolutional layers with Batch Normalization (Ioffe and Szegedy, 2015) and ReLU activation (Nair and Hinton, 2010), maintaining gradient flow through identity shortcut connections. The network concludes with global average pooling, followed by a 2-dimensional fully connected layer that outputs logits optimized for binary classification between lensing systems and astrophysical false positives.

The network is trained using stochastic gradient descent with the Adam optimizer (Kingma and Ba, 2017), configured with a batch size of 32, a learning rate of η = 3 × 1 0 − 6 , and a weight decay of λ = 1 0 − 1 . The cross-entropy loss function (Goodfellow et al., 2016) is minimized: − ∑ i = 1 N y i ⁡ log y i ̂ + 1 − y i log 1 − y i ̂ . (1) In Equation 1, y i ∈ { 0,1 } is label for the i th image (1 for lens and 0 for non-lens), and y i ̂ ∈ [ 0,1 ] is the probability predicted by the network. Training is performed on a Nvidia L40S GPU and required 4 h to complete 200 epochs.

We quantify the detection capability of our gravitational lens identification system using Receiver Operating Characteristic (ROC) analysis. The network assigns a lens probability score p ̂ ∈ [ 0,1 ] to each input image, representing the confidence level of lens detection. By systematically varying the classification threshold τ across the full probabilistic range ( 0 ≤ τ ≤ 1 ) , samples are dynamically classified as positive ( p ≥ τ ) or negative ( p < τ ) , generating four classification categories: True Positives (TPs, correctly identified lenses), False Positives (FPs, non-lens objects misclassified as lenses), True Negatives (TNs, correctly rejected non-lenses) and False Negatives (FNs, undetected lens systems). The diagnostic metrics are mathematically defined as follows: TPR = TP TP + FN , (2) FPR = FP FP + TN , (3) where the True Positive Rate (TPR, in Equation 2) measures the completeness of lens recovery, and the False Positive Rate (FPR, in Equation 3) quantifies contamination from spurious detections. The ROC curve plots TPR against FPR parametrically across all possible τ values, with the Area Under the Curve (AUC) serving as a threshold-independent performance metric.

Figure 4 illustrates the variation of the loss function and AUC with its derivative across training epochs. In addition to the set of hyper-parameter configurations mentioned earlier (Case 1), we also compared and used another set of hyper-parameters configuration (Case 2), which use an initial learning rate of 1 × 1 0 − 4 , a final learning rate of 1 × 1 0 − 5 , a weight decay of 1 × 1 0 − 2 , with Cosine Annealing. It shows that in Case 1, after overcoming an initial local optimum, the AUC increases rapidly to approximately 0.93–0.95 before entering a period of fluctuation. However, the training loss begins to increase abnormally after 100 epochs, rendering the model unreliable beyond this point. In Case 2, the validation loss initially reached a relatively low level, but then gradually increased due to overfitting. Although the AUC in Case 2 reached a value of approximately 0.97, which was higher than that in Case 1, the highly fluctuating loss curve suggests that the model may be unstable. Moreover, the excessively high AUC also poses the risk of overfitting on the mock images. Therefore, we selected a deep learning model with an AUC of 0.94 from those in Case 1. The ROC curve of the best model on the validation set from one of our training sessions is shown in Figure 5. On the simulated testing data, the model achieves an area under the curve (AUC) of 0.939.

We also define Precision (Equation 4) and Recall (Equation 5) as follows: Precision = TP TP + FP , (4) Recall = TPR = TP TP + FN , (5) which are used to determine the probability threshold for human inspection. Figure 5 shows the Precision-Recall Curve (PRC) of the best model in one of our training sessions. In an ideal PRC, the curve would pass through the point (1, 1), where all thresholds would converge. The PRC is consulted to select the probability threshold for inspection, balancing precision (purity) and recall (completeness). For this work, we prioritized a high-precision point to demonstrate the model’s accuracy and minimize the number of false positives, resulting 4,000 positive images.

To improve completeness while maintaining a low false-positive rate, we adopt an iterative search procedure for strong-lensing candidates. In the initial iteration, we manually filter the top 4,000 samples with the highest scores given by the classifier trained by the mock images described in Section 2.2, and the corresponding threshold of the score is 0.99048. We then identified 102 true-positive candidates. After that, we removed them from the data pool for mining lens candidates and incorporated them into the positive sample of the training set for the subsequent training-filtering-updating loops. The details of the performance and outcomes of the classifier in each iteration are shown in Table 2.

Figure 6 shows the number of our candidates with grade A ∼ C (see Section 4.1). We repeat this iterative process until the number of detected strong-lensing candidates no longer shows significant improvement. Based on our experience, nine iterations are typically sufficient to achieve satisfactory results, since the quantity of the highest-quality grade A samples did not increase significantly in the final several iterations. This systematic refinement enriches the training data and enhances the model’s robustness in identifying strong gravitational lenses.

In this section, we present our results. Subsection 4.1 introduces the grading system for the identified candidates and the distribution of their parameters. Subsection 4.2 compares our findings with those of previous studies, and subsection 4.3 highlights specific examples that demonstrate the effectiveness of our scoring system. Finally, subsection 4.4 gives some other discussion.

In this work, we present a systematic search for candidates of galaxy cluster-scale strong gravitational lensing systems based on the galaxy cluster catalog identified by Zou et al. (2021) from the DESI Legacy Imaging Surveys DR8, covering approximately 20,000 d e g 2 of the extragalactic sky. We first designed a 3-step approach, including a strong lens simulator for creating positive samples in the initial training set, a binary classifier built upon Resnet-18 for identifying lens candidates, and human inspection to improve the purity of the outcomes further. Then, we interactively loop the classification and human inspection procedures to achieve stable completeness of the lens candidates. Subsequently, all the authors of this work visually assessed the candidates according to the features of strong lensing and galaxy clusters. As a result, we identified 485 high-confidence cluster lens candidates. The halo mass, log ( M 500 / M ⊙ ) , of these candidates ranges from 1 0 13.67 to 1 0 14.97 M ⊙ , and their brightest cluster galaxy (BCG) photometric redshift spans from 0.04 to 0.89.

The iteration starts with applying a binary classifier trained with the initial training sets (50,000 P + 50,000 N) to galaxy clusters from DESI Legacy Survey DR8. We chose a relatively conservative configuration of hyperparameters (a batch size of 32, a learning rate of 3 × 1 0 − 6 , and a weight decay of 1 × 1 0 − 1 ). The training process continues until the point where the Area Under the ROC Curve (AUC) reaches its maximum value, after which the model’s loss begins to increase abnormally, approximately 0.94. Although the AUC of the model we selected was not as high as that of the more aggressive configuration (a batch size of 32, an adaptive learning rate with an initial value of 1 × 1 0 − 4 and final value of 1 × 1 0 − 5 , and a weight decay of 1 × 1 0 − 2 , with Cosine Annealing), leading to an AUC of 0.97, more robust evaluation beyond AUC is required to avoid the potential overfitting issue and find the optimal configuration of the hyperparameters for the entire process, which is one of our focuses in the furture studies. After that, the first author visually selected candidates from the 4,000 samples with the highest scores given by the model, and added them to the original training set. The classifier is then retrained using this updated dataset. The loop will stop when the classification performance is stable (see Figure 6). We implement nine iterations in this study. At last, through visual inspection using a three-tier grading system based on arc morphology and Einstein radius scale (criteria detailed in Section 4) by all the authors of this paper, we classify 47 candidates as grade A (obvious arcs in cluster scale), 229 as grade B (probable arcs or galaxy scale), and 209 as grade C (possible arcs and galaxy scale).

In comparison to relevant studies, our strategy emphasizes the identification of cluster lens candidates in a Human-in-the-loop-like manner (Fu et al., 2024). First, relative to Storfer et al. (2024) ’s CNN-based galaxy-scale lens search, our training set is optimized based on characteristics of cluster lenses (see Subsection 2.2), which is a key reason why we found hundreds of novel cluster lens candidates. Second, compared with the human-intensive method of COOL-LAMPS (Mork et al., 2025), our CNN-based approach achieves a computational efficiency of < 10   m s per image classification, demonstrating the viability of machine learning for large-area surveys. Furthermore, our strategy achieves a simple human-computer interaction by retraining the model through iterative processes, incorporating manually identified true positive samples into the training set. This approach balances the accuracy of manual search and the speed of machine learning methods. Cross-matching with existing catalogs reveals that 247 candidates (50.9% of the sample) show no counterparts in Storfer et al. (2024) or the COOL-LAMPS project (Mork et al., 2025), which include 16 new grade A systems and 90 grade B, 141 grade C.

However, several aspects of our approach require improvement. First, there are discrepancies between the distributions of the training sets and those of the candidates (see Figure 9), which suggests that we should incorporate more realistic parametric distributions when generating training sets to enhance the classification performance of the deep learning module. Second, the current ResNet-18 achieves an AUC of ≳ 0.94 in validation tests at each step, while advanced architectures like Vision Transformers (Dosovitskiy et al., 2021) have the potential to improve the AUC to 0.999 through better arc morphology capture. Thus, it is sensible to explore novel AI techniques in future investigations. Third, our manual verification process (about 2.5 h per 1000 candidates) in the iterative process could be streamlined using active learning techniques, such as Bayesian active learning (Houlsby et al., 2011), to prioritize ambiguous cases, potentially cutting human effort by 50%.

To summarize, with a 3-step lens-finding approach, we discover 485 high-confidence cluster lens candidates in total cluster lens candidates from DESI Legacy Surveys DR8, and 247 of them are new. This dataset is valuable for examining mass assembly histories through joint strong and weak lensing analyses. Among the findings, 16 Grade A systems are suitable for putting strong constraints on the concentration-mass relation of dark matter halos and the cosmological parameters, such as the matter density, Ω m . In addition, comparing observations taken at different epochs could be beneficial in detecting lensed transients. We encourage observers to explore the best candidates for these applications in the future.