# A New Statistical Model of Star Speckles for Learning to Detect and Characterize Exoplanets in Direct Imaging Observations

Théo Bodrito, Olivier Flasseur, Julien Mairal, Jean Ponce, Maud Langlois, Anne-Marie Lagrange

► To cite this version:

Théo Bodrito, Olivier Flasseur, Julien Mairal, Jean Ponce, Maud Langlois, et al.. A New Statistical Model of Star Speckles for Learning to Detect and Characterize Exoplanets in Direct Imaging Observations. CVPR 2025 - IEEE / CVF Conference on Computer Vision and Pattern Recognition, Jun 2025, Nashville, United States. pp.1-15. hal-05000360

HAL Id: hal-05000360

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

Submitted on 21 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
# A New Statistical Model of Star Speckles for Learning to Detect and Characterize Exoplanets in Direct Imaging Observations

Théo Bodrito *<fnref n="1" />, Olivier Flasseur<fnref n="2" />, Julien Mairal<fnref n="3" />, Jean Ponce<fnref n="1" />, <fnref n="4" />,
Maud Langlois<fnref n="2" />, Anne-Marie Lagrange<fnref n="5" />, <fnref n="6" />

## Abstract

The search for exoplanets is an active field in astronomy, with direct imaging as one of the most challenging methods due to faint exoplanet signals buried within stronger residual starlight. Successful detection requires advanced image processing to separate the exoplanet signal from this nuisance component. This paper presents a novel statistical model that captures nuisance fluctuations using a multi-scale approach, leveraging problem symmetries and a joint spectral channel representation grounded in physical principles. Our model integrates into an interpretable, end-to-end learnable framework for simultaneous exoplanet detection and flux estimation. The proposed algorithm is evaluated against the state of the art using datasets from the SPHERE instrument operating at the Very Large Telescope (VLT). It significantly improves the precision-recall trade-off, notably on challenging datasets that are otherwise unusable by astronomers. The proposed approach is computationally efficient, robust to varying data quality, and well suited for large-scale observational surveys.<fnref n="1" />

## 1. Introduction

Direct imaging [30, 65] is an astronomical technique to probe the vicinity of young, nearby stars, where exoplanets and circumstellar disks—structures of dust and gas from which exoplanets can form—are found [36, 37]. Unlike indirect methods [57], which detect exoplanets via secondary

<figure><img src="image_2.png" /><figcaption>Figure A. Left: typical observations *y* and PSF *h* from the SPHERE instrument in ASDI mode. The synthetic exoplanet is very bright for illustration purposes. Right: temporal slice along the vertical line.</figcaption></figure>

effects like gravitational wobbles or transit dimming, direct imaging captures visual evidence of exoplanets and circumstellar disks by recording a *direct image* of their emitted flux. Analyzing this light across spectral bands provides insights into exoplanet properties (e.g., temperature, gravity), atmospheric compositions, molecular abundances, and formation processes [2, 9, 66]. While existing technology has enabled imaging of young, giant, gaseous exoplanets, next-generation ground- and space-based telescopes may soon image rocky exoplanets in habitable zones, advancing the search for life beyond our Solar System [10, 14].

<footnote n="1">¹Département d'Informatique de l'École normale supérieure (ENS-PSL, CNRS, Inria)</footnote>

<footnote n="2">²Universite Claude Bernard Lyon 1, Centre de Recherche Astrophysique de Lyon UMR 5574,
ENS de Lyon, CNRS, Villeurbanne, F-69622, France</footnote>

<footnote n="3">³Université Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK</footnote>

<footnote n="4">⁴Courant Institute and Center for Data Science, New York University</footnote>

<footnote n="5">⁵Laboratoire d'Études Spatiales et d'Instrumentation en Astrophysique, Observatoire de Paris,
Université PSL, Sorbonne Université, Université Paris Diderot</footnote>

<footnote n="6">⁶Université Grenoble Alpes, Institut de Planétologie et d'Astrophysique de Grenoble</footnote>

<footnote>*Corresponding author: theo.bodrito@inria.fr</footnote>

<footnote n="1">¹Code is available at https://github.com/theobdt/exomild</footnote>
Imaging exoplanets is challenging due to the high contrast and angular resolution required. The primary difficulty lies in the large brightness ratio (or *contrast*) between the host star and exoplanets; in the infrared, where gas giant exoplanets' thermal emission is most detectable, they are typically $10^5-10^6$ times fainter. Additionally, because exoplanets appear close to their host stars from Earth's perspective, high spatial resolution is essential to separate their faint signals. To address these challenges, observatories like the VLT are equipped with specialized instruments for direct imaging (e.g., SPHERE [4]) and cutting-edge optical technologies. Adaptive optics employs a deformable mirror to obtain sharp images by correcting for atmospheric turbulence in real time [20]. Coronagraphs (optical masks) block some starlight, further improving contrast [58].

Despite advanced optical devices, direct imaging remains challenging as exoplanet signals are still dominated by a strong *nuisance component* with approximately $10^3$ times higher contrast. This nuisance is mainly composed of structured speckles [24], caused by imperfections in optical corrections that allow residual starlight, unblocked by the coronagraph, to leak into images as a spatially correlated diffraction pattern. The speckle pattern is both spatially and spectrally correlated, exhibiting quasi-static behavior across exposures with minor fluctuations over time. In addition to speckles, other stochastic noise sources—such as thermal background, detector readout, and photon noise—add further contamination. Together, these factors create a non-stationary nuisance that varies in intensity and structure across the field of view, with higher intensity and correlation near the star. This nuisance, which closely resembles the point-like signals expected from exoplanets (instrumental point-spread function off the optical axis), is the primary limitation to direct imaging. Figure A illustrates these characteristics with a dataset recorded by SPHERE.

In this context, dedicated processing is crucial to separate exoplanet signals from the nuisance component [50]. Observational techniques like angular differential imaging (ADI [41]), spectral differential imaging (SDI [61]), and combined angular and spectral differential imaging (ASDI [66]) introduce diversity that aids in distinguishing exoplanet signals from the nuisance, see Sec. 3 for resulting image formation models. ADI takes advantage of Earth's rotation by keeping the telescope pupil fixed, causing off-axis exoplanets to follow a predictable circular trajectory, while the star remains centered. This apparent motion helps separate exoplanet signals from quasi-static speckles. SDI further improves separation by capturing images across multiple spectral channels, where speckles scale quasi-linearly with wavelength while exoplanet signals stay fixed. ASDI combines ADI and SDI to produce a four-dimensional dataset (spatial, temporal, and spectral dimensions). This diversity in A(S)DI calls for advanced algo-

rithms that can leverage these priors to optimize exoplanet detection and spectral characterization.

We propose a novel hybrid approach that combines the interpretability of a statistical framework with end-to-end learnable components, improving detection sensitivity and characterization accuracy for directly imaged exoplanets. Our main contributions include:

* A learnable architecture based on a statistical model integrating pixel correlations across spatial scales and leveraging the spatial symmetries of the nuisance component, improving both detection and characterization.

* Statistically reliable detection scores and estimates with uncertainty quantification, essential for astrophysics.

* Joint processing of 4-D datasets incorporating the ASDI forward model, shown to significantly boost detection.

## 2. Related work

Various methods have been developed to isolate faint planetary signals from stellar nuisance, broadly categorized as subtraction-, statistical-, and learning-based models [50].

**Subtraction-based methods.** Subtraction-based methods are among the earliest and most common techniques for exoplanet detection in direct imaging, aiming to remove quasi-static speckles that obscure faint signals. The cA(S)DI algorithms [38, 41] subtract a reference model by averaging frames and stacking residuals, enhancing the rotating exoplanet signal. TLOCI [44] and its variants (e.g., [32, 68]) optimize linear combinations of images to model speckles, while KLIP/PCA [3, 60] employs principal component analysis for low-rank subspace projection. Other approaches, such as non-negative matrix factorization [33, 49, 51, 52] and LLSG [34], decompose data into components to isolate sources, while the RSM algorithm [17–19] combines outputs from multiple methods to reduce individual biases. However, these methods often lack statistically rigorous outputs, such as interpretable detection scores and unbiased flux estimates, as they rely on heuristic image combinations rather than fully end-to-end models.

**Statistical approaches.** To address the previous limitations, statistical methods model the nuisance using probabilistic frameworks. The ANDROMEDA [6] and FMMF [54] approaches rely on matched filtering, simplifying the nuisance as uncorrelated Gaussian noise. Similarly, SNAP [64] estimates both the nuisance and exoplanet components jointly through maximum likelihood under the same assumption. The PACO algorithm [25, 26] improves on this by extending beyond the white noise assumption. It uses a patch-based statistical framework to model the spatial and spectral covariances of the nuisance, effectively capturing the local structure of speckles. This approach draws on techniques from computer vision, such as denoising, restoration, super-resolution, collaborative filtering [15], sparse coding [1, 40], or mixture models [71, 72].
**Learning-based approaches.** Machine learning methods are increasingly used in high-contrast imaging, inspired by advances in fields like photography and biomedical imaging. Early applications in exoplanet detection, such as the SVM model by [23], leverage structured high-contrast data, while the (NA)-SODINN algorithms [8, 35] frames detection as a binary classification problem, using KLIP/PCA-processed patches with random forests or CNNs. Despite their effectiveness, SODINN struggles with high false alarm rates and complex hyperparameter tuning [7]. Generative adversarial networks (GANs) have also been applied to simulate nuisance patterns for training deep learning models [70]. Approaches like TRAP [55] and HSR [31] handle unmixing via regularized regression, modeling nuisance evolution with signal-free reference data. Deep PACO [27, 28] combines statistical modeling with deep learning, leveraging PACO's nuisance statistical model and a CNN to detect exoplanets and refine residual mismatches.

All these algorithms are observation-dependent, building nuisance models directly from the dataset of interest. This dependency hinders detection near the host star due to (i) significant temporal fluctuations in the nuisance and (ii) residual self-subtraction, where part of the exoplanet signal is mistakenly removed. Recent observation-independent approaches address these limitations by using archival survey data to model the nuisance. Super-RDI [56] extends KLIP/PCA for large observational databases, while ConStruct [69] uses an auto-encoder to learn typical speckle patterns, and [11] employs a discriminative nuisance model. MODEL&CO [5] improves deep PACO with a deep statistically-modeled nuisance framework. Our proposed approach, ExoMILD (Exoplanet imaging by MIxture of Learnable Distributions), belongs to this new category.

## 3. Image formation models

In direct imaging, a stellar system is observed through an optical system that includes the atmosphere, telescope, and scientific instrument. The response of this system to a point source (e.g., an exoplanet) is defined by its point-spread function off the optical axis (off-axis PSF), describing how light from the source is distributed across the sensor. However, as noted in Sec. 1, residual aberrations uncorrected by adaptive optics produce a quasi-static, structured speckle pattern. To mitigate speckles through numerical processing, several observational strategies are used, as detailed next.

### 3.1. Angular Differential Imaging (ADI)

In ADI, the field derotator of the telescope is turned off causing the field of view (including any exoplanets) to rotate around the target star due to Earth's rotation. Meanwhile, the optical system remains stationary, causing the speckles to remain fixed. This distinction helps separating exoplanet signals from speckle noise. Formally, let $\mathbf{y}$ in $\mathbb{R}^{T \times H \times W}$

denote the sequence of measurements from the observation, where $T$ is the number of exposures and $H, W$ represent the pixel dimensions of each exposure. Each frame $\mathbf{y}_t$ in $\mathbb{R}^{H \times W}$ can then be represented as:

$$ \mathbf{y}_t = \mathbf{s}_t + \sum_{k=0}^{K-1} \alpha^{(k)} \mathbf{h}(\mathbf{x}_t^{(k)}) + \mathbf{\epsilon}_t, \quad (1) $$

where $\mathbf{s}_t$ in $\mathbb{R}^{H \times W}$ is the speckles component, $K$ is the (unknown) number of exoplanets, $\mathbf{x}_t^{(k)}$ is the position of exoplanet $k$ at time $t$, $\mathbf{h}(\mathbf{x})$ in $\mathbb{R}^{H \times W}$ is the PSF model centered on position $\mathbf{x}$ in $\mathbb{R}^2$, $\alpha^{(k)}$ in $\mathbb{R}_+$ is the flux of the exoplanet, and $\mathbf{\epsilon}_t$ in $\mathbb{R}^{H \times W}$ is stochastic noise. The position of exoplanet $k$ at time $t$ can be written as $\mathbf{x}_t^{(k)} = r(\mathbf{x}_0^{(k)}, \phi_t)$, where $\mathbf{x}_0^{(k)}$ in $\mathbb{R}^2$ is its initial position on frame $\mathbf{y}_0$, $\phi_t$ in $\mathbb{R}$ is the predictable *parallactic angle* at time $t$, defined as the cumulated apparent rotation (induced by Earth's rotation) since the start of the sequence, and $r: \mathbb{R}^2 \times \mathbb{R} \to \mathbb{R}^2$ defines a rotation trajectory centered on the star.

### 3.2. Angular Spectral Differential Imaging (ASDI)

The speckles pattern induced by the star can be understood as the PSF on the optical axis (on-axis PSF) of the optical system. According to diffraction laws, this pattern scales homothetically with wavelength, creating a chromatic effect similar to chromatic aberrations in photography. In ADI, this spectral scaling helps further disentangle exoplanet signals from speckles. In this case, each observation is denoted as $\mathbf{y}$ in $\mathbb{R}^{C \times T \times H \times W}$, where $C$ is the number of channels. The forward model (1) becomes

$$ \mathbf{y}_{c,t} = \beta_c \mathbf{D}_{\lambda_c/\lambda_0} \mathbf{s}_{0,t} + \sum_{k=0}^{K-1} \alpha_c^{(k)} \mathbf{h}_c(\mathbf{x}_t^{(k)}) + \mathbf{\epsilon}_{c,t}, \quad (2) $$

where $\beta_c$ in $\mathbb{R}$ is the amplitude of the speckles in channel $c$, $\mathbf{D}_{\lambda_c/\lambda_0}$ is the homothety operator aligning channel 0 on channel $c$ with a dilation coefficient defined as the ratio of their wavelengths, $\mathbf{s}_{0,t}$ in $\mathbb{R}^{H \times W}$ is speckles pattern at wavelength $\lambda_0$, $\mathbf{\epsilon}_{c,t}$ in $\mathbb{R}^{H \times W}$ is the additive thermal noise. In this setting, the exoplanet flux $\alpha_c$ and the off-axis PSF $\mathbf{h}_c$ now depend on channel $c$.

## 4. Proposed method

### 4.1. Convolutional statistical model

**Local Gaussian model of speckles.** Building on the PACO algorithm [25, 26], we propose to capture local spatial correlations between pixels of the nuisance term. We denote by $\mathcal{G}_p$ the grid of spatial pixels, such that $|\mathcal{G}_p| = HW$ and $\forall i \in [0, HW-1]$, $\mathbf{x}_i$ in $\mathbb{R}^2$. Then, we model the statistical distribution of collections of patches positioned on a spatial grid $\mathcal{G}_d$, such that $\mathcal{G}_d \subset \mathcal{G}_p$, with $M := |\mathcal{G}_d|$. We denote by $\mathbf{y}_t^{(j)}$ in $\mathbb{R}^p$, $\forall (j, t) \in [0, M-1] \times [0, T-1]$ the observed patch in spatial location $j$ at time $t$. In absence of exoplanet, i.e., when only the nuisance component
<figure><img src="image_3.png" /><figcaption>Figure B. Workflow of the proposed method: it exploits both the spectral behavior of speckles and the apparent motion of exoplanets to disentangle the exoplanet signal from the nuisance component in the observations $y$. To achieve this, local patches of the nuisance are modeled as Gaussian distributions, leveraging problem symmetries and incorporating multiple scales. These patches are fed to our convolutional statistical model, and combined to form a detection map. Additionally, a learned object prior, represented by a UNet $f_v$, is introduced to denoise this detection map produced by the statistical model. This approach results in an end-to-end learnable architecture.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure C. Proposed convolutional statistical model: spectrally aligned speckles patches, indexed by $j$, with dimensions $N_p$ and $CT$ samples, are first linearly projected into a lower-dimensional space of size $m$. In this space, the parameters of the Gaussian distribution $\hat{m}_j$ and $\hat{C}_j$ are estimated and subsequently combined with the PSF $\hat{h}$ to compute the terms $\mathbf{a}^{(j)}$ and $\mathbf{b}^{(j)}$. As detailed in Appendix A.3, the efficient computation of $\mathbf{a}^{(j)}$ relies on the Cholesky decomposition of the precision matrix $\hat{C}_j^{-1}$.</figcaption></figure>

is present, each collection of patches $\{\mathbf{y}_t^{(j)}\}_{t=0:T-1}$ is modeled by a multivariate Gaussian:

$$ \forall (j, t), \mathbf{y}_t^{(j)} \sim \mathcal{N}(\mathbf{m}_j, \mathbf{C}_j), \quad (3) $$

with $\mathbf{m}_j$ in $\mathbb{R}^p$ and $\mathbf{C}_j$ in $\mathbb{R}^{p \times p}$ the mean and covariance of the Gaussian distribution. These parameters are estimated by maximum likelihood and shrinkage of the covariance matrix. They are denoted as $\hat{\mathbf{m}}_j$ and $\hat{\mathbf{C}}_j$ in the following. Additional details are presented in Appendix A.1.

**Detection criterion at a given position.** This local statistical model estimates the time-invariant flux $\hat{\alpha}$ of an exoplanet at position $\mathbf{x}_0$ in $\mathbb{R}^2$ in the first frame by maximizing the global likelihood:

$$ \hat{\alpha} = \arg\max_{\alpha} \ell(\alpha, \mathbf{x}_0). \quad (4) $$

Assuming that collections of patches on the trajectory of an exoplanet are independent, the likelihood becomes:

$$ \ell(\alpha, \mathbf{x}_0) = \prod_t \mathbb{P} \left( \mathbf{y}_t^{(i_t)} - \alpha \mathbf{h}^{(i_t)}(\mathbf{x}_t) \middle| \hat{\mathbf{m}}_{i_t}, \hat{\mathbf{C}}_{i_t} \right), \quad (5) $$

where $\mathbf{x}_t = r(\mathbf{x}_0, \phi_t)$ is the predictive position of the exoplanet at time $t$, $i_t$ the index of the patch centered on position $\mathbf{x}_t$ in $\mathcal{G}_p$, and $\mathbf{h}^{(i_t)}(\mathbf{x}_t)$ in $\mathbb{R}^p$ the patch $i_t$ of the off-axis PSF model centered on $\mathbf{x}_t$. We propose to consider the modified likelihood:

$$ \ell(\alpha, \mathbf{x}_0) = \prod_t \prod_{j \in S(\mathbf{x}_t)} \mathbb{P} \left( \mathbf{y}_t^{(j)} - \alpha \mathbf{h}^{(j)}(\mathbf{x}_t) \middle| \hat{\mathbf{m}}_j, \hat{\mathbf{C}}_j \right)^{w_j}, \quad (6) $$

where $S(\mathbf{x}_t)$ represents the subset of modeled distributions for patches around location $\mathbf{x}_t$, with each distribution $j$ weighted by $w_j$ in $\mathbb{R}^+$, where $\sum w_j = 1$. Compared to (5), this approach is more robust, as it aggregates overlapping patch contributions at each time-step, modeling the
nuisance component as a convolutional process influenced by multiple, overlapping noise sources. It is also more computationally efficient, as it does not require a dense grid $G_d$ for patch collections, thereby reducing the parameter estimation load. This efficiency is essential for incorporating the model into an end-to-end learning framework (see Sec. 4.3). Finally, this flexible formulation can accommodate additional correlation types (see Sec. 4.2). Solving problem (6) leads to the flux estimator $\hat{\alpha}$ and its standard deviation $\hat{\sigma}_\alpha$:

$$ \hat{\alpha} = \frac{\sum_t b_t(\mathbf{x}_t)}{\sum_t a_t(\mathbf{x}_t)}, \quad \hat{\sigma}_\alpha = \frac{1}{\sqrt{\sum_t a(\mathbf{x}_t)}}, \quad (7) $$

where

$$ b_t(\mathbf{x}_t) = \sum_{j \in S(\mathbf{x}_t)} w_j \mathbf{h}^{(j)}(\mathbf{x}_t)^\top \hat{\mathbf{C}}_j^{-1} (\mathbf{y}_t^{(j)} - \hat{\mathbf{m}}_j), \quad (8) $$

$$ a(\mathbf{x}_t) = \sum_{j \in S(\mathbf{x}_t)} w_j \mathbf{h}^{(j)}(\mathbf{x}_t)^\top \hat{\mathbf{C}}_j^{-1} \mathbf{h}^{(j)}(\mathbf{x}_t). \quad (9) $$

To assess the probability of an exoplanet at $\mathbf{x}_0$, we use the generalized likelihood ratio test (GLRT) to statistically evaluate the parameter $\alpha$. Under the null hypothesis $H_0$ where $\alpha = 0$ (indicating no exoplanet), the statistics

$$ \hat{\gamma} = \hat{\alpha}/\hat{\sigma}_\alpha = (\sum_t b_t(\mathbf{x}_t)) / (\sqrt{\sum_t a(\mathbf{x}_t)}), \quad (10) $$

is controlled and follows a Gaussian distribution $N(0, 1)$. When evaluating the probability of presence of an exoplanet, we test against the alternative hypothesis $H_1: \alpha > 0$. The statistics $\hat{\gamma}$ can be directly mapped to a probability of detection. It can also be interpreted as the output signal-to-noise ratio through the statistical model, and in practice, this is the detection score used by astronomers [63]. To obtain the dense counterparts $\hat{\gamma}, \hat{\alpha}, \hat{\sigma}$ in $\mathbb{R}^{H \times W}$, we adopt a fast approximation similar to [25] detailed in Appendix A.2.

**Iterative procedure for characterization.** Astrometry is the task of precisely estimating the flux $\alpha$ and the sub-pixel position $\mathbf{x}_0$ of an exoplanet. This can be achieved by jointly optimizing the likelihood defined in Eq. (6). Denoting by $\mathbf{z} = [\alpha, \mathbf{x}_0]$ in $\mathbb{R}^3$, the gradient $\mathbf{g}$ and Hessian matrix $\mathbf{H}$ of the neg-log-likelihood are decomposed as follows:

$$ \mathbf{g}(\mathbf{z}) = \sum_t \sum_{j \in S(\mathbf{x}_t)} w_j \mathbf{g}_j(\mathbf{z}), \quad \mathbf{H}(\mathbf{z}) = \sum_t \sum_{j \in S(\mathbf{x}_t)} w_j \mathbf{H}_j(\mathbf{z}) $$

where $\mathbf{g}_j$ and $\mathbf{H}_j$ represent the gradient and Hessian operators for each distribution $j$, respectively. These operators are responsible for updating the statistical parameters $\mathbf{m}_j$ and $\mathbf{C}_j$, which are initially biased due to the presence of the exoplanet. Further details are provided in Sec. A.4. Each iteration indexed by $l$ writes as:

$$ \mathbf{z}^{(l+1)} = \mathbf{z}^{(l)} - \mathbf{H}(\mathbf{z}^{(l)})^{-1} \mathbf{g}(\mathbf{z}^{(l)}). \quad (11) $$

## 4.2. Extensions of the statistical model: mixture of distributions

We now build on the convolutional statistical model presented in Sec. 4.1 to introduce a multi-scale statistical model of the nuisance component.

**Linear projection.** Instead of modeling correlations between the pixels of a patch, we propose to model correlations between projected linear features:

$$ \forall (j, t), \quad \mathbf{A} \mathbf{y}_t^{(j)} \sim N(\mathbf{m}_j, \mathbf{C}_j), \quad (12) $$

where $\mathbf{A}$ in $\mathbb{R}^{m \times p}$ is the projection matrix and $m \le p$. This is a generalization of the model presented in Sec. 4.1, for which $\mathbf{A} = \mathbf{I}_p$. The terms $b_t$ and $a$ introduced in Eqs. (26)-(27) become:

$$ b_{i,t} = \sum_{j \in S(\mathbf{x}_i)} w_j \mathbf{h}_{ji}^\top \mathbf{A}_j^\top \hat{\mathbf{C}}_j^{-1} (\mathbf{A}_j \mathbf{y}_t^{(j)} - \hat{\mathbf{m}}_j), \quad (13) $$

$$ a_i = \sum_{j \in S(\mathbf{x}_i)} w_j \mathbf{h}_{ji}^\top \mathbf{A}_j^\top \hat{\mathbf{C}}_j^{-1} \mathbf{A}_j \mathbf{h}_{ji}. \quad (14) $$

The learnable projection $\mathbf{A}$ decorrelates the feature space, where the statistical distribution is defined, from the pixel space, enabling effective handling of long-range correlations.

**Multi-scale approach.** So far, we have shown how our model can combine multiple neighboring local distributions through weighted averaging of the terms $b_t$ and $a$. Extending this approach, we propose to integrate distributions across multiple spatial scales, by varying the patch sizes $p$. We denote $\mathcal{P}$ the set of patch sizes chosen. We expand the set $S(\mathbf{x}_t)$ of modeled distributions to include different scales, such that:

$$ S(\mathbf{x}_t) = \bigcup_{p \in \mathcal{P}} S_p(\mathbf{x}_t), \quad (15) $$

where $S_p(\mathbf{x}_t)$ is the subset of patches with patch size $p$ containing $\mathbf{x}_t$. The linear projection is particularly useful in controlling the dimensionality of the feature space for larger patches by setting $m < p$, thereby maintaining computational efficiency.

**Leveraging rotational symmetries.** In direct imaging, the speckle pattern exhibits approximate central and rotational symmetries, see Fig. A. These symmetries arise from near-isotropic distortions of the wavefront and from the optical system's near-circular symmetry. Symmetries of the speckles can be theoretically understood through a series expansion of the diffraction pattern, using increasing powers of the Fourier transform of the residual phase error [48, 53, 59]. Leveraging these symmetries is essential for constructing a more robust model of the speckles, particularly to mitigate the effects of self-subtraction. This effect is a well-known challenge in direct imaging that arises
when limited parallactic rotation results in the contamination of speckle parameters by the exoplanet signal, ultimately reducing detection sensitivity [46]. To incorporate these symmetries, we propose modeling the joint distribution of patches extracted from the same location after rotating the observation by $2\pi/N$, also known as N-fold rotational symmetry:

$$ \forall (j, t), \quad \mathbf{A}[\mathbf{R}_{2\pi n/N}(\mathbf{y}_t)^{(j)}]_{n=0:N-1} \sim \mathcal{N}(\mathbf{m}_j, \mathbf{C}_j), \quad (16) $$

where A is in $\mathbb{R}^{m \times N_p}$. The parameters of this joint distribution are less subject to self-subtraction as it is very unlikely to observe exoplanets simultaneously in all jointly modeled patches. We use a mixture model with $N = 1, 2, 4$.

**Joint multi-spectral modeling.** In ASDI, the spectral channels are related by a homothety, as given in Eq. (2). We propose to leverage this relationship, and model:

$$ \forall (j, c, t), \quad \mathbf{A} \beta_{c,j}^{-1} \mathbf{D}_{\lambda_0/\lambda_c}(\mathbf{y}_{c,t})^{(j)} \sim \mathcal{N}(\mathbf{m}_j, \mathbf{C}_j). \quad (17) $$

In this setting, the estimators of statistical parameters $\hat{\mathbf{m}}_j$ and $\hat{\mathbf{C}}_j$ are less prone to self-subtraction, as the spectral diversity reduces the ambiguity between speckles and exoplanet signals. The estimator of the local amplitude $\hat{\beta}_{c,j}$ is computed as the standard deviation of pixel values across $\{\mathbf{D}_{\lambda_0/\lambda_c}(\mathbf{y}_{c,t})^{(j)}\}_t$.

### 4.3. End-to-end trainable approach

**Problem statement.** We propose to combine our statistical model on the nuisance component with a learnable prior on the exoplanet signals. This approach can be formulated as an optimization problem:

$$ \hat{\alpha} = \arg\min_{\alpha \in \mathbb{R}^{H \times W}} \varphi_{\theta}(\alpha, \chi(\mathbf{y})) + \psi_{\nu}(\alpha), \quad (18) $$

where $\varphi$ is the data-fitting term, corresponding to the negative log-likelihood of our statistical model of the nuisance component described in Sec. 4.1, and $\psi$ is a prior on exoplanet signals. We denote by $\theta, \nu$ the learnable parameters associated with these terms, and by $\chi(\mathbf{y})$ the parameters estimated for each observation. In practice $\chi(\mathbf{y}) = \{\hat{\mathbf{m}}_j, \hat{\mathbf{C}}_j, \hat{\beta}_{c,j}\}, \theta = \{\mathbf{A}_j, w_j\}_j$, and $\nu$ corresponds to the parameters of the neural network implementing $\psi$.

**Detection by denoising.** We propose a two-step approach to obtain the detection score corresponding to Eq. (18). First, we compute the detection score corresponding to the statistical model $\varphi$, which admits a fast approximation denoted by $\gamma \in \mathbb{R}^{H \times W}$ in the following, and given by Eq. (28). We recall that under the null hypothesis, i.e., when no exoplanet is present, the elements of $\gamma$ follow a Gaussian distribution $\mathcal{N}(0, 1)$. Therefore, extracting the signals of

exoplanets corresponds to denoising $\gamma$ to remove this background noise. We propose to achieve this step using a neural network, such that:

$$ \tilde{\gamma} = f_{\nu}(\hat{\gamma}), \quad (19) $$

where $f_{\nu}$ is the denoiser implemented by the neural network, and $\tilde{\gamma}$ the final detection score.

**Training objective.** We suppose that the denoised detection score can be decomposed similarly to the GLRT form provided in Eq. (10) for the statistical model. This leads to $\tilde{\gamma} = \tilde{\alpha}/\tilde{\sigma}_e$, where $\tilde{\alpha}, \tilde{\sigma}$ in $\mathbb{R}^{H \times W}$ are the estimated flux of the exoplanet for each pixel, and the standard deviation denoting the uncertainty associated with it. Additionally, we assume that the uncertainty $\tilde{\sigma}$ remains unaffected by the neural network, as it is primarily driven by the high variability of the speckles already captured by the statistical model. Consequently, both components can be expressed as:

$$ \tilde{\alpha} = f_{\nu}(\hat{\alpha}/\hat{\sigma}) \times \hat{\sigma}, \quad \tilde{\sigma} = \hat{\sigma}. \quad (20) $$

This formulation yields a pixel-wise Gaussian distribution: $\mathcal{N}(\tilde{\alpha}, \tilde{\sigma})$. The learnable parameters $\theta$ and $\nu$ of our model are optimized by minimizing the negative log-likelihood between estimates $(\tilde{\alpha}, \tilde{\sigma})$ and ground truth $\alpha_{gt}$:

$$ \mathcal{L}(\tilde{\alpha}, \tilde{\sigma}, \alpha_{gt}) = 0.5 (\tilde{\alpha} - \alpha_{gt})^2 / \tilde{\sigma}^2 + \log \tilde{\sigma}. \quad (21) $$

**Model ensembling.** To improve robustness, we combine outputs from multiple trained models. Specifically, given outputs $\tilde{\alpha}$ and $\tilde{\sigma}$, we can define $\tilde{a}$ and $\tilde{b}$ such that $\tilde{\sigma} = 1/\sqrt{\tilde{a}}$ and $\tilde{\alpha} = \tilde{b}/\tilde{a}$. For $Q$ models indexed by $q$, the outputs are aggregated as follows:

$$ \tilde{a}_Q = \sum_q \tilde{a}_q/Q, \quad \tilde{b}_Q = \sum_q \tilde{b}_q/Q. \quad (22) $$

The combined detection score, $\tilde{\gamma}_Q = \tilde{b}_Q/\sqrt{\tilde{a}_Q}$, represents a weighted average of the individual detection scores.

**Calibration.** The trained model needs to be calibrated in order to relate the detection score to a probability of false alarm. We follow the procedure outlined in [5], and rely on a separate calibration dataset to estimate the cumulative distribution function of $\tilde{\gamma}$ under the null hypothesis. Additional details are provided in the Supplementary Material.

## 5. Experiments

### 5.1. Data, algorithms and evaluation protocol

**Datasets.** We use SPHERE datasets, a cutting-edge exoplanet finder instrument at the VLT [4]. Raw observations were sourced from the public data archive of the European Southern Observatory and calibrated with public tools [21, 47] from the High-Contrast Data Center. This results in
<table><caption>Table A. Impact of multi-scale and N-fold rotational symmetries on detection performance (AUC) of our statistical model.</caption><thead><tr><th></th><th>Spatial scales</th><th>N=1</th><th>N=1,2</th><th>N=1, 2, 4</th></tr></thead><tbody><tr><td rowspan="4">ADI</td><td>8 × 8</td><td>0.554 ± 0.005</td><td>0.571 ± 0.005</td><td>0.575 ± 0.005</td></tr><tr><td>+16 × 16</td><td>0.561 ± 0.005</td><td>0.573 ± 0.004</td><td>0.579 ± 0.005</td></tr><tr><td>+32 × 32</td><td>0.567 ± 0.005</td><td>0.575 ± 0.005</td><td>0.580 ± 0.005</td></tr><tr><td>+64 × 64</td><td>0.569 ± 0.005</td><td>0.577 ± 0.005</td><td><strong>0.581 ± 0.005</strong></td></tr><tr><td rowspan="4">ASDI</td><td>8 × 8</td><td>0.713 ± 0.005</td><td>0.719 ± 0.005</td><td>0.723 ± 0.005</td></tr><tr><td>+16 × 16</td><td>0.720 ± 0.005</td><td>0.723 ± 0.005</td><td>0.725 ± 0.005</td></tr><tr><td>+32 × 32</td><td>0.720 ± 0.005</td><td>0.725 ± 0.005</td><td><strong>0.726 ± 0.005</strong></td></tr><tr><td>+64 × 64</td><td>0.720 ± 0.005</td><td>0.725 ± 0.004</td><td><strong>0.726 ± 0.004</strong></td></tr></tbody></table>

science-ready 4-D datasets $y$ (with $L = 2$ spectral channels, $T \in [15, 300]$ temporal frames, and $H \times W = 1024^2$ pixels per image), along with the off-axis PSF $h$, wavelengths $\lambda_c$, and parallactic angles $\phi_t$ for algorithm input. Many algorithms have achieved optimal detection sensitivity far from the star, limited by photon noise [5, 28]. Thus, our analysis focuses on a smaller, star-centered region of $H = W = 256$ pixels, where detection sensitivity can still be significantly improved [5, 28].

**Training procedure.** For training, we use 220 observations from the SHINE-F150 large survey of SPHERE [39]. For testing, we select 8 datasets representing typical diversity in observing conditions and instrumental settings (e.g., parallactic rotation amplitude $\Delta_{\phi} = |\phi_{T-1} - \phi_0|$). Five test datasets (on stars HD 159911, HD 216803, HD 206860, HD 188228, HD 102647) are used for benchmarking against state of the art methods and conducting a model ablation analysis. To evaluate performance, we simulate synthetic faint point-like sources mimicking exoplanet signatures and inject them into real data. This simulation procedure is common practice in direct imaging to ground the actual performance of algorithms because it is very realistic [7, 8, 16, 17, 28, 35]. Any real source indeed takes the form of the off-axis PSF, which is measured immediately before and after the main observation sequence by offsetting the star from the coronagraph. This simulation procedure is essential due to the lack of ground truth and the limited number of exoplanets (only a few dozen) detected by direct imaging to date. In addition, we rely on simulated sources for training our deep model, which is also a standard practice for learning-based techniques. Three datasets of star HR 8799, hosting three known exoplanets in the field of view [42, 43], are also used as full real data.

**Baselines.** For detection, we benchmark the proposed algorithms against methods from different classes described in Sect. 2. Selection criteria are (i) code availability (often limited in direct imaging), (ii) relevance, and (iii) widespread use. In this context, we include the cADI [41], KLIP/PCA [60] subtraction-based methods, as these are implemented in most data processing pipelines [13, 21, 62],

<table><caption>Table B. Comparison of statistical models for flux estimation.</caption><thead><tr><th>Method</th><th>flux error (ARE)</th><th>position error (RMSE)</th></tr></thead><tbody><tr><td>PACO</td><td>0.56</td><td>0.21</td></tr><tr><td>Proposed</td><td><strong>0.51</strong></td><td><strong>0.11</strong></td></tr></tbody></table>

have been responsible for detecting nearly all imaged exoplanets –including the most recent discoveries [45, 67]– and remain heavily utilized by astronomers. Concerning statistical methods, we focus on PACO [25, 26] that uniquely accounts for data correlations. PACO has consistently outperformed KLIP/PCA in large observational surveys [12, 22] and achieved state of the art performance on SPHERE data in a community benchmark, surpassing various subtraction-based, statistical, and learning-based approaches [7]. We also evaluate the proposed approach against MODEL&CO [5], a recent hybrid method that shown to outperform cADI, KLIP/PCA, and (deep) PACO [5]. For cADI and KLIP/PCA, we use the VIP Python package [13, 33], fine-tuning parameters (e.g., PCA modes) to optimize detection scores guided by the ground truth. PACO and MODEL&CO were processed by their authors using data-driven settings. For flux estimation, we compare only with PACO due to computational constraints. We conduct all analyses in both ADI and ASDI modes (where supported) to assess the gains from joint spectral processing.

**Metrics.** The detection metric used is the area under the receiver operating characteristic curve (AUC), representing the true positive rate against the false discovery rate obtained by varying the detection threshold. Higher AUC values indicate better performance. This standard metric in direct imaging [7, 25, 34] captures the precision-recall trade-off and allows fair algorithm comparisons, as a common threshold does not ensure consistent false alarm rates due to the lack of statistical grounding in some detection maps, see Sec. 1. The primary metric is the absolute relative error (ARE) between the ground truth and estimated flux, with lower values indicating better performance. We also report the root mean square error (RMSE) for sub-pixel localization of exoplanets.

## 5.2. Quantitative and qualitative evaluations

**Statistical model.** We evaluate the impact of leveraging multiple spatial scales and symmetries in our statistical model by testing its detection performance across configurations listed in Table A. Results show that using both symmetries and scales is crucial in ADI, where statistical parameters suffer from self-subtraction without the robustness of added spectral diversity. We then evaluate the statistical model’s performance in estimating flux and sub-pixel position (regression tasks) using the optimization process from
<table><caption>Table C. Comparative detection scores (AUC) in A(S)DI modes. Dataset names (stars) and parallactic amplitude $\Delta_{\phi}$ are reported on top. ASDI mode (which is better than ADI) is not supported by MODEL&CO.</caption><thead><tr><th>Modality</th><th>Method</th><th>HD159911 (54°)</th><th>HD216803 (23°)</th><th>HD206860 (11°)</th><th>HD188228 (6°)</th><th>HD102647 (2°)</th><th>average AUC</th></tr></thead><tbody><tr><td rowspan="5">ADI</td><td>cADI</td><td>0.288 ± 0.014</td><td>0.422 ± 0.007</td><td>0.489 ± 0.010</td><td>0.303 ± 0.014</td><td>0.343 ± 0.009</td><td>0.369 ± 0.005</td></tr><tr><td>PCA</td><td>0.634 ± 0.010</td><td>0.643 ± 0.011</td><td>0.505 ± 0.009</td><td>0.392 ± 0.011</td><td>0.218 ± 0.011</td><td>0.478 ± 0.005</td></tr><tr><td>PACO</td><td>0.629 ± 0.006</td><td>0.669 ± 0.012</td><td>0.579 ± 0.015</td><td>0.517 ± 0.015</td><td>0.207 ± 0.012</td><td>0.520 ± 0.006</td></tr><tr><td>MODEL&amp;CO</td><td>0.653 ± 0.010</td><td>0.731 ± 0.011</td><td>0.646 ± 0.012</td><td><strong>0.638</strong> ± 0.014</td><td><strong>0.554</strong> ± 0.009</td><td><strong>0.645</strong> ± 0.005</td></tr><tr><td>Proposed</td><td><strong>0.673</strong> ± 0.009</td><td><strong>0.740</strong> ± 0.010</td><td><strong>0.661</strong> ± 0.012</td><td>0.631 ± 0.013</td><td>0.518 ± 0.012</td><td><strong>0.645</strong> ± 0.005</td></tr><tr><td rowspan="4">ASDI</td><td>cASDI</td><td>0.448 ± 0.009</td><td>0.537 ± 0.013</td><td>0.451 ± 0.007</td><td>0.352 ± 0.017</td><td>0.294 ± 0.012</td><td>0.417 ± 0.006</td></tr><tr><td>PCA</td><td>0.694 ± 0.013</td><td>0.696 ± 0.007</td><td>0.552 ± 0.011</td><td>0.398 ± 0.011</td><td>0.236 ± 0.009</td><td>0.515 ± 0.005</td></tr><tr><td>PACO</td><td>0.698 ± 0.015</td><td>0.768 ± 0.008</td><td>0.710 ± 0.014</td><td>0.700 ± 0.009</td><td>0.589 ± 0.014</td><td>0.693 ± 0.005</td></tr><tr><td>Proposed</td><td><strong>0.731</strong> ± 0.009</td><td><strong>0.804</strong> ± 0.013</td><td><strong>0.747</strong> ± 0.010</td><td><strong>0.782</strong> ± 0.008</td><td><strong>0.744</strong> ± 0.006</td><td><strong>0.761</strong> ± 0.004</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Figure D. Detection maps on observations of HD 159911 star with synthetic exoplanets. The (calibrated) detection threshold is equivalent for all methods. The proposed approach here detect 1 additional source compared to the second best method (PACO ASDI).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure E. Detection maps on 3 observations (stacked in false RGB colors) of HR 8799 star. The elliptical arcs depict the estimated (projected) orbits of three known exoplanets, with the detection results shown as red, green and blue dots for the corresponding 2016, 2018 and 2021 observations. Squares are for false alarms identified in [5], Fig. 17.</figcaption></figure>

Sec. 4.1. Table B shows average errors across 1,351 synthetic exoplanets detectable by both considered methods. Our approach outperforming others in almost all cases.

**End-to-end learnable model.** We evaluate detection performance on 5 observations using synthetic exoplanet injections through direct models (1)-(2). For each observation, 100 cubes are generated, totaling 1,000 injected exoplanets. This process is repeated 5 times with different random seeds, and AUC scores are reported in Table C. Examples of representative detection maps are shown in Fig. K. The proposed approach matches or outperforms state of the art algorithms in ADI mode. In ASDI mode, it shows a significant boost in detection sensitivity, consistently surpassing comparative methods.

**Validation on real data.** We compare detection methods using 3 observations, spanning several years, of the HR

8799 star. The stacked detection maps, highlighting exoplanets' orbital motion, are shown in Fig. E. All detection maps use a consistent unit (signal-to-noise score) and dynamic range. The proposed method maximizes detection confidence with no false alarms.

## 6. Conclusion

We propose a novel hybrid approach for exoplanet imaging that combines a multi-scale statistical model with deep learning, capturing spatial correlations in the nuisance component for simultaneous detection and flux estimation. This approach provides statistically grounded detection scores, unbiased estimates, and native uncertainty quantification. Tested on VLT/SPHERE data, it outperforms state-of-the-art techniques, demonstrating efficiency and robustness across varied data qualities. Its versatility and reduced computational complexity make it ideal for large-scale surveys. The approach will be extended to handle higher spectral resolution data. Additionally, the nuisance model can be adapted for reconstructing spatially extended objects like circumstellar disks, namely birthplaces of exoplanets.
## Acknowledgments

This work was supported by the French government under management of Agence Nationale de la Recherche as part of the "France 2030" program, PR[AI]RIE-PSAI projet (reference ANR-23-IACL-0008) and MIAI 3IA Institute (reference ANR-19-P3IA-0003), and by the European Research Council (ERC) under grant agreement 101087696 (APHELEIA project). This work was also supported by the ERC under the European Union's Horizon 2020 research and innovation programme (COBREX: grant agreement 885593), the ANR under the France 2030 program (PEPR Origins, reference ANR-22-EXOR-0016), the French National Programs (PNP and PNPS), and the Action Spécifique Haute Résolution Angulaire (ASHRA) of CNRS/INSU co-funded by CNES. This work was granted access to the HPC resources of IDRIS under the allocation 2022-AD011013643 made by GENCI. JP was supported in part by the Louis Vuitton/ENS chair in artificial intelligence and a Global Distinguished Professorship at the Courant Institute of Mathematical Sciences and the Center for Data Science at New York University.

## References

[1] Michal Aharon, Michael Elad, and Alfred Bruckstein. K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation. *IEEE Transactions on signal processing*, 54(11):4311–4322, 2006. 2

[2] France Allard, Nicole F Allard, Derek Homeier, John Kielkopf, Mark J McCaughrean, and Fernand Spiegelman. K-H2 quasi-molecular absorption detected in the T-dwarf Indi Ba. *Astronomy & Astrophysics*, 474(2):L21–L24, 2007. 1

[3] Adam Amara and Sascha P Quanz. Pynpoint: an image processing package for finding exoplanets. *Monthly Notices of the Royal Astronomical Society*, 427(2):948–955, 2012. 2

[4] Jean-Luc Beuzit, Arthur Vigan, David Mouillet, Kjetil Dohlen, Raffaele Gratton, Anthony Boccaletti, J-F Sauvage, Hans Martin Schmid, Maud Langlois, Cyril Petit, et al. SPHERE: the exoplanet imager for the Very Large Telescope. *Astronomy & Astrophysics*, 631:A155, 2019. 2, 6

[5] Théo Bodrito, Olivier Flasseur, Julien Mairal, Jean Ponce, Maud Langlois, and Anne-Marie Lagrange. Model&co: Exoplanet detection in angular differential imaging by learning across multiple observations. *Monthly Notices of the Royal Astronomical Society*, 534(2):1569–1596, 2024. 3, 6, 7, 8, 2

[6] F Cantalloube, D Mouillet, LM Mugnier, J Milli, Olivier Absil, CA Gomez Gonzalez, G Chauvin, J-L Beuzit, and A Corchia. Direct exoplanet detection and characterization using the andromeda method: Performance on vlt/naco data. *Astronomy & Astrophysics*, 582:A89, 2015. 2

[7] Faustine Cantalloube, Carlos Gomez-Gonzalez, Olivier Absil, Carles Cantero, Regis Bacher, MJ Bonse, Michael Botom, C-H Dahlqvist, Célia Desgrange, Olivier Flasseur, et al. Exoplanet imaging data challenge: benchmarking the various image processing methods for exoplanet detection. In

*Adaptive Optics Systems VII*, pages 1027–1062. SPIE, 2020. 3, 7

[8] Carles Cantero, Olivier Absil, C-H Dahlqvist, and Marc Van Droogenbroeck. Na-sodinn: A deep learning algorithm for exoplanet image detection based on residual noise regimes. *Astronomy & Astrophysics*, 680:A86, 2023. 3, 7

[9] Gilles Chabrier, Isabelle Baraffe, France Allard, and P Hauschildt. Evolutionary models for very low-mass stars and brown dwarfs with dusty atmospheres. *The Astrophysical Journal*, 542(1):464, 2000. 1

[10] Gaël Chauvin. Direct imaging of exoplanets at the era of the extremely large telescopes. *arXiv preprint arXiv:1810.02031*, 2018. 1

[11] Pattana Chintarungruangchai, Guey Jiang, Jun Hashimoto, Yu Komatsu, and Mihoko Konishi. A possible converter to denoise the images of exoplanet candidates through machine learning techniques. *New Astronomy*, 100:101997, 2023. 3

[12] A Chomez, A-M Lagrange, P Delorme, M Langlois, G Chauvin, O Flasseur, J Dallant, F Philipot, S Bergeon, D Albert, et al. Preparation for an unsupervised massive analysis of sphere high-contrast data with paco-optimization and benchmarking on 24 solar-type stars. *Astronomy & Astrophysics*, 675:A205, 2023. 7, 2

[13] Valentin Christiaens, Carlos Gonzalez, Ralf Farkas, Henrik Dahlqvist, Evert Nasedkin, Julien Milli, Olivier Absil, Henry Ngo, Carles Cantero, Alan Rainot, et al. Vip: A python package for high-contrast imaging. *Journal of Open Source Software*, 8, 2023. 7

[14] Thayne Currie, Beth Biller, Anne-Marie Lagrange, Christian Marois, Olivier Guyon, Eric Nielsen, Mickael Bonnefoy, and Robert De Rosa. Direct imaging and spectroscopy of extra-solar planets. *arXiv preprint arXiv:2205.05696*, 2022. 1

[15] Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian. Image denoising by sparse 3-d transform-domain collaborative filtering. *IEEE Transactions on image processing*, 16(8):2080–2095, 2007. 2

[16] Hazan Daglayan, Simon Vary, Faustine Cantalloube, P-A Absil, and Olivier Absil. Likelihood ratio map for direct exoplanet detection. In *2022 IEEE 5th International Conference on Image Processing Applications and Systems (IPAS)*, pages 1–5. IEEE, 2022. 7

[17] C-H Dahlqvist, Faustine Cantalloube, and Olivier Absil. Regime-switching model detection map for direct exoplanet detection in adi sequences. *Astronomy & Astrophysics*, 633: A95, 2020. 2, 7

[18] C-H Dahlqvist, Faustine Cantalloube, and Olivier Absil. Auto-rsm: An automated parameter-selection algorithm for the rsm map exoplanet detection algorithm. *Astronomy & Astrophysics*, 656:A54, 2021. 2

[19] C-H Dahlqvist, Gilles Louppe, and Olivier Absil. Improving the rsm map exoplanet detection algorithm-psf forward modelling and optimal selection of psf subtraction techniques. *Astronomy & Astrophysics*, 646:A49, 2021. 2

[20] Richard Davies and Markus Kasper. Adaptive optics for astronomy. *Annual Review of Astronomy and Astrophysics*, 50: 305–351, 2012. 2

[21] Ph Delorme, Nadège Meunier, D Albert, Eric Lagadec, H Le Coroller, R Galicher, D Mouillet, A Boccaletti, DINO Mesa,
J-C Meunier, et al. The SPHERE Data Center: a reference for high contrast imaging processing. *arXiv preprint arXiv:1712.06948*, 2017. 6, 7

[22] P Delorme, A Chomez, V Squicciarini, M Janson, O Flasseur, O Schib, R Gratton, AM Lagrange, M Langlois, L Mayer, et al. Giant planets population around b stars from the first part of the beast survey. *arXiv preprint arXiv:2409.18793*, 2024. 7

[23] Rob Fergus, David W Hogg, Rebecca Oppenheimer, Douglas Brenner, and Laurent Pueyo. S4: A spatial-spectral model for speckle suppression. *The Astrophysical Journal*, 794(2):161, 2014. 3

[24] Michael P Fitzgerald and James R Graham. Speckle statistics in adaptively corrected images. *The Astrophysical Journal*, 637(1):541, 2006. 2

[25] Olivier Flasseur, Loïc Denis, Éric Thiébaut, and Maud Langlois. Exoplanet detection in angular differential imaging by statistical learning of the nonstationary patch covariances-the paco algorithm. *Astronomy & Astrophysics*, 618:A138, 2018. 2, 3, 5, 7

[26] Olivier Flasseur, Loïc Denis, Éric Thiébaut, and Maud Langlois. Paco asdi: an algorithm for exoplanet detection and characterization in direct imaging with integral field spectrographs. *Astronomy & Astrophysics*, 637:A9, 2020. 2, 3, 7

[27] Olivier Flasseur, Théo Bodrito, Julien Mairal, Jean Ponce, Maud Langlois, and Anne-Marie Lagrangev. Combining multi-spectral data with statistical and deep-learning models for improved exoplanet detection in direct imaging at high contrast. In *2023 31st European Signal Processing Conference (EUSIPCO)*, pages 1723–1727. IEEE, 2023. 3

[28] Olivier Flasseur, Théo Bodrito, Julien Mairal, Jean Ponce, Maud Langlois, and Anne-Marie Lagrange. deep paco: Combining statistical models with deep learning for exoplanet detection and characterization in direct imaging at high contrast. *Monthly Notices of the Royal Astronomical Society*, 527(1):1534–1562, 2024. 3, 7

[29] Olivier Flasseur, Eric Thiébaut, Loïc Denis, and Maud Langlois. Shrinkage mmse estimators of covariances beyond the zero-mean and stationary variance assumptions. In *2024 32nd European Signal Processing Conference (EUSIPCO)*, pages 2727–2731, 2024. 1

[30] Katherine B Follette. An introduction to high contrast differential imaging of exoplanets and disks. *Publications of the Astronomical Society of the Pacific*, 135(1051):093001, 2023. 1

[31] Timothy D Gebhard, Markus J Bonse, Sascha P Quanz, and Bernhard Schölkopf. Half-sibling regression meets exoplanet imaging: Psf modeling and subtraction using a flexible, domain knowledge-driven, causal framework. *Astronomy & Astrophysics*, 666:A9, 2022. 3

[32] Benjamin L Gerard and Christian Marois. Planet detection down to a few λ/d: an rsdi/tloci approach to psf subtraction. In *Adaptive Optics Systems V*, pages 1544–1556. SPIE, 2016. 2

[33] Carlos Alberto Gomez Gonzalez, Olivier Wertz, Olivier Absil, Valentin Christiaens, Denis Defrère, Dimitri Mawet, Julien Milli, Pierre-Antoine Absil, Marc Van Droogenbroeck, Faustine Cantalloube, et al. Vip: Vortex image processing package for high-contrast direct imaging. *The Astronomical Journal*, 154(1):7, 2017. 2, 7

[34] CA Gomez Gonzalez, Olivier Absil, P-A Absil, Marc Van Droogenbroeck, Dimitri Mawet, and Jean Surdej. Low-rank plus sparse decomposition for exoplanet detection in direct-imaging adi sequences-the llsg algorithm. *Astronomy & Astrophysics*, 589:A54, 2016. 2, 7

[35] CA Gomez Gonzalez, Olivier Absil, and Marc Van Droogenbroeck. Supervised detection of exoplanets in high-contrast imaging sequences. *Astronomy & Astrophysics*, 613:A71, 2018. 3, 7

[36] SY Haffert, AJ Bohn, J de Boer, IAG Snellen, J Brinchmann, JH Girard, CU Keller, and R Bacon. Two accreting protoplanets around the young star PDS 70. *Nature Astronomy*, 3 (8):749–754, 2019. 1

[37] M Keppler, M Benisty, A Müller, Th Henning, R Van Boekel, F Cantalloube, C Ginski, RG Van Holstein, A-L Maire, A Pohl, et al. Discovery of a planetary-mass companion within the gap of the transition disk around pds 70. *Astronomy & Astrophysics*, 617:A44, 2018. 1

[38] A-M Lagrange, D Gratadour, G Chauvin, T Fusco, D Ehrenreich, D Mouillet, G Rousset, D Rouan, F Allard, É Gendron, et al. A probable giant planet imaged in the β Pictoris disk: VLT/NaCo deep L'-band imaging. *Astronomy & Astrophysics*, 493(2):L21–L25, 2009. 2

[39] Maud Langlois, R Gratton, A-M Lagrange, P Delorme, A Boccaletti, M Bonnefoy, A-L Maire, D Mesa, G Chauvin, S Desidera, et al. The sphere infrared survey for exoplanets (shine)-ii. observations, data reduction and analysis, detection performances, and initial results. *Astronomy & Astrophysics*, 651:A71, 2021. 7

[40] Julien Mairal, Francis Bach, Jean Ponce, and Guillermo Sapiro. Online dictionary learning for sparse coding. In *Proceedings of the 26th annual international conference on machine learning*, pages 689–696, 2009. 2

[41] Christian Marois, David Lafreniere, René Doyon, Bruce Macintosh, and Daniel Nadeau. Angular differential imaging: a powerful high-contrast imaging technique. *The Astrophysical Journal*, 641(1):556, 2006. 2, 7

[42] Christian Marois, Bruce Macintosh, Travis Barman, B Zuckerman, Inseok Song, Jennifer Patience, David Lafrenière, and René Doyon. Direct imaging of multiple planets orbiting the star HR 8799. *Science*, 322(5906):1348–1352, 2008. 7

[43] Christian Marois, B Zuckerman, Quinn M Konopacky, Bruce Macintosh, and Travis Barman. Images of a fourth planet orbiting hr 8799. *Nature*, 468(7327):1080, 2010. 7

[44] Christian Marois, Carlos Correia, Raphael Galicher, Patrick Ingraham, Bruce Macintosh, Thayne Currie, and Rob De Rosa. GPI PSF subtraction with TLOCI: the next evolution in exoplanet/disk high-contrast imaging. In *SPIE Astronomical Instrumentation + Telescopes*, page 91480U. International Society for Optics and Photonics, 2014. 2

[45] D Mesa, R Gratton, P Kervella, M Bonavita, S Desidera, V D’Orazi, S Marino, A Zurlo, and E Rigliaco. Af lep b: The lowest-mass planet detected by coupling astrometric and direct imaging data. *Astronomy & Astrophysics*, 672:A93, 2023. 7
[46] J Milli, D Mouillet, A-M Lagrange, A Boccaletti, D Mawet, G Chauvin, and M Bonnefoy. Impact of angular differential imaging on circumstellar disk images. *Astronomy & Astrophysics*, 545:A111, 2012. 6

[47] Alexey Pavlov, Ole Möller-Nilsson, Markus Feldt, Thomas Henning, Jean-Luc Beuzit, and David Mouillet. Sphere data reduction and handling system: overview, project status, and development. *Advanced Software and Control for Astronomy II*, 7019:1093–1104, 2008. 6

[48] Marshall D Perrin, Anand Sivaramakrishnan, Russell B Makidon, Ben R Oppenheimer, and James R Graham. The structure of high strehl ratio point-spread functions. *The Astrophysical Journal*, 596(1):702, 2003. 5

[49] Sai Krishanth PM, Ewan S Douglas, Justin Hom, Ramya M Anche, John Debes, Isabel Rebolledo, and Bin B Ren. Nmf-based gpu accelerated coronagraphy pipeline. In *Techniques and Instrumentation for Detection of Exoplanets XI*, pages 668–679. SPIE, 2023. 2

[50] Laurent Pueyo. Direct imaging as a detection technique for exoplanets. *Handbook of Exoplanets*, pages 705–765, 2018. 2

[51] Bin Ren, Laurent Pueyo, Guangtun Ben Zhu, John Debes, and Gaspard Duchêne. Non-negative matrix factorization: robust extraction of extended structures. *The Astrophysical Journal*, 852(2):104, 2018. 2

[52] Bin Ren, Laurent Pueyo, Christine Chen, Élodie Choquet, John H Debes, Gaspard Duchêne, François Ménard, and Marshall D Perrin. Using data imputation for signal separation in high-contrast imaging. *The Astrophysical Journal*, 892(2):74, 2020. 2

[53] Erez N Ribak and Szymon Gladysz. Fainter and closer: finding planets by symmetry breaking. *Optics Express*, 16(20): 15553–15562, 2008. 5

[54] Jean-Baptiste Ruffio, Bruce Macintosh, Jason J Wang, Laurent Pueyo, Eric L Nielsen, Robert J De Rosa, Ian Czekala, Mark S Marley, Pauline Arriaga, Vanessa P Bailey, et al. Improving and assessing planet sensitivity of the gpi exoplanet survey with a forward model matched filter. *The Astrophysical Journal*, 842(1):14, 2017. 2

[55] Matthias Samland, J Bouwman, DW Hogg, W Brandner, T Henning, and Markus Janson. Trap: A temporal systematics model for improved direct detection of exoplanets at small angular separations. *Astronomy & Astrophysics*, 646:A24, 2021. 3

[56] Aniket Sanghi, Jerry W Xuan, Jason J Wang, et al. Efficiently searching for close-in companions around young m dwarfs using a multiyear psf library. *The Astronomical Journal*, 168(5):215, 2024. 3

[57] Nuno C Santos. Extra-solar planets: Detection methods and results. *New Astronomy Reviews*, 52(2-5):154–166, 2008. 1

[58] Rémi Soummer. Apodized pupil lyot coronagraphs for arbitrary telescope apertures. *The Astrophysical Journal Letters*, 618(2):L161, 2004. 2

[59] Rémi Soummer, André Ferrari, Claude Aime, and Laurent Jolissaint. Speckle noise and dynamic range in coronagraphic images. *The Astrophysical Journal*, 669(1):642, 2007. 5

[60] Rémi Soummer, Laurent Pueyo, and James Larkin. Detection and characterization of exoplanets and disks using projections on karhunen-loève eigenimages. *The Astrophysical Journal Letters*, 755(2):L28, 2012. 2, 7

[61] William B Sparks and Holland C Ford. Imaging spectroscopy for extrasolar planet detection. *The Astrophysical Journal*, 578(1):543, 2002. 2

[62] Tomas Stolker, Markus J Bonse, Sascha P Quanz, Adam Amara, Gabriele Cugno, Alexander J Bohn, and Anna Boehle. Pynpoint: a modular pipeline architecture for processing and analysis of high-contrast imaging data. *Astronomy & Astrophysics*, 621:A59, 2019. 7

[63] Éric Thiébaut, Loïc Denis, Laurent Mugnier, André Ferrari, David Mary, Maud Langlois, Faustine Cantalloube, and Nicholas Devaney. Fast and robust exo-planet detection in multi-spectral, multi-temporal data. In *Adaptive Optics Systems V*, pages 1534–1543. SPIE, 2016. 5

[64] William Thompson and Christian Marois. Improved contrast in images of exoplanets using direct signal-to-noise ratio optimization. *The Astronomical Journal*, 161(5):236, 2021. 2

[65] Wesley A Traub and Ben R Oppenheimer. Direct imaging of exoplanets. *Exoplanets*, pages 111–156, 2010. 1

[66] Arthur Vigan, Claire Moutou, Maud Langlois, France Alard, Anthony Boccaletti, Marcel Carbillet, David Mouillet, and Isabelle Smith. Photometric characterization of exoplanets using angular and spectral differential imaging. *Monthly Notices of the Royal Astronomical Society*, 407(1):71–82, 2010. 1, 2

[67] Kevin Wagner, Jordan Stone, Andrew Skemer, Steve Ertel, Ruobing Dong, Dániel Apai, Eckhart Spalding, Jarron Leisenring, Michael Sitko, Kaitlin Kratter, et al. Direct images and spectroscopy of a giant protoplanet driving spiral arms in mwc 758. *Nature Astronomy*, 7(10):1208–1217, 2023. 7

[68] Zahed Wahhaj, Lucas A Cieza, Dimitri Mawet, Bin Yang, Hector Canovas, Jozua de Boer, Simon Casassus, François Ménard, Matthias R Schreiber, Michael C Liu, et al. Improving signal-to-noise in the direct imaging of exoplanets and circumstellar disks with MLOCI. *Astronomy & Astrophysics*, 581:A24, 2015. 2

[69] Trevor N Wolf, Brandon A Jones, and Brendan P Bowler. Direct exoplanet detection using convolutional image reconstruction (construct): A new algorithm for post-processing high-contrast images. *The Astronomical Journal*, 167(3):92, 2024. 3

[70] Kai Hou Yip, Nikolaos Nikolaou, Piero Coronica, Angelos Tsiaras, et al. Pushing the limits of exoplanet discovery via direct imaging with deep learning. In *Machine Learning and Knowledge Discovery in Databases: European Conf., ECML PKDD 2019, Würzburg, Germany, September 16–20, 2019, Proceedings, Part III*, pages 322–338. Springer, 2020. 3

[71] Kai Yu, Yuanqing Lin, and John Lafferty. Learning image representations from the pixel level via hierarchical sparse coding. In *CVPR 2011*, pages 1713–1720. IEEE, 2011. 2

[72] Daniel Zoran and Yair Weiss. From learning models of natural image patches to whole image restoration. In *2011 international conference on computer vision*, pages 479–486. IEEE, 2011. 2
# A New Statistical Model of Star Speckles for Learning to Detect and Characterize Exoplanets in Direct Imaging Observations

## Supplementary Material

### A. Technical details

#### A.1. Estimation of the statistical parameters

In this section, we provide additional information regarding the estimation of the statistical parameters of the Gaussian distribution (see Sec. 4.1). These parameters are estimated by maximum likelihood:

$$ \hat{\mathbf{m}}_j = \frac{1}{T} \sum_t \mathbf{y}_t^{(j)}, \quad \hat{\mathbf{S}}_j = \frac{1}{T} \sum_t (\mathbf{y}_t^{(j)} - \hat{\mathbf{m}}_j)(\mathbf{y}_t^{(j)} - \hat{\mathbf{m}}_j)^{\top}, $$

where $\hat{\mathbf{S}}_j$ in $\mathbb{R}^{p \times p}$ is the sample covariance matrix. Given that $\hat{\mathbf{S}}_j$ usually suffers from a large variance due to small sample statistics (i.e., $p \ge T \simeq 70$), we use a shrinkage estimator leading to the final covariance estimator $\hat{\mathbf{C}}_j$:

$$ \hat{\mathbf{C}}_j = (1 - \hat{\rho}_j) \hat{\mathbf{S}}_j + \hat{\rho}_j \hat{\mathbf{D}}_j, \quad (23) $$

where $\hat{\mathbf{D}}_j$ in $\mathbb{R}^{p \times p}$ is the low-variance/high-bias diagonal matrix formed with diagonal elements of $\hat{\mathbf{S}}_j$, and $\hat{\rho}_j$ in [0, 1]. Hyper-parameter $\hat{\rho}_j$ can be estimated optimally in a data-driven fashion by risk minimization between the true (unknown) covariance $\mathbf{C}_j$ and its shrunk counterpart $\hat{\mathbf{C}}_j$. This leads to the closed-form expression [29]:

$$ \hat{\rho}_j = \frac{\mathrm{tr}(\hat{\mathbf{S}}_j^2) + \mathrm{tr}^2(\hat{\mathbf{S}}_j) - 2 \mathrm{tr}(\hat{\mathbf{S}}_j \circ \hat{\mathbf{S}}_j)}{(T+1)(\mathrm{tr}(\hat{\mathbf{S}}_j^2) - \mathrm{tr}(\hat{\mathbf{S}}_j \circ \hat{\mathbf{S}}_j))}, \quad (24) $$

where $\circ$ stands for Hadamard (element-wise) product.

#### A.2. Dense detection map computation details

In this section, we provide additional details on the computation of the dense detection map introduced in Sec. 4.1. To obtain a full-frame detection map, $\hat{\alpha}$ and $\hat{\sigma}_\alpha$ must be computed for all spatial coordinates, which can be formalized by:

$$ \forall \mathbf{x}_i \in \mathcal{G}_p, \quad \hat{\alpha}_i = \arg \max_{\alpha_i} \ell(\alpha_i, \mathbf{x}_i). \quad (25) $$

Note that it corresponds to solve independent optimization problems at each location $i$. To efficiently solve Eq. (25), we employ a fast approximation to compute the terms $a(\mathbf{x}_t)$ and $b_t(\mathbf{x}_t)$ at each time step and for every initial position. First, the terms $a_t$ and $b_t$ are computed for all positions in the grid $\mathcal{G}_p$ of spatial pixels. In the following, we note $b_{ti} = b_t(\mathbf{x}_i)$, $a_i = a(\mathbf{x}_i)$ and $h_{ji} = h^{(j)}(\mathbf{x}_i)$ for all $\mathbf{x}_i$ in $\mathcal{G}_p$, such that:

$$ \forall (i, t), \quad b_{it} = \sum_{j \in \mathcal{S}(\mathbf{x}_i)} w_j \mathbf{h}_{ji}^{\top} \hat{\mathbf{C}}_j^{-1} (\mathbf{y}_t^{(j)} - \hat{\mathbf{m}}_j), \quad (26) $$

$$ \forall i, a_i = \sum_{j \in \mathcal{S}(\mathbf{x}_i)} w_j \mathbf{h}_{ji}^{\top} \hat{\mathbf{C}}_j^{-1} \mathbf{h}_{ji}. \quad (27) $$

With our convolutional approach, these terms can be computed efficiently by Cholesky factorization of the precision matrix, as detailed in Sec. A.3. Let $b_t$ in $\mathbb{R}^{H \times W}$ and $a$ in $\mathbb{R}^{H \times W}$ represent the

resulting quantities. Then, $b_t(\mathbf{x}_t)$ and $a(\mathbf{x}_t)$ can be approximated by spatial interpolation of $b_t$ and $a$. However, the trajectories of all exoplanets in the observation are governed by a common parallactic rotation vector $\phi$ in $\mathbb{R}^T$. So in practice, this interpolation can be computed efficiently for each time-step $t$, by rotating the frames $a$ and $b_t$ by an angle $-\phi_t$. Thus, a dense detection map $\hat{\gamma}$ in $\mathbb{R}^{H \times W}$ can be obtained with:

$$ \hat{\gamma} = \frac{\hat{\alpha}}{\hat{\sigma}} = \frac{\sum_t \mathbf{R}_{-\phi_t}(b_t)}{\sqrt{\sum_t \mathbf{R}_{-\phi_t}(a)}}, \quad (28) $$

where $\mathbf{R}_{-\phi_t}: \mathbb{R}^{H \times W} \rightarrow \mathbb{R}^{H \times W}$ is the frame rotation operator with angle $-\phi_t$, and $\hat{\alpha}$, $\hat{\sigma}$ in $\mathbb{R}^{H \times W}$ such that:

$$ \hat{\alpha} = \frac{\sum_t \mathbf{R}_{-\phi_t}(b_t)}{\sum_t \mathbf{R}_{-\phi_t}(a)} \quad \text{and} \quad \hat{\sigma} = \frac{1}{\sqrt{\sum_t \mathbf{R}_{-\phi_t}(a)}}. \quad (29) $$

#### A.3. Fast computation of $b_t$ and $a$

In this section, we detail an efficient strategy to derive a dense two-dimensional detection map $\hat{\gamma}$ serving as the final quantity to detect exoplanets, see Sec. 4.1.

The terms $b_t$ and $a$ in $\mathbb{R}^{H \times W}$ involved in $\hat{\gamma}$ (see Eqs. (26)–(28)) can be computed efficiently through:

$$ b_t = \sum_j w_j \mathbf{Q}_j (\mathbf{b}_t^{(j)}), \quad (30) $$

$$ a = \sum_j w_j \mathbf{Q}_j (a^{(j)}), \quad (31) $$

where $\mathbf{Q}_j$ is the linear operator that places the patch indexed by $j$ in the correct position in the image, and $b_t^{(j)}$, $a^{(j)}$ are in $\mathbb{R}^p$. The individual terms $a^{(j)}$ and $b_t^{(j)}$ are given by:

$$ b_t^{(j)} = \mathbf{P}^{\top} \hat{\mathbf{C}}_j^{-1} (\mathbf{y}_t^{(j)} - \hat{\mathbf{m}}_j), \quad (32) $$

$$ a^{(j)} = \mathrm{diag} (\mathbf{P}^{\top} \hat{\mathbf{C}}_j^{-1} \mathbf{P}), \quad (33) $$

where $\mathbf{P}$ in $\mathbb{R}^{p \times p}$ is a circulant matrix, with each column representing the flattened PSF centered on a different pixel. Equation (32) can be implemented efficiently with a convolution operation. Implementing Eq. (33) directly is however computationally challenging as it requires storing the matrix $\mathbf{P}$, which is highly inefficient, especially for large patches. To circumvent this issue, we propose to decompose $a^{(j)}$ as:

$$ a^{(j)} = \mathbf{N}(\hat{\mathbf{L}}_j \mathbf{P}), \quad (34) $$

where $\mathbf{N}: \mathbb{R}^{p \times p} \rightarrow \mathbb{R}^p$ is the operator that computes the column-wise squared $\ell_2$-norm, and the matrix $\hat{\mathbf{L}}_j$ in $\mathbb{R}^{p \times p}$ is the lower triangular matrix obtained via Cholesky decomposition of the precision matrix (i.e., $\hat{\mathbf{C}}_j^{-1} = \hat{\mathbf{L}}_j \hat{\mathbf{L}}_j^{\top}$). The formulation of Eq. (34) allows implementing the matrix multiplication with $\mathbf{P}$ efficiently as a convolution operation.
**Algorithm 1** Iterative procedure for flux and position estimation.

**Require:** $y \in \mathbb{R}^{T \times H \times W}, z \in \mathbb{R}^3$
**for** $l = 0 \cdots L-1$ **do**

$\alpha, x_0 \leftarrow z$

**for** $t = 0 \cdots T-1$ **do** ▷ Update nuisance estimate
$x_t \leftarrow r(x_0, \phi_t)$
$s_t \leftarrow y_t - h(x_t)$

**end for**

$g, \mathbf{H} \leftarrow \mathbf{0}_{\mathbb{R}^3}, \mathbf{0}_{\mathbb{R}^{3 \times 3}}$

**for** $t = 0 \cdots T-1$ **do** ▷ Gradient and Hessian
**for** $j \in S(x_t)$ **do**

$\hat{\mathbf{m}}_j, \hat{\mathbf{C}}_j \leftarrow \text{Estim}(\{s_t^{(j)}\}_t)$ ▷ Eqs. (23)–(24)

$g \leftarrow g + w_j \frac{\partial \ell_j}{\partial z}(z, \hat{\mathbf{m}}_j, \hat{\mathbf{C}}_j)$

$\mathbf{H} \leftarrow \mathbf{H} + w_j \frac{\partial^2 \ell_j}{\partial z^2}(z, \hat{\mathbf{m}}_j, \hat{\mathbf{C}}_j)$

**end for**

**end for**

$z \leftarrow \text{Proj}_Z(z - \mathbf{H}^{-1}g)$ ▷ Projection on feasible set
**end for**
**return** $z$

## A.4. Iterative optimization scheme for astrometry

We detail in this section the characterization procedure mentioned in Sec. 4.1 for precisely estimating the flux $\alpha$ and the two-dimensional sub-pixel position $x_0$ of an exoplanet. This procedure is described by Algorithm 1 in pseudo-code form, where $z = [\alpha, x_0]$ in $\mathbb{R}^3$ is the variable to be refined. The position is initialized to the coordinates of the pixel with the highest detection score from the detection step, while the flux is initialized using the value extracted from the flux map $\alpha$ at that location. We decompose the involved loss function $\ell$ as:

$$ \ell(z, \{m_j, C_j\}_j) = \sum_t \sum_{j \in S(x_t)} w_j \ell_j(z, m_j, C_j), \quad (35) $$

where:

$$ \ell_j([\alpha, x_0], m_j, C_j) = \| \mathbf{L}_j (\mathbf{y}_t^{(j)} - \alpha \mathbf{h}^{(j)}(r(x_0, \phi_t)) - m_j) \|_2^2, \quad (36) $$

such that $\mathbf{C}_j^{-1} = \mathbf{L}_j \mathbf{L}_j^\top$. We also define the feasible set $Z$ such that $\alpha > 0$ (enforcing a non-negativity constraint on the exoplanet flux) and the estimated position remains within a fixed radius of the initial estimate. We denote by $\text{Proj}_Z$ the projection associated with it. In Algorithm 1, we also denote by Estim(.) the procedure described in Eqs. (23)–(24) to estimate the statistical parameters of the nuisance component.

## A.5. Prior information on the exoplanet spectrum

The channel-wise detection maps $\tilde{\gamma}_c = \tilde{b}_c/\sqrt{\tilde{a}_c}$ are combined into a single detection map by averaging $\tilde{a}_c$ and $\tilde{b}_c$ with weights $e_c$ in $\mathbb{R}^+$ satisfying $\sum_c e_c = 1$ such that:

$$ \tilde{a}_c = \sum_c e_c \tilde{a}_c \quad \text{and} \quad \tilde{b}_c = \sum_c e_c \tilde{b}_c. \quad (37) $$

<figure><img src="image_7.png" /><figcaption>Figure F. Calibrated PFA versus test PFA: the calibration procedure allows obtaining a reliable approximation of the PFA, except in the extreme low-PFA regime (high detection score under $\mathcal{H}_0$) where samples are exceedingly rare.</figcaption></figure>

The final detection score, spectrally combined over $C$ channels, is $\tilde{\gamma}_C = \tilde{b}_C/\sqrt{\tilde{a}_C}$ in $\mathbb{R}^{H \times W}$. The weights $e_c$ reflect prior assumptions about the exoplanet spectrum, and can be used (if available) to boost the detection of exoplanets having similar spectrum, see e.g., [12, 26]. For all experiments and baselines, we assumed a non-informative flat prior (i.e., $e_c = 1/C$). If additional information about the exoplanet's spectrum is available, the prior can be adjusted accordingly.

## A.6. Calibration procedure

This section outlines the calibration procedure for the full model described in Sec. 4.3, following the approach proposed in [5], which we briefly summarize here.

We aim to calibrate the predicted probability of false alarm (PFA), a critical metric for astronomers to quantify detection uncertainty. This requires characterizing the output distribution of our method under the null hypothesis $\mathcal{H}_0$ (i.e., when no source is present), which we denote $\tilde{\gamma}|\mathcal{H}_0$. For that purpose, we use a calibration dataset $\mathcal{D}_{\text{calib}}$, consisting of 15 observations excluded from training, validation, and testing. Temporal shuffling of the frames and parallactic rotation reversal are applied to augment the dataset and ensure no real source is present. These custom augmentations disrupt the temporal coherence of any (potentially unknown) real exoplanets, preventing them from biasing the training and calibration processes. We apply the model to 10,000 cubes from $\mathcal{D}_{\text{calib}}$, concatenate all output pixels into a single vector $\tilde{\gamma}_{\text{calib}}$ of length $n_{\text{calib}}$, and compute the empirical cumulative distribution function (eCDF):

$$ \forall \tau \in \mathbb{R}, \quad \hat{F}_{\tilde{\gamma}|\mathcal{H}_0}(\tau) := \frac{1}{n_{\text{calib}}} \sum_{i=0}^{n_{\text{calib}}-1} \mathbf{1}_{\tilde{\gamma}_{\text{calib},i} \le \tau}. \quad (38) $$

where 1 is the indicator function. The estimated probability of false alarm is related to the eCDF through:

$$ \forall \tau \in \mathbb{R}, \quad \widehat{\text{PFA}}(\tau) = \hat{P}(\tilde{\gamma} > \tau | \mathcal{H}_0) = 1 - \hat{F}_{\tilde{\gamma}|\mathcal{H}_0}(\tau). \quad (39) $$
<figure><img src="image_8.png" /><figcaption>(a) Distribution of the number of frames per observation</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(b) Repartition of the filters used</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(c) Distribution of the parallactic rotation per observation</figcaption></figure>

<caption>Figure G. Additional details on the F150 dataset</caption>

In Fig. F, we evaluate the PFA across a range of detection thresholds using the calibration set and apply the same thresholds to the test set for comparison. The results show that the calibrated PFA closely matches the test PFA, except in the extremely low-PFA regime. This slight discrepancy is expected, as high detection scores under $H_0$ are exceedingly rare, leading to increased noise in PFA estimation. In this regime, our estimate is nevertheless conservative, i.e., the experienced PFA is slightly smaller than predicted one.

## A.7. Implementation details

Our model implementation and training are based on the PyTorch framework. We optimize the neural network weights using the Adam optimizer with a learning rate of $5 \times 10^{-4}$, a batch size of 16, and 50,000 iterations. We apply data augmentation techniques, including random 90° rotations, vertical/horizontal flipping, and frame shuffling. Bicubic interpolation is used for spectral and exoplanet alignment. For ADI, the model is trained on 64 frames, while for ASDI, data cubes have 32 frames with two channels each. All evaluations and tests are done with 64 frames in both ADI and ASDI. To ensure accurate results, pixels containing real known exoplanets are masked out during training before the temporal summation in Eq. (29). The loss is calculated on the final output for pixels with at least 8 valid frames, unaffected by real exoplanets.

## B. Dataset description

SHINE-F150 is a subset of the SPHERE large survey, with 322 data cubes (20, 462 images in total) from 150 stars. The spatial resolution is 12 milliarcseconds/pixel. The spectral resolution $\lambda/\Delta\lambda$ is 20. Detailed description is in [39]. We followed [5] to exclude corrupted data, leaving 220 datasets, with 165 used for training. The main statistics of the F150 dataset are shown in Fig. G.

## C. Comparison of state-of-the-art methods

<table><caption>Table D. Comparison of advanced exoplanet imaging algorithms based on key desirable properties. P1=observation-independence of the model, P2=physical interpretability, P3=multi-spectral, P4=detection sensitivity, P5=accuracy near the star, P6=position and flux estimation, P7=practicality on large-scale surveys.</caption><thead><tr><th>Method</th><th>nuisance model</th><th>P1</th><th>P2</th><th>P3</th><th>P4</th><th>P5</th><th>P6</th><th>P7</th></tr></thead><tbody><tr><td>PACO</td><td>statistical on data patches</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td><td>✓</td></tr><tr><td>SODINN</td><td>learned on PCA residuals</td><td>✗</td><td>✗</td><td>✗</td><td>✓✓</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>deep PACO</td><td>learned after PACO whitening</td><td>✗</td><td>✗</td><td>✓</td><td>✓✓</td><td>✗</td><td>✓</td><td>✗</td></tr><tr><td>Super-RDI</td><td>PCA</td><td>✓</td><td>✓</td><td>✗</td><td>✓✓</td><td>✗</td><td>✗</td><td>✓</td></tr><tr><td>MWIN5-RB</td><td>learned on PCA residuals</td><td>✓</td><td>✓</td><td>✗</td><td>✓✓</td><td>✓</td><td>✗</td><td>✓✓</td></tr><tr><td>ConStruct</td><td>generative, learned on data</td><td>✓</td><td>✓</td><td>✗</td><td>✓✓</td><td>✓</td><td>✗</td><td>✓✓✓</td></tr><tr><td>MODEL&CO</td><td>statistical on learned features</td><td>✓</td><td>✗</td><td>✗</td><td>✓✓✓</td><td>✓✓</td><td>✗</td><td>✓✓✓</td></tr><tr><td>proposed</td><td>multi-scale stat. on learned features</td><td>✓</td><td>✓</td><td>✓</td><td>✓✓✓</td><td>✓✓✓</td><td>✓</td><td>✓✓✓</td></tr></tbody></table>

Table D compares the key properties of state-of-the-art algorithms for exoplanet detection via direct imaging.

## D. Additional results

In this section we provide additional visual results, and show detection maps obtained with synthetic exoplanets. Consistent with the experimental findings in the main body of the paper, our method performs on par or better than baselines in ADI, and outperforms other approaches in ASDI.
<figure><img src="image_11.png" /><figcaption>Figure H. Detection maps on observations of HD 102647 star with synthetic exoplanets.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure J. Detection maps on observations of HD 206860 star with synthetic exoplanets.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure I. Detection maps on observations of HD 188228 star with synthetic exoplanets.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure K. Detection maps on observations of HD 216803 star with synthetic exoplanets.</figcaption></figure>