# Ghost imaging with non-Gaussian quantum light

Dongyu Liu∗ State Key Laboratory for Mesoscopic Physics and Collaborative
Innovation Center of Quantum Matter, School of Physics, Peking University,
Beijing 10087, China  Mingsheng Tian  These authors contributed equally to
this work. State Key Laboratory for Mesoscopic Physics and Collaborative
Innovation Center of Quantum Matter, School of Physics, Peking University,
Beijing 10087, China  Shuheng Liu  State Key Laboratory for Mesoscopic Physics
and Collaborative Innovation Center of Quantum Matter, School of Physics,
Peking University, Beijing 10087, China  Xiaolong Dong  State Key Laboratory
for Mesoscopic Physics and Collaborative Innovation Center of Quantum Matter,
School of Physics, Peking University, Beijing 10087, China  Jiajie Guo  State
Key Laboratory for Mesoscopic Physics and Collaborative Innovation Center of
Quantum Matter, School of Physics, Peking University, Beijing 10087, China
Qiongyi He<EMAIL_ADDRESS>State Key Laboratory for Mesoscopic Physics
and Collaborative Innovation Center of Quantum Matter, School of Physics,
Peking University, Beijing 10087, China  Collaborative Innovation Center of
Extreme Optics, Shanxi University, Taiyuan, Shanxi 030006, China  Haitan Xu
<EMAIL_ADDRESS>Shenzhen Institute for Quantum Science and Engineering,
Southern University of Science and Technology, Shenzhen 518055, China  School
of Physical Sciences, University of Science and Technology of China, Hefei
230026, China  Zheng Li<EMAIL_ADDRESS>State Key Laboratory for
Mesoscopic Physics and Collaborative Innovation Center of Quantum Matter,
School of Physics, Peking University, Beijing 10087, China  Collaborative
Innovation Center of Extreme Optics, Shanxi University, Taiyuan, Shanxi
030006, China  Peking University Yangtze Delta Institute of Optoelectronics,
Nantong, China

###### Abstract

Non-local point-to-point correlations between two photons have been used to
produce “ghost” images without placing the camera towards the object. Here we
theoretically demonstrated and analyzed the advantage of non-Gaussian quantum
light in improving the image quality of ghost imaging system over traditional
Gaussian light source. For any squeezing degree, the signal-to-noise ratio
(SNR) of the ghost image can be enhanced by the non-Gaussian operations of
photon addition and subtraction on the two-mode squeezed light source. We find
striking evidence that using non-Gaussian coherent operations, the SNR can be
promoted to a high level even within the extremely weak squeezing regime. The
resulting insight provides experimental recipes of quantum imaging using non-
Gaussian light for illumination.

##  I Introduction

Quantum light is the basis for developing technologies of quantum imaging
Lugiato  _et al._ (2002); Gatti  _et al._ (1999); Brambilla  _et al._ (2008);
Delaubert  _et al._ (2008); Brida  _et al._ (2010, 2011a); Santos  _et al._
(2008); D’Angelo  _et al._ (2004), super-resolutionBoto  _et al._ (2000);
Giovannetti  _et al._ (2009); Santos  _et al._ (2003); Sciarrino  _et al._
(2008), etc. Especially, ghost imaging Bussgang  _et al._ (2012) is a
representative technique based on the second-order correlation of quantum
light Klyshko (1988a); Pittman  _et al._ (1995), where entangled photon pairs
are used as the light source. One is able to obtain the image of an object
entirely relying on the correlation between the entangled photon pairs,
without directly viewing the object with a spatially resolving camera. As
conceptual and practical interests are attracted to ghost imaging, the
theoretical description Brida  _et al._ (2011b); Santos  _et al._ (2008);
Erkmen and Shapiro (2009); Valencia  _et al._ (2005); Thiel  _et al._ (2007);
Lopaeva and Chekhova (2010); Thiel  _et al._ (2007) and improvement Brida  _et
al._ (2011b); Ferri  _et al._ (2010); Meyers  _et al._ (2011) of this
technique has been studied. The signal-to-noise ratio (SNR), one of the main
parameters quantifying the quality of a ghost imaging system, characterizes
how well the image is distinguished from the background Brida  _et al._
(2011b). The two-mode squeezed state (TMSS), produced by spontaneous
parametric down-conversion (SPDC) Klyshko (1988b); Shih (2003); Smithers and
Lu (1974); Schumaker and Caves (1985); Barnett and Knight (1987), is commonly
used as the quantum light source in ghost imaging Pittman  _et al._ (1995). It
was shown that the SNR of ghost imaging with TMSS is a monotonically
increasing function of the squeezing parameter of the state, associated with
the degree of entanglement. The advantage of ghost imaging with TMSS has also
been shown over ghost imaging with thermal light Brida  _et al._ (2011b).
However, the SNR approaches saturation with the increase of the squeezing
parameter Brida  _et al._ (2011b), and due to the limitation of experimental
technique, it can be difficult to prepare TMSS in the strong squeezing regime.

Here we propose an experimental scheme to improve the SNR of ghost imaging
with engineered non-Gaussian quantum light by coherently adding and
subtracting photon from the signal mode of TMSS, i.e. $(t\hat{a}_{{\rm
s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\mathrm{TMSS}}\rangle$. While the two-mode
Gaussian states are frequently used as entangled photon source, it has been
demonstrated that non-Gaussian operations of photon addition and subtraction
Lee  _et al._ (2011); Zavatta  _et al._ (2004); Parigi  _et al._ (2007) can
enhance the performance of various applications via distillation of the
entangled states. Although non-Gaussian light has proven to introduce
remarkable advantage to various applications in quantum information
processing, such as secure quantum communications, quantum computing and
continuous variable quantum teleportation Lee  _et al._ (2019); Liao  _et al._
(2018); Ra  _et al._ (2020); Lee  _et al._ (2011), its application beyond the
area of quantum information, such as quantum imaging, remains largely
unexplored. We calculate the SNR of the corresponding ghost imaging system
with coherent operation ($t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger}$),
and give the maximal SNR by optimizing the parameter $r$ of the coherent
operation. We show that the SNR can be significantly enhanced by coherent
superposition operation in the weak squeezing regime, whereas the photon
subtraction or addition operation performs better in the strong squeezing
regime. We also provide an analytical understanding on the mechanism how non-
Gaussian local coherent operation can enhance the SNR of ghost imaging.

##  II Experimental scheme

Figure 1: Schematic of ghost imaging with non-Gaussian light. The upper part
is the diagrammatic sketch of ghost imaging; and the lower part is the
experimental scheme to implement the non-Gaussian coherent operation
$t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger}$ on the signal mode (signal1)
of TMSS. BS1 and BS2 are beam splitters with transmissivities $t_{1}$ and
$t_{2}$, respectively. PD1 and PD2 are photodetectors. The detection of only a
single photon at PD1 or PD2 heralds the success of a coherent operation. The
output mode signal3 after the coherent operation and idler1 form a two-mode
output state, which is then used for ghost imaging.

The diagrammatic sketch of ghost imaging Pittman  _et al._ (1995); Bussgang
_et al._ (2012); Brida  _et al._ (2011b) is shown in the upper part of FIG. 1.
Entangled pairs of signal and idler photons are generated in the process of
spontaneous parametric down conversion (SPDC) in a BBO crystal (BBO1) Klyshko
(1988b); Shih (2003); Smithers and Lu (1974); Barnett and Knight (1987). In
the regime of quantum continuous variables, it is known that the two-mode
squeezed state (TMSS) of squeezing degree $s$, $|{\rm
TMSS}\rangle=\frac{1}{\cosh s}\sum_{k=0}^{\infty}\tanh^{k}s{|{k_{\rm s}k_{\rm
i}}\rangle}$, corresponds to a typical Gaussian light Schumaker and Caves
(1985); Barnett and Knight (1987); Serafini (2017). For ghost imaging, only
the signal photons are sent to the imaged object, and transmitted ones are
collected by a detector without spatial resolution (bucket detector), while
the idler photons are sent to a spatially resolving detector (CCD). The image
of the object is then reconstructed by measuring the function $S(x_{j})$
Bussgang  _et al._ (2012), where $x_{j}$ represents the position of the pixel
$j$ of the CCD. $S$ usually has the form

$S(x_{j})=f(\langle{\mathbb{N}}_{\rm s}^{p}N_{\rm i}^{q}(x_{j})\rangle),\quad
p,q\geq 0,$ (1)

involving the correlation function $\langle{\mathbb{N}}_{\rm s}^{p}N_{\rm
i}^{q}(x_{j})\rangle$ of the total number of photons collected at the bucket
detector, $\mathbb{N}_{\rm s}$, and at the $j\mathrm{th}$ pixel of the CCD,
$N_{\rm i}$ Brida  _et al._ (2011b). Different correlation functions can be
used for ghost imaging. Here we focus on the protocol using the covariance
$\mathrm{cov}(x)=\langle{\mathbb{\mathbb{N}}}_{\rm s}N_{\rm
i}(x)\rangle-\langle{\mathbb{N}}_{\rm s}\rangle\langle N_{\rm i}(x)\rangle$ as
the correlation function Bussgang  _et al._ (2012); Brida  _et al._ (2011b).

To simplify, we consider an object having only binary levels of transmission,
$T=0$ and $T=1$. The SNR of ghost image can be defined as the ratio of the
mean contrast of the correlation functions inside ($T=1$) and outside ($T=0$)
the object profile to the mean relative fluctuation Brida  _et al._ (2011b);
sup :

$\mathrm{SNR}=\frac{|{S_{\mathrm{in}}-S_{\mathrm{out}}}|}{\sqrt{{\delta^{2}(S_{\mathrm{in}}-S_{\mathrm{out}})}}}.$
(2)

In order to improve the SNR of the ghost imaging system, we modify the
Gaussian light source of ghost imaging by applying non-Gaussian operations to
one of the local modes of TMSS. In FIG. 1, the lower part shows the
experimental scheme to implement the non-Gaussian operation, i.e. coherent
operation $t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger}$ with
$|t|^{2}+|r|^{2}=1$ on the signal mode (signal1) Lee and Nha (2010). In
general, both photon subtraction ($\hat{a}$) and addition
($\hat{a}^{\dagger}$) can be implemented in quantum optics experiments. If an
input state, e.g. signal2 in the lower part of FIG. 1, is mixed with an
accessory input of a vacuum state by a beam splitter (BS1) with a
transmissivity $t_{1}\simeq 1$, the detection of a single photon at the
accessory output port (signal4) heralds that a photon is subtracted from the
input state Wenger  _et al._ (2004). On the other hand, when the input state,
e.g. signal1 in the lower part of FIG. 1, is squeezed with an idler mode in a
vacuum state in the process of parametric down conversion with a squeezing
parameter $s_{0}\ll 1$ in a BBO crystal (BBO2), the detection of a single
photon at the idler output port (idler2) heralds that a photon is added to the
input state Zavatta  _et al._ (2004). Combining these two operations, if the
which-path information on the detected single photon is erased by an
additional beam splitter (BS2) with transmissivity $t_{2}$ (FIG. 1), the
operation $t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger}$ can then be
conditionally implemented. The detection of a single photon only at PD1 or PD2
heralds the success of a coherent operation. The parameters $r$ and $t$ of the
coherent operation determined by the squeezing parameter $s_{0}$ and the
transmissivity of BS1 and BS2 are Lee and Nha (2010)

$\displaystyle t=-t_{2}\frac{r_{1}}{t_{1}},\quad r=s_{0}r_{2},\qquad\text{when
a single photon is detected only at PD1}$ (3) $\displaystyle
t=r_{2}^{*}\frac{r_{1}}{t_{1}},\quad r=s_{0}t_{2}^{*},\qquad\text{when a
single photon is detected only at PD2},$ (4)

which is elaborated in the supplementary material sup . The parameters $t$ and
$r$ are then normalized by a factor of
$1/\sqrt{|t_{2}r_{1}/t_{1}|^{2}+|s_{0}r_{2}|^{2}}$ for Eq. 3 and
$1/\sqrt{|r_{2}r_{1}/t_{1}|^{2}+|s_{0}t_{2}|^{2}}$ for Eq. 4 so that
$|t|^{2}+|r|^{2}=1$. For a given BBO2 with a squeezing parameter $s_{0}$
(which is usually much smaller than 1), we can adjust the parameters $t_{1}$
and $t_{2}$ by choosing proper beam splitters to change the parameters $t$ and
$r$, in order to implement an arbitrary coherent operation. When the
transmissivity of BS2 is set to 1, the above operation can be reduced to
photon subtraction or photon addition.

##  III Results

We now demonstrate how the SNR of ghost imaging is enhanced with non-Gaussian
operations. In FIG. 2, we plot the SNR of the ghost imaging with non-Gaussian
state of light $(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ as a function of $r$ for $s=0.01$ and
$s=0.35$, exemplifying the effect of the coherent superposition operation on
SNR. The photon subtraction (addition) is naturally implemented with $r=0\
(r=1)$. We can see that for $s=0.35$, the maximal SNR is obtained at $r=0$
(photon subtraction) and $r=1$ (photon addition). While for $s=0.01$, the
maximal SNR is obtained at $r\simeq s=0.01$, which is much greater than that
at $r=1$ or $r=0$.

In FIG. 2, we plot the SNR of ghost imaging with $(t\hat{a}_{{\rm
s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ in the plane of $s$
and $r$, which reveals the evolution of optimal SNR in the parameter space
(red dashed line). In the large-squeezing regime ($s>0.09$), the optimal
operation is simply photon subtraction or addition. Interestingly, at $s\simeq
0.09$, there appears a critical point where the optimal $(r,t)$ bifurcate
abruptly to $r=0$ and $r=1$. Actually, in the weak squeezing regime, the
optimal SNR is obtained with the coherent superposition operation at $r\simeq
s$.

Figure 2: (a) SNR of ghost imaging as a function of the parameter $r$ of the
coherent operation with (i) $s=0.01$ for the states: $\hat{a}_{{\rm
s}}^{\dagger}|{\mathrm{TMSS}}\rangle$ (blue dashed), $(t\hat{a}_{{\rm
s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ (red solid), and
(ii) $s=0.35$ for the states: $\hat{a}_{{\rm
s}}^{\dagger}|{\mathrm{TMSS}}\rangle$ (green dashed), $(t\hat{a}_{{\rm
s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ (purple solid). (b)
SNR of ghost imaging as a function of $s$ and $r$ for the state
$(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\mathrm{TMSS}}\rangle$.
The SNR can be optimized by choosing proper $r$ for each $s$, and the optimal
$(r,s)$ parameters are shown as red dashed line. Figure 3: SNR of ghost
imaging as a function of $s$ for the states: $|{\mathrm{TMSS}}\rangle$ (yellow
dashed), $\hat{a}_{{\rm s}}|{\mathrm{TMSS}}\rangle$ (or $\hat{a}_{{\rm
s}}^{\dagger}|{\mathrm{TMSS}}\rangle$) (blue dashed), and $(t\hat{a}_{{\rm
s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ (red solid). The
parameter $r$ in the coherent operation is optimized for each $s$.

In FIG. 3, we show the optimized SNR of ghost imaging with non-Gaussian state
of light $(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ as a function of the squeezing
parameter $s$, compared with those using $\hat{a}_{{\rm
s}}|{\mathrm{TMSS}}\rangle$ (or $\hat{a}_{{\rm
s}}^{\dagger}|{\mathrm{TMSS}}\rangle$) and the Gaussian state
$|{\mathrm{TMSS}}\rangle$ as light source. There is a slight improvement of
SNR via photon subtraction (addition). However, the SNR is remarkably enhanced
by the coherent operation $t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger}$ in
the weak squeezing region, and reaches a maximal level even greater than that
of strongly squeezed state. The turning point at $s\simeq 0.091$ comes from
the abrupt change of the optimized operation from coherent operation to photon
subtraction or addition. Given the technical difficulties in preparing TMSS
with large squeezing parameter, such enhancement of SNR can be quite useful
for ghost imaging with quantum light in the weak squeezing regime. The SNR of
ghost imaging with TMSS in the experiment-friendly weak-squeezing regime can
be significantly enhanced via the coherent operation. This circumvents the
difficulty of preparing TMSS with large squeezing parameter, using the
conditional implementation of non-Gaussian operation where strong squeezing
process is not necessary either.

We provide an analytical understanding of the behavior of the SNR with non-
Gaussian light $(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger})|{\mathrm{TMSS}}\rangle$. Assuming the squeezing parameter $s\ll
1$ in the weak-squeezing regime, the TMSS can be approximated by

$\displaystyle|{\rm TMSS}\rangle=\frac{1}{\cosh
s}\sum_{k=0}^{\infty}\tanh^{k}s\ {|{k_{\rm s}k_{\rm i}}\rangle}\simeq{|{0_{\rm
s}0_{\rm i}}\rangle}+{s|{1_{\rm s}1_{\rm i}}\rangle}$ (5)

The non-Gaussian coherent operation $t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger}$ on the TMSS with $|t|^{2}+|r|^{2}=1$ then generates an output
state as

$\displaystyle(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger})|{\mathrm{TMSS}}\rangle\simeq r{|{1_{\rm s}0_{\rm
i}}\rangle}+s(t{|{0_{\rm s}1_{\rm i}}\rangle}+\sqrt{2}r{|{2_{\rm s}1_{\rm
i}}\rangle})$ (6)

According to the optimization strategy in the weak-squeezing regime, the
optimal SNR is obtained at $r\simeq s\ll 1$, with the output state given by

$\displaystyle(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ $\displaystyle\simeq s{|{1_{\rm
s}0_{\rm i}}\rangle}+s\sqrt{1-s^{2}}{|{0_{\rm s}1_{\rm
i}}\rangle}+\sqrt{2}s^{2}{|{2_{\rm s}1_{\rm i}}\rangle}$
$\displaystyle\rightarrow\frac{1}{\sqrt{2}}({|{1_{\rm s}0_{\rm
i}}\rangle}+{|{0_{\rm s}1_{\rm i}}\rangle})$ (7)

Thus we obtain a Bell state $\frac{1}{\sqrt{2}}({|{1_{\rm s}0_{\rm
i}}\rangle}+{|{0_{\rm s}1_{\rm i}}\rangle})$ as an approximation of the non-
Gaussian light in the weak squeezing regime. The calculated SNR given by this
Bell state is about 0.07, which coincides well with the optimal SNR of
$(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\rm TMSS}\rangle$ in the
weak squeezing regime. On the other hand, the output states cannot be reduced
to a Bell state by photon subtraction or addition alone, resulting in less
enhancement of SNR than by the coherent $t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger}$ operation. The optimized $(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ greatly differs from $\hat{a}_{{\rm
s}}|{\mathrm{TMSS}}\rangle$ (or $\hat{a}_{{\rm
s}}^{\dagger}|{\mathrm{TMSS}}\rangle$) in the weak-squeezing regime. As the
squeezing parameter $s$ increases, the higher order terms of $s$ in Eq.
(5),(6) can no longer be neglected, and the approximation to Bell state in Eq.
(7) becomes invalid. Due to the contribution of the component states related
to the higher order terms of $s$ in $(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm
s}}^{\dagger})|{\rm TMSS}\rangle$, such as ${|{1_{\rm s}2_{\rm i}}\rangle}$
and ${|{2_{\rm s}1_{\rm i}}\rangle}$, the SNR decreases relevantly. When
$(t\hat{a}_{{\rm s}}+r\hat{a}_{{\rm s}}^{\dagger})|{\mathrm{TMSS}}\rangle$ no
longer shows superiority to $\hat{a}_{{\rm s}}|{\mathrm{TMSS}}\rangle$ (or
$\hat{a}_{{\rm s}}^{\dagger}|{\mathrm{TMSS}}\rangle$) in ghost imaging, the
turning point appears, and the enhancement of SNR is then optimized by the
photon subtraction or addition operation.

##  IV Conclusion

In summary, we have proposed an experimental scheme to improve the SNR of
ghost imaging by non-Gaussian coherent operation, instead of increasing the
squeezing parameter of TMSS. We have achieved a remarkable increase of SNR of
ghost imaging in the weak squeezing regime, demonstrated by numerical
simulation. We have also given an analytical understanding of the enhancement
of SNR via coherent operation, which is related to the EPR correlation. The
enhancement of SNR can be regarded as a significant advance for ghost imaging
with quantum light, overcoming the limitation of SNR due to the technical
difficulty in preparing states in the strong squeezing regime. In addition, it
can be considered as a feasible application of the coherent superposition
operation of photon subtraction and addition.

##  V Acknowledgements

This work is supported by National Natural Science Foundation of China (Grant
Nos. 12174009, 11974031).

## References

  * Lugiato  _et al._ (2002) L A Lugiato, A Gatti, and E Brambilla, “Quantum imaging,” J. Optics B: Quantum and Semiclassical Optics 4, S176–S183 (2002).
  * Gatti  _et al._ (1999) A. Gatti, E. Brambilla, L. A. Lugiato, and M. I. Kolobov, “Quantum entangled images,” Phys. Rev. Lett. 83, 1763–1766 (1999).
  * Brambilla  _et al._ (2008) E. Brambilla, L. Caspani, O. Jedrkiewicz, L. A. Lugiato, and A. Gatti, “High-sensitivity imaging with multi-mode twin beams,” Phys. Rev. A 77, 053807 (2008).
  * Delaubert  _et al._ (2008) V. Delaubert, N. Treps, C. Fabre, H. A. Bachor, and P. Réfrégier, “Quantum limits in image processing,” EPL (Europhys. Lett.) 81, 44001 (2008).
  * Brida  _et al._ (2010) G. Brida, M. Genovese, and I. Ruo Berchera, “Experimental realization of sub-shot-noise quantum imaging,” Nature Photonics 4, 227–230 (2010).
  * Brida  _et al._ (2011a) Giorgio Brida, Marco Genovese, Alice Meda, and Ivano Ruo Berchera, “Experimental quantum imaging exploiting multimode spatial correlation of twin beams,” Phys. Rev. A 83, 033811 (2011a).
  * Santos  _et al._ (2008) Ivan F. Santos, J. G. Aguirre-Gómez, and S. Pádua, “Comparing quantum imaging with classical second-order incoherent imaging,” Phys. Rev. A 77, 043832 (2008).
  * D’Angelo  _et al._ (2004) Milena D’Angelo, Yoon-Ho Kim, Sergei P. Kulik, and Yanhua Shih, “Identifying entanglement using quantum ghost interference and imaging,” Phys. Rev. Lett. 92, 233601 (2004).
  * Boto  _et al._ (2000) Agedi N. Boto, Pieter Kok, Daniel S. Abrams, Samuel L. Braunstein, Colin P. Williams, and Jonathan P. Dowling, “Quantum interferometric optical lithography: Exploiting entanglement to beat the diffraction limit,” Phys. Rev. Lett. 85, 2733–2736 (2000).
  * Giovannetti  _et al._ (2009) Vittorio Giovannetti, Seth Lloyd, Lorenzo Maccone, and Jeffrey H. Shapiro, “Sub-rayleigh-diffraction-bound quantum imaging,” Phys. Rev. A 79, 013827 (2009).
  * Santos  _et al._ (2003) Ivan F. Santos, M. A. Sagioro, C. H. Monken, and S. Pádua, “Resolution and apodization in images generated by twin photons,” Phys. Rev. A 67, 033812 (2003).
  * Sciarrino  _et al._ (2008) Fabio Sciarrino, Chiara Vitelli, Francesco De Martini, Ryan Glasser, Hugo Cable, and Jonathan P. Dowling, “Experimental sub-rayleigh resolution by an unseeded high-gain optical parametric amplifier for quantum lithography,” Phys. Rev. A 77, 012324 (2008).
  * Bussgang  _et al._ (2012) Julian J. Bussgang, Julian J. Bussgang, Leon Cohen, H. Vincent Poor, and Marlan O. Scully,  _Classical, Semi-classical and Quantum Noise_ , 1st ed. (Springer Verlag, 2012).
  * Klyshko (1988a) D. N. Klyshko, “A simple method of preparing pure states of an optical field, of implementing the einstein–podolsky–rosen experiment, and of demonstrating the complementarity principle,” Soviet Physics Uspekhi 31, 74–85 (1988a).
  * Pittman  _et al._ (1995) T. B. Pittman, Y. H. Shih, D. V. Strekalov, and A. V. Sergienko, “Optical imaging by means of two-photon quantum entanglement,” Phys. Rev. A 52, R3429–R3432 (1995).
  * Brida  _et al._ (2011b) G. Brida, M. V. Chekhova, G. A. Fornaro, M. Genovese, E. D. Lopaeva, and I. Ruo Berchera, “Systematic analysis of signal-to-noise ratio in bipartite ghost imaging with classical and quantum light,” Phys. Rev. A 83, 063807 (2011b).
  * Erkmen and Shapiro (2009) Baris I. Erkmen and Jeffrey H. Shapiro, “Signal-to-noise ratio of gaussian-state ghost imaging,” Phys. Rev. A 79, 023833 (2009).
  * Valencia  _et al._ (2005) Alejandra Valencia, Giuliano Scarcelli, Milena D’Angelo, and Yanhua Shih, “Two-photon imaging with thermal light,” Phys. Rev. Lett. 94, 063601 (2005).
  * Thiel  _et al._ (2007) C. Thiel, T. Bastin, J. Martin, E. Solano, J. von Zanthier, and G. S. Agarwal, “Quantum imaging with incoherent photons,” Phys. Rev. Lett. 99, 133603 (2007).
  * Lopaeva and Chekhova (2010) E. D. Lopaeva and M. V. Chekhova, “Ghost imaging with the use of the variance of the difference photocurrent,” JETP Lett. 91, 447–451 (2010).
  * Ferri  _et al._ (2010) F. Ferri, D. Magatti, L. A. Lugiato, and A. Gatti, “Differential ghost imaging,” Phys. Rev. Lett. 104, 253603 (2010).
  * Meyers  _et al._ (2011) Ronald E. Meyers, Keith S. Deacon, and Yanhua Shih, “Turbulence-free ghost imaging,” Appl. Phys. Lett. 98, 111115 (2011).
  * Klyshko (1988b) D.N. Klyshko,  _Photons and nonlinear optics_ , rev. and enl. ed ed. (Gordon and Breach, 1988).
  * Shih (2003) Y. Shih, “Entangled photons,” IEEE Journal of Selected Topics in Quantum Electronics 9, 1455–1467 (2003).
  * Smithers and Lu (1974) Martin E. Smithers and Eugene Y. C. Lu, “Quantum theory of coupled parametric down-conversion and up-conversion with simultaneous phase matching,” Phys. Rev. A 10, 1874–1880 (1974).
  * Schumaker and Caves (1985) Bonny L. Schumaker and Carlton M. Caves, “New formalism for two-photon quantum optics. ii. mathematical foundation and compact notation,” Phys. Rev. A 31, 3093–3111 (1985).
  * Barnett and Knight (1987) S.M. Barnett and P.L. Knight, “Squeezing in correlated quantum systems,” J. Mod. Optics 34, 841–853 (1987).
  * Lee  _et al._ (2011) Su-Yong Lee, Se-Wan Ji, Ho-Joon Kim, and Hyunchul Nha, “Enhancing quantum entanglement for continuous variables by a coherent superposition of photon subtraction and addition,” Phys. Rev. A 84, 012302 (2011).
  * Zavatta  _et al._ (2004) Alessandro Zavatta, Silvia Viciani, and Marco Bellini, “Quantum-to-classical transition with single-photon-added coherent states of light,” Science 306, 660–662 (2004).
  * Parigi  _et al._ (2007) Valentina Parigi, Alessandro Zavatta, Myungshik Kim, and Marco Bellini, “Probing quantum commutation rules by addition and subtraction of single photons to/from a light field,” Science 317, 1890–1893 (2007).
  * Lee  _et al._ (2019) Jaehak Lee, Jiyong Park, and Hyunchul Nha, “Quantum non-gaussianity and secure quantum communication,” npj Quantum Information 5, 49 (2019).
  * Liao  _et al._ (2018) Qin Liao, Ying Guo, Duan Huang, Peng Huang, and Guihua Zeng, “Long-distance continuous-variable quantum key distribution using non-gaussian state-discrimination detection,” New Journal of Physics 20, 023015 (2018).
  * Ra  _et al._ (2020) Young-Sik Ra, Adrien Dufour, Mattia Walschaers, Clément Jacquard, Thibault Michel, Claude Fabre, and Nicolas Treps, “Non-gaussian quantum states of a multimode light field,” Nature Phys. 16, 144–147 (2020).
  * Serafini (2017) Alessio Serafini,  _Quantum Continuous Variables: A Primer of Theoretical Methods_ , 1st ed. (CRC Press, Taylor and Francis Group, 2017).
  * (35) See the Supplementary Material at URL for details of calculation of SNR and non-Gaussian coherent operation.
  * Lee and Nha (2010) Su-Yong Lee and Hyunchul Nha, “Quantum state engineering by a coherent superposition of photon subtraction and addition,” Phys. Rev. A 82, 053812 (2010).
  * Wenger  _et al._ (2004) Jérôme Wenger, Rosa Tualle-Brouri, and Philippe Grangier, “Non-gaussian statistics from individual pulses of squeezed light,” Phys. Rev. Lett. 92, 153601 (2004).


# A Bayesian Approach to Strong Lens Finding in the Era of Wide-area Surveys

Philip Holloway1, Philip J. Marshall2,3, Aprajita Verma1, Anupreeta More4,5,
Raoul Cañameras6,7, Anton T. Jaelani8,9, Yuichiro Ishida10, Kenneth C.
Wong11,12
1 Sub-department of Astrophysics, University of Oxford, Denys Wilkinson
Building, Keble Road, Oxford OX1 3RH, UK
2 Kavli Institute for Particle Astrophysics and Cosmology, Department of
Physics, Stanford University, Stanford, CA 94305, USA
3 SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
4 The Inter-University Centre for Astronomy and Astrophysics (IUCAA), Post Bag
4, Ganeshkhind, Pune 411007, India
5 Kavli Institute for the Physics and Mathematics of the Universe (IPMU),
5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8583, Japan
6 Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748
Garching, Germany
7 Technical University of Munich, TUM School of Natural Sciences, Department
of Physics, James-Franck-Str. 1, D-85748 Garching, Germany
8 Astronomy Research Group and Bosscha Observatory, FMIPA, Institut Teknologi
Bandung, Jl. Ganesha 10, Bandung 40132, Indonesia
9 U-CoE AI-VLB, Institut Teknologi Bandung, Jl. Ganesha 10, Bandung 40132,
Indonesia
10 Department of Earth and Planetary Sciences, School of Science, Kyushu
University, 744 Motooka, Nishi-ku, Fukuoka 819-0395, Japan
11 Research Center for the Early Universe, Graduate School of Science, The
University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
12 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo
181-8588, Japan
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The arrival of the Vera C. Rubin Observatory’s Legacy Survey of Space and Time
(LSST), Euclid-Wide and Roman wide area sensitive surveys will herald a new
era in strong lens science in which the number of strong lenses known is
expected to rise from $\mathcal{O}(10^{3})$ to $\mathcal{O}(10^{5})$. However,
current lens-finding methods still require time-consuming follow-up visual
inspection by strong-lens experts to remove false positives which is only set
to increase with these surveys. In this work we demonstrate a range of methods
to produce calibrated probabilities to help determine the veracity of any
given lens candidate. To do this we use the classifications from citizen
science and multiple neural networks for galaxies selected from the Hyper
Suprime-Cam (HSC) survey. Our methodology is not restricted to particular
classifier types and could be applied to any strong lens classifier which
produces quantitative scores. Using these calibrated probabilities, we
generate an ensemble classifier, combining citizen science and neural network
lens finders. We find such an ensemble can provide improved classification
over the individual classifiers. We find a false positive rate of $10^{-3}$
can be achieved with a completeness of $46\%$, compared to $34\%$ for the best
individual classifier. Given the large number of galaxy-galaxy strong lenses
anticipated in LSST, such improvement would still produce significant numbers
of false positives, in which case using calibrated probabilities will be
essential for population analysis of large populations of lenses.

###### keywords:

gravitational lensing: strong – methods: statistical – methods: data analysis
– surveys

††pubyear: 2023††pagerange: A Bayesian Approach to Strong Lens Finding in the
Era of Wide-area Surveys–C

##  1 Introduction

Gravitational lensing is the deflection of light due to gravity as it passes a
massive object. This can occur at a wide range of mass scales, ranging from
microlensing of single stars by a compact source (e.g. a black hole, Sahu et
al. (2022)), to lensing of galaxies by galaxy clusters (e.g. SMACS 0723,
Pascale et al. (2022)). Strong lensing occurs when, from the perspective of
the observer, the source appears as multiple images surrounding the lens.
Strong lensing is a rare occurrence: only the most massive galaxies can
deflect the source light sufficiently to be resolvable from the lensing
galaxy. The number of galaxy-scale strong lenses identified to date is of
order 1000 (hundreds of which have been confirmed, e.g. Bolton et al. (2008);
Brownstein et al. (2012)). These have been found across a range of surveys and
telescopes, such as DES (e.g. Jacobs et al. (2019a, b)), CFHTLS (Marshall et
al., 2016; More et al., 2016; Jacobs et al., 2017), CFIS (Savary et al.,
2022), SDSS (e.g. Bolton et al. (2006, 2008)), HSC (Sonnenfeld et al., 2020;
Shu et al., 2022; Cañameras et al., 2021), HST (Faure et al., 2008; Jackson,
2008; Pawase et al., 2014), KiDS (Li et al., 2021; He et al., 2020; Petrillo
et al., 2017, 2019) including in the radio (e.g. the VLA, Myers et al. (2003);
Browne et al. (2003)). The known lenses therefore span a wide resolution range
and all come with different selection functions. The number of lenses useful
for a particular science case is typically much smaller than the number of
known lenses, particularly if that science case requires spectroscopic
redshifts or data at particular wavelengths. Such cases include constraining
the stellar IMF over a lens population (Sonnenfeld et al., 2019), the
redshift-evolution of the velocity dispersion function of ETG’s (Geng et al.,
2021) (both of which required spectroscopic redshifts), or the mass of the
fuzzy-dark matter particle through high-resolution VLBI observations (Powell
et al., 2022, 2023). The $\sim 100$ SLACS lenses (Bolton et al., 2006, 2008)
have been a go-to lens sample owing to the availability of redshift and high-
resolution data. With the arrival of wide-area surveys such as LSST (Ivezić et
al., 2019), the Euclid Wide Survey (Euclid Collaboration et al., 2022) and
Roman High Latitude Wide Area survey (Akeson et al., 2019; Spergel et al.,
2015), strong lensing science will move from small samples of strong lenses to
_population-level_ analysis of the hundreds of thousands of lenses (Collett,
2015; Holloway et al., 2023) identified by these surveys. Although the 4MOST
Strong Lensing Spectroscopic Legacy Survey (4SLSLS) (Collett et al., 2023)
will provide spectroscopic redshifts for $\sim 10,000$ lens-source pairs, this
will leave the majority without spectroscopic confirmation and follow-up
visual inspection of the >>100,000 strong lens candidates put forward by lens
finding methods will be an unenviable task.
Strong lenses have been found using a wide variety of methods. These include
machine learning (typically CNN’s, e.g. Lanusse et al. (2018); Schaefer et al.
(2018), though other methods including Support Vector Machines (e.g. Hartley
et al. (2017)) and self-attention-based encoding models (Thuruthipilly et al.,
2022) have also been used), citizen science (Marshall et al., 2016; More et
al., 2016; Sonnenfeld et al., 2020; Garvin et al., 2022; Geach et al., 2015),
arc-finding algorithms (e.g YATTALENS (Sonnenfeld et al., 2018) & Arcfinder
(Seidel & Bartelmann, 2007)) and expert visual inspection (Faure et al., 2008;
Jackson, 2008; Pawase et al., 2014). Although the automated methods e.g. CNN’s
are continually improving over time, they still require visual follow-up
inspection of high-scoring candidates by strong lens experts to remove false
positives. Even with small-medium sized surveys, this is time-consuming and
difficult: strong lens identification is not clear-cut, particularly with
ground-based imaging or when using single-band imaging. From recent lens
searches (Sonnenfeld et al., 2020; Cañameras et al., 2021; Rojas et al., 2022)
the number of high-scoring candidates which underwent visual inspection
exceeded the resulting number of A-B grade lenses by factors $\geq 10\times$.
Furthermore, Rojas et al. (2023) found that $\sim$ 6 expert graders are
required to reduce the classification error below 0.1 (on a 0-1 grading
system). With wide-area surveys, lens finding methods which scale-up easily
will be key to identifying most lenses however these will still leave a large
number of high-scoring false positives to remove. This paper aims to address
this problem. The goals of this paper are two-fold:

  * •

To produce calibrated probabilities that a given galaxy system is a lens. We
define a calibrated score X as one for which systems with this score are
indeed lenses X% of the time, as measured on a distinct test set. For wide
area surveys, follow-up spectroscopy or visual inspection of all lens
candidates will not be possible. It will therefore be important to have
accurate probabilities for lens candidates in order to perform unbiased
population-level analysis, the alternative being significantly reducing the
sample-size to those which have been spectroscopically or visually confirmed.
Furthermore, having accurate probabilities allows direct comparison of lens
candidates across different lens finders.

  * •

To create an ensemble strong lens classifier. A given survey may be targeted
by multiple lens finding methods, each with their own strengths and
weaknesses. Neural network ensembles have been investigated previously (e.g.
Canameras et al. (2023); Schaefer et al. (2018) for galaxy-galaxy lensing and
Andika et al. (2023) for lensed quasars), simply averaging over the individual
network scores. We wish to take this method further, by incorporating the
calibrated probabilities generated above in a Bayesian framework, as well as
including a citizen science classifier to diversify our ensemble beyond neural
networks.

In the process of reaching these goals, we aim to answer the following
questions:

  1. 1.

Is there a significant improvement in performance using an ensemble method
compared to a single classifier?

  2. 2.

Are there regions of parameter space which suit certain methods best?

  3. 3.

How does the degree of improvement change when combining only neural network
classifiers, compared to combining neural networks with a citizen science
classifier?

This paper is structured as follows. We detail the data used in this work in
Section 2. In Sections 3.1 and 3.2 we summarise and apply different classifier
calibration methods. In Section 3.3 we then combine these calibrated
probabilities into a single ensemble classifier. Our results are presented in
Section 4 and discussed further in Section 5, including implications for
forthcoming lens searches in LSST and other wide-area surveys in Section 5.4.
We conclude in Section 6.

##  2 Data

To develop our ensemble classifier, we used the classification outputs from
six strong lens finders applied to Hyper-Suprime Cam Subaru Strategic Program
(HSC SSP) data (Aihara et al., 2022). This, being a ground based survey over
hundreds of square degrees, was a useful proxy for LSST-like data. The wide-
area survey is being conducted in _grizy_ bands, covering $\sim 1200\deg^{2}$.
The HSC S17A and PDR2 data releases (Aihara et al., 2018, 2019) covered
$1026\deg^{2}$ (225) and $1114\deg^{2}$ (305) respectively in at least one
band (all bands) to a $5\sigma$ depth of $i\sim 26.2$. Only the $gri$ bands
were used in this work. The strong lens classifiers used are described below:

  * •

Citizen Science (Sonnenfeld et al., 2020): A citizen science search by the
Space Warps project, which used the HSC S17A data release (Aihara et al.,
2018, 2019) and provided citizen classifications for $\sim 300,000$ objects.
The galaxies selected for this search had photometric redshifts $0.2<z<1.2$
and inferred stellar masses of log $M_{*}$ > 11.2.

  * •

Neural Network 1 (Cañameras et al., 2021): A ResNet classifier, which used HSC
PDR2 data (Aihara et al., 2019) over a much larger sample of objects (network
scores and positions were available for $5.4\times 10^{7}$ objects). The
objects chosen for this search were required to have i-band Kron radius $\geq
0.8\arcsec$.

  * •

Neural Network 2: (Shu et al., 2022): A ResNet classifier from Lanusse et al.
(2018), using HSC PDR2 data and in particular targeting high redshift lens
systems. The primary selections for this search were colour cuts:
$0.6<g-r<3.0$ and $2.0<g-i<5.0$ to identify red, high-redshift galaxies.

  * •

Neural Network 3: The second classifier presented by Shu et al. (2022),
differing from the first in having a higher fraction of z>0.6 lenses in the
training set, and image pre-processing.

  * •

Neural Network 4 (Ishida et al., in prep.): A convolutional neural network
applied to HSC PDR2 data over the galaxy-sample with classifier outputs
available for the above 4 classifiers.

  * •

Neural Network 5 (Jaelani et al., in prep.): A convolutional neural network
applied to HSC PDR2 and in particular targeting smaller Einstein radii > 0.5”.
The galaxies were originally selected based on pre-selection criteria, such as
stellar mass > $5\times 10^{10}M_{\odot}$, specific star formation rate <
$10^{-10}$, and redshift range $0.2<z<1.2$. However, in this work, the CNN was
applied to the sample of galaxies with classifier outputs available from all
of the first 4 classifiers listed.

The objects in each catalogue were cross-matched, with a maximum separation of
1″, to produce a sample of 126,312 galaxies with both citizen science and
neural network outputs. A subset (109,128 objects), with $\geq 5$ citizen
classifications, was used in the following analysis. In order to calibrate the
output scores of the classifiers, we required a ‘ground-truth’, i.e. a list of
classified objects of known lenses/non-lenses. We collated such a list using
expert grades from four sources: Sonnenfeld et al. (2020), Cañameras et al.
(2021), Shu et al. (2022) and the SuGOHI lens database111http://www-
utap.phys.s.u-tokyo.ac.jp/~oguri/sugohi/. We then assigned each cross-matched
object with an overall grade according to the following order: HOLISMOKES VI,
SuGOHI VI, SuGOHI database and HOLISMOKES VIII. There were 34 objects whose
grades differed by $\geq 2$ which were inspected by PH and re-assigned an
appropriate grade (typically taking the average of the relevant grades except
in clear-cut cases). In total, 3,744 cross-matched objects had assigned
grades, of which 189 were grade A-B and were treated as true lenses.

##  3 Method

Our first aim is to produce calibrated scores for each classifier
individually. The graded objects were split into three sets: a training
(‘calibration’) set, a validation set and a test set, in ratios 2:1:1. For
this work, we considered A+B grade lenses as true lenses, and the remainder
(graded or otherwise) as non-lenses. Excluding the ungraded objects would lead
to unrealistic values for the receiver operating characteristic (ROC) curve
and purity-completeness summary statistics and it is likely the vast majority
of true lenses were identified by at least one of our sources of expert
grades.

###  3.1 Summary of Calibration Methods

Figure 1: Calibration curves applied to a toy model for the range of
calibration methods considered. The perfect calibration curve is shown on the
far left. Platt scaling is the worst performing as the class distribution does
not follow a sigmoid shape and thus is not considered further in this work. A
subset of the bins used to generate the calibration curve in (d) are shown;
the bin widths become smaller for larger classifier outputs as the fraction of
class-1 objects increases.

Here we discuss the various common methods available to map a classifier score
to a probability. The calibration methods considered are described below and a
demonstration of these methods applied to toy data is shown in Figure 1.

  * •

Isotonic regression (Zadrozny & Elkan, 2002) is a non-parametric, discontinous
fit of a monotonically increasing curve to the classifier output distribution.

  * •

Platt scaling (Platt, 2000) is a parametric calibration method of the form:

$P(C=1|x)=\frac{1}{1+\exp(Ax+B)}$ (1)

where $A$ and $B$ are fitted to the data. Although commonly used, this
requires the calibration curve to closely match a sigmoid function, which _a
priori_ may not be true.

  * •

Kullback Leibler importance estimation procedure (KLIEP, Sugiyama et al.
(2008)): This minimises the KL divergence of the ratio $f(S|C_{1})/f(S)$ with
probability density functions $f(S|C_{1})$ and $f(S)$ of a classifier score
$S$ (given the subject is in class 1 in the first case). From Bayes’ Theorem,
weighting this by $P(C_{1})$ gives the probability a subject is in class 1
given a classifier score. KLIEP uses a Gaussian mixture model, with Gaussians
of fixed width placed at the positions of the lenses in classifier-output
space. Their respective weights are tuned to minimise the KL divergence.

  * •

Variable Bin Fitting: This was designed to account for the small number of
lenses in our sample, compared to the much larger number of non-lenses. This
generates overlapping bins containing equal numbers of lenses (we chose 5),
rather than bins of equal width. The calibration curve is then formed from the
fraction of lenses in each bin compared to the total bin size. This provides
higher resolution calibration in regions where there are lots of lenses (i.e.
for high classifier scores), while averaging out the calibration where lenses
are sparse (i.e. for low classifier scores).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Calibration mapping for a range of calibration methods for each of
the classifiers used in this work. For clarity, only the top $20\%$ of
classifier scores are shown in each panel. The red and blue histograms refer
to the fraction of lenses (blue) and non-lenses (red) in each bin. The
mappings should be expected to roughly follow these histograms but not always,
for example at the highest classifier scores, where the fraction of true
lenses increases rapidly within a single bin. The bars are shaded by the
number of objects in each bin. The dashed line is the $y=x$ line, which the
calibration curves would follow if the outputs of the classifiers were already
calibrated.

###  3.2 Application of Calibration Methods

We applied a subset of calibration methods to our data. We did not apply Platt
Scaling as we did not expect our calibration curves to match a sigmoid
function. A key feature of the data was the significant class imbalance: only
$\sim 0.2\%$ of all objects used in this work were lenses, and these are more
concentrated at high classifier scores. We used this to our advantage however,
to achieve higher-resolution calibration for high scoring subjects.
For each calibration method, we applied the calibration algorithm to the rank
values of the classifier outputs, rather than the outputs themselves. We then
interpolated back from rank to classifier score to produce the calibration
map. We wished our calibration method to be independent of the type of
classifier used and the score distribution of its output. The distribution of
outputs can change significantly between different classifiers acting on the
same objects, however an excess of high (raw) scoring objects from a
particular classifier should not be treated as all of these objects having
high calibrated probabilities. Using the ranks removed this effect and simply
assumes that a higher classifier score (qualitatively) implies higher
confidence of a lens.
The KLIEP and variable binning methods have hyperparameters we fixed prior to
calibration. We used a gaussian kernel of width $\sigma=40$; a balance between
overfitting to the training data (occurring from a smaller kernel) and
preventing over-smoothing at high scores (from a larger kernel). The GMM model
had the same number of kernels as lenses. As shown by the calibration curve,
the fraction of objects which are lenses increases significantly for high
scores: using the rank values allowed finer tuning to these regions while not
overfitting at lower classifier scores. Figure 2 shows the calibration curves
for each of the methods applied. In all cases the curves steepen for high
scores: choice of score threshold will have a significant impact on purity in
this region. The calibration curves do not follow the y=x line (dotted),
indicating that (as expected) the original classifier outputs were not already
calibrated. The mappings from the variable bin and isotonic regression method
are relatively similar to each other, across all classifiers; the KLIEP method
differs more significantly as the PDF only becomes non-zero close to the
positions of the lenses in the training set. This effect would be reduced by
increasing the kernel size, but would lead to underfitting at the highest
classifier scores. Since we were primarily focussed on high-probability
candidates, we prioritised reducing this underfitting.
We validated these calibration curves against a separate test set, shown in
Figure 3. These show that spanning a wide range on log-p, the calibration
mapping can produce accurate probabilities. The isotonic calibration produced
the smallest variation upon bootstrapping, so is used for the subsequent
analysis.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Validation of the calibration curves, applied to a separate test set
of graded images. The y-axis is determined using the variable-binning method,
and the error-bars are determined from bootstrapping. Each method shows good
calibration (perfect calibration would be along the $y$ = $x$ line), but
occasionally err on the side of underconfident, particularly for the highest
scores.

###  3.3 Summary of Combination Methods

Given now-calibrated classifier outputs, we considered methods to combine them
to produce a single score. The combined score will not necessarily itself be
calibrated, thus further calibration stage using the methods above may be
required. The aim of this classifier combination was to maximise the purity of
the resultant sample. We first considered a simple generalised mean of the
form:

$P_{combo}=\left(\frac{1}{N}\sum_{i}^{N}{p_{i}^{\alpha}}\right)^{\frac{1}{\alpha}}$
(2)

This takes on a variety of useful functions for different values of $\alpha$,
in particular the arithmetic mean, harmonic mean, minimum and maximum of
$\\{P_{i}\\}$ for $\alpha=1,-1,-\infty$ and $+\infty$ respectively.

####  3.3.1 Dependent Bayesian Classifier Combination

Figure 4: Distributions of the difference in ranking of the classifier scores
for the cross-matched objects. The ranks were first normalised to take values
in the range [0,1]. The black and blue histograms show the overall and lens
distributions respectively. The respective curves show the projections of the
best-fit multivariate gaussian to the data. The dashed curve shows a triangle
distribution (which hasn’t been fitted to the data) which results from taking
the difference between two uniformly distributed variables.

Although we have generated calibrated probabilities for each classifier, the
dependence of each classifier on another is not quantified. For example if two
identical classifiers both gave probabilities of 0.9 for the same object, our
combined posterior should be 0.9 (the second, identical, classifier adds no
new information), but if those classifiers were independent, we would expect
the posterior would be $>0.9$. We here outline a Bayesian approach to model
this dependence as follows. We wish to find $P(L|\\{R_{i}\\})$, where $R_{i}$
are the set of score rankings for a given object from each classifier. We
denote the corresponding set of calibrated classifier scores as $\\{C_{i}\\}$.
From Bayes’ Theorem:

$P(L|\\{R_{i}\\})=\frac{f(\\{R_{1},...R_{N}\\}|L,R_{0})P(L|R_{0})}{f(\\{R_{1},...R_{N}\\}|R_{0})}$
(3)

where $f$ denotes a probability density function, and P denotes the
probability of a discrete random variable. $P(L|R_{0})$ is known, given a 1-1
relation between rank and calibrated score (i.e. the calibration mapping is
strictly monotonic) already as it is given by the calibrated output:
$P(L|R_{0})=C_{0}$. We re-write Equation 3, in terms of the difference between
the classifier outputs, with respect to classifier $C_{0}$:
$\\{\Delta_{i}\\}=\\{R_{i}\\}-R_{0}$:

$P(L|\\{R_{i}\\})=\frac{f(\\{\Delta_{1},...\Delta_{N}\\}|L,R_{0})}{f(\\{\Delta_{1},...\Delta_{N}|R_{0})\\}}\cdot
C_{0}$ (4)

We now model both numerator and denominator PDFs as gaussian distributions,
constant with respect to ${R_{i}}$. We use output ranking, rather than the
classifier output itself as this provides a better fit to this distribution
but this does not change the result. We now have:

$P(L|\\{R_{i}\\})=\frac{n(\\{\Delta_{1},...\Delta_{N}\\}|\bm{\mu_{lens}},\bm{\Sigma_{lens}})}{n(\\{\Delta_{1},...\Delta_{N}\\}|\bm{\mu_{full}},\bm{\Sigma_{full}})}\cdot
C_{0}$ (5)

where $n$ denotes the Gaussian function. For 6 classifiers, the
multidimensional gaussians in Eqn 5 are 5 dimensional. These multivariate
normal distributions are independent to permutation of classifiers (i.e. which
classifier is chosen to be $C_{0}$ does not change the best fitting gaussian
distribution). However, this would change the value of $C_{0}$ (assuming all
classifiers aren’t in exact agreement). With no reason, _a priori_ , to favour
one classifier over another, it is reasonable to average over the classifiers.
Equation 5 becomes:

$P(L|\\{R_{i}\\})=\left(\frac{n(\\{\Delta_{1},\Delta_{2},\Delta_{3},\Delta_{4},\Delta_{5}\\}|\bm{\mu_{lens}},\bm{\Sigma_{lens}})}{n(\\{\Delta_{1},\Delta_{2},\Delta_{3},\Delta_{4},\Delta_{5}\\}|\bm{\mu_{full}},\bm{\Sigma_{full}})}\right)\cdot{\langle
C_{i}\rangle}$ (6)

Note, this stems from our (strong) assumption that the ratio of PDFs in
Equation 4 can be modelled as a ratio of two multivariate gaussians; more
complex functions would not suffer this problem. More flexible models (2 and
3-component mixture models of multivariate gaussian distributions) were also
trialled, but didn’t provide significant improvements in the resultant
combined calibrated score. The cases where the bracketed term in Eqn. 6
exceeds 1 refer to when there is greater agreement between the individual
classifiers than would be expected from the overall distribution (i.e. when
the blue curve exceeds the black in Figure 4) so the subject should therefore
belong to the lens class.

####  3.3.2 Independent Bayesian Classifier Combination

We also considered a case where the results from each classifier were entirely
independent of each other. From Bayes’ Theorem, for a single classifer:

$P(L|C_{1})=\frac{f(C_{1}|L)\cdot P(L)}{f(C_{1})}=\frac{f(C_{1}|L)\cdot
P_{0}}{f(C_{1}|L)\cdot P_{0}+f(C_{1}|\bar{L})\cdot(1-P_{0})}$ (7)

For a set of calibrated probabilities, $\\{C_{i}\\}:$

$P(L|\\{C_{i}\\})=\frac{P_{0}\cdot\prod_{i}f(C_{i}|L)}{P_{0}\cdot\prod_{i}f(C_{i}|L)+(1-P_{0})\cdot\prod_{i}f(C_{i}|\bar{L})}$
(8)

where $\bar{L}$ denotes ‘non-lens’. Since $\\{C_{i}\\}$ are in fact calibrated
probabilities, we know:

$C_{i}=\frac{N_{L}\cdot f(C_{i}|L)}{(N_{L}+N_{NL})\cdot f(C_{i})}$ (9)

where $N_{L}$ and $N_{NL}$ refer to the number of lenses and non-lenses in the
(training) sample. We can therefore simplify Eqn 8:

$P(L|\\{C_{i}\\})=\frac{P_{0}\cdot\prod_{i}\frac{C_{i}}{N_{L}}}{P_{0}\cdot\prod_{i}\frac{C_{i}}{N_{L}}+(1-P_{0})\cdot\prod_{i}\frac{(1-C_{i})}{N_{NL}}}$
(10)

For an accurate prior: $P_{0}$ = $N_{L}/(N_{L}+N_{NL})$:

$P(L|\\{C_{i}\\})=\frac{N_{L}^{1-N_{cl}}\cdot\prod_{i}C_{i}}{\left(N_{L}^{1-N_{cl}}\cdot\prod_{i}C_{i}\right)+\left(N_{NL}^{1-N_{cl}}\cdot\prod_{i}(1-C_{i})\right)}$
(11)

where $N_{cl}$ denotes the number of classifiers in the ensemble.

##  4 Results

###  4.1 Testing the Bayesian Combination Approaches

We generated a toy model to test the differences between the dependent and
independent combination methods described above. We trialled 3 different sets
of 6 classifiers, with varying levels of inter-dependency; we found that while
both combination methods showed an improvement in the ROC curve compared to
the individual classifiers, which of the dependent and independent methods
showed the greatest improvement varied depending on the level of dependency
set by the toy model. For sets of real classifiers, we would therefore
recommend applying both to verify which performs best.

###  4.2 Applying the Bayesian Combination methods

We investigated the dependence of the six classifiers used in this analysis.
Figure 5 shows the distribution of classifier rank (normalised to 1) for all
cross-matched subjects. The greatest correlation is seen between the two
HOLISMOKES VIII neural network classifiers; this is perhaps expected as these
networks have the same architecture and similar non-lens training data. The
distributions against the citizen science classifier are near-uniform - this
suggests that the networks and citizen classifiers find different objects
easier/more difficult to classify (otherwise, the same objects would receive
similar rankings from each). Figure 4 shows the binned distributions of these
data, along with a (single) multivariate gaussian fit. There is greater
agreement between classifiers when presented with a true lens, than there is
with non-lenses. It is this property which is used by the dependent Bayesian
classifier combination method described in Section 3.3.1.

Figure 5: Distribution of the ranking of the classifier scores for the cross-
matched objects. The ranks have been normalised to take values in the range
[0,1]. The points have been shaded according to their expert grade (Red:
C-Grade, Orange: B-Grade, Green: A-Grade, Black: Ungraded).

Upon inspecting the ensemble probabilities for a small number of
spectroscopically confirmed lenses, we observed that in some cases, while the
citizen science classifier would correctly identify the lens, the scores from
the networks would not be sufficiently high to map to high probabilities. The
networks could then effectively ‘outvote’ the citizen science classifier in
the independent Bayesian ensemble. We therefore trialled generating a network-
only ensemble, recalibrating this, then further combining this ensemble with
the citizen science classifier. We show the ROC curve for this method as
‘Network Ensemble + CS’ in Figure 6 which shows the ROC curves for the
individual classifiers along with the best performing ensemble methods. A
wider comparison of the ensemble methods described in Section 3.3 is given in
Appendix A. When considering the A-B lenses (Figure 6(a)), the ensemble
methods show improved classification over their individual constituent
classifiers. We find $46\%$ completeness can be achieved with a false positive
rate of $10^{-3}$; by comparison, the best individual classifier achieved
$34\%$ completeness on the same dataset. How this relates to the sample purity
is discussed in Section 5.4.

(a) Including A-B grade lenses as ‘true lenses’, without excluding cluster-
scale candidates.

(b) Including A-C grade lenses as ‘true lenses’ and excluding cluster-scale
candidates.

Figure 6: Receiver operating characteristic (ROC) curve for the individual
lens classifiers (dashed) and combined methods, applied to a separate test set
to that which the calibration methods were tuned. Ungraded subjects were
treated as non-lenses. The dashed guidelines show the $50\%$ completeness and
$0.1\%$ false positive rate. Note the y-axes are linear, while the x-axes are
logarithmic.

##  5 Discussion

###  5.1 Comparison with previous work

The overlap in lenses found by machine learning, citizen science and
spectroscopy was investigated by Knabel et al. (2020). They found very little
overlap between the 3 methods: out of 107 lenses identified, only two were
identified by more than one method (ML + CS). They attributed this two
differences in the parent sample (e.g. different redshift cuts) and particular
behaviours of each method (e.g. ML typically finding lenses similar to its
training set). Our results show citizen science and machine learning can be in
much greater agreement. There are however two significant differences between
our methods: firstly, Knabel et al. (2020) employed GalaxyZoo (Lintott et al.
(2008), Marshall et al. (2016), Holwerda et al. (2019), Kelvin et al. in prep)
as their citizen science classifier, which uses a question tree to identify
the overall galaxy morphology, including the presence of lensing whereas
Sonnenfeld et al. (2020) only looked for strong lenses. Furthermore, in this
work we only compare cross-matched objects which both techniques classified as
opposed to objects simply in the same field, removing the effect of differing
selection functions for each sample. The difference in our results highlights
the importance of object selection when conducting lens searches; too narrow a
selection could significantly reduce the number of lenses identified.

###  5.2 Comparison of Citizen Science vs a Network Ensemble

Figure 7: Example cutouts of high and low scoring objects from the citizen
science search and the ensemble of 5 neural network classifiers. Objects with
expert grades $1.5<G<2$ are outlined in orange, and those with $2\leq G$ are
outlined in green. The outer and inner dashed lines show the respective cutout
sizes shown to the networks in Cañameras et al. (2021) and Shu et al. (2022)
respectively, while the whole cutout was shown to the citizen scientists.

We investigated the qualitative properties of the galaxy systems identified
and rejected by the citizen science classifier compared to those receiving
high/low scores from an ensemble of 5 neural networks. Figure 7 shows a
selection of such cutouts. Firstly, the vast majority of those which received
very high (low) probabilities from both the ensemble and citizen science were
correctly identified as (non) lenses. Those which received high citizen scores
but low ensemble scores contained a mix of true lenses and interlopers
(according to the expert grades used). There is no clear trend in these
objects, but the presence of many bright bulges suggest that the network
classifiers may have learnt to reject these; it should be noted that the
networks did not have access to lens-subtracted images (unlike the citizen
scientists) which would make some images, for example those with bright bulges
but small Einstein radii, harder for the networks to classify. Furthermore, a
small number of candidates were groups or had lensing features outside of the
cutouts provided to the network, so it is unsurprising the networks rejected
these - we tested the effect of this below (Section 5.3). The top-left panel
of Figure 7 shows a selection of objects which received high ensemble
probabilities but low citizen science probabilities. These contained a number
of face-on spiral galaxies where some networks may have been misidentified the
spiral arms as lensed arcs.

###  5.3 Effect of ground-truth selection on classifier performance

As demonstrated in Section 5.2, some of the objects which the network ensemble
assigned a low probability to were galaxy clusters, which would not have been
included in their training sets. These would, however, have been identifiable
in a citizen science search. We identified galaxy clusters which had
previously been assigned as ‘true lenses’ as follows. We cross-matched the A-C
grade lenses in the masterlens database222https://test.masterlens.org/
(Moustakas, 2012) with our object sample and retrieved those with a ‘CLUST-
GAL’ flag from the database. Since not all the objects in our object sample
were also in the database, we conducted a further visual inspection of the A-C
graded candidates in our sample, flagging the cluster-scale lenses. Any
objects identified by either method which had received a grade greater than
the relevant cutoff ($\geq$ B grade, or $\geq$ C grade, as specified) were
removed from the sample. These accounted for $\sim 25\%$ of the lens systems
in our sample. We show the corresponding ROC curves for the individual
classifiers in this work with/without the clusters removed in Appendix B. We
find there is a small narrowing in the performance difference between the
citizen science and neural network classifiers when the cluster-scale lenses
are excluded, however; a combination of the classification method and use of
lens subtraction means the citizen science classifier still outperforms the
networks for our object sample.
Since removing cluster-scale lenses reduced our ground-truth sample size, we
also investigated the inclusion of C-grade lenses as ‘true lenses’ in our
ensemble. While a sizeable fraction may not be lenses in reality, this mimics
the effect of increasing the sample size of lenses for the calibration, as may
be available for a wider survey. We found this improved the calibration for
our classifiers; validation plots for this set-up are shown in Appendix C. In
turn, we found this improved the performance of the classifier ensemble
compared to the individual classifiers. Figure 6(b) shows the ROC curves using
A-C grade lenses as true lenses and with cluster-scale candidates removed. We
found the ensemble of networks only (no CS) provided substantial improvement
above the best individual network and performs better (in AUROC, FPR at $50\%$
completeness and TPR at FPR$=10^{-3}$) than the individual citizen science
classifier. When citizen science was included, the ensemble further improved,
increasing the completeness from $27\%$ (Network-only) to $40\%$ at
FPR=$10^{-3}$.

###  5.4 Expectations and Implications for LSST

In Figure 8, we compare the number of true and false positives expected for an
LSST-like survey for an ensemble including and excluding the citizen science
classifier. In these plots, we have scaled the total number of strong lenses
to $10^{5}$ and used the true and false positive rate functions stemming from
the ensemble and individual classifiers applied to the HSC data in this work.
We find that, when including citizen science in our sample, a $40\%$ complete
sample can be achieved with $49\%$ purity, compared to $32\%$ purity for the
best individual classifier. However, for higher completeness (towards the left
of each plot), the sample would remain overwhelmingly false positives. We will
further discuss the implications and mitigation for population-level analysis
in a future work (Holloway et al. 2024 in prep). The difference is more
substantial when the citizen science classifier is excluded, nearly doubling
the purity to $28\%$ for a $40\%$ complete sample.
For such large samples, expert grading of all but the highest ranked
candidates becomes intractable. It thus becomes even more important that the
outputs of lens-finders are calibrated to allow statistical analysis of large
samples of strong lenses including a known proportion of false positives.
Figure 9 shows the effective sample size of both the individual calibrated
classifiers and the ensemble. It demonstrates the ensemble can retrieve a
larger effective sample of lenses from the data, as well as a clear plateau,
the ‘knee’ of which would be a useful initial starting point as a statistical
sample of uncertain lenses. Having calibrated probabilities also allows
comparison of objects inspected by different classifiers. This would allow
rank-ordering across different samples of objects, not necessarily all seen by
the same classifier which could be used for identifying candidates for follow-
up and may be useful for selecting the forthcoming 4MOST sample (Collett et
al., 2023).

(a) An ensemble of 6 classifiers (solid lines) compared to only the citizen
science classifier (dashed)

(b) An ensemble of the 5 neural networks (solid lines) compared to a single
network (Neural Network 2, dashed)

Figure 8: A plot of the expected number of true and false positives for an
LSST; the total number of lenses has been fixed to $10^{5}$. These plots use
the false-positive and true-positive rates as a function of the score
threshold for the HSC classifiers used in this work (with A-B grade lenses as
‘true’ lenses), assuming the classifier performance in LSST matches those
applied to HSC here. The curves have been cutoff when there were $<5$ lenses
or non-lenses remaining in the test set, to reduce the effect of small-number
statistics. Figure 9: A plot of the effective sample size of the calibrated
individual and ensemble classifiers presented above, as a function of total
sample size. The effective sample size was determined by rank-ordering the
objects in order of calibrated probability, then cumulatively summing them.
The ensemble classifier used here is comprised of all 6 classifiers, and has
been further recalibrated by isotonic regression after forming the ensemble.
The ensemble classifier can produce the largest effective sample size out of
all the classifiers.

Given that at for high completeness, the number of expected false positives
from LSST is still large, we discuss here further avenues for improving lens
searches. Due to the huge data volume anticipated with forthcoming surveys,
such alternatives must still minimise human intervention. Individual lens
finding algorithms will continue to improve: the best network in Canameras et
al. (2023) further improves upon that of Cañameras et al. (2021) used here
(the former was not applied to the whole HSC survey) and suggest
TPR${}_{0}\sim 60\%$ could be achievable. With respect to ensemble
classifiers, the large differences in the rank ordering of objects between
classifiers (Figure 5) suggest classifiers trained on a diverse range of
training data find different types of non-lenses easier/more difficult to
classify. This suggests that larger ensembles could offer further improvement.
We investigated this with the classifiers used in this work. We measured the
AUROC, $FPR_{50}$ (False positive rate at $50\%$ completeness) and $TPR_{-3}$
(Completeness at an FPR of $10^{-3}$) as a function of number of classifiers
in the ensemble, averaged over the combinations of available classifiers. We
find the primary benefit in these metrics is achieved when adding $\geq 3$
classifiers, though the ensemble continued to improve up to the 6 used in this
work. While very diverse, some aspects of the lens classifiers were similar
(e.g. the lens model and non-lens sample used in Network 2 and 3) causing
correlation between their scores. Although beyond the scope of this current
work, it would be interesting to compare the performance of an ensemble of
classifiers trained with non-overlapping samples of non-lenses with a single
network trained on the whole dataset, the benefit of an ensemble being that it
should be able to mitigate the effects of any biases in an individual
classifier.
The citizen science classifier (which used lens subtracted imaging) performed
the best out of the classifiers used here however without significant
increases in citizen numbers or object selection, it is likely unfeasible to
apply citizen science to the whole LSST field (Marshall et al., 2016).
Nevertheless, citizen science will still remain a key tool for lens-finding in
this new era. The use of citizen-informed networks (i.e. active learning) such
as applied to galaxy morphology classification (Walmsley et al., 2020) could
be used instead to further improve automated classification. Similarly, using
citizen science classifications as labelled training data for lens-finding
neural networks would increase the size of a training set which used
exclusively real images. Currently the number of known strong lenses is too
small to train a network, which therefore typically rely on simulated lenses,
however once data begins to arrive from wide-area surveys, citizen
classifications of the real data could provide large volumes of labelled
training data. Simultaneous fast lens modelling such as developed by Park et
al. (2021); Wagner-Carena et al. (2021) could also help improve current lens
finding algorithms.
A key component of the methodology described above is the presence of a
ground-truth from which to base the calibration. In this work, we compiled
expert grades from a range of sources, in order to maximise the number of
objects with assigned grades. While it may initially appear that having a
known ground-truth for a random sample of objects would be optimal for
calibration, in practice the high-score regime is where calibration is most
important and the calibration mapping changes the most rapidly. Citizen
science could be used here to provide high quality lens candidates to use as a
ground truth for network calibration. For LSST, this ground-truth could then
be used for calibrating automated methods such as neural networks across the
whole survey area.

##  6 Conclusion

This work had two fundamental aims: 1) to provide calibrated probabilities for
a sample of galaxy cutouts that a given galaxy is a strong lens system and 2)
to combine neural network and citizen science classifiers into an ensemble
classifier, to maximise the purity of the resulting sample, without
compromising on completeness. We used 6 classifiers (1 citizen science search
and 5 neural networks) previously applied to HSC data, chosen as a proxy for
the forthcoming LSST. Having achieved these aims, our conclusions are as
follows:

  1. 1.

It is possible to post-process the outputs of a given lens classifier to
produce accurate calibrated probabilities that a given classified object is a
lens. The original outputs of such classifiers are not _a priori_ calibrated
probabilities, and should not be treated as such in a statistical analysis.

  2. 2.

There is very little correlation between classifiers of the scores of non-
lenses. This can be used advantageously to help remove false positives, as
different classifiers find certain systems easier/more difficult to classify.

  3. 3.

An ensemble classifier can provide improved classification above its
constituent components. For an FPR of $10^{-3}$, an ensemble classifier can
increase the completeness of the resultant sample from $34\%$ to $46\%$.

  4. 4.

Given $\sim 10^{5}$ strong lenses in LSST, further improvements in scalable
lens finding methods will be needed in order to achieve completenesses $>50\%$
without significant contamination by false positives.

## Acknowledgements

PH would like to thank the KIPAC Strong Lensing group for very useful
discussions and Alessandro Sonnenfeld and Yiping Shu for providing data for
this project.
PH acknowledges funding from the Science and Technology Facilities Council,
Grant Code ST/W507726/1. AV acknowledges support from the Science and
Technology Facilities Council, grant code ST/S006168/1 and ST/X00127X/1. This
research is supported in part by the Max Planck Society, and by the Excellence
Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy – EXC- 2094 –
390783311. This project has received funding from the European Research
Council (ERC) under the European Unions Horizon 2020 research and innovation
programme (LENSNOVA: grant agreement No 771776). A.T.J. is supported by the
Program Riset Unggulan Pusat dan Pusat Penelitian (RU3P) of LPIT Insitut
Teknologi Bandung 2023. This work is supported by JSPS KAKENHI Grant Number
JP20K14511.

## Data Availability

Data underlying this article were provided by the sources listed in Section 2.
Requests for such data should be made to the corresponding authors.

## References

  * Aihara et al. (2018) Aihara H., et al., 2018, PASJ, 70, S8
  * Aihara et al. (2019) Aihara H., et al., 2019, PASJ, 71, 114
  * Aihara et al. (2022) Aihara H., et al., 2022, PASJ, 74, 247
  * Akeson et al. (2019) Akeson R., et al., 2019, arXiv e-prints, p. arXiv:1902.05569
  * Andika et al. (2023) Andika I. T., et al., 2023, arXiv e-prints, p. arXiv:2307.01090
  * Bolton et al. (2006) Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Moustakas L. A., 2006, ApJ, 638, 703
  * Bolton et al. (2008) Bolton A. S., Burles S., Koopmans L. V. E., Treu T., Gavazzi R., Moustakas L. A., Wayth R., Schlegel D. J., 2008, ApJ, 682, 964
  * Browne et al. (2003) Browne I. W. A., et al., 2003, MNRAS, 341, 13
  * Brownstein et al. (2012) Brownstein J. R., et al., 2012, ApJ, 744, 41
  * Cañameras et al. (2021) Cañameras R., et al., 2021, A&A, 653, L6
  * Canameras et al. (2023) Canameras R., et al., 2023, arXiv e-prints, p. arXiv:2306.03136
  * Collett (2015) Collett T. E., 2015, ApJ, 811, 20
  * Collett et al. (2023) Collett T. E., et al., 2023, The Messenger, 190, 49
  * Euclid Collaboration et al. (2022) Euclid Collaboration et al., 2022, A&A, 662, A112
  * Faure et al. (2008) Faure C., et al., 2008, ApJS, 176, 19
  * Garvin et al. (2022) Garvin E. O., Kruk S., Cornen C., Bhatawdekar R., Cañameras R., Merín B., 2022, A&A, 667, A141
  * Geach et al. (2015) Geach J. E., et al., 2015, MNRAS, 452, 502
  * Geng et al. (2021) Geng S., Cao S., Liu Y., Liu T., Biesiada M., Lian Y., 2021, MNRAS, 503, 1319
  * Hartley et al. (2017) Hartley P., Flamary R., Jackson N., Tagore A. S., Metcalf R. B., 2017, MNRAS, 471, 3378
  * He et al. (2020) He Z., et al., 2020, MNRAS, 497, 556
  * Holloway et al. (2023) Holloway P., Verma A., Marshall P. J., More A., Tecza M., 2023, MNRAS, 525, 2341
  * Holwerda et al. (2019) Holwerda B. W., et al., 2019, AJ, 158, 103
  * Ivezić et al. (2019) Ivezić Ž., et al., 2019, ApJ, 873, 111
  * Jackson (2008) Jackson N., 2008, MNRAS, 389, 1311
  * Jacobs et al. (2017) Jacobs C., Glazebrook K., Collett T., More A., McCarthy C., 2017, MNRAS, 471, 167
  * Jacobs et al. (2019a) Jacobs C., et al., 2019a, ApJS, 243, 17
  * Jacobs et al. (2019b) Jacobs C., et al., 2019b, MNRAS, 484, 5330
  * Knabel et al. (2020) Knabel S., et al., 2020, AJ, 160, 223
  * Lanusse et al. (2018) Lanusse F., Ma Q., Li N., Collett T. E., Li C.-L., Ravanbakhsh S., Mandelbaum R., Póczos B., 2018, MNRAS, 473, 3895
  * Li et al. (2021) Li R., et al., 2021, ApJ, 923, 16
  * Lintott et al. (2008) Lintott C. J., et al., 2008, MNRAS, 389, 1179
  * Marshall et al. (2016) Marshall P. J., et al., 2016, MNRAS, 455, 1171
  * More et al. (2016) More A., et al., 2016, MNRAS, 455, 1191
  * Moustakas (2012) Moustakas L., 2012, The Master Lens Database and The Orphan Lenses Project, HST Proposal ID 12833. Cycle 20
  * Myers et al. (2003) Myers S. T., et al., 2003, MNRAS, 341, 1
  * Park et al. (2021) Park J. W., Wagner-Carena S., Birrer S., Marshall P. J., Lin J. Y.-Y., Roodman A., LSST Dark Energy Science Collaboration 2021, ApJ, 910, 39
  * Pascale et al. (2022) Pascale M., et al., 2022, ApJ, 938, L6
  * Pawase et al. (2014) Pawase R. S., Courbin F., Faure C., Kokotanekova R., Meylan G., 2014, MNRAS, 439, 3392
  * Petrillo et al. (2017) Petrillo C. E., et al., 2017, MNRAS, 472, 1129
  * Petrillo et al. (2019) Petrillo C. E., et al., 2019, MNRAS, 484, 3879
  * Platt (2000) Platt J., 2000, Adv. Large Margin Classif., 10
  * Powell et al. (2022) Powell D. M., Vegetti S., McKean J. P., Spingola C., Stacey H. R., Fassnacht C. D., 2022, MNRAS, 516, 1808
  * Powell et al. (2023) Powell D. M., Vegetti S., McKean J. P., White S. D. M., Ferreira E. G. M., May S., Spingola C., 2023, MNRAS, 524, L84
  * Rojas et al. (2022) Rojas K., et al., 2022, A&A, 668, A73
  * Rojas et al. (2023) Rojas K., et al., 2023, MNRAS, 523, 4413
  * Sahu et al. (2022) Sahu K. C., et al., 2022, ApJ, 933, 83
  * Savary et al. (2022) Savary E., et al., 2022, A&A, 666, A1
  * Schaefer et al. (2018) Schaefer C., Geiger M., Kuntzer T., Kneib J. P., 2018, A&A, 611, A2
  * Seidel & Bartelmann (2007) Seidel G., Bartelmann M., 2007, A&A, 472, 341
  * Shu et al. (2022) Shu Y., Cañameras R., Schuldt S., Suyu S. H., Taubenberger S., Inoue K. T., Jaelani A. T., 2022, A&A, 662, A4
  * Sonnenfeld et al. (2018) Sonnenfeld A., et al., 2018, PASJ, 70, S29
  * Sonnenfeld et al. (2019) Sonnenfeld A., Jaelani A. T., Chan J., More A., Suyu S. H., Wong K. C., Oguri M., Lee C.-H., 2019, A&A, 630, A71
  * Sonnenfeld et al. (2020) Sonnenfeld A., et al., 2020, A&A, 642, A148
  * Spergel et al. (2015) Spergel D., et al., 2015, arXiv e-prints, p. arXiv:1503.03757
  * Sugiyama et al. (2008) Sugiyama M., Suzuki T., Nakajima S., Kashima H., Bünau P. V., Kawanabe M., 2008, Annals of the Institute of Statistical Mathematics, 60, 699
  * Thuruthipilly et al. (2022) Thuruthipilly H., Zadrozny A., Pollo A., Biesiada M., 2022, A&A, 664, A4
  * Wagner-Carena et al. (2021) Wagner-Carena S., Park J. W., Birrer S., Marshall P. J., Roodman A., Wechsler R. H., LSST Dark Energy Science Collaboration 2021, ApJ, 909, 187
  * Walmsley et al. (2020) Walmsley M., et al., 2020, MNRAS, 491, 1554
  * Zadrozny & Elkan (2002) Zadrozny B., Elkan C., 2002, Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining

##  Appendix A ROC Curves for Different Ensemble Methods

In Figure 10, we show the ROC curves for an ensemble classifier, combining the
classifiers using the methods described in Section 3.3. We find that, in the
low false-positive region, the Independent Bayesian methods provide the
greatest completeness.

Figure 10: ROC Curves for a selection of different classifier-combination
methods. The best performing methods at the low false-positive end were the
Independent Bayesian methods (with/without generating a network-only ensemble
first), which are used in the rest of the paper.

##  Appendix B ROC Curves using A-B grade candidates with/without clusters

Here we show the ROC curves for the individual lens classifiers when altering
the ground-truth for ‘true lenses’. We observe a small narrowing between the
performance of the neural networks and that of the citizen science classifier
when cluster-scale lenses are excluded.

(a) Including A-B grade lenses as ‘true lenses’, without excluding cluster-
scale candidates.

(b) Including A-B grade lenses as ‘true lenses’, and excluding cluster-scale
candidates.

Figure 11: Receiver operating characteristic (ROC) curve for the individual
lens classifiers, applied to the whole cross-matched dataset. Note the y-axes
are linear, while the x-axes are logarithmic.

##  Appendix C Calibration Validation curves with A-C grade candidates,
excluding clusters

Here we show the validation plots for each calibration method, when treating
A-C grade lenses (but not cluster-scale systems) as ‘true lenses’. Overall we
see improved calibration accurately reaching higher calibrated probabilities.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12: Validation of the calibration curves, applied to a separate test
set of graded images. A-C grade lenses are counted as ‘true lenses’ for these
plots, and cluster-scale candidates receiving A-C grades have been excluded.