# Use the 4S (Signal-Safe Speckle Subtraction): Explainable Machine Learning
reveals the Giant Exoplanet AF Lep b in High-Contrast Imaging Data from 2011

Markus J. Bonse ETH Zurich, Institute for Particle Physics & Astrophysics,
Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland  Max Planck Institute for
Intelligent Systems, Max-Planck-Ring 4, 72076 Tübingen, Germany  Timothy D.
Gebhard ETH Zurich, Institute for Particle Physics & Astrophysics, Wolfgang-
Pauli-Str. 27, 8093 Zurich, Switzerland  Max Planck Institute for Intelligent
Systems, Max-Planck-Ring 4, 72076 Tübingen, Germany  Max Planck ETH Center for
Learning Systems, Max-Planck-Ring 4, 72076 Tübingen, Germany  Felix A. Dannert
ETH Zurich, Institute for Particle Physics & Astrophysics, Wolfgang-Pauli-Str.
27, 8093 Zurich, Switzerland  National Center of Competence in Research
PlanetS 111www.nccr-planets.ch Olivier Absil University of Liège, STAR
Institute, Allée du Six Août 19C, 4000 Liège, Belgium  Faustine Cantalloube
Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France  Valentin
Christiaens University of Liège, STAR Institute, Allée du Six Août 19C, 4000
Liège, Belgium  KU Leuven, Institute for Astronomy, Celestijnenlaan 200D,
Leuven, Belgium  Gabriele Cugno ETH Zurich, Institute for Particle Physics &
Astrophysics, Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland  University of
Michigan, Department of Astronomy, Ann Arbor, MI 48109, USA  Emily O. Garvin
ETH Zurich, Institute for Particle Physics & Astrophysics, Wolfgang-Pauli-Str.
27, 8093 Zurich, Switzerland  Jean Hayoz ETH Zurich, Institute for Particle
Physics & Astrophysics, Wolfgang-Pauli-Str. 27, 8093 Zurich, Switzerland
Markus Kasper European Southern Observatory, Garching bei München, Germany
Elisabeth Matthews Max Planck Institute for Astronomy, Königstuhl 17, 69117
Heidelberg, Germany  Bernhard Schölkopf Max Planck Institute for Intelligent
Systems, Max-Planck-Ring 4, 72076 Tübingen, Germany  Department of Computer
Science, ETH Zurich, 8092 Zurich, Switzerland  Sascha P. Quanz ETH Zurich,
Institute for Particle Physics & Astrophysics, Wolfgang-Pauli-Str. 27, 8093
Zurich, Switzerland  National Center of Competence in Research PlanetS
222www.nccr-planets.ch

###### Abstract

The main challenge of exoplanet high-contrast imaging (HCI) is to separate the
signal of exoplanets from their host stars, which are many orders of magnitude
brighter. HCI for ground-based observations is further exacerbated by speckle
noise originating from perturbations in the Earth’s atmosphere and
imperfections in the telescope optics. Various data post-processing techniques
are used to remove this speckle noise and reveal the faint planet signal.
Often, however, a significant part of the planet signal is accidentally
subtracted together with the noise. In the present work, we use explainable
machine learning to investigate the reason for the loss of the planet signal
for one of the most used post-processing methods: Principal Component Analysis
(PCA). We find that PCA learns the shape of the telescope point spread
function for high numbers of PCA components. This representation of the noise
captures not only the speckle noise, but also the characteristic shape of the
planet signal. Building upon these insights, we develop a new post-processing
method (4S) that constrains the noise model to minimize this signal loss. We
apply our model to 11 archival HCI datasets from the VLT-NACO instrument in
the L’-band and find that our model consistently outperforms PCA. The
improvement is largest at close separations to the star ($\leq\qty{4}{}$)
providing up to 1.5 magnitudes deeper contrast. This enhancement enables us to
detect the exoplanet AF Lep b in data from 2011, 11 years before its
subsequent discovery. We present updated orbital parameters for this object.

Direct imaging (387), Astronomy data reduction (1861), Exoplanets (498),
Interdisciplinary astronomy (804), High angular resolution (2167), Astronomy
image processing (2306)

tablenum

##  1 Introduction

The primary challenge of exoplanet high contrast imaging (HCI) is to separate
and detect the light of faint exoplanets from their host stars, which are many
orders of magnitude brighter. During an observation, the light of the star
gets scattered in the image plane, creating a halo of highly structured
speckle noise superimposed on the planet signal. Speckles tend to mimic the
expected shape of the planet, which makes it difficult to distinguish them
from a real planet signal. There are two types of speckle noise (Males et al.,
2021): quasi-static speckles, caused by nanometer scale imperfections in the
telescope optics, and atmospheric speckles, resulting from the turbulence of
the Earth’s atmosphere. Over the last decade, considerable efforts have been
spent to push the limits of HCI observations and reduce speckle noise. In this
context, one key aspect is the development of faster and more accurate eXtreme
Adaptive Optics (XAO) systems. Examples of such systems are the Gemini Planet
Imager (GPI), the Spectro-Polarimetric High-contrast Exoplanet REsearch
instrument (SPHERE), the Subaru Coronagraphic Extreme Adaptive Optics (SCExAO)
and MagAO-X (Macintosh et al., 2014; Jovanovic et al., 2015; Males et al.,
2018; Beuzit et al., 2019). XAO uses new detectors, advanced wavefront sensors
and deformable mirrors with higher actuator counts, along with faster
computing hardware and control software. For HCI, they are often coupled with
different types of coronagraphs to efficiently suppress the stellar point
spread function (PSF) (Kenworthy et al., 2007; Mawet et al., 2009; Snik et
al., 2012; N’Diaye et al., 2015; Otten et al., 2017).

Equally crucial for achieving deeper contrast and for removing speckle noise
is the development of post-processing techniques (see Cantalloube et al.,
2021, for an overview). Post-processing techniques are combined with dedicated
observing strategies like angular differential imaging (ADI) (Marois et al.,
2006), spectral differential imaging (SDI) (Racine et al., 1999; Sparks &
Ford, 2002), reference star differential imaging (RDI) (Mawet et al., 2012;
Ruane et al., 2019) or polarimetric differential imaging (PDI) (Kuhn et al.,
2001; Quanz et al., 2011) which are designed to enable planet signals to be
distinguished from speckle noise. A major breakthrough for the processing of
HCI data was the introduction of principal component analysis (PCA / KLIP;
Amara & Quanz, 2012; Soummer et al., 2012). Since its adoption to HCI, PCA has
formed the basis for the analysis of many HCI surveys (Nielsen et al., 2019;
Launhardt et al., 2020; Janson et al., 2021; Langlois et al., 2021; Cugno et
al., 2023) and allowed the detection of several new planets (Macintosh et al.,
2015; Chauvin et al., 2017).

However, after more than a decade of HCI observations, only a relatively small
number of objects have been found (see Currie et al. 2023 for an overview and
table of sub-stellar companions) and the results of several surveys have
consistently shown that super-Jovian planets with masses between 513 at large
separations (10100) are extremely rare (between $\sim$ 0.32.2 depending on
stellar type; see Vigan et al. 2021). Only recently, taking advantage of the
latest generation of instruments, it was possible to search for less massive
planets at closer separations. While the sample size is still small, new
planets have been detected in this regime (e.g., 51 Eri b; Macintosh et al.
2015, AF Lep b; Mesa et al. 2023; Franson et al. 2023; De Rosa et al. 2023)
hinting at an elevated occurrence rate towards smaller masses and separations.
On the other hand, radial velocity (RV) studies have found that the occurrence
rate of giant planets peaks around 3.6(2.0:1.8), decreasing closer to the star
and further out (Fulton et al., 2021). However, the decrease beyond 58 has not
been confirmed and is still the subject of ongoing research (Lagrange et al.,
2023). It is therefore crucial to push down the limits of HCI observations
close to the star to bridge the gap and create an overlap between RV and HCI
surveys. This will provide empirical access to an important parameter space
for exoplanet population and formation models and allow for direct imaging
follow-up of exoplanets detected by RV surveys.

Unfortunately, close to the star, it becomes increasingly difficult to
suppress speckle noise using data post-processing techniques. If we use
methods like PCA, a significant fraction of the planet signal is accidentally
subtracted together with the noise. While at large separations the loss of
planet signal is moderate, close to the star more than 95 of the planet signal
might get lost (e.g., Gebhard et al., 2022; Bonse et al., 2023). This effect
is referred to as planet over- and self-subtraction (Pueyo, 2016) and
represents the main bottleneck for post-processing of HCI data.

#### Contributions:

In this paper, we gain a deep understanding of the reason for the loss of
planet signal in PCA and related methods and use it to develop a new post-
processing algorithm. To this end, we revisit the mathematical background of
PCA and give an overview of more recent post-processing methods in Section 2.
By using methods from the field of explainable machine learning, we explore
the reason for planet over- and self-subtraction (Section 3). We identify two
main limitations of PCA, which we use as a basis to develop our new post-
processing algorithm 4S (Section 4). A detailed analysis of the contrast
performance of 4S on archival data from the VLT-NACO instrument is given in
Section 5. We find that 4S reaches significantly deeper contrast at close
separations to the star, and learns a model of the noise that matches our
understanding of the speckle noise. Using 4S we can reveal the faint exoplanet
AF Lep b in archival VLT-NACO data from 2011 (Section 6) which is not
detectable with PCA. In Sections 7 and 8, we outline directions for future
work and provide concluding remarks.

#### Reproducibility:

Along with this paper, we are releasing our algorithm 4S as a Python package
called fours. The code is available on GitHub and comes with a documentation
page on ReadtheDocs. All raw data (Bonse, 2024a) and intermediate results
(Bonse, 2024b) to reproduce our plots are available on Zenodo.

##  2 PCA Revisited

Figure 1:  _Top row:_ Illustration of how PCA is applied to ADI data in HCI.
The planetary signal is shown in blue, while the stellar speckle noise is
represented by the red splash. Details about the steps ( 1 to 5) are given in
the text. _Bottom row:_ Examples of the first 5 principal components are
shown, estimated for the dataset #5 (HD 22049) in Table 1. Further details on
the data are given in Section 5. The first components of PCA are often
interpretable. For example, $\mathbf{u}_{1}$ in the given example models the
overall brightness variations of the post-coronagraphic PSF. The components
$\mathbf{u}_{2}$ and $\mathbf{u}_{3}$ represent horizontal and vertical
variations caused by small errors in the centering of the star behind the
coronagraph.

The concept of using PCA (or Karhunen–Loève transformation) in HCI is to
estimate the stellar PSF and speckle noise to subtract it from the science
data. For this purpose, PCA utilizes the temporal variability and brightness
of the speckles. A typical HCI dataset consists of a sequence of individual
exposures $\mathbf{x}_{t,\text{raw}}\in\mathbb{R}^{D}$ with pixel resolution
$D=n_{i}\cdot n_{j}$ taken as a function of time $t=1,\ldots,N$. If we ignore
background and detector noise, each image is composed of light from the star
$\mathbf{s}_{t}$ and possibly light from the planet $\mathbf{p}_{t}$:

$\mathbf{x}_{t,\text{raw}}=\mathbf{s}_{t}+\mathbf{p}_{t}\,.$ (1)

The stellar light $\mathbf{s}_{t}$ outshines the planet signal
$\mathbf{p}_{t}$ by several orders of magnitude. As we go through the temporal
stack of science images, the PSF and speckle pattern changes. This means that
the strongest brightness variations in the data are due to the light of the
star and not due to the planet. This data characteristic is the foundation of
speckle-subtraction with PCA: PCA is a linear method which tries to find a
low-dimensional representation of the data with dimension $K\ll D$, which only
captures the strong covariances of the data. If the planet signal is weak
compared to the speckle noise, it should thus not be part of this low-
dimensional representation. If we reconstruct the data from the low-
dimensional space, we get an estimate of the speckle noise that can be
subtracted from the science data. Applying PCA to HCI data involves the
following steps (see Figure 1):

  1. 1

Data preparation: The temporal average of the science data
$\overline{\mathbf{x}}$ is subtracted from each individual science frame:

$\mathbf{x}_{t}=\mathbf{x}_{t,\text{raw}}-\overline{\mathbf{x}}\,.$

The mean subtracted 2D science frames are then flattened into one dimensional
vectors and stacked into a data matrix $\mathbf{X}\in\mathbb{R}^{N\times D}$:

$\mathbf{X}=\begin{bmatrix}\mathbf{x}_{1}^{T}\\\ \mathbf{x}_{2}^{T}\\\
\vdots\\\
\mathbf{x}_{N}^{T}\end{bmatrix}=\begin{bmatrix}x_{1,1}&x_{1,2}&\ldots&x_{1,D}\\\
x_{2,1}&x_{2,2}&\ldots&x_{2,D}\\\ \vdots\\\
x_{N,1}&x_{N,2}&\ldots&x_{N,D}\end{bmatrix}\,.$

Each row in $\mathbf{X}$ corresponds to one mean subtracted science
frame.333We choose the notation to be consistent with our Python
implementation. Our notation might deviate from the conventions used in other
papers (e.g., Lewis et al., 2023).

  2. 2

Calculation of the projection matrix: We choose the dimensionality of the low-
dimensional space $K$ and calculate the PCA basis (i.e., the projection matrix
$\mathbf{U}\in\mathbb{R}^{D\times K}$) which transforms our data into the low-
dimensional space:

$\mathbf{U}=\begin{bmatrix}\mathbf{u}_{1}&\ldots&\mathbf{u}_{K}\end{bmatrix}\,.$

The column vectors of $\mathbf{U}$ are the principal components
$\mathbf{u}_{1},\ldots,\mathbf{u}_{K},\forall\mathbf{u}\in\mathbb{R}^{D}$.
Details on how to calculate these components are given in the next section.

  3. 3

Dimensionality reduction: We project our data into the low dimensional space
$\mathbf{x}_{t}\rightarrow\mathbf{a}_{t},\quad\mathbf{a}_{t}\in\mathbb{R}^{K}$:

$\mathbf{X}\cdot\mathbf{U}=\begin{bmatrix}\mathbf{a}_{1}^{T}\\\
\mathbf{a}_{2}^{T}\\\ \vdots\\\
\mathbf{a}_{N}^{T}\end{bmatrix}=\mathbf{A}\in\mathbb{R}^{N\times K}\,.$ (2)

If the dimensionality $K$ is chosen carefully, the low-dimensional
representation $\mathbf{A}$ should mostly capture the variations caused by the
star.

  4. 4

Data reconstruction: We reconstruct our data from the low-dimensional space
$\mathbf{a}_{t}\rightarrow\hat{\mathbf{x}}_{t}$,
$\hat{\mathbf{x}}_{t}\in\mathbb{R}^{D}$:

$\mathbf{A}\cdot\mathbf{U}^{T}=\begin{bmatrix}\mathbf{\hat{x}}_{1}^{T}\\\
\mathbf{\hat{x}}_{2}^{T}\\\ \vdots\\\
\mathbf{\hat{x}}_{N}^{T}\end{bmatrix}=\mathbf{\hat{X}}\,.$

Since $\mathbf{a}_{t}$ mostly captures the strong covariances of the data, we
can use its reconstruction to approximate the speckle noise:

$\hat{\mathbf{x}}_{t}\approx\mathbf{s}_{t}\,.$ (3)

  5. 5

Calculation of the residual: We subtract the noise estimate
$\hat{\mathbf{x}}_{t}$ for each mean-subtracted science frame
$\mathbf{x}_{t}$. In the case of angular differential imaging (ADI; Marois et
al. 2006), we have to de-rotate all frames such that North is up in all
images. The result is the so-called residual sequence $\mathbf{r}_{t}$:

$\mathbf{r}_{t}=\phi_{t}(\mathbf{x}_{t}-\hat{\mathbf{x}}_{t})\,,$ (4)

where $\phi_{t}$ is the de-rotation function, which uses the parallactic
angles to re-align North in all images. The final residual image
$\mathbf{r}_{\text{final}}$ can then be obtained by averaging along the time
axis:444Some implementations take the median along time instead of the mean,
which can be advantageous in the case of detector artifacts such as bad
pixels.

$\mathbf{r}_{\text{final}}=\frac{1}{N}\sum_{t=1}^{N}\mathbf{r}_{t}\,.$ (5)

Figure 2:  Variation in planet visibility as a function of principal
components ($K=20$, $K=60$, $K=120$, $K=300$). The shown dataset is dataset #5
(HD 22049) in Table 1. More details on data pre-processing are given in
Section 5. Two fake planets have been inserted into the data: one positioned
at $\sim$ 3 separation from the star with a contrast brightness of 9 (white
circle), and one at $\sim$ 10 with a contrast brightness of 12.5 (white dashed
circle). The inner planet reaches the highest S/N at $\sim K=60$ components
(left plot, solid line) and the outer planet at $\sim K=300$ components (left
plot, dashed line). This example highlights the challenge of choosing the
right number of principal components within a single dataset.

###  2.1 How to find the optimal projection?

Depending on the type of the dataset, the projection matrix $\mathbf{U}$ can
be calculated based on the science data directly (e.g., in case of ADI, Amara
& Quanz, 2012; Soummer et al., 2012) or based on reference stars (RDI, Ruane
et al., 2019). In this work, we focus on ADI, but the same concepts apply for
RDI. As motivated in the previous section, we want to choose $\mathbf{U}$ such
that the strongest variations of the data, that is, the speckle noise, are
captured in the low-dimensional space. This objective is equivalent to the
maximum variance derivation of PCA (see, e.g., chapter 12.1 of Bishop 2006 for
the details).

For simplicity, let us consider the projection onto the first component
$\mathbf{u}_{1}$. We want to choose $\mathbf{u}_{1}$ such that the variance
along the first axis of the low-dimensional space $\sigma^{2}(a_{t,1})$ is
large:

$\displaystyle\sigma^{2}(a_{t,1})$
$\displaystyle=\mathbf{u}_{1}^{T}\mathbf{C}\mathbf{u}_{1}\,,$

where $\mathbf{C}\in\mathbb{R}^{D\times D}$ is the covariance matrix of the
data

$\displaystyle\mathbf{C}=\frac{1}{T}\mathbf{X}^{T}\mathbf{X}\,.$

We add another constraint
$||\mathbf{u}_{1}||=\mathbf{u}_{1}^{T}\mathbf{u}_{1}=1$ to avoid
$||\mathbf{u}_{1}||\rightarrow\infty$. To enforce this constraint, we
introduce a Lagrange multiplier $\lambda_{1}$ which gives us the following
optimization objective:

$\operatorname*{arg\,max}_{\mathbf{u}_{1}}=\mathbf{u}_{1}^{T}\mathbf{C}\mathbf{u}_{1}+\lambda_{1}\left(1-\mathbf{u}_{1}^{T}\mathbf{u}_{1}\right)\,.$

If we set the derivative w. r. t. $\mathbf{u}_{1}$ equal to zero, we get:

$\mathbf{C}\mathbf{u}_{1}=\lambda_{1}\mathbf{u}_{1}\,.$

This means that $\mathbf{u}_{1}$ should be chosen to be an eigenvector of
$\mathbf{C}$. Using $\mathbf{u}_{1}^{T}\mathbf{u}_{1}=1$ we see that:

$\sigma^{2}(a_{t,1})=\mathbf{u}_{1}^{T}\mathbf{C}\mathbf{u}_{1}=\lambda_{1}\,.$

To maximize the variance of the projection, we should therefore choose the
eigenvector $\mathbf{u}$ with the largest eigenvalue. We sort all eigenvectors
of $\mathbf{C}$ by their eigenvalues:

$\displaystyle\begin{array}[]{ccccccccccc}\lambda_{1}&>&\lambda_{2}&>&\ldots&>&\lambda_{K}&>&\ldots&>&\lambda_{D}\\\
\mathbf{u}_{1}&&\mathbf{u}_{2}&&\ldots&&\mathbf{u}_{K}&&\ldots&&\mathbf{u}_{D}\end{array}\,,$

and select the first $K$ eigenvectors with the largest eigenvalues. In
practice, $\mathbf{U}$ is usually obtained by computing a singular value
decomposition (SVD) of the data matrix $\mathbf{X}$. This approach ensures
numerical stability and reduces computation time.

A big advantage of PCA is that the first components of the matrix $\mathbf{U}$
can be interpreted visually (see Figure 1). For example, the first component
$\mathbf{u}_{1}$ often models the general brightness variations of the stellar
PSF. The data used to calculate the components in Figure 1 was taken with a
vortex coronagraph (Mawet et al., 2013). Thus, the first component shows the
typical donut-shape of the post-coronagraphic PSF. The second and third
components show the response of the coronagraphic PSF to tip-tilt errors
(i.e., the first modes of the AO).

###  2.2 Limitations of PCA

PCA is widely used by the exoplanet imaging community and available in most
HCI packages such as VIP (Christiaens et al., 2023), PynPoint (Stolker et al.,
2019) or pyKLIP (Wang et al., 2015). Its appeal is mainly due to its
simplicity and to the fact that it comes with only a single hyperparameter:
the number of principal components $K$. However, this apparent simplicity
hides a complex interplay of several factors which makes it difficult to tune
$K$ in practice. The optimal choice is highly dataset-specific, influenced by
variables such as the instrument, observing mode, wavelength, and observing
conditions. Even within a single dataset, the optimal number of components can
vary as a function of the separation from the star (Meshkat et al., 2013).
This is to be expected, as low order wavefront modes, which affect the closest
separations to the star, vary the most. The further away from the center, the
higher the modal order and the larger the required $K$. An example of this
effect is shown in Figure 2. Two artificial planets are inserted into the same
dataset: one close to the star ($\sim$ 3) and another farther out ($\sim$ 10).
The closer planet is best seen if 60 principal components are selected,
whereas the farther one requires about 300 components to be fully visible.
Notably, at 300 components, the signal of the inner planet is completely lost.
This demonstrates a critical limitation of PCA: as the number of components
increases, PCA does not only model the noise but starts to include parts of
the planet signal, leading to what is known as planet self- and over-
subtraction (Pueyo, 2016). As discussed further in Section 5, the choice of
the number of components can affect the final contrast by several magnitudes.
Hence, selecting the right number of components is crucial; it can be the
difference between a detection and a non-detection.

###  2.3 Beyond PCA

The loss of planet signal in PCA has encouraged the development of more
advanced post-processing techniques. Nevertheless, many of these methods
either enhance PCA or incorporate its principles. Annular PCA, for instance,
applies PCA in concentric annular segments around the star (Absil et al.,
2013; Gomez Gonzalez et al., 2017). A specific exclusion criterion for ADI
prevents the signal of the planet from entering the PCA basis. Lewis et al.
(2023) extends PCA to utilize the temporal and spatial axis simultaneously.

After the speckle-subtraction the residual images are commonly analyzed with
meta-techniques which quantify detections and non-detections; for example
using $t$-tests (Mawet et al., 2014), STIM maps (Pairet et al., 2019) or
performance maps (Jensen-Clem et al., 2017). A recent overview of these
methods can be found in Bonse et al. (2023). Beyond the analysis of a single
residual, supervised machine learning techniques have been proposed (Gomez
Gonzalez et al., 2018; Cantero et al., 2023). These methods use several PCA
residuals with different numbers of components as input and train a neural
network to detect the planet. The RSM map by Dahlqvist et al. (2020) and
Dahlqvist & Absil (2021) replaces the temporal averaging in the last step of
PCA. It allows the combination of multiple PCA outputs alongside other post-
processing techniques.

Conceptually similar to PCA are other subtraction-based techniques which
replace the bottleneck of PCA (e.g., NNMF Arcidiacono & Simoncini 2018; Ren et
al. 2018 or LLSG Gomez Gonzalez et al. 2016). For example, LLSG adds an extra
sparse term to the model with the goal of limiting the inclusion of the planet
in the noise model. Daglayan et al. (2023) propose to change the objective
function to a L1 loss to account for the heavy-tailed distribution of the
speckle noise. The recently proposed ConStruct method uses non-linear auto-
encoders to locally predict the speckle noise (Wolf et al., 2023). Further
examples of subtraction-based techniques are LOCI and its variations
(Lafrenière et al., 2007; Marois et al., 2014; Wahhaj et al., 2015; Thompson &
Marois, 2021), or the more recently proposed HSR (Gebhard et al., 2022). Bonse
et al. (2018) proposed to use wavelets to reduce speckle noise in time. Most
of these methods discuss different exclusion criteria to avoid subtraction of
the planet signal.

Complementary to subtraction-based techniques are inverse-problem based
methods. Instead of modeling and subtracting the noise, these methods rely on
forward-modeling to search for the expected signature of planets in the data.
Examples are ANDROMEDA (Cantalloube et al., 2015), FMMF (Pueyo, 2016; Ruffio
et al., 2017), PACO (Flasseur et al., 2018, 2024) and TRAP (Samland et al.,
2021). In this context, FMMF is the direct extension of PCA to forward-
modeling.

Although the concepts behind these modern post-processing techniques are
diverse, they share the common idea of better using the knowledge about the
shape and behavior of the planet signal. However, despite all these efforts,
PCA still represents the main benchmark in practice: For example, Franson et
al. (2023) used it for the discovery of the exoplanet AF Lep b. Further, none
of these advanced algorithms have managed to completely avoid the problem of
planet signal loss, which is the main bottleneck for achieving deeper
detection limits.

##  3 Explaining signal loss in PCA

In the previous section, we observed that with an increasing number of
principal components, progressively more signal of the planet is lost. Is
there a specific component responsible for this effect? While the first
components of PCA can be interpreted visually (Figure 1), it is difficult to
understand the higher order modes. To better understand why we lose signal in
PCA, we propose to use methods from the field of eXplainable Artificial
Intelligence (XAI). With the recent boom in machine learning, there is a trend
towards increasingly complex deep network architectures. Unfortunately, the
decisions made by such models are hard to understand and therefore hard to
trust. XAI seeks to gain a better understanding of why and how machine
learning models make their decisions, thus giving us a glimpse into the black
box. A large number of XAI methods have been proposed in the literature (see
Ali et al. 2023, for a recent overview). Some of them are model-specific
(e.g., tailored to neural networks), while other methods are more generally
applicable (e.g., LIME by Ribeiro et al. 2016). In the following paragraph, we
will keep our explanation of the XAI method used in the present paper general,
to provide the basis for a better understanding of other post-processing
methods beyond PCA.

###  3.1 Saliency maps

A widely used approach to explain the predictions of machine learning methods
is the use of input gradients. The most basic approach was proposed by
Simonyan et al. (2013). It is often referred to as Vanilla Gradient and
straightforward to apply to subtraction-based post-processing methods in HCI:

Let us consider a single science image $\mathbf{x}_{t}$ at time $t$. Given a
speckle-subtraction method $S$, we estimate the speckle noise for this
particular image:

$S(\mathbf{x}_{t})=\hat{\mathbf{x}}_{t}.$

For each pixel position in the noise estimate $\hat{\mathbf{x}}_{t}$, we can
use input gradients to calculate so-called _saliency maps_. Saliency maps are
heatmaps which highlight areas in the input image $\mathbf{x}_{t}$ which
contribute the most to the noise estimate at the selected pixel position.
Given a pixel position, we calculate the derivative of $S_{(i,j)}$ with
respect to the input image $\mathbf{x}$ at the point $\mathbf{x}_{t}$:

$\mathbf{w}=\frac{\partial
S_{(i,j)}}{\partial\mathbf{x}}\bigg{|}_{\mathbf{x}_{t}}.$ (6)

The magnitude of the derivative $\mathbf{w}$ points to pixels in $\mathbf{x}$
that cause the largest difference on the noise estimate at the position
$(i,j)$. For each pixel position $(i,j)$ in the noise estimate, we get the
saliency map $M_{(i,j)}$ by first calculating the input gradient using
Equation 6, then rearranging the one-dimensional vector $\mathbf{w}$ back into
a two-dimensional image, and, finally, taking the absolute value:

$M_{(i,j)}=\left|\mathbf{w}_{(i,j)}\right|$

Saliency maps can be computed for any post-processing method $S$ for which we
can calculate the gradient in Equation 6. For non-linear methods such as
ConStruct (Wolf et al., 2023) or SODINN (Gomez Gonzalez et al., 2018), more
advanced methods such as SmoothGrad (Smilkov et al., 2017) or Integrated
Gradients (Sundararajan et al., 2017) are likely to give better results than
Vanilla Gradients.

###  3.2 Right for the wrong reasons

Figure 3:  Illustration of how to compute saliency maps for PCA. For each
pixel position in the noise estimate, we can compute one saliency map. Each
saliency map highlights the areas in the science frame that are used to obtain
the noise estimate at the respective position. Since the noise model of PCA is
linear, every pixel in the noise estimate is a weighted sum of the pixels in
the input image. The saliency map shows the absolute value of these weights.
Two examples of such saliency maps are shown in the middle row (300 principal
components). The example used in the text is marked in orange. For large
numbers of principal components, the saliency map converges to the shape of
the PSF. A gallery of saliency maps for different separations and number of
principal components is given in Figure 14.

In the case of PCA, the noise subtraction method $S$ is just the projection
onto the low-dimensional space and its reconstruction:

$S_{\text{PCA}}(\mathbf{x}_{t})=\mathbf{x}^{T}_{t}\mathbf{U}\mathbf{U}^{T}=\mathbf{\hat{x}}_{t}$
(7)

Given the number of principal components $K$, we compute $\mathbf{U}$ using
the procedure explained in Section 2.1. Assuming that $\mathbf{U}$ is fixed,
the saliency map of PCA is then given by:

$M_{(i,j),\text{PCA}}=\left|\mathbf{U}\mathbf{U}^{T}\right|_{[:,n_{i}j+i]}$
(8)

Each pixel value in the noise estimate of PCA is a weighted sum of the pixel
values in the science frame. The weights of this sum are the column vectors of
the matrix $\mathbf{U}\mathbf{U}^{T}$. This means that _each column_ in
$\mathbf{U}\mathbf{U}^{T}$ gives us an explanation for which information is
used to estimate the noise at _one output pixel position_. To get the saliency
map for the pixel $(i,j)$, we have to select the ($n_{i}j+i$)-th column vector
of $\mathbf{U}\mathbf{U}^{T}$ and take its absolute value. Note that since PCA
is a linear model, the saliency map for a given position is the same for all
science frames along time.

Figure 3 illustrates the concept and provides examples of saliency maps for
PCA based on dataset #4 (HD 22049) (see Section 5 for more details on the
dataset). To better understand the loss of planet signal in PCA, let us assume
that the planet signal in the science frame (Figure 3, top left image) is at
the position marked by the blue circle. If parts of the planet signal are lost
in the subtraction step, it means that some signal leaks into the noise
estimate at the same position (top right image, orange $+$). We can use
Equation 8 to compute the saliency map for this position (Figure 3,middle row,
left image). It turns out that most of the information used to estimate the
noise at the planet’s position is only taken from the input image at the
planet’s position. This means that PCA just copies the planet signal from the
input image to the noise estimate. The same pattern can be observed for other
pixel positions which do not contain the planet signal (Figure 3, middle row,
right image). This means that even if the planet’s position changes as a
function of time (e.g., in the case of ADI), it is still subtracted.

It is interesting to note that the shape of the information used corresponds
to the shape of the unsaturated PSF. We know that, to first order, the
speckles follow the shape of the PSF (i.e., discarding wavelength-smearing
effects). This means that PCA successfully learns what the speckles look like.
Since, not only the speckles follow the shape of the PSF, but also the signal
of the planet, the latter is partially subtracted.

It should be noted that this problem occurs regardless of whether the planet
is actually in the data. That is, even if we build the PCA basis using data
from a different star (RDI), we still get a model that can subtract the planet
signal. As the number of components increases, the projection
$\mathbf{U}\mathbf{U}^{T}$ converges to a de facto identity by learning the
shape of the PSF. This convergence prevents us from learning more complex
noise patterns.

#### Quantification of the effect

The convergence rate of PCA to the shape of the PSF is very dataset- and
instrument-specific, and depends on whether the noise is dominated by
speckles. For reasons explained later in Section 5, we confine our analysis to
the VLT-NACO data taken in the $L^{\prime}$ filter, leaving research on other
instruments for future work.

Figure 4: Convergence of the saliency map to the shape of the PSF as a
function of the number of principal components used. At closer separations
from the star, the convergence is faster, resulting in a greater loss of the
planet signal with a smaller number of components.

To quantify how many components are needed to converge to the shape of the
PSF, we compute $\mathbf{U}\mathbf{U}^{T}$ for different numbers of principal
components $K\in[1,...,500]$. We calculate saliency maps for pixel positions
at 1.5, 4 and 8, and estimate the match between the saliency map and the
unsaturated PSF template. More precisely, we compute a scalar product between
the local saliency map and the unsaturated PSF and normalize it by the maximum
value over all separations and components. The result of this experiment is
shown in Figure 4. For separations close to the star, PCA converges much
faster to the shape of the PSF than further out. This observation is in line
with the fake planet experiment shown in Figure 2. While the signal of the
fake planet close to the star is lost at 300 components, the planet further
out is still clearly visible. The speckle intensity gradually decreases with
separation from the star. As PCA models the strongest covariances first, the
convergence to the shape of the PSF is faster close to the star. As a result,
the risk of losing parts of the planet signal is also higher close to the
star.

###  3.3 Misguided loss function

The origin of the planet signal loss described in the previous section occurs
regardless of whether a planet signal is present in the data used to compute
the PCA basis. It is therefore often referred to as planet over-subtraction.
In the case of ADI datasets, in addition to over-subtraction, the planet
signal can also directly affect the PCA basis. This second effect is known as
planet self-subtraction and can lead to even greater signal loss. In Section
2.1, we explained how to derive PCA by maximizing the variance in the low-
dimensional space. This derivation is well motivated for HCI, since the
speckles are usually much brighter and more variable than the planet signal.
An alternative derivation of PCA is the minimum reconstruction error
formulation. We want to find a linear dimensionality reduction that minimizes
the following reconstruction loss:

$\mathcal{L}_{\text{PCA}}=\frac{1}{N}\sum_{t=1}^{N}\underbrace{\left\|\mathbf{x}_{t}-\mathbf{\hat{x}}_{t}\right\|^{2}}_{\mathbf{s}_{t}+\mathbf{p}_{t}-\mathbf{\hat{x}}_{t}}.$
(9)

That is, we want our reconstruction $\mathbf{\hat{x}}_{t}$ to be close to our
science data $\mathbf{x}_{t}$. Both the maximum variance derivation and the
minimum reconstruction error formulation yield the same algorithm (cf. chapter
12.1.2 in Bishop 2006). In the context of HCI, minimizing Equation 9 can be
problematic, as $\mathbf{x}_{t}$ contains both speckle noise and planet
signal. This means that $\mathcal{L}_{\text{PCA}}$ can only become small if
the noise estimate $\mathbf{\hat{x}}_{t}$ also fits the signal of the planet
$\mathbf{p}_{t}$. For large numbers of principal components, and for bright
companions, PCA tends to include the companion signal in its basis.

Figure 5:  Retained flux after post-processing with PCA as a function of
planet brightness and separation from the star. The grid was computed using
applefy. Each value in the grid is the result of 3 fake planet experiments
with planets inserted at different position angles. A significant fraction of
the planet signal is lost, especially close to the star.

To quantify this effect, we compute a throughput grid using the Python package
applefy (Bonse et al., 2023). A throughput grid illustrates the attenuation of
the planet signal due to the data post-processing. The grid is computed by
inserting fake companions with different known brightnesses and distances to
the star, and then measuring the fraction of the flux that remains after PSF
subtraction. The result of this experiment is shown in Figure 5. We observe
that for faint companions the retained flux is only a function of the
separation and does not depend on the brightness of the companion. This
gradual loss of planet signal closer to the star is caused by the convergence
of PCA to the shape of the PSF (see Section 3.2). For bright companions,
however, PCA begins to include parts of the signal into its basis, causing an
even greater loss of planet signal. The brightness at which planets are
affected depends on the number of principal components. If more components are
used, fainter companions start to be affected.

##  4 4S: Signal Safe Speckle Subtraction

In the previous section, we identified two main reasons for the loss of planet
signal in PCA: first, the convergence of the noise model to the telescope’s
PSF, and second, a problematic loss function that accidentally leads the noise
estimate to fit the planet signal. Based on these findings, we propose a novel
post-processing method that adopts the main steps of PCA but overcomes its
limitations. We have named this method 4S, short for Signal Safe Speckle
Subtraction. 4S is characterized by three innovations: First, a linear noise
model with a _right reason_ constraint, second, a new loss function that is
invariant to the planet signal, and third, a domain knowledge-informed
regularization.

###  4.1 The right reason constraint

In Section 3.2, we found that the noise model of PCA converges to a de facto
identity that reproduces the shape of the telescope’s PSF. This is problematic
as both the speckle noise and the planet signal follow this shape. Ideally, we
would like our noise model to specifically recognize the speckle noise, that
is, we want to find explanations for the speckle noise that are unique to the
noise and do not describe the planet signal. In the context of XAI, Ross et
al. (2017) proposed to add a _right reason_ penalty to guide the predictions
of neural networks using domain knowledge. The domain knowledge is represented
as a mask that penalizes _wrong_ explanations. The method has been
successfully applied in practice (Schramowski et al., 2020). Inspired by this
work, we propose a _right reason_ constraint to prevent the noise model of 4S
from learning the shape of the telescope’s PSF.

Like PCA (cf. Equation 7), 4S uses a linear model, but its projection matrix
$\mathbf{B}$ is modified to circumvent the bottleneck of PCA:

$\text{4S}(\mathbf{x}_{t})=\mathbf{x}_{t}^{T}\mathbf{B}=\mathbf{\hat{x}}_{t}\quad.$
(10)

For PCA, each column of the projection matrix $\mathbf{U}\mathbf{U}^{T}$ gives
us an explanation of which information of the input $\mathbf{x}_{t}$ is used
to estimate the noise at one pixel position in $\mathbf{\hat{x}}_{t}$.
Similarly, the columns of $\mathbf{B}\in\mathbb{R}^{D\times D}$ contain the
explanations of 4S. Following the idea of Ross et al. (2017), we limit the
information of $\mathbf{x}_{t}$ that can be used to predict the noise. For
each pixel position in the noise estimate $\mathbf{\hat{x}}_{(i,j)}$, we have
one parameter vector $\beta_{(i,j)}\in\mathbb{R}^{D}$ as a column of
$\mathbf{B}$. We design a circular _right reason_ mask for each pixel position
which sets the values near the de facto identity to zero:

$\mathbf{B}_{[:,n_{i}j+i]}=\beta_{(i,j)}\circ\mathbf{m}_{(i,j)},$

where $\circ$ is the element-wise multiplication. The right reason mask
$\mathbf{m}_{(i,j)}$ is defined as:555Note: The mask does not have to be
binary, which is beneficial for datasets with low resolution. Our
implementation makes use of circular apertures based on the photutils package.

$\mathbf{m}_{(i,j)}(i^{\prime},j^{\prime})=\begin{cases}0,&\text{if
}\sqrt{(i-i^{\prime})^{2}+(j-j^{\prime})^{2}}\leq d\\\
1,&\text{otherwise}\end{cases},$ (11)

where $d$ sets the radius of the mask. We choose $d$ relative to the size of
the PSF: $d=0.75\cdot\text{FWHM}$. This mask size ensures that the core of the
PSF is protected from the de facto identity, while local information about the
speckle noise is still preserved. A larger mask size allows us to protect even
more signal, but also weakens the noise model. The value of
$d=0.75\cdot\text{FWHM}$ is a compromise for which we have obtained good
results in practice.

The right reason mask should not be confused with the exclusion criteria used
by other post-processing techniques (e.g., Absil et al., 2013; Gomez Gonzalez
et al., 2017; Lafrenière et al., 2007; Wahhaj et al., 2015). The mask does not
block the signal of the planet, but instead forces the noise model to explore
other areas of the input image that explain the speckle noise. Our approach is
similar to the post-processing methods based on the HSR (Samland et al., 2021;
Gebhard et al., 2022). However, while Gebhard et al. (2022) select areas in
the input images that promise to contain useful information about the speckle
noise, we use the entire science frame, but exclude areas that correspond to
the shape of the telescope’s PSF near the de facto identity.

###  4.2 Signal invariant loss function

Figure 6:  Illustration of our new post-processing technique 4S. To overcome
the limitations of PCA, 4S makes use of three innovations: I. a right reason
mask to prevent over-subtraction (bottom left box), II. a loss function which
is invariant to the signal of the planet (bottom right box), and III. a
regularization based on domain knowledge to prevent overfitting. The main
steps 1 to 3 are explained in detail in Section 4.4.

A second limitation of PCA is its loss function (see Equation 9), which can
only become small if the noise estimate includes the planet signal. This is
problematic in the case of ADI, for which the planet signal is part of the
dataset used to compute the PCA basis. Different types of exclusion regions
(Lafrenière et al., 2007; Cantalloube et al., 2015; Thompson & Marois, 2021)
or explicit planet models (Samland et al., 2021) have been proposed in the
literature to protect the planet’s signal from becoming part of the noise
estimate. These methods, however, often come with additional hyperparameters
which are hard to optimize in practice. Therefore, instead of minimizing the
reconstruction loss, as done in Equation 9, we propose to minimize the
variance along the temporal axis in the de-rotated residual sequence:

$\mathcal{L}_{\text{4S}}=\frac{1}{N}\sum_{t=1}^{N}\left\|\mathbf{r}_{t}-\overline{\mathbf{r}}\right\|^{2},\quad\overline{\mathbf{r}}=\frac{1}{N}\sum_{t=1}^{N}\mathbf{r}_{t}\,.$
(12)

An illustration of why $\mathcal{L}_{\text{4S}}$ is sensitive to the speckle
noise but does not affect the planet signal is shown in the bottom right box
of Figure 6. In the de-rotated residual sequence $\mathbf{r}_{t}$, the signal
of the planet $\mathbf{p}_{t}$ is always at the same position and has
approximately constant brightness666We do not account for changing conditions
during the night. However, if information about the variability of the planet
signal as a function of time is available (e.g., for data taken with a vAPP
coronagraph such as in Otten et al. 2017), one could normalize for these
variations as part of the pre-processing.. Consequently, it only contributes
to the temporal mean $\overline{\mathbf{r}}$ but does not influence the
temporal variance. Since the temporal mean $\overline{\mathbf{r}}$ is
subtracted in Equation 12, the loss function becomes invariant to the presence
of the planet signal. Quasi-static and atmospheric speckles, however, vary
with time, which is why they affect the loss function (red dashed line in
Figure 6). Note that the noise model of 4S is still applied to the science
images before de-rotation (cf. Equation 10). Only the loss function is
evaluated on the de-rotated residual sequence.

The signature of the noise learned by 4S is encoded in the projection matrix
$\mathbf{B}\in\mathbb{R}^{D\times D}$, which contains the parameters $\beta$
of our noise model. Finding good values for $\beta$ that minimize
$\mathcal{L}_{\text{4S}}$ is difficult in practice, since the image rotation
$\phi_{t}$ must be taken into account during optimization. We tackle this
problem by using the Python package PyTorch for gradient-based optimization
(Paszke et al., 2019), leveraging its automatic gradient computation (Baydin
et al., 2018). Gradients are calculated through sampler grids, initially
developed for spatial transformer networks (Jaderberg et al., 2016), allowing
effective back-propagation through the image rotations. We chose the L-BFGS
optimizer (Liu & Nocedal, 1989), which converged faster for our application
than other optimizers such as SGD with momentum (Sutskever et al., 2013) or
Adam (Kingma & Ba, 2017). We also find that the use of mini-batches, as often
used in neural networks, leads to significantly worse results. Therefore, the
whole data set is used during each gradient step. This is likely due to noisy
values of $\overline{\mathbf{r}}$, which contains the signal of the planet for
each mini-batch. For large datasets that do not fit into GPU memory, gradient
accumulation can be used.

###  4.3 Regularization

So far, the number of free model parameters of 4S is huge: Given a frame
resolution of $100\times 100$ pixels, $\mathbf{B}$ contains 100 million
parameters. To mitigate the risk of overfitting, we use two types of
regularization to reduce the expressiveness of the model: First, a classical
L2 penalty term on the model parameters (similar to ridge regression; see
chapter 3.1.4 of Bishop 2006), which we add to the loss function:

$\mathcal{L}_{\text{4S}}=\underbrace{\frac{1}{N}\sum_{t=1}^{N}\left\|\mathbf{r}_{t}-\overline{\mathbf{r}}\right\|^{2}}_{\text{temporal
variance}}+\underbrace{\lambda\cdot\left\|\beta\right\|^{2}\vphantom{\sum_{t=1}^{N}}}_{\text{L2
penalty}}\,,$ (13)

where $\lambda$ is the only hyperparameter of 4S. Larger values of $\lambda$
force the parameters $\beta$ to be closer to zero, reducing the expressiveness
of the model. This means that a smaller value of $\lambda$ results in a
stronger reduction of the speckle noise while coming at the risk of
overfitting (i.e., memorizing the data). In contrast to PCA, we notice that
the best choice of $\lambda$ does not depend on the separation from the star.
Instead, it depends on the temporal and spatial resolution of the science
frames. If only a few frames are available in time, a stronger regularization
(i.e., larger values of $\lambda$) is needed and the reduction of speckle
noise is less effective. However, a detailed analysis of how to tune $\lambda$
is beyond the scope of this paper.

In addition to the L2 regularization, we also convolve the column vectors of
$\mathbf{B}$ with the unsaturated PSF:

$\mathbf{B}_{[:,n_{i}j+i]}=\left(\beta_{(i,j)}\circ\mathbf{m}_{(i,j)}\right)\circledast\text{PSF}\,.$
(14)

The convolution ties spatially adjacent parameters together, reducing the
effective number of model parameters. From our domain knowledge, we know that
speckles approximately follow the shape of the PSF. All features in the data
that are smaller are not due to the speckle noise, but are caused, for
example, by detector noise. By convolving the parameters $\beta$ with the PSF,
we prevent our model from overfitting pixel level artifacts.

###  4.4 Main steps of 4S

Applying 4S to ADI data in HCI consists of the following steps, which we
illustrate in Figure 6:

  1. 1

Data preparation: Similar to PCA, 4S works on normalized data. We first
subtract the temporal mean $\overline{\mathbf{x}}$ and divide by the standard
deviation along the time axis:

$\mathbf{x}_{t}=\left(\mathbf{x}_{t,\text{raw}}-\overline{\mathbf{x}}\right)/\text{std}(\mathbf{x}_{t,\text{raw}})$
(15)

The additional division by $\text{std}(\mathbf{x}_{t,\text{raw}})$ is required
to ensure that the L2 regularization affects all areas in the image equally.
The 2D science frames are flattened into one-dimensional vectors and stacked
into the data matrix $\mathbf{X}\in\mathbb{R}^{N\times D}$ (cf. the PCA
procedure in Section 2).

  2. 2

Optimization of model parameters: We initialize all parameters $\beta$ with
zeros and create the right reason mask $\mathbf{m}_{(i,j)}$ for all output
pixel positions. Afterward, we optimize $\beta$ to minimize the loss function
in Equation 13 using an L-BFGS optimizer.

  3. 3

Calculation of the residual: After the optimization is completed, we calculate
the residual sequence using Equation 4. As for PCA, the final residual is
obtained by averaging along the time dimension (Equation 5).

We note that the subtraction of the temporal average in step 1 is not free
from planet self-subtraction. More advanced ideas to overcome this step should
be explored in future work. Our Python implementation takes advantage of
modern NVIDIA A100 and H100 GPUs to speed up the optimization. Depending on
the spatial and temporal resolution of the dataset, the optimization converges
in about 10 to 30 minutes using one GPU.

Table 1:  VLT-NACO datasets used for performance evaluation.  # | Target | Program ID | Observing date |  $r_{0}$ () | Field rotation () | DIT () | total () | AGPM
---|---|---|---|---|---|---|---|---
1 | HD 2262 | 199.C-0065 (C)${}^{\textit{a}}$ | 2017-11-01 | $0.56$ | $73.7$ | $0.35$ | $5043$ | ✓
2 | HD 7570 | 1101.C-0092 (C)${}^{\textit{c}}$ | 2018-11-27 | $0.65$ | $68.9$ | $0.35$ | $6895$ | ✓
3 | HD 11171 | 1101.C-0092 (C)${}^{\textit{c}}$ | 2018-11-28 | $0.63$ | $108.7$ | $0.35$ | $7350$ | ✓
4 | HD 22049 | 096.C-0679 (A)${}^{\textit{c}}$ | 2015-12-17 | $0.57$ | $90.1$ | $0.08$ | $4608$ | ✓
5 | HD 22049 | 199.C-0065 (C)${}^{\textit{a}}$ | 2017-10-30 | $0.61$ | $104.7$ | $0.08$ | $5304$ | ✓
6 | HD 38678 | 084.C-0396 (A)${}^{\textit{b}}$ | 2009-11-27 | $1.47$ | $39.2$ | $0.2$ | $2460$ |
7 | HD 40136 | 1101.C-0092(C)${}^{\textit{c}}$ | 2018-11-29 | $0.87$ | $114.7$ | $0.3$ | $7500$ | ✓
8 | HD 115892 | 1101.C-0092 (E)${}^{\textit{c}}$ | 2019-05-23 | $0.48$ | $95.2$ | $0.35$ | $6895$ | ✓
9 | HD 169022 | 1101.C-0092 (E)${}^{\textit{c}}$ | 2019-05-20 | $1.28$ | $113.8$ | $0.35$ | $6195$ | ✓
10 | HD 177724 | 091.C-0654 (A)${}^{\textit{c}}$ | 2013-08-19 | $0.65$ | $24.9$ | $0.2$ | $3200$ |
11 | HD 209952 | 089.C-0149 (A)${}^{\textit{b}}$ | 2012-07-14 | $1.12$ | $40.4$ | $0.2$ | $2880$ |
12 | HD 35850 | 088.C-0085(A)${}^{\textit{b}}$ | 2011-10-21 | $1.1$ | $69.6$ | $0.2$ | $3120$ |

##  5 Quantitative analysis

Comparing the contrast performance of post-processing techniques can be
difficult for two main reasons: First, different post-processing methods
follow fundamentally different concepts. PCA and 4S are both subtraction-based
techniques. They return a residual image, which is analyzed in a second step
with statistical methods like t-tests or STIM maps (Mawet et al., 2014; Pairet
et al., 2019). Methods like FMMF, PACO or TRAP (Ruffio et al., 2017; Flasseur
et al., 2018; Samland et al., 2021) belong to the family of inverse problem
approaches. They produce detection maps that show the models’ belief in the
presence of a planet. The detection process for subtraction-based and inverse-
problem approaches rely on very different statistics. Therefore, their results
are not directly comparable. In addition, the residual noise characteristics
are method-specific and usually unknown. This makes a fair comparison even
more difficult. Especially if we want to compare different methods given the
same False Positive Fraction (FPF). For this reason, methods should only be
compared within the same category (see, e.g., Cantalloube et al., 2021).

Secondly, the contrast performance is influenced by the instrument and the
observing mode. Depending on the wavelength, different noise sources dominate,
such that the best choice of post-processing method is likely to be
instrument-dependent. To keep our analysis consistent, we decided to focus on
archival data from the VLT-NACO instrument in the $L^{\prime}$-band
($\lambda=\qty{3.8}{\micro}$). The exposure times of these data are typically
very short, which gives us a high temporal resolution. This is crucial for 4S
to minimize the risk of overfitting. It also allows us to correct not only for
quasi-static speckles, but also for atmospheric speckles on shorter time
scales. To ensure the robustness of our results under changing observing
conditions and targets, we base our analysis on 11 archival datasets (see
Table 1). This guarantees that any observed improvements are not random, but
consistent across multiple datasets. All targets are bright
($L^{\prime}<\qty{4}{}$) and were chosen such that no obvious companions or
disc signals are present. We ensure that each dataset has at least 30 minutes
of observations in seeing conditions no worse than $r_{0}<1.5$. The final list
is a random selection of datasets in the archive that satisfy these
conditions.

We process all our datasets with the state-of-the-art pipeline PynPoint. The
preprocessing steps include: 1\. a simple dark and flat calibration, 2\. an
interpolation of bad pixels, 3\. a simple background subtraction using sky
frames, and 4\. a centering of the star using cross-correlations. To keep the
computation time manageable, we combine every 5 frames of the science sequence
and crop the frames to a field of view of radius 1.2 (equivalent to 12). Even
after temporal binning, we still have between $2400$ and $13\,000$ frames
along the time dimension. The final pre-processed data cubes are publicly
available on Zenodo (Bonse, 2024a).

###  5.1 Right for the right reasons

Figure 7:  Saliency maps for 4S, computed for the datasets #4 (HD 22049) (left
column) and #7 (HD 40136) (right column). The saliency maps show which
information from the science frames $\mathbf{x}_{t}$ is used by 4S to estimate
the noise at the position marked by the white cross. In contrast to PCA (cf.
Figure 3), 4S successfully learns the point-symmetric signatures of the
speckle noise. The red star marks the center of the frame.

In Section 3.2, we used saliency maps to better understand which information
is used by PCA to estimate the noise at a given position. This approach
revealed that PCA converges to the telescope’s PSF for high numbers of
components, causing significant signal loss. 4S is designed to overcome this
limitation and forces the model to explore other explanations for the speckle
noise using a right reason constraint. But is the noise estimate of 4S right
for the right reasons?

Similar to PCA, we can compute saliency maps for 4S which are just the column
vectors of our noise model:

$M_{(i,j)}=\left|\mathbf{B}\right|_{[:,n_{i}j+i]}.$ (16)

We run 4S on the datasets #4 (HD 22049) and #7 (HD 40136) with a
regularization strength $\lambda=10^{3}$ and compute saliency maps for three
different positions. The two datasets were selected randomly, with about three
years between observations. The output saliency maps are shown in Figure 7. We
notice that 4S no longer learns the shape of the telescope’s PSF, but instead
explores other areas in the science frames. The most dominant feature in the
saliency maps of 4S are strong point symmetries. This is a strong indication
that 4S indeed learns the characteristic behavior of speckles, which are known
to exhibit symmetric or antisymmetric correlations (Bloemhof et al., 2001;
Perrin et al., 2003; Bloemhof, 2007; Ribak & Gladysz, 2008). Another pattern,
characterized by four points of equal distance from the star, can be seen in
the lower left saliency map of Figure 7. This pattern is due to waffle mode, a
distortion of the wavefront that is not visible to the wavefront sensor. It is
an artifact in the data that we want to remove in the data post-processing.
Due to the width of the $L^{\prime}$-band filter, the speckles are slightly
elongated. 4S finds this characteristic of the speckle noise and uses
information on a line pointing towards the center of the frame (bottom two
frames of Figure 7).

The two datasets, #4 (HD 22049) and #7 (HD 40136), were acquired with the same
observing mode. However, there is a significant gap of three years between the
two observations. Despite this time difference, it is remarkable to observe
the similarity in the saliency maps produced by 4S. This suggests that 4S can
learn signatures of the noise that are linked to the instrument and not the
individual observation. This result opens up several possibilities for future
work, which we discuss in Section 7.

###  5.2 Fake planet experiments

Figure 8:  Results of our fake planet experiments for dataset #5 (HD 22049)
(see Table 1). On the left, example residuals are shown for fake planets at
$\blacktriangle=(\qty{1.5}{},\qty{7.5}{})$,
$\blacksquare=(\qty{2.5}{},\qty{9.5}{})$ and
$\bigstar=(\qty{11}{},\qty{13}{})$. In the case of PCA, we tune the number of
principal components to get the deepest contrast (_optimal PCA_). For 4S, we
use the same setup with $\lambda=100$ for all separations. An interactive
version of the plot is available on our documentation page. The plot on the
right shows contrast curves of 4S and PCA with different numbers of principal
components. Each contrast curve was obtained by computing a contrast grid
using the Python package applefy and is based on 945 fake planet experiments.
The dotted black line gives the cADI contrast, which is equivalent to a simple
subtraction of the temporal median. Note that for large numbers of components
and small separations, none of the inserted fake planets reaches a
$\text{S/N}>5$.

To reliably calculate the contrast performance of 4S, we perform extensive
fake planet experiments. We use contrast grids implemented in the python
package applefy (Bonse et al., 2023). Contrast grids systematically calculate
the S/N of inserted fake companions as a function of distance from the star
and planet brightness. We insert fake planets every 0.5 mag for contrasts
$[5,15]$ mag, and for separations $[1.5,12.5]$ $\lambda/D$ every
$0.5\lambda/D$. For each combination of planet brightness and separation,
three planets are inserted at different position angles. We insert only one
fake planet at a time, resulting in 945 fake planet experiments per dataset.
Each dataset with an inserted fake planet is processed with PCA and 4S to
obtain a residual image. Afterward, the S/N of the companion is calculated
using the method proposed by Mawet et al. (2014). Since our limits are in the
speckle-dominated regime, we use spaced pixels instead of apertures to extract
the photometry, as suggested by Bonse et al. (2023). The final contrast curves
are obtained by thresholding the contrast grid using all 945 fake planet
experiments.

The results of PCA depend strongly on the number of components used.
Therefore, 33 different setups with $K\in[1,500]$ were calculated for each
fake planet experiment. For 4S, the regularization strength was set to
$\lambda\in[10^{2},10^{3},10^{4},10^{5}]$. We compute contrast grids for each
dataset (see Table 1) and all algorithm setups. In total, we compute more than
$300\,000$ PCA and more than $30\,000$ 4S residuals. This experiment required
substantial compute and was thus realized on a compute cluster, allowing us to
use up to 200 NVIDIA A100 / H100 GPUs simultaneously. To further speed up the
computations, we reimplemented the PynPoint version of PCA using PyTorch to
enable GPU computations. Our new GPU implementation of PCA uses the fast SVD
from Halko et al. (2011), resulting in a 50-fold speedup on one NVIDIA H100
GPU compared to 64 CPUs. The code for our new implementation is available
along with our Python package fours. Example residuals and contrast curves for
one of our datasets #5 (HD 22049) are shown in Figure 8. As shown in the plot,
4S reaches significantly deeper contrast compared to PCA, especially at
separations close to the star. The two close fake companions at $1.5\lambda/D$
and $2.5\lambda/D$ are undetectable with PCA, regardless of the number of
principal components. For separations $<$ 2.5, we reach about 1.4 deeper
contrast with 4S compared to PCA. Even at larger separations, a decent
improvement of about $\sim 0.3$ mag can be achieved. Although 4S is designed
to prevent planet self-subtraction, we still observe negative _wings_ next to
the planet’s signal at close separations. These negative wings are caused by
the subtraction of the temporal average in the first step of 4S. If we were
able to subtract the temporal mean without the planet signal, the negative
wings would disappear and even better results could be achieved. This path
should be explored in future work within the context of RDI.

To achieve the deepest contrast of PCA, we need to carefully tune the number
of components for each separation. While 16 components give the best results
at $\sim$ 3, 400 components are required at $\sim$ 12\. For 4S, on the other
hand, $\lambda=100$ consistently outperforms all PCA setups. Since tuning the
number of principal components is difficult in practice and could suffer from
human bias, this is a clear advantage of 4S.

We repeat the same experiment shown in Figure 8 for all datasets summarized in
Table 1. Depending on the dataset, we select the best
$\lambda\in[10^{2},10^{3},10^{4},10^{5}]$ for 4S and compare it to the optimal
PCA limits. For datasets #10 and #11 we had to exclude the result for the
innermost 2.5 due to strong over-saturation from the host star. The relative
improvement of 4S over PCA accumulated over all datasets is shown in Figure 9.

Figure 9:  Relative improvement in contrast of 4S over PCA, accumulated over
the 11 datasets shown in Table 1. The blue line shows the median improvement,
the gray boxes give the interquartile range (IQR), and the whiskers extend to
points that lie within 1.5 IQR. The points marked as circles are outliers.

The results are consistent with the improvement achieved on dataset #5, with
about 0.5 variation depending on the dataset. The improvement is largest at
close separations to the star. We know that the loss of the planet signal in
PCA is at its maximum at these separations (see Section 3). 4S successfully
mitigates this limitation, allowing it to learn a better representation of the
speckle noise.

The improvement of 4S over PCA is probably larger in practice, mainly for
three reasons:

  * •

If we go through the stack of PCA residuals for different numbers of
components, the speckles pattern changes considerably. For 4S the speckle
pattern in the residual images is very static for different values of
$\lambda$. Since we try 33 different setups for PCA, the chance of a speckle
co-aligning with the inserted fake companion increases. This artificially
increases the S/N in the PCA residuals.

  * •

It is likely that the optimal PCA limit is not reached in practice. This is
because it is difficult to find the best number of principal components. For
dataset #4 (HD 22049), Launhardt et al. (2020) report a contrast of 10.5 at
0.5 and 11.6 at 0.75. Our optimal PCA limits reach down to 11.5 at 0.5 and
12.3 at 0.75, which is almost one magnitude deeper. This difference is likely
due to the large number of setups tried in our PCA analysis. In addition,
better pre-processing of the data with PynPoint could explain part of this
improvement.

  * •

The noise distribution of the PCA residuals often deviates from Gaussian noise
(Pairet et al., 2019; Bonse et al., 2023). This is especially true for small
numbers of principal components. Heavy tailed non-Gaussian residual noise
increases the risk of false positives. In Appendix A, we study the residual
noise distribution of 4S and PCA. We find that the noise of 4S is not
perfectly Gaussian either, but it is more Gaussian than the residual noise of
PCA. This means that the chance of observing a $\text{S/N}>5$ that is caused
by noise is greater for PCA than for 4S. However, since the true noise
distribution is unknown for both methods, the limits cannot be calculated at
the same FPF level.

Over all datasets, the median contrast performance of 4S is consistently
better than limits reached with _optimal_ PCA. This improvement opens a new
discovery space on archival data and for future observations.

##  6 Scientific Demonstration

Figure 10:  Residual images of the AF Lep dataset (#12 in Table 1) taken with
VLT-NACO in 2011. The two images on the left show the residuals obtained with
PCA with $K=20$ and $K=60$ components, respectively. The image on the right
shows the residual of 4S, which reveals the giant exoplanet AF Lep b at
$\text{S/N}=6.8$, 11 years before its discovery. For better visibility, we
filter the residual images with a Gaussian kernel ($\sigma=\qty{0.8}{}$). The
S/N is calculated on the raw residuals before applying the Gaussian filter. An
interactive version of the plot for different setups of PCA and 4S is
available in our online documentation. A gallery of all PCA and 4S residuals
is given in Figure 15.

To demonstrate that 4S is not only capable of improving the contrast
performance for artificial companions, we have searched the ESO archive for
companions that have been found with VLT-SPHERE but not with VLT-NACO. VLT-
SPHERE has found several new companions, but unfortunately most of them have
never been observed with VLT-NACO.

AF Leporis is a young ($\sim$28\text{\,}\mathrm{Myr}$$ based on moving group
age estimate, Bell et al. 2015) F8V-type star (Gray et al., 2006) belonging to
the $\beta$ Pictoris moving group (Zúñiga-Fernández et al., 2021). Since AF
Lep presents a significant astrometric Hipparcos \- Gaia acceleration, it was
observed by multiple direct imaging surveys targeted at accelerating stars.
Three direct detections of the companion AF Lep b obtained with VLT-SPHERE and
Keck-NICR2 were claimed almost simultaneously in 2023 (De Rosa et al., 2023;
Mesa et al., 2023; Franson et al., 2023). While these observations firmly
place AF Lep b as a faint L-type dwarf in the planetary regime, estimations of
its dynamical mass via the astrometric acceleration are conflicting between
instruments and with age-luminosity model-based estimates. While the latter
disagreement is often observed in recent analyses (Cheetham et al., 2018;
Brandt et al., 2021a; Tobin et al., 2024), the former is likely due to the
very short time-baseline on which the orbit of AF Lep b could be sampled
(December 2021 to February 2023; see Figure 11).

Franson et al. (2023) discovered the planet in the $L^{\prime}$-band at a
separation of $\rho_{2023.090342}=\qty{342\pm 8}{}$ and at a contrast of
$\Delta L^{\prime}=\qty{9.94\pm 0.14}{}$. This is precisely the regime where
we expect an improvement with 4S over PCA. The star AF Lep was imaged with
VLT-NACO in 2011 at $L^{\prime}$-band using the normal imaging mode (see Table
1). The dataset was taken in medium observing conditions (seeing 1.1), but
provides $\qty{69.6}{}$ of field rotation. We process this archival dataset
with PynPoint (Stolker et al., 2019) using the same preprocessing routine
explained in Section 5, and then subtract the stellar PSF using PCA and 4S.
Our 4S reduction reveals the giant exoplanet AF Lep b on the opposite side of
the star, 11 years before its discovery in 2022 (see Figure 10). The signal of
the planet is also present in the residuals of PCA, but it is
indistinguishable from the speckle noise.

This precovery of AF Lep b presents a unique opportunity to extend the time-
baseline of its orbit sampling from about $1\text{\,}\mathrm{yr}$ to more than
$11\text{\,}\mathrm{yr}$. Furthermore, it enables an independent confirmation
of the existing $L^{\prime}$-band magnitude estimate (Franson et al., 2023).

#### Astrometry and Photometry

Because of the attenuation of the planet signal due to the post-processing, we
cannot extract the photometry and astrometry from the residual images
directly. In the case of PCA, it is common practice to insert a negative fake
planet into the raw data at the position of the detection. The photometry and
astrometry, as well as their uncertainties, are then calculated using MCMC
(see, e.g., Stolker et al., 2019). Each step in the MCMC requires a complete
recalculation of the PCA basis and the residual. A full MCMC often relies on
several thousand re-reductions, which becomes computationally expensive. For
4S, a reduction of the AF Lep dataset takes about 10 minutes, resulting in a
total computation time of about one month. It is likely possible to
significantly reduce the computation time required for an MCMC with 4S. This
could be achieved by avoiding the recalculation of the projection matrix
$\mathbf{B}$ for each iteration. Instead, $\mathbf{B}$ could be calculated
once and then updated continuously during the MCMC optimization. However, to
verify the accuracy of such an approach, further tests are needed that are
beyond the scope of this paper. Further details for future work in this
direction are discussed in Section 7. Instead, we decided to estimate the
astrometry and photometry based on the PCA residuals using the well-
established and tested MCMC code presented in Stolker et al. (2019). For this,
we first calculated PCA residuals for $K\in[1,200]$ every 10 components. The
peak $\text{S/N}=4.7$ is reached at 65 components. We fix $K=65$ and use the
MCMC procedure in PynPoint to estimate the astrometry and photometry. As in
Stolker et al. (2019), we first use a simplex minimization to find good
initialization parameters. Afterward, we run the MCMC for 500 steps using 100
walkers. Our results are summarized in Table 2. The estimated photometry of
$\Delta L^{\prime}=\qty{10.03(13:12)}{}$ is in perfect agreement (i.e., within
one standard deviation) with the results of Franson et al. (2023), who
reported $\Delta L^{\prime}=\qty{9.94\pm 0.14}{}$.

Table 2: Properties of AF Lep b. $\rho_{2011.800}$ () | ${323.24}_{-6.44}^{+6.71}$
---|---
$\theta_{2011.800}$ () | ${258.81}_{-0.59}^{+0.53}$
$\Delta L^{\prime}$ () | ${10.03}_{-0.12}^{+0.13}$
$L^{\prime}$ () | ${14.96}_{-0.13}^{+0.14}$
$M_{L^{\prime}}$ () | ${12.81}_{-0.13}^{+0.14}$
Primary mass () | ${1.219}_{-0.054}^{+0.052}$
Secondary mass () | ${3.74}_{-0.50}^{+0.53}$
Semi-major axis () | ${9.01}_{-0.19}^{+0.20}$
Inclination () | ${55.8}_{-7.2}^{+6.2}$
Ascending node () | ${69.7}_{-5.4}^{+5.6}$
Mean longitude at 2010.0 () | ${170}_{-13}^{+11}$
Eccentricity | ${0.031}_{-0.020}^{+0.027}$
Argument of periastron () | ${-10}_{-52}^{+66}$

#### Orbit Fit

Figure 11:  Results from the updated orbit fit of AF Lep b. The corner diagram
reports the posterior distribution of an MCMC-based orbit fit using orvara.
The median with the 16 and 84 quantiles is reported. The top-right inset shows
all available relative astrometry measurements with the best orbital fit. To
illustrate the confidence of the fit, 100 orbits are drawn from the fit
posterior (gray). For a summary of all results, see Table 2.

We combine the new relative astrometry data point obtained in this work with
the relative astrometry points from the three AF Lep b detection papers (De
Rosa et al., 2023; Mesa et al., 2023; Franson et al., 2023), Hipparcos (see
van Leeuwen, 2007) and Gaia (Prusti et al., 2016; Vallenari et al., 2023)
acceleration as well as archival RV data published by Butler et al. (2017). To
combine the data into a single orbit fit, we use the MCMC-based orbit fitting
package orvara (Brandt et al., 2021b). We fit for all orbital parameters and
the primary and secondary mass, setting the default priors as defined in Table
4 of Brandt et al. (2021b). To stay comparable to the results in Franson et
al. (2023), we also set a highly constrained Gaussian prior on the mass of AF
Lep: $M=\qty{1.2\pm 0.06}{}$ (Kervella et al., 2022). We run the Monte Carlo
chain to $500\,000$ samples over $100$ walkers, discarding the first 20 of
each chain as burn-in. The parameter distributions and resulting fit values
are reported in Figure 11 and Table 2. We report a dynamical mass estimate of
$M=\qty{3.74(53:50)}{}$, confirming its status well within the planetary-mass
regime. This estimate offers a significant improvement over and is compatible
with the masses reported in Franson et al. (2023) and De Rosa et al. (2023).

As the additional VLT-NACO data point is almost precisely opposite to the more
recent observations, it offers an ideal opportunity to constrain the orbit’s
eccentricity. Indeed, we find AF Lep b to move on a circular orbit
($\epsilon=\qty{0.031(27:20)}{}$). Considering that due to the low mass of the
companion it likely formed in-situ in the proto-planetary disk, our
measurement agrees much better with tentative evidence of planet-mass
companions occurring preferentially on circular orbit (Bowler et al., 2020).
We note, that due to the relative astrometry point being opposite, we do not
improve the fit of the inclination. This circular orbit agrees with recent
interferometric observations of the planet with VLTI-GRAVITY (Balmer et al.,
subm.).

In addition to the 2011 observation, AF Lep was observed a second time with
VLT-NACO in 2015 (prog. ID 096.C-0679(A); Launhardt et al. 2020). Our re-
analysis of this dataset did not yield a detection. In 2015, the projected
separation of the planet from the star is much smaller ($\lesssim$ 2) than in
2011, resulting in higher contrast requirements. Therefore, the non-detection
in 2015 is consistent with our contrast limits and orbit fit.

##  7 Future Work

The detection of AF Lep b in archival VLT-NACO data demonstrates the potential
of 4S to find previously missed companions. The archive of VLT-NACO
observations comprises more than 800 datasets from more than 450 stars
(considering only $L^{\prime}$-band data taken in pupil tracking mode). It is
possible that a re-analysis of these datasets will reveal additional
companions. Future work should consider a re-analysis not only for VLT-NACO
data, but also for other instruments. In this way, a better understanding of
the contrast performance of 4S for different observing modes and instruments
can be obtained. The use of saliency maps together with 4S may help us to
better understand the optics and speckle pattern. If some difficult-to-
interpret patterns are systematically visible in the saliency maps of an
instrument, they may indicate overlooked effects in the instrument optics.
This could open the door to new developments on the hardware side.

The primary drawback of 4S is its computational cost. Processing one dataset
takes about 10–30 minutes on a single NVIDIA A100/H100 GPU. While this
computation time is unproblematic for detection purposes, it quickly becomes
infeasible for scientific analysis. Extracting astrometry and photometry, or
calculating detection limits, requires re-processing the same dataset hundreds
of times. In fact, the most computationally expensive step in 4S is the
optimization of the projection matrix in step 2 (see Figure 6). We recompute
the projection matrix from scratch for all our fake planet experiments. This
might not be necessary. Future work should explore strategies to compute the
matrix once and update it as the scientific analysis progresses. Fast GPU-
based MCMC implementations compatible with our code are available, for
example, in pyro (Bingham et al., 2019). The L-BFGS optimizer used in this
work requires a lot of GPU memory, especially for large science frames. This
sets the maximum resolution which can still be processed with one NVIDIA
A100/H100 GPU to about $150\times 150$ pixels, depending on the temporal
resolution of the dataset. We recommend cropping the science frames to benefit
from the better contrast performance of 4S at close separations.
Unfortunately, this memory requirement currently prevents us from processing
the datasets of the exoplanet imaging data challenge (Cantalloube et al.,
2015). Future work should focus on making our implementation more memory
efficient, for example, by using gradient accumulation during optimization.

#### Further development of the method

The noise model of 4S is a simple linear model. Especially at separations
close to the star and in data taken with a coronagraph, non-linear artifacts
can occur. For this purpose, the linear model of 4S in Equation 10 could be
replaced by a neural network. Neural networks have already been proposed for
speckle subtraction in HCI (see Wolf et al. 2023). In Section 5, we observed
that a substantial part of the information used by 4S to explain the speckle
noise are symmetries across the entire frame. Therefore, neural networks
should not be used locally as in Wolf et al. (2023), but globally. A possible
direction for this could be masked auto-encoders (He et al., 2022).

A number of new post-processing methods based on PCA have been proposed in the
literature. One example is the RSM map of Dahlqvist et al. (2020), which
replaces the last step in PCA, the averaging along time and allows combining
multiple post-processing algorithms. Another example is the combination of PCA
with forward modeling (FMMF; Ruffio et al. 2017). The basic steps of 4S are
identical to those of PCA, which is why methods such as RSM and FMMF can also
be used with 4S. For this reason, 4S should be seen as a building block for
new post-processing methods.

Our analysis of the saliency maps in Section 5 has shown that the noise model
learned by 4S hardly changes over years. This is a strong indicator that 4S
could also be used with RDI. For RDI, the noise model of 4S needs to be
learned only once and can then be applied to many datasets. This procedure
also avoids the high computational costs of 4S for ADI. The parameters of 4S
could be conditioned on observation conditions, as suggested by Gebhard et al.
(2022) for half-sibling regression. In combination with RDI, we might be able
to obtain a universal model of the speckle noise that utilizes all previously
collected data.

##  8 Summary & Conclusions

In this paper, we have gained a deep understanding for the loss of planet
signal in PCA-based post-processing. Using saliency maps, we identified two
main reasons for this loss: first, the convergence of the PCA noise model to
the telescope PSF; and second, a problematic loss function that misleads the
noise estimate to fit the planet signal. Building upon these insights, we
developed a new post-processing algorithm that overcomes the limitations of
PCA. We named this method 4S (Signal Safe Speckle Subtraction). 4S is
characterized by three innovations: First, a linear noise model with a right
reason constraint, second, a new loss function that is invariant to the planet
signal, and third, a domain knowledge-informed regularization.

We applied 4S to archival data from the VLT-NACO instrument to compare its
contrast performance to PCA. On a sample of 11 datasets, we observed a
contrast improvement of up to 1.5 magnitudes, especially at close separations
to the star. In addition to achieving deeper contrast, 4S produces residual
noise that is more Gaussian compared to PCA (see Appendix A), reducing the
risk of false positives. The main reason for this improvement is the ability
of 4S to learn a better representation of the speckle noise. A detailed
analysis revealed that 4S can recover speckle pattern that we would expect
from the theory. The increase in contrast performance allowed us to detect the
giant exoplanet AF Lep b in archival data from 2011, 11 years before its
discovery. Using this new astrometric point, we significantly improve the
constraints on the orbital parameters of the planet.

The introduction of 4S marks a leap for the post-processing of high-contrast
imaging (HCI) data. Looking ahead, our findings have the potential to boost
exoplanet detection on archival data and for future observations. We believe
that the findings from our work will provide inspiration for the development
of new post-processing methods beyond 4S.

##  9 Acknowledgments

This work was supported by an ETH Zurich Research Grant. M.J.B. and S.P.Q.
gratefully acknowledge the financial support from ETH Zurich. Parts of this
work has been carried out within the framework of the National Centre of
Competence in Research PlanetS supported by the Swiss National Science
Foundation (SNSF). J.H. acknowledges the financial support from the Swiss
National Science Foundation (SNSF) under project grant number 200020_200399.

Author contributions: M.J.B. carried out the main analyses, developed the 4S
algorithm and wrote the manuscript. He programmed the publicly available
python package fours and wrote the documentation page. T.D.G. helped with the
analysis and the development of 4S. His previous work on half-sibling
regression algorithm has set the basis for this paper. F.A.D. carried out the
orbit fits of AF Lep b and helped writing the corresponding section. S.P.Q.
and B.S. provided access to the compute resources needed for this project. All
authors discussed the algorithm, analysis and results. All authors commented
on the manuscript.

Used software: This work has made use of a number of open-source Python
packages, including astropy (The Astropy Collaboration et al., 2013, 2018,
2022), matplotlib (Hunter, 2007), numpy (Harris et al., 2020), pandas
(McKinney, 2010), photutils (Bradley et al., 2024), scikit-image (van der Walt
et al., 2014), scikit-learn (Pedregosa et al., 2011), scipy (Virtanen et al.,
2020), seaborn (Waskom, 2021), torch (Paszke et al., 2019), and tqdm (da
Costa-Luis, 2019).

## References

  * Absil et al. (2013) Absil, O., Milli, J., Mawet, D., et al. 2013, Astronomy & Astrophysics, 559, L12, doi: 10.1051/0004-6361/201322748
  * Ali et al. (2023) Ali, S., Abuhmed, T., El-Sappagh, S., et al. 2023, Information Fusion, 99, 101805, doi: 10.1016/j.inffus.2023.101805
  * Amara & Quanz (2012) Amara, A., & Quanz, S. P. 2012, Monthly Notices of the Royal Astronomical Society, 427, 948, doi: 10.1111/j.1365-2966.2012.21918.x
  * Arcidiacono & Simoncini (2018) Arcidiacono, C., & Simoncini, V. 2018, in Adaptive Optics Systems VI, Vol. 10703 (International Society for Optics and Photonics), 1070331
  * Balmer et al. (subm.) Balmer, W., et al. subm.
  * Baydin et al. (2018) Baydin, A. G., Pearlmutter, B. A., Radul, A. A., & Siskind, J. M. 2018, Journal of Machine Learning Research, 18, 1
  * Bell et al. (2015) Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, Monthly Notices of the Royal Astronomical Society, 454, 593, doi: 10.1093/mnras/stv1981
  * Beuzit et al. (2019) Beuzit, J. L., Vigan, A., Mouillet, D., et al. 2019, Astronomy and Astrophysics, 631, A155, doi: 10.1051/0004-6361/201935251
  * Bingham et al. (2019) Bingham, E., Chen, J. P., Jankowiak, M., et al. 2019, Journal of Machine Learning Research, 20, 1
  * Bishop (2006) Bishop, C. M. 2006, Pattern Recognition and Machine Learning, Information Science and Statistics (New York: Springer)
  * Bloemhof (2007) Bloemhof, E. E. 2007, Optics Express, 15, 4705, doi: 10.1364/OE.15.004705
  * Bloemhof et al. (2001) Bloemhof, E. E., Dekany, R. G., Troy, M., & Oppenheimer, B. R. 2001, The Astrophysical Journal, 558, L71, doi: 10.1086/323494
  * Bonse (2024a) Bonse, M. J. 2024a, Raw data for: Use the 4S (Signal-Safe Speckle Subtraction), 0.1, Zenodo, doi: 10.5281/zenodo.11456704
  * Bonse (2024b) —. 2024b, Intermediate results for: Use the 4S (Signal-Safe Speckle Subtraction), 0.1, Zenodo, doi: 10.5281/zenodo.11457071
  * Bonse et al. (2018) Bonse, M. J., Quanz, S. P., & Amara, A. 2018, Wavelet Based Speckle Suppression for Exoplanet Imaging - Application of a de-Noising Technique in the Time Domain, doi: 10.48550/arXiv.1804.05063
  * Bonse et al. (2023) Bonse, M. J., Garvin, E. O., Gebhard, T. D., et al. 2023, The Astronomical Journal, 166, 71, doi: 10.3847/1538-3881/acc93c
  * Bowler et al. (2020) Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, The Astronomical Journal, 159, 63, doi: 10.3847/1538-3881/ab5b11
  * Bradley et al. (2024) Bradley, L., Sipőcz, B., Robitaille, T., et al. 2024, astropy/photutils: 1.12.0, Zenodo, doi: 10.5281/zenodo.10967176
  * Brandt et al. (2021a) Brandt, G. M., Dupuy, T. J., Li, Y., et al. 2021a, The Astronomical Journal, 162, 301, doi: 10.3847/1538-3881/ac273e
  * Brandt et al. (2021b) Brandt, T. D., Dupuy, T. J., Li, Y., et al. 2021b, The Astronomical Journal, 162, 186, doi: 10.3847/1538-3881/ac042e
  * Butler et al. (2017) Butler, R. P., Vogt, S. S., Laughlin, G., et al. 2017, The Astronomical Journal, 153, 208, doi: 10.3847/1538-3881/aa66ca
  * Cantalloube et al. (2015) Cantalloube, F., Mouillet, D., Mugnier, L. M., et al. 2015, Astronomy & Astrophysics, 582, A89, doi: 10.1051/0004-6361/201425571
  * Cantalloube et al. (2021) Cantalloube, F., Gomez-Gonzalez, C., Absil, O., et al. 2021, arXiv:2101.05080 [astro-ph], doi: 10.1117/12.2574803
  * Cantero et al. (2023) Cantero, C., Absil, O., Dahlqvist, C. H., & Van Droogenbroeck, M. 2023, Astronomy and Astrophysics, 680, A86, doi: 10.1051/0004-6361/202346085
  * Chauvin et al. (2017) Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, Astronomy and Astrophysics, 605, L9, doi: 10.1051/0004-6361/201731152
  * Cheetham et al. (2018) Cheetham, A., Ségransan, D., Peretti, S., et al. 2018, Astronomy & Astrophysics, 614, A16, doi: 10.1051/0004-6361/201630136
  * Christiaens et al. (2023) Christiaens, V., Gonzalez, C. A. G., Farkas, R., et al. 2023, Journal of Open Source Software, 8, 4774, doi: 10.21105/joss.04774
  * Cugno et al. (2023) Cugno, G., Pearce, T. D., Launhardt, R., et al. 2023, Astronomy & Astrophysics, 669, A145, doi: 10.1051/0004-6361/202244891
  * Currie et al. (2023) Currie, T., Biller, B., Lagrange, A.-M., et al. 2023, Direct Imaging and Spectroscopy of Extrasolar Planets, arXiv.  http://ascl.net/2205.05696
  * da Costa-Luis (2019) da Costa-Luis, C. O. 2019, JOSS, 4, 1277, doi: 10.21105/joss.01277
  * Daglayan et al. (2023) Daglayan, H., Vary, S., Leplat, V., Gillis, N., & Absil, P.-A. 2023, Direct Exoplanet Detection Using L1 Norm Low-Rank Approximation, arXiv.  http://ascl.net/2304.03619
  * Dahlqvist & Absil (2021) Dahlqvist, C.-H., & Absil, O. 2021, Astronomy & Astrophysics, 646, A49, doi: 10.1051/0004-6361/202039597
  * Dahlqvist et al. (2020) Dahlqvist, C.-H., Cantalloube, F., & Absil, O. 2020, Astronomy & Astrophysics, 633, A95, doi: 10.1051/0004-6361/201936421
  * De Rosa et al. (2023) De Rosa, R. J., Nielsen, E. L., Wahhaj, Z., et al. 2023, Astronomy and Astrophysics, 672, A94, doi: 10.1051/0004-6361/202345877
  * Flasseur et al. (2024) Flasseur, O., Bodrito, T., Mairal, J., et al. 2024, Monthly Notices of the Royal Astronomical Society, 527, 1534, doi: 10.1093/mnras/stad3143
  * Flasseur et al. (2018) Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, Astronomy & Astrophysics, 618, A138, doi: 10.1051/0004-6361/201832745
  * Franson et al. (2023) Franson, K., Bowler, B. P., Zhou, Y., et al. 2023, The Astrophysical Journal, 950, L19, doi: 10.3847/2041-8213/acd6f6
  * Fulton et al. (2021) Fulton, B. J., Rosenthal, L. J., Hirsch, L. A., et al. 2021, The Astrophysical Journal Supplement Series, 255, 14, doi: 10.3847/1538-4365/abfcc1
  * Gebhard et al. (2022) Gebhard, T. D., Bonse, M. J., Quanz, S. P., & Schölkopf, B. 2022, Astronomy & Astrophysics, 666, A9, doi: 10.1051/0004-6361/202142529
  * Gomez Gonzalez et al. (2016) Gomez Gonzalez, C. A., Absil, O., Absil, P.-A., et al. 2016, Astronomy & Astrophysics, 589, A54, doi: 10.1051/0004-6361/201527387
  * Gomez Gonzalez et al. (2018) Gomez Gonzalez, C. A., Absil, O., & Van Droogenbroeck, M. 2018, Astronomy & Astrophysics, 613, A71, doi: 10.1051/0004-6361/201731961
  * Gomez Gonzalez et al. (2017) Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, The Astronomical Journal, 154, 7, doi: 10.3847/1538-3881/aa73d7
  * Gray et al. (2006) Gray, R. O., Corbally, C. J., Garrison, R. F., et al. 2006, The Astronomical Journal, 132, 161, doi: 10.1086/504637
  * Halko et al. (2011) Halko, N., Martinsson, P. G., & Tropp, J. A. 2011, SIAM Review, 53, 217, doi: 10.1137/090771806
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * He et al. (2022) He, K., Chen, X., Xie, S., et al. 2022, in 2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) (New Orleans, LA, USA: IEEE), 15979–15988, doi: 10.1109/CVPR52688.2022.01553
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/mcse.2007.55
  * Jaderberg et al. (2016) Jaderberg, M., Simonyan, K., Zisserman, A., & Kavukcuoglu, K. 2016, Spatial Transformer Networks, arXiv, doi: 10.48550/arXiv.1506.02025
  * Janson et al. (2021) Janson, M., Squicciarini, V., Delorme, P., et al. 2021, Astronomy and Astrophysics, 646, A164, doi: 10.1051/0004-6361/202039683
  * Jensen-Clem et al. (2017) Jensen-Clem, R., Mawet, D., Gomez Gonzalez, C. A., et al. 2017, The Astronomical Journal, 155, 19, doi: 10.3847/1538-3881/aa97e4
  * Jovanovic et al. (2015) Jovanovic, N., Martinache, F., Guyon, O., et al. 2015, Publications of the Astronomical Society of the Pacific, 127, 890, doi: 10.1086/682989
  * Kenworthy et al. (2007) Kenworthy, M. A., Codona, J. L., Hinz, P. M., et al. 2007, The Astrophysical Journal, 660, 762, doi: 10.1086/513596
  * Kervella et al. (2022) Kervella, P., Arenou, F., & Thévenin, F. 2022, Astronomy & Astrophysics, 657, A7, doi: 10.1051/0004-6361/202142146
  * Kingma & Ba (2017) Kingma, D. P., & Ba, J. 2017, Adam: A Method for Stochastic Optimization, arXiv, doi: 10.48550/arXiv.1412.6980
  * Kuhn et al. (2001) Kuhn, J. R., Potter, D., & Parise, B. 2001, The Astrophysical Journal, 553, L189, doi: 10.1086/320686
  * Lafrenière et al. (2007) Lafrenière, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, The Astrophysical Journal, 660, 770, doi: 10.1086/513180
  * Lagrange et al. (2023) Lagrange, A. M., Philipot, F., Rubini, P., et al. 2023, Astronomy and Astrophysics, 677, A71, doi: 10.1051/0004-6361/202346165
  * Langlois et al. (2021) Langlois, M., Gratton, R., Lagrange, A.-M., et al. 2021, arXiv:2103.03976 [astro-ph].  http://ascl.net/2103.03976
  * Launhardt et al. (2020) Launhardt, R., Henning, Th., Quirrenbach, A., et al. 2020, Astronomy & Astrophysics, 635, A162, doi: 10.1051/0004-6361/201937000
  * Lewis et al. (2023) Lewis, B., Fitzgerald, M. P., Dodkins, R. H., Davis, K. K., & Lin, J. 2023, The Astronomical Journal, 165, 59, doi: 10.3847/1538-3881/aca761
  * Liu & Nocedal (1989) Liu, D. C., & Nocedal, J. 1989, Mathematical Programming, 45, 503, doi: 10.1007/BF01589116
  * Macintosh et al. (2014) Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, Proceedings of the National Academy of Science, 111, 12661, doi: 10.1073/pnas.1304215111
  * Macintosh et al. (2015) Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64, doi: 10.1126/science.aac5891
  * Males et al. (2021) Males, J. R., Fitzgerald, M. P., Belikov, R., & Guyon, O. 2021, arXiv:2107.04604 [astro-ph].  http://ascl.net/2107.04604
  * Males et al. (2018) Males, J. R., Close, L. M., Miller, K., et al. 2018, in Adaptive Optics Systems VI, Vol. 10703 (SPIE), 76–89, doi: 10.1117/12.2312992
  * Marois et al. (2014) Marois, C., Correia, C., Galicher, R., et al. 2014, in Adaptive Optics Systems IV, Vol. 9148 (International Society for Optics and Photonics), 91480U
  * Marois et al. (2006) Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, The Astrophysical Journal, 641, 556, doi: 10.1086/500401
  * Mawet et al. (2009) Mawet, D., Serabyn, E., Liewer, K., et al. 2009, The Astrophysical Journal, 709, 53, doi: 10.1088/0004-637X/709/1/53
  * Mawet et al. (2012) Mawet, D., Absil, O., Montagnier, G., et al. 2012, Astronomy & Astrophysics, 544, A131, doi: 10.1051/0004-6361/201219662
  * Mawet et al. (2013) Mawet, D., Absil, O., Delacroix, C., et al. 2013, Astronomy & Astrophysics, 552, L13, doi: 10.1051/0004-6361/201321315
  * Mawet et al. (2014) Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, The Astrophysical Journal, 792, 97, doi: 10.1088/0004-637X/792/2/97
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 56–61, doi: 10.25080/Majora-92bf1922-00a
  * Mesa et al. (2023) Mesa, D., Gratton, R., Kervella, P., et al. 2023, Astronomy and Astrophysics, 672, A93, doi: 10.1051/0004-6361/202345865
  * Meshkat et al. (2013) Meshkat, T., Kenworthy, M. A., Quanz, S. P., & Amara, A. 2013, The Astrophysical Journal, 780, 17, doi: 10.1088/0004-637X/780/1/17
  * N’Diaye et al. (2015) N’Diaye, M., Pueyo, L., & Soummer, R. 2015, The Astrophysical Journal, 799, 225, doi: 10.1088/0004-637X/799/2/225
  * Nielsen et al. (2019) Nielsen, E. L., De Rosa, R. J., Macintosh, B., et al. 2019, The Astronomical Journal, 158, 13, doi: 10.3847/1538-3881/ab16e9
  * Otten et al. (2017) Otten, G. P. P. L., Snik, F., Kenworthy, M. A., et al. 2017, The Astrophysical Journal, 834, 175, doi: 10.3847/1538-4357/834/2/175
  * Pairet et al. (2019) Pairet, B., Cantalloube, F., Gomez Gonzalez, C. A., Absil, O., & Jacques, L. 2019, Monthly Notices of the Royal Astronomical Society, 487, 2262, doi: 10.1093/mnras/stz1350
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., et al. 2019, in Advances in Neural Information Processing Systems, Vol. 32 (Curran Associates, Inc.)
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., et al. 2019, arXiv preprints
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Perrin et al. (2003) Perrin, M. D., Sivaramakrishnan, A., Makidon, R. B., Oppenheimer, B. R., & Graham, J. R. 2003, The Astrophysical Journal, 596, 702, doi: 10.1086/377689
  * Prusti et al. (2016) Prusti, T., de Bruijne, J. H. J., Brown, A. G. A., et al. 2016, Astronomy & Astrophysics, 595, A1, doi: 10.1051/0004-6361/201629272
  * Pueyo (2016) Pueyo, L. 2016, The Astrophysical Journal, 29
  * Quanz et al. (2011) Quanz, S. P., Schmid, H. M., Geissler, K., et al. 2011, The Astrophysical Journal, 738, 23, doi: 10.1088/0004-637X/738/1/23
  * Racine et al. (1999) Racine, R., Walker, G. A. H., Nadeau, D., Doyon, R., & Marois, C. 1999, Publications of the Astronomical Society of the Pacific, 111, 587, doi: 10.1086/316367
  * Rameau et al. (2013) Rameau, J., Chauvin, G., Lagrange, A.-M., et al. 2013, Astronomy & Astrophysics, 553, A60, doi: 10.1051/0004-6361/201220984
  * Ren et al. (2018) Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, The Astrophysical Journal, 852, 104, doi: 10.3847/1538-4357/aaa1f2
  * Ribak & Gladysz (2008) Ribak, E. N., & Gladysz, S. 2008, Optics Express, 16, 15553, doi: 10.1364/OE.16.015553
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., & Guestrin, C. 2016, in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (San Francisco California USA: ACM), 1135–1144, doi: 10.1145/2939672.2939778
  * Ross et al. (2017) Ross, A. S., Hughes, M. C., & Doshi-Velez, F. 2017, arXiv:1703.03717 [cs, stat].  http://ascl.net/1703.03717
  * Ruane et al. (2019) Ruane, G., Ngo, H., Mawet, D., et al. 2019, The Astronomical Journal, 157, 118, doi: 10.3847/1538-3881/aafee2
  * Ruffio et al. (2017) Ruffio, J.-B., Macintosh, B., Wang, J. J., et al. 2017, The Astrophysical Journal, 842, 14, doi: 10.3847/1538-4357/aa72dd
  * Samland et al. (2021) Samland, M., Bouwman, J., Hogg, D. W., et al. 2021, Astronomy & Astrophysics, 646, A24, doi: 10.1051/0004-6361/201937308
  * Schramowski et al. (2020) Schramowski, P., Stammer, W., Teso, S., et al. 2020, arXiv:2001.05371 [cs, stat].  http://ascl.net/2001.05371
  * Simonyan et al. (2013) Simonyan, K., Vedaldi, A., & Zisserman, A. 2013, arXiv preprint arXiv:1312.6034.  http://ascl.net/1312.6034
  * Smilkov et al. (2017) Smilkov, D., Thorat, N., Kim, B., Viégas, F., & Wattenberg, M. 2017, SmoothGrad: Removing Noise by Adding Noise, arXiv.  http://ascl.net/1706.03825
  * Snik et al. (2012) Snik, F., Otten, G., Kenworthy, M., et al. 2012, in Modern Technologies in Space- and Ground-based Telescopes and Instrumentation II, Vol. 8450 (SPIE), 224–234, doi: 10.1117/12.926222
  * Soummer et al. (2012) Soummer, R., Pueyo, L., & Larkin, J. 2012, The Astrophysical Journal, 755, L28, doi: 10.1088/2041-8205/755/2/L28
  * Sparks & Ford (2002) Sparks, W. B., & Ford, H. C. 2002, The Astrophysical Journal, 578, 543, doi: 10.1086/342401
  * Stolker et al. (2019) Stolker, T., Bonse, M. J., Quanz, S. P., et al. 2019, Astronomy & Astrophysics, 621, A59, doi: 10.1051/0004-6361/201834136
  * Sundararajan et al. (2017) Sundararajan, M., Taly, A., & Yan, Q. 2017, in Proceedings of the 34th International Conference on Machine Learning (PMLR), 3319–3328
  * Sutskever et al. (2013) Sutskever, I., Martens, J., Dahl, G., & Hinton, G. 2013, in Proceedings of the 30th International Conference on Machine Learning (PMLR), 1139–1147
  * The Astropy Collaboration et al. (2013) The Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * The Astropy Collaboration et al. (2022) The Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Thompson & Marois (2021) Thompson, W., & Marois, C. 2021, arXiv:2103.09252 [astro-ph].  http://ascl.net/2103.09252
  * Tobin et al. (2024) Tobin, T. L., Currie, T., Li, Y., et al. 2024, The Astronomical Journal, 167, 205, doi: 10.3847/1538-3881/ad3077
  * Vallenari et al. (2023) Vallenari, A., Brown, A. G. A., Prusti, T., et al. 2023, Astronomy & Astrophysics, 674, A1, doi: 10.1051/0004-6361/202243940
  * van der Walt et al. (2014) van der Walt, S., Schönberger, J. L., Nunez-Iglesias, J., et al. 2014, PeerJ, 2, e453, doi: 10.7717/peerj.453
  * van Leeuwen (2007) van Leeuwen, F. 2007, Astronomy & Astrophysics, 474, 653, doi: 10.1051/0004-6361:20078357
  * Vigan et al. (2021) Vigan, A., Fontanive, C., Meyer, M., et al. 2021, Astronomy & Astrophysics, 651, A72, doi: 10.1051/0004-6361/202038107
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wahhaj et al. (2015) Wahhaj, Z., Cieza, L. A., Mawet, D., et al. 2015, Astronomy & Astrophysics, 581, A24, doi: 10.1051/0004-6361/201525837
  * Wang et al. (2015) Wang, J. J., Ruffio, J.-B., De Rosa, R. J., et al. 2015, Astrophysics Source Code Library, ascl:1506.001
  * Waskom (2021) Waskom, M. 2021, JOSS, 6, 3021, doi: 10.21105/joss.03021
  * Wolf et al. (2023) Wolf, T. N., Jones, B. A., & Bowler, B. P. 2023, Direct Exoplanet Detection Using Deep Convolutional Image Reconstruction (ConStruct): A New Algorithm for Post-Processing High-Contrast Images, arXiv.  http://ascl.net/2312.03671
  * Zúñiga-Fernández et al. (2021) Zúñiga-Fernández, S., Bayo, A., Elliott, P., et al. 2021, Astronomy & Astrophysics, 645, A30, doi: 10.1051/0004-6361/202037830

Figure 12:  Two examples of the Gaussianity test performed in Figure 13. The
top row shows examples of the residuals obtained with PCA and 4S. The pixels
outside 2.54.5 were masked to extract the noise. Q-Q plots are shown below. If
the noise were perfectly Gaussian, we would expect all values to be on the
gray diagonal line and the value of $R^{2}$ to be close to 1.

##  Appendix A Residual noise distribution

A common assumption for the quantification of detections and non-detections in
HCI is that the residual noise after post-processing with PCA is Gaussian
distributed (see, e.g., Mawet et al., 2014). This assumption is generally
justified by the central limit theorem (we average many frames along time) and
the whitening effect of PCA. Nevertheless, the residual noise after post-
processing often deviates from Gaussian noise (Pairet et al., 2019; Bonse et
al., 2023). This is problematic because the risk of false positives (we claim
a detection that is just noise) depends on the type of noise. More
specifically, the probability of observing a value of $\text{S/N}=5$ that is
just noise is higher if the noise is heavy-tailed, that is, the frequency of
bright values such as speckles is higher than for Gaussian noise. Bonse et al.
(2023) found that if the residual noise follows a Laplacian distribution, the
final contrast could be overestimated by up to one magnitude if the test of
Mawet et al. (2014) is used.

In practice, the true noise distribution is usually unknown, making it
difficult to compare limits across datasets. This is also true if we want to
compare different post-processing techniques. We study the residuals of PCA
and 4S for the datasets #1–#9 (see Table 1) to get a better understanding of
their noise distributions. The datasets #10 and #11 are excluded due to strong
over-saturation at close separations. We consider two regions: 1. within
2.54.5, and 2. within 5.07.0 distance to the star. For each dataset and region
we, plot Q-Q plots and calculate the coefficient of determination $R^{2}$
(see, e.g., Pairet et al. 2019 for a detailed explanation). A comparison of
the residual noise for datasets #5 and #9 is shown in Figure 12. The Q-Q plots
shown in the figure compare the pixel noise distribution with Gaussian noise.
If the residual noise was perfectly Gaussian, we would expect all values to
lie on the gray diagonal line. The better the observed noise can be explained
by Gaussian noise, the closer the coefficient of determination $R^{2}$ is to
1. For PCA, the large values are above the diagonal line, indicating that the
noise is heavy-tailed. In the case of 4S, the residual noise is much closer to
Gaussian, although it is still not perfectly Gaussian. The quantification
based on the $R^{2}$ values confirms this observation. A summary of the
results for all datasets is shown in Figure 13. In general, the residual noise
at close separations tends to be more affected by speckle and is therefore
often more non-Gaussian. For most datasets, especially those for which PCA
deviates the most from Gaussian noise, 4S improves the Gaussianity of the
residuals. This means that 4S does not only improve the contrast limits as
shown in Section 5, but also reduces the risk of false positives.

Figure 13:  Comparison of the Gaussianity of the noise in the residuals of 4S
and PCA for the datasets #1–#9. If the coefficient of determination $R^{2}$ is
closer to 1, the noise is better explained by Gaussian noise. We choose the
number of principal components $K$ and the regularization strength $\lambda$
of 4S that maximizes the contrast (_optimal_ PCA).  Figure 14: PCA saliency
maps for different numbers of principal components $K$ and separations from
the star. The saliency maps show which information from the science frames is
used by PCA to estimate the noise at the position marked by the white circle.
Especially for large $K$ and separations close to the star, only the
information within the white circle is used, resulting in a substantial loss
of the planet signal. The red star marks the center of the frame.  Figure 15:
Residual images of the AF Lep dataset shown in Figure 10 for several different
settings of the algorithm’s hyperparameters (principal components $K$ for PCA
and regularization strength $\lambda$ for 4S). The planet AF Lep b is barely
visible in the PCA residuals. For 4S, any choice of parameters yields a clear
detection, highlighting the robustness of the method.