© The Authors 2024

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 Sub- aru/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 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 maximumlikelihood 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 ${ℝ}^{t\times n^2}$ containing observations of t unfolded frames, where each row represents a single vectorized frame of size n × n. Assuming there is a single planet located at position 𝑔 є [n] × [n] within the first frame, with [n] = {1,2,…, n}, the model for M can be written as $$M=L+a_g{P}_g+E,\mathrm{rank}(L)\le k,P_g\in \mathcal{P},$$(1)

where L denotes the low-rank model for the stellar diffraction pattern; E stands for the noise; a_𝑔 is the intensity of the planet referred to as the flux; $P_g\in \mathcal{P}\subset {ℝ}^{t\times n^2}$ is the planet signature along the trajectory, illustrated in Fig. 1; and 𝒫 is the set of all feasible planet signatures. To construct P_𝑔, we started from a time-by-pixel matrix with zero entries, and in each unfolded frame (i.e., row of the matrix), we placed a copy of the normalized reference PSF at the location occupied by a planet that is at position 𝑔 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_𝑔 is a constant matrix for one single planet position 𝑔.

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 contrast, 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 maximum 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: $$\widehat{L}=\underset{L}{\arg \mathrm{\hspace{0.17em}min}}\Vert M-L{\Vert }_2\text{\hspace{1em} subject to \hspace{1em}}\mathrm{rank}(L)\le k,$$(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: $$\widehat{L}=\underset{L}{\arg \mathrm{\hspace{0.17em}min}}\Vert M-L{\Vert }_1\text{\hspace{1em} subject to \hspace{1em}}\mathrm{rank}(L)\le k.$$(3)

This approach allows us to maintain a consistent noise assumption 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 alternating 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 $\widehat{U}\text{ and }\widehat{V}$ of the low-rank approximation $\widehat{L}=\widehat{U}{\widehat{V}}^{\top }$. This algorithm relies on the fact that the columns of the matrix $\widehat{V}\in {ℝ}^{n^2\times k}$ form a basis of a subspace of ${ℝ}^{n^2}$ that best approximates the data frames in the L1 sense, and for each data frame, the corresponding row of $\widehat{U}\in {ℝ}^{t\times k}$ contains the weights of the L1-best approximation of the data frame in this basis. Here, A_j denotes the j-th column of the matrix A.

To solve the minimization problem in steps 5–6 of Algorithm 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 assumptions, we fit Gaussian and Laplacian distributions to the residual 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λ/D separation and one more distant from the star at 8λ/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).

Fig. 1 Cube of the planet signature constructed by rotating the position 𝑔 of the PSF function along the trajectory.

Fig. 2 Histograms of the residual cube after low-rank approximation is applied for small and large separations for a 51 Eri dataset.

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_𝑔, we propose estimating them simultaneously with the following optimization problem: $$\underset{L\in {ℝ}^{t\times n^2},a_g\in ℝ}{\min }\Vert M-L-a_g{P}_g{\Vert }_{\text{\#}}\text{\hspace{1em}s}.\text{t}.\text{\hspace{1em}rank}(L)\le k,P_g\in \mathcal{P},$$(4)

where || · ||_# 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 alternatingly solving the following two sub-problems until a stopping criterion is met: $$L^{(i)}=\underset{L\in {ℝ}^{t\times n^2}}{\arg \mathrm{\hspace{0.17em}min}}\Vert M-L-a_g^{(i-1)}{P}_g{\Vert }_{\text{\#}},$$(5a) $$a_g^{(i)}=\underset{a\in ℝ}{\arg \mathrm{\hspace{0.17em}min}}\Vert M-L^{(i)}-a{P}_g\Vert \text{\#}.$$(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⁽ⁱ⁾ in Eq. (5a) amounts to a k-truncated SVD, denoted by ${\text{H}}_k^{\text{SVD}}(\cdot )$. In practice, in order to speed up the computations, we compute L⁽ⁱ⁾ 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_𝑔 = M − L⁽ⁱ⁾ with the planet signature P_𝑔: $$a_g^{(i)}=\frac{{{\displaystyle \sum }}_{(\theta ,r)\in {\text{Ω}}_g}{R}_g(\theta ,r)P_g(\theta ,r){\sigma }_{R_g(r)}^{-2}}{{{\displaystyle \sum }}_{(\theta ,r)\in {\text{Ω}}_g}{\left(P_g(\theta ,r)\right)}^2{\sigma }_{R_g(r)}^{-2}},$$(6)

where ${\sigma }_R^2$ is the empirical variance of the residual frames computed along the time dimension and we define the set Ω_𝑔 as the indices (θ, r) of pixels whose distance from the trajectory is smaller than half the chosen aperture diameter ρλ/D, with ρ > 0. This set can be expressed as: $${\text{Ω}}_g=\{(\theta ,r)\in [t]\times {[n]}^2\mid {\Vert r-g_t\Vert }_2<\frac{1}{2}\frac{\lambda }{D}\rho \}.$$(7)

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)}=\underset{a}{\text{arg\hspace{0.17em}min}}{\displaystyle \sum _{(\theta ,r)\in {\text{Ω}}_g}}\frac{\left|R_g(\theta ,r)-a{P}_g(\theta ,r)\right|}{{\sigma }_{R(r)}}.$$(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 ℝ → ℝ with intervals between the points R(θ, r)/P_𝑔(θ, r), (θ, r) є Ω_𝑔, its minimum is attained at one of the (tn²) 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_𝑔 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_𝑔 for the P_𝑔 in the planet pixels, whereas the change in the flux a_𝑔 for the P_𝑔 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.

3. Apply the AMAT algorithm: Iteratively estimate a low- rank matrix that encapsulates the background, including quasi-static speckles, and the flux a_𝑔 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_𝑔 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 𝑔 in G, the flux a_𝑔 can be interpreted as follows: If there is a single planet located at 𝑔 in the first frame (and assuming that Eq. (5b) is solved exactly), then its flux that best explains the data is a_𝑔 . Since we apply this algorithm for each position 𝑔, we calculate ag 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¹. 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_𝑔 multiplied by the intensity a_𝑔 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.

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𝑔 in the location of the planet. The orange plots show how the intensity changes in each iteration when we choose P𝑔 in a location without a planet.

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.

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_𝑔 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 underlying reason for this is traced back to AMAT_L1’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_L1 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_L2 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.

Fig. 5 AMATL2 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.

Fig. 6 AMATL1 for the flux estimation of an exoplanet. Same inputs as Fig. 5.

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_L1 and AMAT_L2 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 AMATL1 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_L2, and 20 for AnnL1-LRA and AMAT_L1. 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 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_L2 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_L1 algorithm demonstrates increased efficacy, outperforming other algorithms across all separations.

Fig. 7 Signal-to-noise ratio maps after (from left to right) AnnPCA, AnnL1-LRA, AMATL2 , and AMATL1, 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.

Fig. 8 Contrast curves for full-frame PCA, AnnPCA and AnnL1-LRA, AMATL2, and AMATL1 for the 51 Eri dataset.

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λ/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 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.

Fig. 9 Receiver operating characteristic curve of S/N maps. We compare full-frame PCA, AnnPCA and AnnL1-LRA, AMATL2, and AMATL1 for the 51 Eri dataset. We used the square root to scale the axes in order to better see the low FPR regime.

Fig. 10 Signal-to-noise ratio maps after AMATL1. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.

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 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{F}1\text{ score }=\frac{2\text{TP}}{2\text{TP}+\text{FP}+\text{FN}}.$$(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_L1, we could detect all planets without false positives, and this allowed us to achieve the highest F1 score.

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 AMATL1 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_L2 algorithm exhibits a performance between AMAT_L1 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.

Fig. 11 Signal-to-noise ratio maps after AMATL2. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.

Fig. 12 Box plot for flux estimation for planets that are injected in 4λ/D (left) and 10λ/D (right) separation. The dashed red line represents the injected flux value.

3.4 Flux estimation performance

For the purpose of flux estimation, we conducted a comparative analysis of our AMAT_L2 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.

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λ/D) and the other at a large separation (10λ/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λ/D and the 10λ/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 $\left|a_g^{(i)}-a_g^{(i-1)}\right|/\left|a_g^{(i)}\right|<{10}^{-3}$), we simply stopped the algorithm.

In Fig. 12, the AMAT_L2 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.

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.

Fig. 13 Receiver operating characteristic curve of LRMs. We compare full-frame PCA, AnnPCA and AnnL1-LRA, AMATL2, and AMATL1.

Fig. 14 Likelihood ratio maps after AMATL1 for EIDC Datasets. In these maps, white circles represent TP, red square denotes FP, and red circles signify FN.

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 ℒ_𝑔(â_𝑔|R) to the likelihood of the null hypothesis ℒ_𝑔(0|R), which corresponds to the absence of a planet: $$\begin{array}{l}\log {\text{\hspace{0.17em}Λ}}_g(R)=\log \left(\frac{{ℒ}_g\left({\widehat{a}}_g\mid R\right)}{{ℒ}_g(0\mid R)}\right)\\ =-{\displaystyle \sum _{(\theta ,r)\in {\text{Ω}}_g}}\frac{\left|R(\theta ,r)-{\widehat{a}}_g{P}_g(\theta ,r)\right|-|R(\theta ,r)|}{{\sigma }_R(r)},\end{array}$$(10)

where R is the residual cube and â_𝑔 is the estimated flux derived from solving the following optimization problem: $$\begin{array}{l}{\widehat{a}}_g={\mathrm{argmax}}_a\log {ℒ}_g(a\mid R)\\ ={\text{argmin}}_a{\displaystyle \sum _{(\theta ,r)\in {\text{Ω}}_g}}\frac{\left|R(\theta ,r)-a{P}_g(\theta ,r)\right|}{{\sigma }_R(r)}.\end{array}$$(11)

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_L2 demonstrated a higher performance compared to the S/N map results in Fig. 9. We believe this discrepancy arises from AMAT_L2 being based on a Gaussian distribution, while LRM relies on a Laplacian distribution. Among the algorithms, AMAT_L1 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 AMATL1, 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_L1 is indeed in this range in order to maximize the F1 score on the EIDC datasets.

Figure 14 shows that the version AMAT_L1 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 AMATL1 method with both S/N map and LRM comes out as the most effective within the category of classical speckle subtraction methods. Furthermore, the AMATL1 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_L1 in the realm of exoplanet detection.

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.

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_L1 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_L1 takes approximately a day to complete the same task. The AMAT_L1 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 Scien- tifique – 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).

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.

Appendix B F1 score versus threshold for AMAT_L1

Fig. B.1 F1 score of AMATL1 algorithm with S/N map using various thresholds for 51 Eri datasets (left) and EIDC datasets (right).

Fig. B.2 F1 score of AMATL1 algorithm with LRM using various thresholds for 51 Eri datasets (left) and EIDC datasets (right).

Appendix C EIDC results after AnnPCA

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.

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.

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.

Appendix D The location of injected planets

Fig. D.1 Signal-to-noise ratio maps for 51 Eri datasets with injected planets. In these S/N maps, white circles represent TP.

Fig. D.2 Likelihood ratio maps for 51 Eri datasets with injected planets. In these maps, white circles represent TP.

References

All Tables

All Figures

	Fig. 1 Cube of the planet signature constructed by rotating the position 𝑔 of the PSF function along the trajectory.
In the text

	Fig. 2 Histograms of the residual cube after low-rank approximation is applied for small and large separations for a 51 Eri dataset.
In the text

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𝑔 in the location of the planet. The orange plots show how the intensity changes in each iteration when we choose P𝑔 in a location without a planet.

In the text

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.

In the text

	Fig. 5 AMATL2 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.
In the text

	Fig. 6 AMATL1 for the flux estimation of an exoplanet. Same inputs as Fig. 5.
In the text

	Fig. 7 Signal-to-noise ratio maps after (from left to right) AnnPCA, AnnL1-LRA, AMATL2 , and AMATL1, 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.
In the text

	Fig. 8 Contrast curves for full-frame PCA, AnnPCA and AnnL1-LRA, AMATL2, and AMATL1 for the 51 Eri dataset.
In the text

	Fig. 9 Receiver operating characteristic curve of S/N maps. We compare full-frame PCA, AnnPCA and AnnL1-LRA, AMATL2, and AMATL1 for the 51 Eri dataset. We used the square root to scale the axes in order to better see the low FPR regime.
In the text

	Fig. 10 Signal-to-noise ratio maps after AMATL1. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.
In the text

	Fig. 11 Signal-to-noise ratio maps after AMATL2. In these S/N maps, white circles represent TP, red squares denote FP, and red circles signify FN.
In the text

	Fig. 12 Box plot for flux estimation for planets that are injected in 4λ/D (left) and 10λ/D (right) separation. The dashed red line represents the injected flux value.
In the text

	Fig. 13 Receiver operating characteristic curve of LRMs. We compare full-frame PCA, AnnPCA and AnnL1-LRA, AMATL2, and AMATL1.
In the text

	Fig. 14 Likelihood ratio maps after AMATL1 for EIDC Datasets. In these maps, white circles represent TP, red square denotes FP, and red circles signify FN.
In the text

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.

In the text

	Fig. B.1 F1 score of AMATL1 algorithm with S/N map using various thresholds for 51 Eri datasets (left) and EIDC datasets (right).
In the text

	Fig. B.2 F1 score of AMATL1 algorithm with LRM using various thresholds for 51 Eri datasets (left) and EIDC datasets (right).
In the text

	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.
In the text

	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.
In the text

	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.
In the text

	Fig. D.1 Signal-to-noise ratio maps for 51 Eri datasets with injected planets. In these S/N maps, white circles represent TP.
In the text

	Fig. D.2 Likelihood ratio maps for 51 Eri datasets with injected planets. In these maps, white circles represent TP.
In the text