{'fraction_non_alphanumeric': 0.0267758961952936, 'fraction_numerical': 0.007220878234733524, 'mean_word_length': 4.923252279635259, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Finding quadruply imaged quasars with machine learning. I. Methods

Strongly lensed quadruply imaged quasars (quads) are extraordinary objects. They are very rare in the sky -- only a few tens are known to date -- and yet they provide unique information about a wide range of topics, including the expansion history and the composition of the Universe, the distribution of stars and dark matter in galaxies, the host galaxies of quasars, and the stellar initial mass function. Finding them in astronomical images is a classic"needle in a haystack"problem, as they are outnumbered by other (contaminant) sources by many orders of magnitude. To solve this problem, we develop state-of-the-art deep learning methods and train them on realistic simulated quads based on real images of galaxies taken from the Dark Energy Survey, with realistic source and deflector models, including the chromatic effects of microlensing. The performance of the best methods on a mixture of simulated and real objects is excellent, yielding area under the receiver operating curve in the range 0.86 to 0.89. Recall is close to 100% down to total magnitude i~21 indicating high completeness, while precision declines from 85% to 70% in the range i~17-21. The methods are extremely fast: training on 2 million samples takes 20 hours on a GPU machine, and 10^8 multi-band cutouts can be evaluated per GPU-hour. The speed and performance of the method pave the way to apply it to large samples of astronomical sources, bypassing the need for photometric pre-selection that is likely to be a major cause of incompleteness in current samples of known quads.


INTRODUCTION
Strong gravitational lenses are extremely valuable sources of information about the Universe. They provide unique information about the expansion rate of the Universe, the properties of distant sources that would be too faint (compact) to be detected (resolved), and about the distribution of mass in the Universe (Treu 2010, and references therein). Unfortunately, they are very rare on the sky, because the phenomenon requires the almost perfect alignment of a background source with a foreground deflector.
Quadruply imaged quasars are a very special case of strong lensing. They are especially valuable because of the wealth of information they provide, including, for example, three independent time delays and flux ratios. At the same time, they are especially rare because they require an intrinsically rare source (quasar) to be within the inner caustic of a foreground massive galaxy. Based on the model by , the density of quads in the sky is expected to be of order 10 −2 deg −2 with total flux brighter than i ∼ 20 (i.e. ∼ 400 in the full sky), but only a fraction of those will be resolved and identifiable in ground based wide-field imaging of the kind obtained by the Dark Energy Survey (DES, e.g., . Even though the numbers have improved considerably in the past few years, only ∼60 quads are known across the entire sky at the time of this writing.
There are two main challenges in identifying quads from ground based optical imaging data. The first challenge is the sheer volume of data one has to inspect, considering that there are about ∼ 10 8 stars and galaxies in the sky brighter than i ∼ 20 . The second challenge is that many of the quads are only partially resolved in ground based images and thus difficult to identify and separate from astronomical contaminants. In order to overcome the first challenge, many search teams rely on color preselection to reduce the number of astronomical sources. The second challenge then becomes more manageable with a combination of algorithms applied to the image pixels and visual inspection. In the end, even in the most successful cases, confirmation via spectroscopy or higher resolution imaging (from space or from the ground with adaptive optics) is needed. Considering the cost of spectroscopy and high resolution imaging, most searches so far have focused on obtaining high purity candidate lists with high confirmation rates (e.g. . The drawback of this process is that many lenses are lost along the way, as evidenced by the low completeness of searches so far.
We present a new machine learning based approach to finding quadruply imaged quasars. Machine learning techniques have been applied with success to lens searches before . While building on the work in this area, our effort differs in two main ways. First, we focus exclusively on quadruply imaged quasars, developing a realistic training set using real astronomical images from the Dark Energy Survey coupled with macro and millilensing models. Second, we avoid any need for image preselection with the goal of running our algorithm on complete flux-limited samples, which in turns requires our method to be extremely fast. This may allow us to recover quads that would have otherwise been lost in preselection steps, while retaining sufficient purity to be cost-effective for follow up. To achieve these goals we apply several new techniques and methodological improvements over previous astrophysical work, such as polar convolutions, the use of multiple networks, and attention masking. In addition to focusing exclusively on quads, tests on validation datasets suggest our method outperforms machine learners used for previous searches for lensed quasars .
In this first paper of a series we describe the method, the training set, and the results on validating datasets. A follow-up paper will present the results on a search on actual Dark Energy Survey data. The paper is organized as follows. Section 2 provides some background on the machine learning methods. Section 3 describes the training set. Section 4 describes the machine learning methods used. Section 5 describes the application of our machine learning algorithms to the problem and evaluates the performance on validating datasets. A summary is given in Section 6.

ELEMENTS OF DEEP LEARNING
In recent decades, machine learning, and particularly deep learning, have demonstrated extraordinary success in tackling a wide range of tasks related to computer vision and natural language processing, benefiting fields ranging from healthcare to the development of self-driving cars, among many others . In this section, We briefly review the basic elements of deep learning that we use in our models. A more thorough introduction to deep learning can be found in .
Deep learning builds on simple artificial neural network (ANN) models dating back to the perceptron algorithm . Loosely inspired by biological neural networks, ANNs employ computational units referred to as neurons. A single neuron receives a set of inputs, represented by vector X, either directly from the stimulus (data), or from the outputs of a preceding set of neurons. Each neuron takes a weighted sum of its inputs, w, X , adds a offset (aka "intercept" or "bias") term b, then passes this weighted sum through a function g to arrive at that neuron's activation level, g( w, X + b). The term b is often absorbed into the weight vector by simply appending a constant 1 to the dimensions of X and letting b be the corresponding weight, allowing the activation level to be written more simply as g( w, X ). When g(·) is chosen to be non-linear, this allows for more complicated models to be built than those represented by conventional (linear) models, particularly as layers are added to the neural network.
The arrangement of neurons in layers is specified by an architecture that, for each group of neurons (called layers), determines from which other groups they receive their inputs and to which groups they send their outputs. Layers in which every neuron is connected to every neuron of the next layer are called "dense" (see Figure 1). The first layer of neurons are activated directly by the data. The final layer is the output layer that generates the quantity of interest, for example, the predicted class of the input in a classification task. An ANN with more than a few layers is often called deep, in contrast to early-generation ANNs that typically employed only one or two hidden layers located between the input and output layers. Adding layers increases the learning capacity of the network, but exacerbates the difficulties of fitting or "training" the model. Figure 1. An example of a feed-forward neural network comprised of the input layer, two hidden layers, and the output layer. The information flows in one direction from each layer to the next, and the connections are called dense since every neuron in one layer is connected to every neuron in the next one.

Training the model
Neural networks are trained by adjusting the weights so that they optimally convert input data X (e.g. an astronomical image) into the desired output Y (e.g. a prediction of whether the image contains a lensed quasar). During the training, the ANN takes input data X and generates an output,Ŷ , for each training example, by propagating computations forward through the network: each neuron computes g( w, X ) using the inputs (X) given to it and the current value of the weights (inclusive of the bias term), w. The difference between the generated output,Ŷ , and the correct answer, Y , is assessed using a loss function, L(Y,Ŷ ), chosen according to the nature of the problem. Various algorithms are available to determine how the weights of the network should be adjusted in light of each error. Typically these approaches assess how each weight contributes to the error, in a process that works backwards through the model computing the gradient of the loss using the chain rule. The weights can then be updated by some fraction of the value needed to correct the response, a procedure generally referred to as backpropagation as it propagates the error signal backwards through the network. Numerous optimization methods have been proposed with varying performance. One popular choice, used here, is the adaptive momentum (Adam) algorithm , which computes individual, adaptive learning rates for the model parameters considering the first and second moments of the gradients.
Training is usually iterative, with each step using a portion of the training data (called a "mini-batch") and adjusting the weights according to learning rates that control the speed of convergence. The end goal of the optimization is to reach a global minimum, a set of values of the parameters where the loss function is minimized. However, as the number of parameters is typically very large (in the millions for all models here), we generally expect multiple local minima. Stochastic gradient descent avoids getting "stuck" in suboptimal local minima by randomly selecting mini-batches.
Because these models are so flexible, they can lead to "overfitting", where the model learns specific characteristics of the training dataset associated with particular outcomes, but which do not generalize well to unseen data. This results in large errors when the model is applied to new data. Numerous methods help to prevent this. One important approach is to add terms to the loss function that effectively penalize models with more extreme weights. This constrains the model and avoids results that depend heavily on very specific features but that might have turned out very differently in different samples. Such an penalization scheme is known as regularization and can be represented as seeking to minimize where the scalar λ controls the strength of regularization and Ω defines a regularizing functional. In the past, the first and second norms were frequently chosen as Ω and referred as L1 and L2 regularization respectively. Later, orthogonality of the weights was argued to be a desirable property since multiplication by an orthogonal matrix leaves the norm of the input unchanged. This led to orthogonal regularization, Ω(W ) = W, W T − I 1, where I is the identity matrix.

Choice of the activation function
In many practical problems the model must be made to fit a non-linear function between the inputs and outputs, calling for a non-linear choice of activation function. Until recent years, the logistic (or sigmoid) function g(x) = e x /(1 + e x ), hyperbolic tangent tanh(x) = e x −e −x e x +e −x , softmax [s(x)]i = e x i k e x k , and linear rectifier ReLU(x) = max(0, x) were among the most popular activation functions. However, deeper models require a choice of activation function that protects against the risk of having an error gradient equal to zero for many weights. To this end, functions such as parametric ReLU (PReLU) (He et al. 2015), exponential linear unit (ELU) , scaled exponential linear unit (SELU) Swish (Ramachandran et al. 2018) have become widely used. These are defined as: (2)

Convolutional neural networks
While dense ANNs with at least one hidden layer and an appropriate activation function can approximate any function, in practice ANNs can be made far more powerful, with less training data and fewer parameters to tune when they can be designed to extract the more relevant and informative features from the data. To this end, convolutional neural networks (CNNs) (see Figure 2) have proven extremely powerful in image processing applications. Figure 2. A convolutional neural network consisting of two convolutional layers and three dense layers. The first convolutional layer produces a feature map by sliding a convolution window over the input image. Within each window the layer takes a weighted sum of the corresponding pixels to produce a single output pixel. The output feature map of the second convolutional layer is flattened into a one-dimensional vector, used as an input the following fully-connected ANN.
These networks have one or more layers that, like layers of visual cortex, contain neurons whose activation's summarize key features in designated patches of the input image. Specifically, neurons in a given layer are receptive to a particular area of the input/image, known as their "receptive fields". The convolution or weighted average they perform over their receptive fields summarizes the information in it. The weights used in this convolution, which constitute a kernel or filter, can be preset or learned during training. In CNNs with multiple convolutional layers, each layer takes the previous layer's activation as its input, with wider kernels in subsequent layers, so that the receptive field corresponding to a neuron grows to cover larger receptive fields over the image.
Stacking these convolutional layers, therefore, enables learning representations of the data at different levels of abstraction. This architecture has made deep CNNs very effective in computer vision, sometimes achieving superhuman performance in classification and discrimination tasks. Units within a given layer typically share the same kernels (weights) to reduce the number of learned parameters. Nontrainable convolutional layers are often referred as pooling layers and can be used to apply basic operations such as max(·) or mean(·) that corresponds to "MaxPooling" and "AvgPooling" layers respectively.

Blocks
Blocks of neurons can be designed and connected together to develop powerful network architectures. Here we discuss three types of blocks that have proven valuable in related computer vision problems and that we consider in our own architecture.
The "inception" block, is best known from InceptionNet (also known as GoogLeNet) , and later ResNet ). This block is described in Figure 3 A). It concatenates four versions of the processed input data, each processed in parallel: the original image, two convolutional layers that use a sequence of increasingly large receptive fields, and one simpler local-averaged/smoothed version of the image.
The next block type is the residual block. A key problem with deeper networks-and part of the reason they did not emerge in earlier decades of ANNs-is the "vanishing gradient" problem: the updates to the weights computed by backpropogation factor in the gradient at each step, and a long sequence of such factors produces update values close to zero. One way to address this, employed in the ResNet architecture ) is an "identity shortcut connection", also known as a residual or skip-connection. In this configuration, the input of each learning block is added to the output before propagating to the next one (see Figure 3 B). This makes it easier to propagate information forward and backward without significant alterations and simplifies training of deep models.
The third block type we consider is the dense block in the convolutional setting, as in the DenseNet architecture . While ResNet and its residual block use element-wise addition, dense convolutional blocks combine processed inputs of one layer (here, a convolutional layer) with lower level data by concatenation instead of addition. Each layer thus receives feature maps from all preceding convolutions, within the same block (see Figure 3 C). Later blocks may then use pooling or other approaches to reduce dimensionality. DenseNet has become widely used in various computer vision problem such as image classification, object detection, and image segmentation due to superior computational efficiency and quality of the learned features.

Generative modeling and data representation
Some tasks in computer vision are "image-to-image" problems, in which we have one input image and desire to create another related image of the same size. These including denoising (removing artifacts from an input image), segmentation (estimating a set of binary masks that encode different regions on the input image), and detection (locating objects on the input image) tasks. The U-Net is one important architecture for such problems. The model was first popularized in biomedical image segmentation . In this context it takes input images (e.g. CT scans) and outputs segmentation masks that show regions of interest (e.g. malignant tumors), based on information from a labeled dataset. The U-net architecture is so named because it contains contracting and then expansive paths (see Figure 4 A). The contracting path employs a series of feature extraction blocks followed by one or more "scaling down" layers, such as a pooling layer. The expansive path concatenates the features of the same resolution, fuses the extracted features, and up-scales the image representation. The up-scaling method could be anything from a non-trainable nearest-neighbor interpolation to a trainable deep CNN. Additional skip connections between individual blocks of the contracting and expansive pathways enable easier gradient propagation. Importantly, the symmetry between two paths and the U-shaped architectures lets the network propagate context information to higher resolution layers.
Another valuable class of image-to-image architectures is the autoencoder (AE), which operate on unlabeled data and provide efficient, latent space representations. The architecture has two major parts connected sequentially: encoder and decoder (see Figure 4 B). The learning objective is to reconstruct the original image as effectively as possible, as judged by a loss function, while forcing the information The inception block applies multiple convolutional operations of different window size (shown in orange) and a pooling operation (shown in green) in parallel, instead of being restricted to a single window size, and then concatenates the extracted features together. As concatenation leads to a quickly growing size of the output tensor, 1-by-1 convolutional layers (shown in blue) are used to reduce the dimensionality. The residual block employs a skip-connection that adds the input tensor to the output tensor, which mitigates the "vanishing gradient" problem. The dense block similarly uses skip-connections, but concatenates the tensors instead of adding them. The U-Net model contains a contracting path and an expansive path that down-scales and up-scales the input image to allow simultaneous feature extraction at higher and lower image resolutions. The autoencoder model has two major parts: an encoder that maps the input image to the latent features, and a decoder that maps the latent features to a reconstruction of the input image. As the dimensionality of the latent features is often smaller than the dimensionality of the input image, the model learns a compressed representation of the data and removes noisy components.
to pass through a lower-dimensional middle block, labeled "latent features" above. The activations of neurons in this layer thus offer a lower dimensional representation of the input image, sufficient to recover the best reconstructed image possible. It is thus valuable as an unsupervised means of learning the important features of images.
One limitation of the lower dimensional representations learned by AE is that they can follow an arbitrary distribution, which may lead to situations where samples of the same class have drastically different latent features. Variational AE (VAE) is often used to alleviate this deficiency . VAE differs from AE in two ways: its latent features are selected pseudo-randomly and its loss function is extended by a penalty term. During training, each latent feature is independently drawn from a Gaussian distribution N (µ, σ), where the parameters µ and σ are taken from the output vector of the encoder {µ, log (σ)}. This exposes the decoder to a range of encoding vectors (as opposed to a single vector in AE) forcing it to map neighbouring feature vectors to the same image. The loss function is penalized by a Kullback-Leibler divergence of the latent features for a given input sample and the standard normal distribution. The penalty term ensures that the encoder refrains from producing extreme values and encourages it to evenly distribute around the center of the latent space. This leads to a continuous and orthogonal latent space, a highly desired property for data representation.

METHODOLOGY: GENERATION OF TRAINING DATASETS
We generate a training sample of lensed quasars and known sets of contaminants for training of the network. Since only a few tens of true lensed quasars have been confirmed, we cannot use them alone to construct a training set, and must generate simulated lensed quasars based on their well understood physics. Ensuring the realism of these simulated observations is essential to the ultimate generalization of our model to real data. We use a version of simct (More et al. 2016) modified for this purpose. Leaving the details to More et al. (2016), we begin by using the redMaGiC galaxy catalog ) as our parent galaxy sample. All galaxies from this sample are considered potential lenses. By using the measured redshift and magnitudes and known scaling relations, we estimate the lens mass assuming that the mass density profile follows a singular isothermal ellipsoid. We assume mass follows light to determine the centroid, ellipticity and position angle of the lens. We also include external shear to account for effects due to objects in the immediate environments of the lens galaxy. We draw sources from known luminosity functions with a certain i-band magnitude and redshift. Colors are then extracted from the quasar catalog of  by cross-matching the source i-band magnitude and redshift. Given the lens parameters and source parameters, we calculate the lensing cross-section and determine if a source would be lensed by the potential lensing galaxy such that the multiple images can be well-resolved and above the limiting magnitude. We further implemented the microlensing magnification effect by stars within the lensing galaxy which can affect the fluxes of the lensed quasars. For a given lens and source, we calculate the positions and fluxes of the lensed quasar images. The microlensing effect increases or decreases the flux of the lensed images as determined by the local convergence (κ), shear (γ) and smooth matter fraction (s = 1 − κ * /κtot) as described, e.g., by . In order to optimize computing resources, we compute microlensing magnification maps for a large number of fixed values of κ, γ and s and interpolate from this grid to real cases. Stars are assumed to have masses 1 M and the stellar density profile is assumed to follow the de Vaucouleurs profile (de Vaucouleurs 1948). We determine the convergence due to compact (stellar) population κ * in the image plane following Vernardos (2019). The resulting sample consists of about 28500 simulated lensed images of background quasar which are then added on top of the redmagic galaxies from the DES-Y3 data ). We show a few simulated lenses in the top row of Figure 5 with fainter systems on the left end and brighter on the right. There are two examples for each of the faint (i > 19.0), intermediate (18.5 < i < 18.) and bright (i < 17.5) magnitude bins. The image cutouts are 6.75 arcsec on the side where each pixel is 0.27 arcsec wide matching the DES pixel resolution.
For the training set of non-lenses, we use the spectroscopically confirmed stars from the Sloan Digital Sky Survey data and photometrically selected quasars  and blue cloud galaxies ). As we are interested in quads, we randomly draw objects from these catalogs and place them randomly around a massive galaxy which could mimic a lensed quad. About 2000 such systems are generated and this sample size is increased by a factor of five by applying rotations. We show examples of these non-lenses in the bottom row of Figure 5. As before, fainter systems are on the left end and brighter on the right. These examples correspond to the same magnitude bins as the simulated lenses in the top row. As part of the non-lens sample, we also include the same redMaGiC galaxies that were used to generate simulated lenses but without any lensing features around them. This resulted in 28,500 of simulated positive examples (lenses) and 28,500 of negative examples (contaminant galaxies).

METHODOLOGY: MODEL ARCHITECTURE AND TRAINING
Our modeling methodology involves several steps. First, we pre-process, augment, and split the data for purposes of model training and tuning. Then, we employ unsupervised learning methods to explore the data and aid in feature extraction. Next, we train a series of supervised learning models with a variety of architectures. Finally, we develop an ensemble over the resulting models. In this section we describe each of these steps in turn.

Data pre-processing and splitting
We first standardize each image so that the pixels in each image (pooled together across all four griz-bands) have zero mean and unit standard deviation. More specifically, let Iipg denote the intensity in griz-band g ∈ {1, 2, 3, 4} of pixel p in image i. Let µi and σi be the sample mean and standard deviation of {Iipg}p,g, respectively. We normalize each pixel by computing across all i, g and p, where ε > 0 is a small perturbation introduced for numerical stability. We refer to this procedure as instance normalization.
The resulting images are augmented with random rotations to ensure a rotation invariant result. We then split the data into training, validation, and testing subsets. This produces roughly 10 6 images in each of the training and validation sets and 5 × 10 5 in the test set.
The models described below are trained on the training set, with hyperparameters optimized by minimizing the prediction loss on the validation dataset. After training and L NL Faint Intermediate Bright  choosing hyperparameters, the training and validation data are combined to be re-used for training. The resulting model is then evaluated for its performance on the testing set.

Unsupervised learning
To aid in constructing features that would facilitate the supervised learning process, we begin with an unsupervised exploratory process on a subset of 100, 000 images sampled from augmented data. We first decompose the dataset using a Gaussian kernel principal components analysis (kPCA). Explaining 95% of the variance in the dataset requires 8 components. Plotting just the first two principal components ( Figure 6A) shows little evidence of the separability of these two classes. Fortunately, however, the higher-dimensional information does suggest the classes are reasonably separable. This is visualized using the two-dimensional t-distributed stochastic neighbor embedding (tSNE) of the images in Figure 6B, which uses the relative similarity of lens images with each other as compared to non-lens images to infer the possibility of successfully distinguishing these objects in a sufficiently rich non-linear feature space. Next, Figure 7A depicts the architecture of VAE which we use to generate an orthogonal, lower-dimensional latent space characterization. The encoding CNN (left side of the VAE) contains two blocks, each with two convolutional layers followed by normalization and dropout layers. The pur-pose of the first convolutional block is to extract low-level features using 5 × 5 kernels. The second convolutional block further expands the effective receptive field using 9 × 9 kernels to extract higher order features. Note that this configuration makes the effective receptive field of the last convolutional layer of 25 × 25, covering the entire input image. The drop-out layer is used to reduce overfitting by dropping a portion of the randomly selected inputs.
The subsequent block consists of two parallel dense layers that process the input tensor and map it to the latent space of dimension p. The upper layer applies the sigmoid activation to its output to estimate a feature-selecting vector sigmoid Wupperz for each training sample. By multiplying it with the latent space representation of an image z produced by the output of the lower layer, W lower z, we locally suppress irrelevant latent features of each sample: where W lower and Wupper are p × q matrices and z is a flattened image of size q × 1. The output dense layer, shown in pink in Figure 7A, linearly combines the remaining features to estimate the means and standard deviations. Finally, as shown on the right half, the process is repeated in reverse to reconstruct an image from the compressed features.
To explore the result, Figure 7B shows how the two classes are distributed along four selected latent features. As VAE gives a mean (µ) and a standard deviation (σ) for Figure 6. Visualization of a subset of the data consisting of galaxies (blue) and lensed quasars (red) after dimensionality reduction: (A) the first two components of the Gaussian Principal Component Analysis (PCA) and (B) two-dimensional t-distributed stochastic neighbor embedding (tSNE). PCA shows little evidence of the separability of lenses from non-lenses. The two-dimensional tSNE, however, suggest that the classes are reasonably separable in a higher dimensional space. each latent feature, the aforementioned figure depicts empirical distribution of the means of the latent features. We observed that although many features such as feature 1 do not discriminate contaminant galaxies from lensed quasars, some of the features such as feature 4 do separate the two classes well. To better understand the meaning of the latent features, we plot some of the weights of the lower dense layer in Figure 7C. The plots depict the underlying modes that naturally span the dataset including those that prevail among lensed quasars. To see how dimensionality of the latent space impacts the explained variance along individual bands (g, r, i, z) or their combination (griz), we apply linear PCA on the features of 1024-dimensional VAE models that were trained to reconstruct a single band or the entire image. The results in Figure 7 D suggest that 16-dimensional space covers 95% of the variance in griz-bands (around 12 latent variables for R), while 128 features correspond to 99% of the variance (around 64 latent features in R). This shows that VAE models can learn meaningful features and describe the dataset with a 128 latent variables. We later employ an encoding portion of the VAE models for feature extraction and subsequent image classification.

Supervised learning
This section provides details for the construction and fitting of several neural network architectures, before describing how they are combined in an ensemble.
The unsupervised learning explored in Section 4.2 indicates that lensed quasars have distinctive geometric features that could be utilized in their detection. We frame lens detection as a binary classification problem in which "one" corresponds to a sample with a lensed quasar and "zero" corresponds to anything else. A classifier takes an image as an input and outputs a binary label. In practice, the output layer has two units, whose activations we can refer to as z0 and z1. To these we apply the softmax function to obtain quantities we interpret as probabilities, i.e., y1 = P(Y = "lens") = exp(z1)/(exp(z0) + exp(z1)), witĥ y0 = P(Y = "not lens") = 1 − P(Y = "lens"). We found empirically that training with the softmax activation results in a better generalization compared to training with the sigmoid activation.
For the loss function, we employ the binary crossentropy, also known as the negative log-likelihood, between a class label y ∈ Y and the predicted probabilities (ŷ1,ŷ0), L(y,ŷ1,ŷ0) = −y · log (ŷ1) − (1 − y) · log (ŷ0). Moreover, binary cross-entropy leads to a more consistent gradient propagation as the log-terms mitigate exponential behavior due to gradient saturation for extreme values.
Instead of the conventional batch normalization, we use instance normalization introduced in (3). While instance normalization makes the samples of both classes statistically indistinguishable in terms of the pixelwise mean and variance, it has several important advantages. First, it limits the range of the values that tensors can take, preventing saturation of hyperbolic activation functions, and improving gradient propagation. Second, it improves generalization by avoiding overfitting on synthetic data, as the gap between the simulated and the real objects is often the major issue in problems that rely on statistical models trained on synthetically-generated data. In particular, it reduces the effect of the outliers that could cause significant covariate shift in hidden layers of the ANN.
To penalize overfitting, we extensively use dropout layers. The exact order of the hidden layers, activation functions, normalization layer, and dropout layers is typically optimized for every problem. We empirically found that the combination of activation layers followed by normalization and dropout layers achieve a better bias-variance trade-off. that despite nearly identical statistical properties in galaxies and lensed quasars across the majority of learned latent variables, some of them describe features that can be used to distinguish between these two classes; (C) a visualization of a few learned features of the lower dense layer of the VAE model generated by putting the learned weights into square shapes of the same size as the output of the previous convolutional layer; (D) empirical relationship between the explained variance and the number of principal components of the latent variables of 1024-dimensional VAE suggests that a 128-dimensional latent space can be sufficient to describe the dataset.
To additionally limit overfitting, we introduce zero-mean Gaussian noise at the input of each model, which adds random noise to the input images at every training step. We set the noise covariance matrix to σ 2 I, where I is the identity matrix and σ 2 = 0.062, with the value determined by hyper-parameter optimization. Note that this is equivalent to additional data augmentation performed simultaneously with training.
Models are trained using the Adam algorithm . Adam uses adaptive learning rates for every model parameter, enabling faster convergence. We also employ the scheduled (staircase exponential) and triggered decrease of the learning rate. The latter changes the step size when the loss rate of decay is below a certain threshold.
To further combat overfitting, we employ early stopping based on the validation loss. More specifically, after every iteration on the training set, we evaluate the loss function on the validation set which includes data not used for optimizing model parameters. We stop once the validation, rather than the training, loss has converged.
Besides trainable parameters, each model has hyperparameters, a small set of values that define model architecture and training dynamics. These include model depth (number of hidden layers), model width (number of convolutional filters and dense units), model resolution (size of convolutional kernels), and the choice of the activation function, regularization strength, and learning rate. We use the recently proposed Hyperband algorithm , an efficient bandit-based method for hyperparameter optimization. It allocates computational resources to as many configurations as possible and throws out those that show poor performance over time until a single configuration remains. This method maximizes the number of tested configurations and results in a more efficient resource utilization compared to the grid search, random search, or Bayesian optimization.

Prior methods
Numerous prior works proposed using the traditional computer vision models such as InceptionNet, ResNet, and DenseNet for lens searches. These architectures have demonstrated exceptional performance on traditional computer vision problems, but the vast number of trainable parameters they require can result in either poor generalization (overfitting) or weak convergence (underfitting). Pre-training these models on the ImageNet, a large dataset of ordinary images collected for object recognition, degrades the performance further. We hypothesize that this is due to drastic differences in visual features between the lensed quasars and ImageNet samples. Fortunately, several earlier projects have sought to detect lenses using various deep CNN architectures. We describe these here before proposing two novel architectures. Previously proposed CNN models include CMU DeepLens , LensFlow , CNNS , and several others . We reproduce and retrain LensFlow and CNNS on our data for comparison. LensFlow is based on a classical architecture comprised of three convolutional layers, maximum pooling layer and four dense layers all connected in series. Drop-out layers in-between the dense layers mitigate over-fitting. Unusually, LensFlow applies hyperbolic tangent to the outputs of the convolutional layers and ReLu activation function to the output of the dense layers. CNNS combines four convolutional layers into blocks of two layers each with three subsequent dense layers. Each convolutional block is followed by pooling operations, halving the dimensionality of the feature maps the output. The activation function is ReLu, and drop-outs are again used to mitigate over-fitting.

NaiveNet
Inspired by these early examples, we first propose a CNN involving several convolutional layers followed by dense layers, which we refer to as NaiveNet, illustrated in Figure 8. The number of layers, number of dense units (or convolutional kernels) within each layer, and size of the convolutional kernels respectively, are optimized together with other hyperparameters. We intentionally avoid specific constraints on the architecture and, instead, optimize performance through a hyperparameter search. Next, we introduce NaiveNetV2, replacing convolutional and dense blocks with their residual alternatives. Each residual path has two consecutive layers: the first one with a chosen non-linearity and the second one with a linear activation function.

Polar convolution
Next, we propose an alternative approach to feature extraction. Since we can ensure positioning of the objects at the center of the griz images and given that the gravitational field originates from the center of mass of the objects, a polar coordinate system is a natural choice. We define the corresponding 2D polar convolution similarly to its rectangular counterpart: Figure 9. Visualization of the proposed polar convolution operation: each feature map gets converted to higher resolution polar coordinates where a regular rectangular convolution is applied and then translated back to Cartesian coordinates. Effectively, the operation has a convolutional window of an angular sector of annuli of different radius (shown in red) with smaller window size at the center of the image.
where T (·) translates an input tensor from Cartesian to a polar coordinate system. It first transforms a fixed size rectangular tensor z into a scatter field of points z * with polar coordinates {ρ, θ} = { (x 2 + y 2 ), arctan (y/x)}. Then, it estimates a smooth rectangular tensor s using biharmonic spline interpolation: where g is Green's function, pi = [xi, yi] is a vector pointing at the corresponding scatter point with value z * i , and the weight vector α can be found by solving g · α = z * . An illustration of the proposed operation is shown in Figure 9. Beyond fitting the underlying geometrical structure on the images, polar convolution yields a heterogeneous receptive field size and a particularly sparse connectivity matrix by focusing at the central part of the image, which is useful as the corners of the images are expected to be less informative, with more noise.

Visual attention masking
We train an additional U-Net-like model that produces a binary mask representing regions of lensed quasars (named as AttnCNN in Figure 10 D). The goal of this auxiliary masking network is to guide feature extraction by bringing attention to regions likely to have important features. We found empirically that thresholding image pixels with 22.5% of their maximum intensity creates meaningful binary masks as the brightest regions are the most informative. We tried several loss functions for training the AttnCNN and found that mean absolute error leads to the most stable results. It was beneficial to train the entire assembly simultaneously in an end-to-end fashion. At each optimization step i (batch size of n), we update the parameters of the AttnCNN model Ψ(·) to minimize the discrepancy between its outputs Ψ(Xi) and the binary masks Mi: At the same time, we independently update parameters of the feature extracting model Λ(·) (named as FeatureCNN in the figure) and the binary classification model Π(·) by backpropagating the gradient of the cross-entropy loss Li Π Ψ(Xi) Λ(Xi) , Yi . This way AttnCNN directly affects the performance of the feature extracting CNN and the classifier, which forces FeatureCNN to adjust its learned parameters.

Argus: ensemble of models
One of the challenges we encountered was overfitting stemming from the gap between the simulated and real data. Particularly, we observed a significantly lower performance during testing on real samples generated from the DES, despite the imposed regularization and optimal synthetic testing results. Given the limited size of the real dataset, we found it impractical to directly bridge the gap by finding a projection between real and synthetic samples. Moreover, since the confirmed objects do not represent the entire range of possible lenses, additional training on them would further amplify the bias. Instead, a more practical and efficient method to improve generalization is model stacking (sometimes called blending), which combines the outputs of pre-trained models. We employ this technique by blending classifier models and VAEs (see Figure 11). The upper branch maps the hidden feature spaces of the AttnNet and NaiveNetV2 based on rectangular and polar convolutions to another 128-dimensional space. It employs two dense layers with ReLu non-linearity to fuse a 512-dimensional vector from the At-tnNet with a pair of 128-dimensional vectors from the rectangular and polar NaiveNetV2: Φs : R 512+128+128 → R 128 . The lower branch mixes the latent spaces of VAE models that were pre-trained to reconstruct specific components on the images including lensed quasars (VAELQ), contaminant galaxies (VAEG), underlying source quasars (VAEQ), and deflectors (VAED) that represent the strength of the gravitational field of the lensing galaxy. We found that reconstruction on the quasars and deflectors achieved the best results when conditioned on the output of the VAELQ and VAEG respectively. This is implemented by stacking the second order VAE model on top of the decoder and using the mean component of the latent features of the first encoder. We choose a 128-dimensional latent space since it is able to converge to a global minimum of the loss function and it was earlier suggested by the results of PCA to keep the explained variance ratio above 99%. However, because galaxies and lensed quasars have an identical distribution along the majority of the latent space axes, we apply a stronger dimensionality reduction: Φu : R 512+512 → R 64 . Note that the standard deviations estimated by the encoding CNN models carry information about the noise in the images and the associated confidence levels. Since µ and σ often have different scales and distributions we first blend them separately and then combine new intermediate feature spaces. The output layer linearly combines the features of the upper and lower branches to classify the input samples.
To train the entire end-to-end model we split the training dataset into two halves. We train each individual component (supervised and unsupervised) on the first half. Then, after stacking, we train the blending dense layers using the second half of the data. This scheme allows us to minimize additional overfitting as the blending layers learn how to combine the feature spaces of the pre-trained models on previously unseen samples. Moreover, to foster generalization we use normalization and dropout layers as well as combination of L1 and L2 ("elastic net") regularization. The corresponding hyperparameters of the model and its components are selected using the Hyperband search described earlier.

Training and testing setup
After training these models we tested their performance on three datasets that were not used during the training. First, the simulated objects is a conventional "testing" dataset that was sampled from the data generated for model training as described in Section 3. It contains around 250, 000 lensed quasars and 250, 000 contaminant galaxies. Second, the past candidates dataset contains 20 confirmed lensed quasars and 108 confirmed non-lenses that were previously proposed as candidates for spectroscopic follow-up within the STRIDES collaboration (as collated by one of us, C.L.). Finally, the mixed dataset contains 1216 objects randomly sampled from the simulated objects and past candidates. We augmented the samples from past candidates via random rotations and mirroring to achieve a roughly 1:1 ratio. Figure 11. Argus: an end-to-end ensemble of deep neural networks. The proposed method includes feature extraction elements of pretrained models: a 512-dimensional output from the first dense layer of AttnNet, and a 128-dimensional vector from the second dense layer of NaiveNetV2. These features are concatenated together and then reduced to 128-dimensional space through two dense layers. In addition, the model takes 128-dimensional latent features of four serially-connected VAE models trained on reconstructing images of lensed quasars (LQ), galaxies (G), delensed quasars (Q), and deflector fields (D) from the blended input image. As the encoders produce both means and standard deviations of these features, they are first processed independently and then combined into a 64-dimensional feature vector. The final output is produced by a dense layer that classifies samples based on 192 features. We trained each model on a computer with a single GPU (Nvidia GeForce GTX 1080 Ti). A full training cycle of a single model with hyperparameter optimization took between 20 to 180 hours depending on the model. In addition to our proposed models, we reproduced two previously reported models, CNNS and LensFlow, retraining these on our data. In doing so we followed the design choices as reported, but adjusted learning rates and the corresponding parameters (decay rate, learning rate schedule, etc.) to maximize performance.

Performance
As a first performance metric we consider the area under the receiving operating curve (AUROC). This can be interpreted as the probability that the model ranks a random positive example more highly than a random negative example.  Attention masking (AttnNet), however, led to a substantial improvement in performance with an AUROC of 0.940 on the simulated objects and 0.774 on the past candidates. AttnNet was the most conservative in detecting lenses and had the highest AUROC on the past candidates dataset, which is driven by the built-in feature selection and stronger orthogonal regularization. Finally, the proposed ensemble model (Argus) reached a nearly perfect AUROC score of 0.997 on the simulated data, but had a slightly lower score on the past candidates compared to AttnNet. On the mixed data, Argus achieved an AUROC of 0.894, the highest among all models.
ROC curves of the AttnNet and Argus models on the "mixed dataset" are shown in Figure 12. The Argus ensemble achieves a nearly perfect true positive rate while keeping the false positive rate under 0.3. The AttnNet demonstrates true positive rate of about 0.8 for the same false positive rate under 0.3. Figure 13 depicts how precision (the fraction of true lenses among the predicted lenses) and recall (the fraction of correctly identified lenses among true lenses) of two best performing models change with AB magnitude of the whole system (obtained by summing the flux over the entire cutout; for comparison, the 10-σ limiting magnitude of the data for point sources is 23.34). To produce these plots we grouped the samples of the mixed dataset into five bins of equal width in AB magnitude and evaluated the metrics for each group. We found that as AB magnitude increases the proposed models start suffering from a noticeable drop in precision. This is expected because of the declining signal to noise ratio of the images and false positives become more difficult to identify as such. Remarkably, however, recall is uniformly close to 100% suggesting that our method should achieve high completeness even when the signal to noise ratio is low. Figure 14 shows objects misclassified as lenses (false positives) by the Argus Ensemble model. Those in the top row were also selected by previous searches and were considered credible enough to warrant telescope time for investigation, so it is not catastrophic that they mislead Argus. Those in the bottom row were simulated false positives and indeed look like plausible lenses to a human classifier, so it does not seem catastrophic that they mislead Argus as well. At some level, a small number of false positive is inevitable due to the finite amount of information in the imaging data used for the search. Those false positives will need additional data to be resolved (spectroscopy and/or higher resolution images).  We discuss the practical implications of this and other limitations to be kept in mind when applying this method to search for quads in the next section.

Practical Considerations
Considering that telescope time is costly and limited, we wish to estimate how effectively lists of candidates produced by Argus or AttnNet can prioritize search time. The two key aspects of performance in this context are precision-the proportion of candidates examined that prove to be lensesand recall, the fraction of true lenses that are identified in a given catalog. Figure 15 shows the estimated precision and recall as we adjust the number of selected candidates ranked  by the models. Both Argus and AttnNet look promising in the sense that over a wide range of the number of candidates that might feasibly be considered, the precision rate remains quite high, well over 80% until the number of candidates begins to approach the actual number of lenses in this sample (400). Naturally, the recall can at best grow linearly with a slope of 1, and the observed recall is not far below this. If these results generalize to images outside of those used in this paper, this plot would suggest that a search using our methods would be very efficient, yielding a complete sample of quadruply imaged quasars with a relatively modest fraction of contaminants. A follow-up paper will put our methods to further test by applying them to new DES datasets.

SUMMARY
In this work, we proposed and trained novel deep learningbased methods for detecting quadruply lensed quasars. We investigated several ANN architectures and compared them with previously proposed LensFlow and CNNS models. Overall, our newly developed tools and in particular the ensemble model improve on the previously proposed algorithms ) in terms of precision and recall, although prior work was not tailored to the detection of quadruply imaged quasars, and therefore the comparison is not direct.
A major strength of our method is its ability to process 10 8 objects in a single day on an ordinary hardware setup with a single GPU. This makes it suitable for discovering lenses across wide areas from surveys such as DES with minimal pre-screening of candidates, reducing the risk of discarding lenses. In a follow-up paper we will carry out such a search.
Two particularly novel and powerful elements of our approach are the use of polar convolution, which is well tuned to detecting circularly-arranged features at different scales, and attention masking. In the future, yet other ideas for feature extraction approaches, likely motivated by expert knowledge, may further improve performance.
Finally, the challenges that remain are not only algorithmic, but in the acquisition of training data, real or simulated. That said, one area in which further progress could be made is in taking account of additional data, such as bands beyond griz, or using time series rather than single epoch images.

DATA AVAILABILITY STATEMENT
All data and code required to reproduce these results and apply the trained model will be made publicly available upon publication at https://des.ncsa.illinois. edu/releases/other/paper-data.

ACKNOWLEDGMENTS
TT acknowledges support from the National Science Foundation through grant NSF-AST-1906976, from NASA through grants HST-GO-15320 and HST-GO-15652, and