# An Alternating Minimization Algorithm with Trajectory for Direct Exoplanet Detection - The AMAT Algorithm

Hazan Daglayan, Simon Vary, Olivier Absil, Faustine Cantalloube, Valentin Christiaens, Nicolas Gillis, Laurent Jacques, Valentin Leplat, P. -A. Absil

## ► To cite this version:

Hazan Daglayan, Simon Vary, Olivier Absil, Faustine Cantalloube, Valentin Christiaens, et al., An Alternating Minimization Algorithm with Trajectory for Direct Exoplanet Detection - The AMAT Algorithm. Astronomy & Astrophysics - A&A, 2024, 10.1051/0004-6361/202451242. doi:10.1007/s00176-024-1051-8

HAL Id: hal-04749785

https://hal.science/hal-04749785v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An alternating minimization algorithm with trajectory for direct exoplanet detection

## The AMAT algorithm

H. Daglayan<fnref n="1" />,<fnref n="2" />,*, S. Vary<fnref n="3" />, O. Absil<fnref n="4" />,**, F. Cantalloube<fnref n="5" />, V. Christiaens<fnref n="4" />, N. Gillis<fnref n="6" />, L. Jacques<fnref n="1" />, V. Leplat<fnref n="7" />, and P.-A. Absil<fnref n="1" />

Received 25 June 2024 / Accepted 8 October 2024

## ABSTRACT

*Context.* Effective image post-processing algorithms are vital for the successful direct imaging of exoplanets. Standard point spread function (PSF) subtraction methods use techniques based on a low-rank approximation to separate the rotating planet signal from the quasi-static speckles and rely on signal-to-noise ratio maps to detect the planet. These steps do not interact or feed each other, leading to potential limitations in the accuracy and efficiency of exoplanet detection.

*Aims.* We aim to develop a novel approach that iteratively finds the flux of the planet and the low-rank approximation of quasi-static signals in an attempt to improve upon current PSF subtraction techniques.

*Methods.* In this study, we extend the standard L2 norm minimization paradigm to an L1 norm minimization framework in order to better account for noise statistics in the high contrast images. Then, we propose a new method, referred to as the alternating minimization algorithm with trajectory (AMAT), that makes more advanced use of estimating the low-rank approximation of the speckle field and the planet flux by alternating between them and utilizing both L1 and L2 norms. For the L1 norm minimization, we propose using L1 norm low-rank approximation (L1-LRA), a low-rank approximation computed using an exact block-cyclic coordinate descent method, while we use randomized singular value decomposition for the L2 norm minimization. Additionally, we enhance the visibility of the planet signal using a likelihood ratio as a post-processing step.

*Results.* Numerical experiments performed on a VLT/SPHERE-IRDIS dataset show the potential of AMAT to improve upon the existing approaches in terms of higher S/N, sensitivity limits (contrast curves), and receiver operating characteristic curves. Moreover, for a systematic comparison, we used datasets from the exoplanet data challenge to compare our algorithm with other algorithms in the challenge, and we find AMAT with a likelihood ratio map performs better than most algorithms tested on the exoplanet data challenge.

**Key words.** methods: data analysis – techniques: image processing – planets and satellites: detection

## 1. Introduction

High-contrast imaging (HCI) is an essential observing technique for exoplanet discovery, especially because of its ability to directly observe young, massive planets orbiting their host stars at large distances (Galicher & Mazoyer 2023). This direct observation method is complementary to indirect techniques and also provides crucial information such as bolometric luminosity, effective temperature, surface gravity, composition of the planets, and – assuming a formation and evolutionary model – an estimate of their mass (Bowler 2016). Progress in this field is supported by remarkable technological advances in extreme adaptive optics (AO) systems (Guyon 2018) and coronagraphy, which provide better atmospheric turbulence correction and raw stellar light suppression, respectively. These technologies equip

new generation HCI instruments, such as VLT/SPHERE (Beuzit et al. 2019), Gemini/GPI (Macintosh et al. 2014), and Subaru/SCExAO (Martinache et al. 2009), with the ability to achieve unprecedented raw contrasts. Despite these important steps, direct imaging of exoplanets remains a challenging task, and only 1% of known exoplanets have been discovered using this method (NASA Exoplanet Archive 2024). This limitation underscores the inherent difficulty in distinguishing these faint planets from the overwhelming brightness of their host stars. Residuals of star light in the form of speckles, which are caused by atmospheric turbulence and imperfections in telescopes and instruments, are very similar in shape and contrast to planets, posing a major challenge to direct imaging. To address this issue, angular differential imaging (ADI) has emerged as a common strategy (Marois et al. 2006). ADI involves capturing a sequence of frames in pupil-stabilized mode, wherein the telescope tracks the star's motion over time, keeping it centered in the image. This approach results in the star and the speckles appearing static

<footnote n="1">¹ ICTEAM Institute, UCLouvain, 1348 Louvain-la-Neuve, Belgium</footnote>

<footnote n="2">² Vlaamse Instelling voor Technologisch Onderzoek (VITO), 2400 Mol, Belgium</footnote>

<footnote n="3">³ Department of Statistics, University of Oxford, Oxford, UK</footnote>

<footnote n="4">⁴ STAR Institute, Université de Liège, 4000 Liège, Belgium</footnote>

<footnote n="5">⁵ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="6">⁶ Department of Mathematics and Operational Research, University of Mons, Mons, Belgium</footnote>

<footnote n="7">⁷ Innopolis University, Innopolis, Russia</footnote>

<footnote>* Corresponding author; hazan.daglayan@uclouvain.be</footnote>

<footnote>** F.R.S.-FNRS Senior Research Associate.</footnote>
or quasi-static, while planets exhibit movement as a function of the parallactic angle due to Earth's rotation. By exploiting this differential motion, ADI enables the isolation and detection of exoplanetary signals from the surrounding noise, increasing the chances of successful direct imaging.

Using ADI sequences, several post-processing methods have been proposed. Among the most common are those that build a model for the stellar PSF (including the static and quasi-static speckle field) and subtract it to detect planets. To build a PSF model, methods such as principal component analysis (PCA; Amara & Quanz 2012; Soummer et al. 2012) and non-negative matrix factorization (NMF, Gomez Gonzalez et al. 2017; Ren et al. 2018) aim to obtain a low-rank approximation of the time-by-pixel matrix containing the ADI observing sequence. After subtracting the PSF model from each frame, the residual matrix consists of both planetary signals and noise, which includes some residual stellar signal. Some studies suggest modeling the residual matrix as a sparse matrix and noise (LLSG and LRPT, Gomez Gonzalez et al. 2016; Vary et al. 2023). The locally optimized combination of images (LOCI) algorithm (Lafreniere et al. 2007) employs a least-squares approach to construct a model PSF and has variants such as TLOCI (Marois et al. 2010) and MLOCI (Wahhaj et al. 2015). Once the PSF model is obtained and the residual matrix is derived, one must still use a method to extract the planetary signal. The most popular way to do so is to build a signal-to-noise ratio map (Mawet et al. 2014), with an alternative being the use of a standardized trajectory intensity mean (STIM) map (Pairet et al. 2019). More recent methods to extract the planetary signal from residual data cubes include the regime-switching model (RSM, Dahlqvist et al. 2020), which attempts planetary detection using multiple PSF subtraction techniques at the same time, and the likelihood ratio map (LRM; Daglayan et al. 2022), which proposes a map consisting of likelihood ratios based on maximum likelihood estimation. Besides PSF model subtraction, other post-processing methods are based on inverse problem approaches to estimate the speckle field and planetary signal simultaneously in a maximum-likelihood approach, such as ANDROMEDA (Cantalloube et al. 2015), FMMF (Ruffio et al. 2017), PACO (Flasseur et al. 2018), and SNAP (Thompson & Marois 2021). Additionally, another type of algorithm utilizes unsupervised machine learning for planet detection, such as SODIRF and SODINN (Gonzalez et al. 2018) as well as their refined variant NA-SODINN (Cantero et al. 2023). As an aside, most of these post-processing methods are designed or optimized to detect point sources around the target star, and they generally struggle to reconstruct extended sources such as circumstellar disks. Iterative versions of PSF-subtraction methods, including iterative PCA, have recently been proposed to address this shortcoming (Pairet et al. 2021; Stapper & Ginski 2022; Juillard et al. 2023).

One recurrent assumption made by the post-processing methods described above is that residual noise after PSF subtraction is Gaussian. However, recent studies (Ruffio et al. 2017; Pairet et al. 2019; Dahlqvist et al. 2020; Daglayan et al. 2022; Cantero et al. 2023) have shown that the noise in the residual cube and/or processed frame, particularly in the tails of the distribution, tends to be non-Gaussian. In light of this, in Sect. 2, our paper proposes a low-rank approximation based on the Laplacian distribution, termed L1-LRA, that leverages the L1 norm. We demonstrate that this approach more accurately fits the data. On the other hand, the minimization of the L2 norm is computationally easier due to the smoothness of the objective function. Subsequently, we present an iterative method named alternating minimization with trajectory (AMAT) that is designed to

enhance algorithmic performance and more effectively differentiate between planetary signals and static or quasi-static signals using both L1 and L2 norms. Additionally, we establish that this iterative method, in its L2 norm version, slightly outperforms state-of-the-art methods for determining the planetary flux. To assess the performance of our proposed algorithms, we employ various benchmarks, such as S/N maps, contrast curves, and receiver operating characteristic (ROC) curves in Sect. 3. For these empirical analyses, we use an ADI cube obtained on 51 Eri with the VLT/SPHERE-IRDIS instrument, as used in the Samland et al. (2017) publication, containing 256 frames obtained in the K1 (2.11 µm) band over a parallactic angle range of 42° (Samland et al. 2017). Additionally, we leverage datasets from the exoplanet data challenge for comparative analysis against other state-of-the-art methods (Cantalloube et al. 2020). In Sect. 4, we propose to further improve the performance of our method by computing an LRM based on the residual cube processed with the AMAT algorithm, instead of using a standard S/N map. Section 5 concludes the paper.

Preliminary results related to the present work appeared in the proceedings of two machine learning conferences (Daglayan et al. 2023a,b). Daglayan et al. (2023b) suggested the potential use of the L1 norm low-rank approximation for exoplanet detection, while Daglayan et al. (2023a) provided a brief description of the AMAT algorithm and performed ablation studies to evaluate its performance. Both studies evaluated the algorithms using a single dataset. This paper expands on these descriptions, offering a comprehensive explanation of the algorithms and comparing them using various metrics. We demonstrate the performance of the algorithm across different datasets to illustrate its data-independent capabilities. Additionally, we explore the application of the AMAT algorithm for flux estimation and compare its efficacy with existing methods.

## 2. The AMAT algorithm

In this section, we present a novel method that employs an iterative technique for exoplanet detection that distinguishes planetary signals from the star and its associated speckles as well as from the sky background. This method aims to find a low-rank matrix that better fits the quasi-static signal and reveals the planetary signal more clearly and proposes the use of the L1 norm as a solution to mitigate the effects of Laplacian noise. Firstly, we describe the L1 low-rank approximation (L1-LRA) in detail, explaining its reliance on the L1 norm. We then present the AMAT algorithm, which is designed to accommodate both L1 and L2 norm scenarios.

### 2.1. L1 low-rank approximation

We let $M$ be a matrix in $\mathbb{R}^{t \times n^2}$ containing observations of $t$ unfolded frames, where each row represents a single vectorized frame of size $n \times n$. Assuming there is a single planet located at position $g \in [n] \times [n]$ within the first frame, with $[n] = \{1, 2, \dots, n\}$, the model for $M$ can be written as

$$M = L + a_g P_g + E, \quad \mathrm{rank}(L) \le k, \quad P_g \in \mathcal{P}, \quad (1)$$

where $L$ denotes the low-rank model for the stellar diffraction pattern; $E$ stands for the noise; $a_g$ is the intensity of the planet referred to as the flux; $P_g \in \mathcal{P} \subset \mathbb{R}^{t \times n^2}$ is the planet signature along the trajectory, illustrated in Fig. 1; and $\mathcal{P}$ is the set of all feasible planet signatures. To construct $P_g$, we started from a time-by-pixel matrix with zero entries, and in each unfolded
<figure><img src="image_2.png" /><figcaption>Fig. 1. Cube of the planet signature constructed by rotating the position g of the PSF function along the trajectory.</figcaption></figure>

frame (i.e., row of the matrix), we placed a copy of the nor-
malized reference PSF at the location occupied by a planet that
is at position g in the first frame. We normalized the reference
PSF using the method described in VIP (Gomez Gonzalez et al.
2017). This involves dividing the pixel values by the sum of
pixel intensities measured within a full-width half maximum
(FWHM) aperture. As a result, P<sub>g</sub> is a constant matrix for one
single planet position g.

In models employing low-rank approximations, selecting the
appropriate rank value is critical. If the rank is too small, speckle
signals may persist in the residual *M* − *L*, making it challenging
to distinguish the signal of the planet from the speckles. In con-
trast, the signal of the planet may be absorbed by the low-rank
matrix if the rank is too large, making it more difficult to locate
the planet signal in the residual.

When the error *E* follows a Gaussian distribution, the maxi-
mum likelihood estimator for *L* is obtained by minimizing the L2
norm. Classical methods such as PCA and LLSG fit the low-rank
component as follows:

$$
\hat{L} = \arg \min_{L} \|M - L\|_{2} \quad \text{subject to} \quad \text{rank}(L) \le k, \qquad (2)
$$

where ||A||₂ denotes the entry-wise L2 norm of A (the Frobenius
norm), which can be solved using the truncated singular value
decomposition (SVD), according to the Eckart–Young–Mirsky
theorem (Eckart & Young 1936). Such an approach to estimate L
has been used in Amara & Quanz (2012); Soummer et al. (2012);
Gomez Gonzalez et al. (2016). Recent findings indicate that the

error term *E* exhibits heavy tails, aligning more closely with the
Laplacian distribution (Pairet et al. 2019; Cantero et al. 2023).
Consequently, we propose fitting the low-rank component using
the component-wise L1 norm:

$$
\hat{L} = \arg \min_{L} \|M - L\|_{1} \quad \text{subject to} \quad \text{rank}(L) \le k. \qquad (3)
$$

This approach allows us to maintain a consistent noise assump-
tion across the low-rank speckle subtraction. Moreover, in PCA,
a common issue is the sensitivity to outliers when using the L2
norm. In contrast, the L1 norm demonstrates a robust approach to
outliers, leading to a better fit (Ke & Kanade 2003, 2005b; Song
et al. 2017). This makes the L1 norm a more suitable choice for
data with potential outliers.

The L1 low-rank approximation in Eq. (3), however, is an
NP-hard problem, even in the rank-one case (Gillis & Vavasis
2018). Hence, most algorithms to tackle Eq. (3), such as alter-
nating convex optimization (Ke & Kanade 2005a), the Wiberg
algorithm (Eriksson & Van Den Hengel 2010), and augmented
Lagrangian approaches (Zheng et al. 2012), do not guarantee
that a global optimal solution will be found, unlike in the case
of PCA. Moreover, the computed solutions are sensitive to the
initialization of the algorithms. We used Algorithm 1 (L1-LRA)
suggested by Gillis & Vavasis (2018) to solve Eq. (3). It solves the
problem using an exact block-cyclic coordinate descent method,
where the blocks of variables are the columns of Ŭ and Ŵ of
the low-rank approximation Ĭ = ŬW̄. This algorithm relies on
the fact that the columns of the matrix Ŵ ∈ ℝⁿ²×ᵏ form a basis
of a subspace of ℝⁿ² that best approximates the data frames in
the L1 sense, and for each data frame, the corresponding row of
Ū ∈ ℝ²×ᵏ contains the weights of the L1-best approximation of
the data frame in this basis. Here, Aⱼ denotes the j-th column of
the matrix A.

**Algorithm 1 L1-LRA (Gillis & Vavasis 2018)**

**Input:** Image sequence *M* ∈ ℝ<sup>*t*×*n*<sup>2</sup></sup>, the initial components
*Û* ∈ ℝ<sup>*t*×*k*</sup> and *V̂* ∈ ℝ<sup>*n*<sup>2</sup>×*k*</sup> of *M* (default initialization with the
randomized SVD), rank *k*, maximum number of iteration
*ℓ*<sub>max</sub>.

**Output:** The components Û and V̂.

1: **for** *ℓ* = 1 : *ℓ*<sub>max</sub> **do**

2: R = M - ŬV̂<sup>T</sup>

3: **for** *j* = 1 : *k* **do**

4: R ← R + Ū<sub>j</sub>V̂<sub>j</sub><sup>T</sup>

5: Ū<sub>j</sub> ← arg min<sub>u∈R<sup>t</sup></sub> ||R - uV̂<sub>j</sub><sup>T</sup>||<sub>1</sub>

6: V̂<sub>j</sub> ← arg min<sub>v∈R<sup>n<sup>2</sup></sup></sub> ||R<sup>T</sup> - vV̂<sub>j</sub><sup>T</sup>||<sub>1</sub>

7: R ← R - Ū<sub>j</sub>V̂<sub>j</sub><sup>T</sup>

8: end for

9: end for

To solve the minimization problem in steps 5–6 of Algo-
rithm 1, we used the exact method from Gillis & Plemmons
(2011); these sub-problems are weighted median problems that
can be solved in closed form. In our experiments, we applied an
annular version, similar to annular PCA (AnnPCA, Absil et al.
2013; Gomez Gonzalez et al. 2017), that selects only the pixels
of *M* in a certain annulus. Indeed, as the intensities of pixels
decrease away from the star, it is usually better to calculate the
low-rank approximation of each annulus separately.

In order to analyze the suitability of different noise assump-
tions, we fit Gaussian and Laplacian distributions to the residual
<figure><img src="image_3.png" /><figcaption>Fig. 2. Histograms of the residual cube after low-rank approximation is applied for small and large separations for a 51 Eri dataset.</figcaption></figure>

data, that is, the data after subtracting the low-rank component using PCA or L1-LRA. We looked at two different annuli separately, one that is close to the star at $4\lambda/D$ separation and one more distant from the star at $8\lambda/D$, and we measured the goodness of fit visually. In Fig. 2, we observed that the residual data follows somewhere between Gaussian and Laplacian on the peak and Laplacian on the tails of the distribution after applying PCA. However, after applying L1-LRA, the Laplace distribution provides a better fit for the residual cube distribution in general for both small and large separations. This supports that L1 is the indicated norm in a noise model where the error follows a Laplacian distribution (Gao et al. 2009).

## 2.2. The AMAT algorithm: Planet detection

Instead of traditional approaches of estimating first the low-rank component $L$ and then estimating the flux $a_g$, we propose estimating them simultaneously with the following optimization problem:

$$ \min_{L \in \mathbb{R}^{R \times n^2}, a_g \in \mathbb{R}} \|M - L - a_g P_g\|_{\#} \quad \text{s.t.} \quad \text{rank}(L) \le k, \quad P_g \in \mathcal{P}, \quad (4) $$

where $\| \cdot \|_{\#}$ denotes either the L1 or L2 norm depending on the assumed distribution of the error $E$ in Eq. (1). The optimization problem Eq. (4) is addressed by alternately solving the following two sub-problems until a stopping criterion is met:

$$ L^{(i)} = \arg \min_{L \in \mathbb{R}^{R \times n^2}} \|M - L - a_g^{(i-1)} P_g\|_{\#}, \quad (5a) $$

$$ a_g^{(i)} = \arg \min_{a \in \mathbb{R}} \|M - L^{(i)} - aP_g\|_{\#}. \quad (5b) $$

The stopping criterion is defined as either reaching a maximum number of iterations or ensuring that the relative changes in the intensity $a_g^{(i)}$ are less than a specified threshold.

When we select the norm L2 in (4) as given in Algorithm 2, computing $L^{(i)}$ in Eq. (5a) amounts to a k-truncated SVD, denoted by $H_k^{\text{SVD}}(\cdot)$. In practice, in order to speed up the computations, we compute $L^{(i)}$ by a randomized SVD of $M - a_g^{(i-1)} P_g$ (Halko et al. 2011). We checked that the resulting approximation does not affect the planet detection performance. The optimal value of $a_g^{(i)}$ can be computed by cross correlating the residual cube $R_g = M - L^{(i)}$ with the planet signature $P_g$:

$$ a_g^{(i)} = \frac{\sum_{(\theta,r) \in \Omega_g} R_g(\theta, r) P_g(\theta, r) \sigma_{R_g(r)}^{-2}}{\sum_{(\theta,r) \in \Omega_g} (P_g(\theta, r))^2 \sigma_{R_g(r)}^{-2}}, \quad (6) $$

where $\sigma_R^2$ is the empirical variance of the residual frames computed along the time dimension and we define the set $\Omega_g$ as the indices $(\theta, r)$ of pixels whose distance from the trajectory is smaller than half the chosen aperture diameter $\rho\lambda/D$, with $\rho > 0$. This set can be expressed as:

$$ \Omega_g = \{(\theta, r) \in [t] \times [n]^2 \mid \|r - g_t\|_2 < \frac{1}{2} \frac{\lambda}{D} \rho\}. \quad (7) $$
**Algorithm 2** *AMAT<sub>L2</sub>*

**Input:** Image sequence $M \in \mathbb{R}^{\times n^2}$, possible trajectories $\mathcal{P}$, rank $k$, maximum number of iteration for AMAT $i_{\max}$, threshold for relative change $\epsilon$

**Output:** low-rank component $L^{(i)}$ and flux $a_g^{(i)}$ for each trajectory.

1: **for** all trajectories $P_g \in \mathcal{P}$ **do**

2: $a_g^{(0)} = 0$

3: **for** $i = 1 : i_{\max}$ **do**

4: $\dot{U}, \dot{S}, \dot{V}^{\top} = H_k^{\text{SVD}}(M - a_g^{(i-1)}P_g)$

5: $L^{(i)} = \dot{U}\dot{S}\dot{V}^{\top}$

6: Compute $a_g^{(i)}$ using (6)

7: **if** $|a_g^{(i)} - a_g^{(i-1)}|/|a_g^{(i)}| < \epsilon$ **then**

8: **break**

9: **end if**

10: **end for**

11: **end for**

If (4) is set with the norm L1, we solve the problem Eq. (5a) with the algorithm suggested in Sect. 2.1. We initialize the algorithm with the randomized SVD solution in Algorithm 3. Then, we tackle the problem Eq. (5b) by

$$a_g^{(i)} = \arg \min_a \sum_{(\theta,r) \in \Omega_g} \frac{|R_g(\theta, r) - aP_g(\theta, r)|}{\sigma_{R(r)}}. \quad (8)$$

Solving Eq. (8) is an instance of the weighted least absolute deviation (LAD) problem, which unlike least squares does not have a closed form solution. In general, L1 minimization can be solved by a number of efficient iterative methods. However, in our specific case, it is possible to compute the solution even more efficiently. Since the objective function is a convex piecewise linear function $\mathbb{R} \to \mathbb{R}$ with intervals between the points $R(\theta, r)/P_g(\theta, r)$, $(\theta, r) \in \Omega_g$, its minimum is attained at one of the $(\tau n^2)$ kink points. These kink points are the boundary points where different linear segments of the piecewise function meet, and they can be easily exhaustively searched.

As part of our evaluation of the performance of the AMAT algorithms, we delve into the impact of iteration counts on its results. Iterative processes are integral to the algorithm, making it essential to investigate how its outcomes evolve over successive iterations. To illustrate the behavior of the algorithm, we show in Fig. 3 how the estimated flux $a_g$ evolves in iterations when the trajectory corresponds to the correct location of the planet and when it does not. We simulated this scenario by injecting a fake planet into the 51 Eri dataset in order to observe the changes in the flux values at the planet pixels. The results for both norms show that there is a considerable amount of change in the flux $a_g$ for the $P_g$ in the planet pixels, whereas the change in the flux $a_g$ for the $P_g$ in the pixels without a planet is very small.

We defined the set $G$ of positions of the planet as a collection of all points excluding the pixels of the host star and the pixels in the corners and edges of the images. In our algorithm, we apply the following steps to construct the residual cube and the flux map:

1. Select a pixel, denoted as $(x, y)$ from the set $G$.

2. Take the pixels of the annulus centered on the star with an inner radius of $r$-FWHM and an outer radius of $r$+FWHM, where $r$ represents the distance from the center of the star to the point $(x, y)$ resulting in an annulus width of 2FWHM.

**Algorithm 3** *AMAT<sub>L1</sub>*

**Input:** Image sequence $M \in \mathbb{R}^{\times n^2}$, possible trajectories $\mathcal{P}$, rank $k$, maximum number of iteration for AMAT $i_{\max}$, maximum number of iteration for L1-LRA $\ell_{\max}$, threshold for relative change $\epsilon$

**Output:** Low-rank component $L^{(i)}$ and flux $a_g^{(i)}$ for each trajectory.

1: **for** all trajectories $P_g \in \mathcal{P}$ **do**

2: $a_g^{(0)} = 0$

3: $\dot{U}, \dot{S}, \dot{V}^{\top} = H_k^{\text{SVD}}(M - a_g^{(0)}P_g)$

4: $U^{(0)} = \dot{U}\dot{S}$; $V^{(0)} = \dot{V}$

5: **for** $i = 1 : i_{\max}$ **do**

6: $U^{(i)}, V^{(i)} =$

7: L1-LRA($M - a_g^{(i-1)}P_g$, $U^{(i-1)}$, $V^{(i-1)}$, $k$, $\ell_{\max}$)

8: $L^{(i)} = U^{(i)}V^{(i)\top}$

9: Compute $a_g^{(i)}$ using (8)

10: **if** $|a_g^{(i)} - a_g^{(i-1)}|/|a_g^{(i)}| < \epsilon$ **then**

11: **break**

12: **end if**

13: **end for**

14: **end for**

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Intensity $a$ of the planet against the number of iterations of Eqs. (5a)–(5b). Top: intensity $a_g^{(i)}$ obtained using the L1 norm. Bottom: intensity $a_g^{(i)}$ obtained using the L2 (Frobenius) norm. The blue plots show how the intensity changes in each iteration when we choose $P_g$ in the location of the planet. The orange plots show how the intensity changes in each iteration when we choose $P_g$ in a location without a planet.</figcaption></figure>

3. Apply the AMAT algorithm: Iteratively estimate a low-rank matrix that encapsulates the background, including
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Representation of the residual cube and flux map construction. Example pixels, visualized in blue and green, have been selected from possible trajectories. For each selected pixel, a low-rank matrix is obtained, which is used to construct a residual cube and calculate a flux value. After de-rotating the residual cube, an array is produced. Combining these arrays from all possible trajectories forms the de-rotated residual cube. Similarly, flux values from all points along the possible trajectories are combined to form the flux map.</figcaption></figure>

quasi-static speckles, and the flux $a_g$ associated with the exoplanet.

4. For residual cube:

(a) Subtract the low-rank matrix from the annulus pixels of the original data matrix.

(b) Assign the residual values from each frame corresponding to (x, y) into the same positions in an empty cube referred to as the de-rotated residual cube.

(c) Repeat the processes 1–4(b) for all pixels in set G to fill the de-rotated residual cube.

5. For flux map:

(a) Assign the flux $a_g$ to (x, y) into an empty frame while ensuring it matches the dimensions of any frame in the data cube.

(b) Repeat the processes 1–3 and 5(a) for all pixels in set G in order to fill the flux map.

For every $g$ in G, the flux $a_g$ can be interpreted as follows: If there is a single planet located at $g$ in the first frame (and assuming that Eq. (5b) is solved exactly), then its flux that best explains the data is $a_g$. Since we apply this algorithm for each position $g$, we calculate $a_g$ for each trajectory, regardless of whether there is more than one planet or no planets at all. The entire process is illustrated in Fig. 4 and implemented in the AMAT Python package<fnref n="1" />. From the flux map, detection is then performed by producing an S/N map and applying a threshold, as described in Sect. 3.1. As for the residual cube, it is made use of in Sect. 4.

In existing studies that employ iterative PCA for disk and exoplanet detection (Pairet et al. 2021; Stapper & Ginski 2022; Juillard et al. 2023), the process typically involves applying the full process of PCA steps to generate a median frame. This median frame is then rotated according to parallactic angles and subtracted from each frame of the cube, leading to the creation of a new residual cube and a new median frame. This cycle of

rotating, subtracting, and then creating new residual and median frames is the part that is iteratively repeated. In contrast, when using the AMAT algorithm, which is tailored for point source detection, instead of subtracting the median frame, we subtract a matricized cube $P_g$ multiplied by the intensity $a_g$ identified at each step, effectively isolating only the planetary signature because the pixels outside the trajectory are zero. This approach significantly enhances the separation performance between the background signal and the planetary signature, improving overall detection efficiency.

2.3. *Flux estimation*

The detection of exoplanets is followed by the characterization of the planets, which encompasses estimating their positions and the intensity relative to the host star. Different algorithms, such as the negative fake companion (NEGFC) method (Lagrange et al. 2010; Marois et al. 2010; Wertz et al. 2017), ANDROMEDA via maximum likelihood estimation (Cantalloube et al. 2015), PACO estimation (Flasseur et al. 2018) are employed for this purpose. In our work, when using the AMAT algorithm, we produced a flux map consisting of the intensities $a_g$ corresponding to each trajectory. Following planet detection, we used this flux map to estimate the intensity of the planet at its detected location. We obtained these intensity values directly from the flux map without accounting for sub-pixel precision.

To test this method, we injected a fake planet into the 51 Eri dataset. In this scenario, we already know the intensity of the injected planet. In Fig. 5, we apply the L2 norm version of the AMAT algorithm ($AMAT_{L2}$) for different ranks from two to 20. Although we observed that the algorithm requires more iterations at large ranks, the algorithm approximates the injected intensity for many rank possibilities. Similarly, we applied the L1 norm version of the AMAT algorithm ($AMAT_{L1}$) in Fig. 6. However, unlike $AMAT_{L2}$, $AMAT_{L1}$ only gets close to the injected intensity at small ranks and not at all when using large ranks. The

<footnote n="1">¹ The AMAT algorithm Python package is available on GitHub: https://github.com/hazandaglayan/AMAT</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** *AMAT*<sub>L2</sub> for the flux estimation of an exoplanet. Black dashed line represents the intensity of the injected planet. The algorithm is applied with different ranks ranging from 2 to 20. The iteration number for each rank varies (terminated before reaching the maximum iteration), as the changes in intensity become smaller than the given threshold.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** *AMAT*<sub>L1</sub> for the flux estimation of an exoplanet. Same inputs as Fig. 5.</figcaption></figure>

underlying reason for this is traced back to *AMAT*<sub>L1</sub>'s initialization process; it begins with a randomized SVD for its initial approximation but subsequently relies on the outcomes of preceding steps for further initialization within its iterative process. This dependency introduces a sensitivity to the initial conditions, which coupled with the inability of the algorithm to guarantee recovery of the global optimum, restricts the effectiveness of *AMAT*<sub>L1</sub> to accurately capture the injected intensity.

Our aim was also to assess the effect of injecting a planet at different positions within the pixel. Initially, we injected a planet with an intensity of 100 at the center of the pixel, which resulted in an intensity of 100.07 after the *AMAT*<sub>L2</sub> algorithm ran. Next, we injected planets with the same intensity at each of the four corners of the pixel in four separate instances. Upon analysis of the trajectories originating from this pixel in the algorithm, the resulting intensities were observed to be 97.23, 96.31, 97.19, and 96.36, respectively. We attribute this flux loss to the intra-pixel variation of flux in the injected PSF, considering the fact that the flux map is only estimated at the center of each pixel.

### 3. Performance evaluation

In this section, we perform a comprehensive performance evaluation of the proposed algorithms by using S/N maps for visual comparison, contrast curves to assess sensitivity limits, and receiver operating characteristic (ROC) curves derived

from datasets into which synthetic planets were injected. These approaches allow for a detailed comparison of algorithmic effectiveness. We first made use of the SPHERE-IRDIS 51 Eri dataset described in Sect. 1, cropping frames to 100-by-100 pixels to reduce computation time for S/N maps and ROC curves. For contrast curves, we used a larger frame size of 200-by-200 pixels to evaluate the performance at angular separations up to the edge of the SPHERE well-corrected field. Our analysis then extends to the Exoplanet Imaging Data Challenge (EIDC) datasets, using S/N maps as a comparison metric. This evaluation strategy was designed to rigorously assess the performance of algorithms for distinguishing planetary signals from noise in different datasets, allowing us to compare them on different well-characterized observational datasets. Finally, we present a comparative study to showcase the efficacy of the AMAT algorithm in flux estimation in Sect. 3.4.

#### 3.1. Signal-to-noise ratio and sensitivity limits

To begin our comparison, we aimed to assess the performance of various algorithms in detecting the real planet within the 51 Eri dataset. To do so, we compared the S/N map obtained after AnnPCA and after the annular version of the L1-LRA algorithm (AnnL1-LRA) described in Sect. 2.1, with the S/N maps built from the output of the annular version of the *AMAT*<sub>L1</sub> and *AMAT*<sub>L2</sub> algorithms. To build these S/N maps, we paved the annulus containing the pixel of interest with non-overlapping apertures, and we extracted the central pixel of each aperture, as proposed by Bonse et al. (2023). This ensured that the signal and noise were defined in the same way (pixel-wise) for AMAT algorithms and that the pixels used to build the noise estimation were not correlated. We also obtained S/N maps using the VIP package Gomez Gonzalez et al. (2016); Christiaens et al. (2023), which is a more common version using the concept of S/N computation of Mawet et al. (2014). Figure 7 shows the S/N maps generated using the median frame obtained from AnnPCA and AnnL1-LRA as well as the flux map obtained from AMAT algorithms. The S/N map after *AMAT*<sub>L1</sub> outperforms the other three algorithms on both S/N map versions. For each algorithm, the S/N maps proposed by Bonse et al. (2023) yield higher values than the VIP version at the location of the planet. However, in AnnPCA and AnnL1-LRA, this version causes more false positives, which can be seen in Appendix C.1–C.2. Therefore, we used the VIP package for S/N map after (Ann)PCA and AnnL1-LRA in the following sections of this paper.

A more comprehensive way to assess the sensitivity of various algorithms is to build contrast curves (Mawet et al. 2014). To do so, we relied on the VIP package. In Fig. 8, we compare full-frame PCA, AnnPCA, AnnL1-LRA, and our AMAT algorithm using both norms. To generate the AnnPCA contrast curves, we used VIP with its default values except for the rank. We used the same intensities, which are used to obtain the contrast curve of AnnPCA, for the injected planets. Because the different ranks might yield varying results, we applied each algorithm at various ranks, increasing from five to 30 in five steps. For each algorithm, we selected the deepest contrast curve to represent the best outcome. The optimal ranks were found to be 25 for full-frame PCA, 15 for AnnPCA and *AMAT*<sub>L2</sub>, and 20 for AnnL1-LRA and *AMAT*<sub>L1</sub>. The intensities of the injected planets used to build the contrast curve for each algorithm are based on the noise values in each annulus obtained when applying AnnPCA with default values other than rank. This ensured that the contrast curves were obtained using the dataset with the same injected planets. As can be observed in Fig. 8, the
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Signal-to-noise ratio maps after (from left to right) AnnPCA, AnnL1-LRA, AMAT<sub>L2</sub>, and AMAT<sub>L1</sub>, respectively. Top: S/N maps calculated using VIP package. Bottom: S/N maps proposed by Bonse et al. (2023). The white circles highlight the location of the planet.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Contrast curves for full-frame PCA, AnnPCA and AnnL1-LRA, AMAT<sub>L2</sub>, and AMAT<sub>L1</sub> for the 51 Eri dataset.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Receiver operating characteristic curve of S/N maps. We compare full-frame PCA, AnnPCA and AnnL1-LRA, AMAT<sub>L2</sub>, and AMAT<sub>L1</sub> for the 51 Eri dataset. We used the square root to scale the axes in order to better see the low FPR regime.</figcaption></figure>

performance using full-frame PCA tends to be the worst among the compared methods. It is followed by AnnPCA, which shows a slightly better performance. AMAT<sub>L2</sub> outperforms AnnPCA by a small margin. The performance of AnnL1-LRA fluctuates based on the separation, indicating a performance where it sometimes outperforms or underperforms compared to the other methods. Finally, the AMAT<sub>L1</sub> algorithm demonstrates increased efficacy, outperforming other algorithms across all separations.

## 3.2. ROC curves

While S/N and contrast curves are useful to illustrate the gain provided by AMAT, they do not explore the behavior of the algorithms as a function of the detection threshold, which can be done through ROC curves. Building ROC curves relies on injecting a large quantity of synthetic planets into the chosen dataset. Our process began with the removal of the real planet present in the dataset using the VIP package (Gomez Gonzalez et al. 2017; Christiaens et al. 2023). Subsequently, we injected

synthetic planets with an intensity of 1.5 times the standard deviation of the values in the cube at a distance of $2\lambda/D$ from the star. The injections were placed methodically, starting from 0 to 360 degrees in increments of 3.6 degrees and placing a synthetic planet per scenario. This approach resulted in a total of 100 different cases for evaluation, effectively covering the entire 360-degree span around the star. For each scenario, we applied the algorithms and then examined the location where the planet was injected. A detection within a specified aperture exceeding a predefined threshold was counted as a true positive (TP). Then, we checked the other apertures for the presence of the signal. If a signal above the threshold was found within these apertures, it was classified as a false positive (FP), and the absence of such a signal resulted in a true negative (TN) classification. This examination across all apertures facilitated the construction of a ROC curve. Given the critical importance of maintaining a low number of FPs in exoplanet detection, our analysis primarily focuses on achieving a high true positive rate (TPR) without incurring
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Signal-to-noise ratio maps after $AMAT_{L1}$. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.</figcaption></figure>

false positives. To better visualize and compare the ROC curves, especially to highlight the performance at the minimal false positive rate (FPR), we employed a transformed plot of the square root of TPR versus FPR, allowing for an enhanced representation of algorithmic efficiency.

We compared our AMAT algorithms with the results of the full-frame PCA, AnnPCA, and AnnL1-LRA algorithms. We used the ranks where we obtained the best contrast curves for each algorithm. The results displayed in Fig. 9 show that our method consistently outperforms the results of full-frame PCA, AnnPCA, and AnnL1-LRA in terms of ROC curves. Moreover, similar to the findings for the S/N map, the results obtained using the L1 norm are better than those obtained with the L2 norm.

Furthermore, AnnL1-LRA shows a slightly better performance compared to both AnnPCA and the full-frame PCA.

### 3.3. EIDC results

As part of our evaluation process, we utilized the datasets from EIDC (Cantalloube et al. 2020), as specified in Table A.1, to provide varying datasets for comparing and assessing our algorithm alongside state-of-the-art HCI algorithms. This challenge encompasses a diverse array of ADI sequences, including nine ADI datasets, with a total of 20 injected planet signals. These signals exhibited varying contrasts and positional
<figure><img src="image_12.png" /><figcaption>Fig. 11. Signal-to-noise ratio maps after AMAT<sub>L2</sub>. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.</figcaption></figure>

coordinates to create an evaluation from three different instruments: VLT/SPHERE-IRDIS, Keck/NIRC2, and LBT/LMIRCam, with each providing three datasets.

As with the EIDC, our assessment involves the generation of detection maps by each algorithm for every ADI sequence. A detection map refers to applying a threshold for detection to the S/N maps or LRMs detailed in Sect. 4. In the result report of EIDC, a set of standard metrics were employed to compare the performance of these detection maps. Therefore, to apply a common metric for comparing the EIDC algorithms, we calculated the F1 scores using the same definition:

$$ \text{F1 score} = \frac{2\text{TP}}{2\text{TP} + \text{FP} + \text{FN}} \qquad (9) $$

To determine the values of TP, FP, and FN, we needed to decide on a threshold value, as this can significantly impact performance. To determine the appropriate threshold, we injected synthetic planets into the 51 Eri dataset. We created three distinct datasets, illustrated in Appendix D, and each set was injected with three, two, or four planets, respectively, placed at different locations. Subsequently, we ran our algorithm. We experimented with varying the thresholds in order to identify the threshold yielding the highest F1 score. Fig. B.1 illustrates the plot of the F1 score against the threshold for the 51 Eri dataset with synthetic planets. Based on this analysis, if we selected a threshold between 6.6 and 7.7 for the S/N maps after AMAT<sub>L1</sub>, we could detect all planets without false positives, and this allowed us to achieve the highest F1 score.
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Box plot for flux estimation for planets that are injected in $4\lambda/D$ (left) and $10\lambda/D$ (right) separation. The dashed red line represents the injected flux value.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

In the published results report, we have knowledge of the locations of the planets, and we evaluated various threshold values as evidenced by the F1 scores presented in Fig. B.2. Since we also established our threshold range by testing injections into the 51 Eri dataset, we conducted a posteriori verification to ensure that the thresholds also yield a high F1 score when applied to the EIDC datasets. Participating in the EIDC, where we had multiple opportunities to submit and observe the F1 score, allowed us to undergo a similar process. Selecting the value that yields the highest F1 score represents a fair and comparable approach. In our comparison, we used the threshold of 7.5 within the 6.6–7.7 range in which we obtained the highest F1 score for the 51 Eri dataset.

The effectiveness of the *AMAT*<sub>L1</sub> algorithm is illustrated in Fig. 10 through S/N maps. The analysis of the VLT/SPHERE-IRDIS datasets showcases the proficiency of the algorithm, where it successfully identified all six planets, with only one false positive. In the Keck/NIRC2 dataset, the algorithm accurately detected four out of seven planets with three false positives, and in the LBT/LMIRCam dataset five out of seven planets without any false positives. The results of our experiments using AnnPCA on the EIDC dataset are shown in Fig. C.1 for reference, to illustrate the significant gain provided by AMAT in terms of TPR. Additionally, the *AMAT*<sub>L2</sub> algorithm exhibits a performance between *AMAT*<sub>L1</sub> and AnnPCA in Fig. 11. It detects five out of six planets with one false positive in the VLT/SPHERE-IRDIS datasets, three out of seven planets in the Keck/NIRC2 dataset with one false positive, and five out of seven planets with two false positives in the LBT/LMIRCam dataset. We compute F1 scores based on these maps and compare them with the best results of the published results report in Fig. B.2.

### 3.4. Flux estimation performance

For the purpose of flux estimation, we conducted a comparative analysis of our *AMAT*<sub>L2</sub> algorithm against the NEGFC method, which uses AnnPCA, implemented in VIP, where the companion parameters are estimated either through a Nelder Mead minimization method or through a Markov chain Monte Carlo (MCMC) method (Wertz et al. 2017). Some other forward modeling or inverse problem approaches (Cantalloube et al. 2015; Wang et al. 2016) also have the potential to provide accurate flux measurements and may outperform both NegFC and AMAT.

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Receiver operating characteristic curve of LRMs. We compare full-frame PCA, AnnPCA and AnnL1-LRA, *AMAT*<sub>L2</sub>, and *AMAT*<sub>L1</sub>.</figcaption></figure>

However, a full comparison of these methods is beyond the scope of this paper. We designed two distinct sets of comparisons, one with the synthetic planets injected at a small separation ($4\lambda/D$) and the other at a large separation ($10\lambda/D$). In both scenarios, we created various cases by rotating the position of the injected planet from 0 to 350 degrees in increments of 10 degrees. This resulted in 36 different cases for each scenario. This systematic approach was applied to both the $4\lambda/D$ and the $10\lambda/D$ separations, thereby ensuring a comprehensive evaluation. Planets were injected at specified radii and angles without consideration for whether the location falls precisely in the center or on the edge of the pixel, ensuring consistency across all cases. To decide on the rank for each algorithm, we used the approach suggested in VIP (Christiaens et al. 2023; Gomez Gonzalez et al. 2017), which finds the optimal number of principal components in terms of S/N using PCA. In the AMAT algorithm, after a maximum of 50 iterations or when the relative change of the intensity $a_g^{(i)}$ is minor (i.e., when $|a_g^{(i)} - a_g^{(i-1)}|/|a_g^{(i)}| < 10^{-3}$), we simply stopped the algorithm.

In Fig. 12, the *AMAT*<sub>L2</sub> algorithm demonstrates excellent accuracy in flux estimations for both small and large separations, as evidenced by the median line of the boxplots precisely aligning with the injected intensities. In contrast, the boxplots for the Nelder-Mead and MCMC methods exhibit the presence of outliers, indicating lower precision in their flux estimations.
<figure><img src="image_16.png" /><figcaption>Fig. 14. Likelihood ratio maps after AMAT<sub>L1</sub> for EIDC Datasets. In these maps, white circles represent TP, red square denotes FP, and red circles signify FN.</figcaption></figure>

Specifically, the Nelder-Mead method tends to produce higher flux values, often deviating significantly from the injected values, while the MCMC method frequently results in lower flux estimations.

4. Improving detection performance with likelihood ratio maps

In Daglayan et al. (2022), we presented an enhanced approach for exoplanet detection utilizing LRMs, which offers improvements over traditional S/N maps. The LRM is derived through a maximum likelihood estimation that employs Laplacian distributions. It consists of the ratio of the maximum likelihood $\mathcal{L}_g(\hat{a}_g|R)$ to the likelihood of the null hypothesis $\mathcal{L}_g(0|R)$, which corresponds

to the absence of a planet:

$$ 
\begin{aligned}
\log \Lambda_g(R) &= \log \left( \frac{\mathcal{L}_g(\hat{a}_g|R)}{\mathcal{L}_g(0|R)} \right) \\
&= -\sum_{(\theta, r) \in \Omega_g} \frac{|R(\theta, r) - \hat{a}_g P_g(\theta, r)| - |R(\theta, r)|}{\sigma_R(r)},
\end{aligned} 
\quad (10) $$

where R is the residual cube and $\hat{a}_g$ is the estimated flux derived from solving the following optimization problem:

$$ 
\begin{aligned}
\hat{a}_g &= \underset{a}{\operatorname{argmax}} \log \mathcal{L}_g(a|R) \\
&= \underset{a}{\operatorname{argmin}} \sum_{(\theta, r) \in \Omega_g} \frac{|R(\theta, r) - a P_g(\theta, r)|}{\sigma_R(r)}.
\end{aligned} 
\quad (11) $$
<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Ranking based on the F1 score of the different algorithms in the EIDC results report. The algorithms are classified as AMAT algorithms with S/N maps and LRMs (purple), classical speckle subtraction providing residual maps (red), advanced speckle subtraction building detection maps (orange), inverse problems (blue), and supervised machine learning (green). The light, medium, and dark colors correspond to the three VLT/SPHERE-IRDIS, Keck/NIRC2, and LBT/LMIRCam datasets, respectively.</figcaption></figure>

In light of this information, we combined the AMAT algorithm with the LRM. After obtaining the residual cube, we applied the LRM algorithm to our results.

We obtained the ROC curve using the residual cubes presented in Sect. 3.2 and applied the LRM instead of the S/N map in Fig. 13. Each algorithm except AMAT<sub>L2</sub> demonstrated a higher performance compared to the S/N map results in Fig. 9. We believe this discrepancy arises from AMAT<sub>L2</sub> being based on a Gaussian distribution, while LRM relies on a Laplacian distribution. Among the algorithms, AMAT<sub>L1</sub> showed the best performance, while AnnL1-LRA also exhibited strong results. Both algorithms are Laplacian-based, which likely contributes to their superior performance.

To decide on the threshold, we used the same method as described in Sect. 3.2 for finding the threshold of S/N maps. Based on this analysis, if we choose a value between 67 and 77 for the LRMs after AMAT<sub>L1</sub>, we get the highest F1 score for the 51 Eri dataset, and as for the case of AMAT-S/N, we checked that the threshold for LRMs after AMAT<sub>L1</sub> is indeed in this range in order to maximize the F1 score on the EIDC datasets.

Figure 14 shows that the version AMAT<sub>L1</sub> using LRM proficiently identifies all exoplanets within the VLT/SPHERE-IRDIS datasets and accurately detects four out of seven planets and six out of seven planets in the LBT/LMIRCam datasets without any false positive in any dataset. These findings showcase the algorithm's high success rate, especially when compared with other algorithms reported in the EIDC results.

In the paper reporting on the EIDC results (Cantalloube et al. 2020), the algorithms are classified into four categories: classical speckle subtraction providing residual maps, advanced speckle subtraction building detection maps, inverse problems, and supervised machine learning. A comparison of the most successful of these methods in their categories, focusing on the F1 score, is presented in Fig. 15. Our AMAT<sub>L1</sub> method with both S/N map and LRM comes out as the most effective within

the category of classical speckle subtraction methods. Furthermore, the AMAT<sub>L1</sub> method with LRM exhibits a level of success comparable to that of advanced algorithms such as RSM and FMMF. This underlines the significant potential and robustness of AMAT<sub>L1</sub> in the realm of exoplanet detection.

## 5. Conclusions

In this paper, we have investigated AMAT, a new exoplanet detection method designed to improve the separation of planetary flux from static and quasi-static signals using iterative techniques. The method enhances models based on low-rank approximations such as PCA. Current approaches typically assume Gaussian noise, but recent studies suggest that residuals often follow a Laplacian distribution. Most low-rank approximation techniques still rely on PCA, which assumes Gaussian noise. To address this inconsistency, we proposed L1-LRA, which is based on assuming Laplace-distributed noise, and integrated it within the AMAT algorithm.

The AMAT algorithm, which has an open-source implementation, was thoroughly tested with various approaches. First, we compared sensitivity limits using S/N comparisons and contrast curves on the 51 Eri dataset acquired with the VLT/SPHERE-IRDIS instrument. We evaluated the performance of AMAT alongside AnnPCA and AnnL1-LRA, and our results demonstrated significant performance improvements. This was supported by ROC curve comparisons. We then benchmarked AMAT against state-of-the-art algorithms using EIDC datasets. When comparing similar speckle subtraction algorithms, AMAT delivered a competitive performance, achieving results comparable to the most successful algorithms in its category. By utilizing the LRM based on Laplacian noise, we further enhanced the algorithm's performance, allowing AMAT<sub>L1</sub> with LRM to achieve the highest F1 score among all categories. Additionally, the flux map generated by the AMAT algorithm provided accurate planetary flux measurements, even for faint planets.

One limitation of this study is the high computational cost due to the iterative nature of the algorithm and the need to apply it to each pixel, given the possibility of the planet appearing in any pixel. The computations were performed on a server equipped with an Intel CPU with 18 cores and 125 GB of RAM. While AnnPCA takes only a few minutes to process the 51 Eri dataset, AMAT<sub>L1</sub> takes approximately a day to complete the same task. The AMAT<sub>L1</sub> algorithm, while producing better results, relies on L1-LRA, which is computationally more demanding than PCA. Future studies should explore faster implementations of L1-LRA in order to improve the algorithm's speed and applicability. Another limitation involves determining the threshold for the LRM, which should be standardized rather than observation based, ensuring a more automated and reliable method.

*Acknowledgement.* This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 819155); VC thanks the Belgian Federal Science Policy Office (BELSPO) for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000142531; the Fonds de la Recherche Scientifique - FNRS and the Fonds Wetenschappelijk Onderzoek - Vlaanderen under EOS Project no. 30468160; the Fonds de la Recherche Scientifique-FNRS under Grant no. T.0001.23; NG acknowledges the support by the European Union (ERC consolidator, eLinoR, no. 101085607); LJ thanks the FRS-FNRS for the funding related to the PDR project QuadSense (T.0160.24).
References

Absil, O., Milli, J., Mawet, D., et al. 2013, A&A, 559, L12

Amara, A., & Quanz, S. P. 2012, MNRASS, 427, 948

Beuzit, J.-L., Vigan, A., Mouillet, D., et al. 2019, A&A, 631, A155

Bonde, M. J., Garvin, E. O., Gebhard, T. D., et al. 2023, AJ, 166, 71

Bowler, B. P. 2016, PASP, 128, 102001

Cantalloube, F., Mouillet, D., Mugnier, L., et al. 2015, A&A, 582, A89

Cantalloube, F., Gomez-Gonzalez, C., Absil, O., et al. 2020, in Adaptive Optics Systems VII, 11448, SPIE, 1027

Cantero, C., Absil, O., Dahlqvist, C.-H., & Van Droogenbroeck, M. 2023, A&A, 680, A86

Christiaens, V., Gonzalez, C., Farkas, R., et al. 2023, J. Open Source Softw., 8, 4774

Daglayan, H., Vary, S., Cantalloube, F., Absil, P.-A., & Absil, O. 2022, in 2022 IEEE 5th International Conference on Image Processing Applications and Systems (IPAS), IEEE, 1

Daglayan, H., Vary, S., & Absil, P.-A. 2023a, in ESANN 2023-European Symposium on Artificial Neural Networks, Computational Intelligence and Machine Learning

Daglayan, H., Vary, S., Leplat, V., Gillis, N., & Absil, P.-A. 2023b, Proceedings of BNAIC/BeNeLearn 2023, 1

Dahlqvist, C.-H., Cantalloube, F., & Absil, O. 2020, A&A, 633, A95

Eckart, C., & Young, G. 1936, Psychometrika, 1, 211

Eriksson, A., & Van Den Hengel, A. 2010, in 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, IEEE, 771

Flasseur, O., Denis, L., Thiébaut, É., & Langlois, M. 2018, A&A, 618, A138

Galicher, R., & Mazoyer, J. 2023, Comptes Rend. Phys., 24, 1

Gao, J., Kwan, P. W., & Guo, Y. 2009, Neurocomputing, 72, 1242

Gillis, N., & Plemmons, R. J. 2011, Opt. Eng., 50, 027001

Gillis, N., & Vavasis, S. A. 2018, Math. Oper. Res., 43, 1072

Gomez Gonzalez, C. A., Absil, O., Absil, P.-A., et al. 2016, A&A, 589, A54

Gomez Gonzalez, C. A., Wertz, O., Absil, O., et al. 2017, AJ, 154, 7

Gonzalez, C. G., Absil, O., & Van Droogenbroeck, M. 2018, A&A, 613, A71

Guyon, O. 2018, ARA&A, 56, 315

Halko, N., Martinsson, P.-G., & Tropp, J. A. 2011, SIAM Rev., 53, 217

Juillard, S., Christiaens, V., & Absil, O. 2023, A&A, 679, A52

Ke, Q., & Kanade, T. 2003, Robust subspace computation using L1 norm, Tech. Rep. CMUCS-03-172, Carnegie Mellon University, Pittsburgh, PA

Ke, Q., & Kanade, T. 2005a, in 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'05), 1, 739

Ke, Q., & Kanade, T. 2005b, in 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'05), 1, IEEE, 739

Lafreniere, D., Marois, C., Doyon, R., Nadeau, D., & Artigau, É. 2007, ApJ, 660, 770

Lagrange, A.-M., Bonnefoy, M., Chauvin, G., et al. 2010, Science, 329, 57

Macintosh, B., Graham, J. R., Ingraham, P., et al. 2014, PNAS, 111, 12661

Marois, C., Lafreniere, D., Doyon, R., Macintosh, B., & Nadeau, D. 2006, ApJ, 641, 556

Marois, C., Macintosh, B., & Véran, J.-P. 2010, in Adaptive Optics Systems II, 7736, SPIE, 595

Martinache, F., Guyon, O., Lozi, J., et al. 2009, The Subaru coronagraphic extreme AO project

Mawet, D., Milli, J., Wahhaj, Z., et al. 2014, ApJ, 792, 97

NASA Exoplanet Archive 2024, Planetary Systems

Pairet, B., Cantalloube, F., Gomez Gonzalez, C. A., Absil, O., & Jacques, L. 2019, MNRAS, 487, 2262

Pairet, B., Cantalloube, F., & Jacques, L. 2021, MNRAS, 503, 3724

Ren, B., Pueyo, L., Zhu, G. B., Debes, J., & Duchêne, G. 2018, ApJ, 852, 104

Ruffio, J.-B., Macintosh, B., Wang, J. J., et al. 2017, ApJ, 842, 14

Samland, M., Mollière, P., Bonnefoy, M., et al. 2017, A&A, 603, A57

Song, Z., Woodruff, D. P., & Zhong, P. 2017, in Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing, 688

Soummer, R., Pueyo, L., & Larkin, J. 2012, ApJ, 755, L28

Stapper, L. M., & Ginski, C. 2022, A&A, 668, A50

Thompson, W., & Marois, C. 2021, AJ, 161, 236

Vary, S., Daglayan, H., Jacques, L., & Absil, P.-A. 2023, in ICASSP 2023 – 2023 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 1

Wahhaj, Z., Cieza, L. A., Mawet, D., et al. 2015, A&A, 581, A24

Wang, J. J., Graham, J. R., Pueyo, L., et al. 2016, AJ, 152, 97

Wertz, O., Absil, O., González, C. G., et al. 2017, A&A, 598, A83

Zheng, Y., Liu, G., Sugimoto, S., Yan, S., & Okutomi, M. 2012, in 2012 IEEE Conference on Computer Vision and Pattern Recognition, IEEE, 1410
## Appendix A: Datasets

The properties of the datasets which we used in this paper are given in Table A.1. Nine of these datasets are from the EIDC, while the final dataset, 51 Eri, is used in all experiments except those specific to the EIDC results.

<table><caption>Table A.1: Properties of the datasets.</caption><thead><tr><th>ID</th><th>Telescope/Instr</th><th>N<sub>img</sub><br>(px × px)</th><th>N<sub>t</sub></th><th>Δ<sub>field</sub><br>(°)</th><th>N<sub>p</sub></th></tr></thead><tbody><tr><td>sph1</td><td>VLT/SPHERE-IRDIS</td><td>160×160</td><td>252</td><td>40.3</td><td>1</td></tr><tr><td>sph2</td><td>VLT/SPHERE-IRDIS</td><td>160×160</td><td>80</td><td>31.5</td><td>0</td></tr><tr><td>sph3</td><td>VLT/SPHERE-IRDIS</td><td>160×160</td><td>228</td><td>80.5</td><td>5</td></tr><tr><td>nrc1</td><td>Keck/NIRC2</td><td>321×321</td><td>29</td><td>53.0</td><td>4</td></tr><tr><td>nrc2</td><td>Keck/NIRC2</td><td>321×321</td><td>40</td><td>37.3</td><td>3</td></tr><tr><td>nrc3</td><td>Keck/NIRC2</td><td>321×321</td><td>50</td><td>166.9</td><td>0</td></tr><tr><td>lmr1</td><td>LBT/LMIRCam</td><td>200×200</td><td>4838</td><td>153.4</td><td>2</td></tr><tr><td>lmr2</td><td>LBT/LMIRCam</td><td>200×200</td><td>3219</td><td>60.6</td><td>2</td></tr><tr><td>lmr3</td><td>LBT/LMIRCam</td><td>200×200</td><td>4620</td><td>91.0</td><td>3</td></tr><tr><td>51 Eri</td><td>VLT/SPHERE-IRDIS</td><td>200×200</td><td>256</td><td>42.0</td><td>1</td></tr></tbody></table>

**Notes.** N<sub>img</sub> is the size of the images, N<sub>t</sub> is the number of frames, Δ<sub>field</sub> is the total field rotation of the planets, and N<sub>p</sub> is the number of the injected planet for EIDC datasets and the number of the real planet for 51 Eri dataset.
Appendix B: F1 score versus threshold for AMAT<sub>L1</sub>

<figure><img src="image_18.png" /><figcaption>Fig. B.1: F1 score of AMAT<sub>L1</sub> algorithm with S/N map using various thresholds for 51 Eri datasets (left) and EIDC datasets (right).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. B.2: F1 score of AMAT<sub>L1</sub> algorithm with LRM using various thresholds for 51 Eri datasets (left) and EIDC datasets (right).</figcaption></figure>
Appendix C: EIDC results after AnnPCA

<figure><img src="image_20.png" /><figcaption>Fig. C.1: Signal-to-noise ratio maps after AnnPCA using VIP package. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. C.2: Signal-to-noise ratio maps proposed by Bonse et al. (2023), after AnnPCA. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. C.3: Likelihood ratio maps after AnnPCA for EIDC datasets. In these maps, white circles represent TP, red squares denote FP, and red circles signify FN.</figcaption></figure>
**Appendix D: The location of injected planets**

<figure><img src="image_23.png" /><figcaption>Fig. D.1: Signal-to-noise ratio maps for 51 Eri datasets with injected planets. In these S/N maps, white circles represent TP.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. D.2: Likelihood ratio maps for 51 Eri datasets with injected planets. In these maps, white circles represent TP.</figcaption></figure>