# Joint unmixing and deconvolution
for angular and spectral differential imaging

Olivier Flasseur, Loïc Denis, Éric Thiébaut and Maud Langlois O. Flasseur is
with LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université,
Univ. Paris Diderot, Sorbonne Paris Cité, 5 Place Jules Janssen, 92195 Meudon,
France and was previously at Université de Lyon, Université Lyon1, ENS de
Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR 5574, F-69230,
Saint-Genis-Laval, France where most of this work was performed.L. Denis is
with Université de Lyon, UJM-Saint-Etienne, CNRS, Institut d Optique Graduate
School, Laboratoire Hubert Curien UMR 5516, F-42023, Saint-Étienne, France.É.
Thiébaut and M. Langlois are with Université de Lyon, Université Lyon1, ENS de
Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR 5574, F-69230,
Saint-Genis-Laval, France.

###### Abstract

Angular and spectral differential imaging is an observational technique used
in astronomy to study the close environment of stars. The relative angular
motion and spectral scaling between on-axis and off-axis sources are exploited
by processing techniques to separate two components: the residual light from
the star and the light coming from surrounding objects such as circumstellar
disks. This paper describes a method to unmix these two components and
deconvolve the image of disks. It is based on a statistical modeling of the
residual star light, in particular its spatial and spectral correlations.
These correlations have so far been neglected by standard reconstruction
techniques. We show on several datasets from the SPHERE instrument at the Very
Large Telescope, Chile, that accounting for correlations strongly improves the
reconstructions.

###### Index Terms:

Inverse problems, unmixing, deconvolution, Gaussian scale mixture, covariance
matrix.

##  I Introduction

Our understanding of the formation and evolution of planetary systems has
radically evolved these last two decades thanks to unprecedented progress in
observational techniques. In particular, several methods have emerged to
detect exoplanets [1]. Most of them are classified as _indirect_ observational
methods since the detections comes from the analysis of perturbations of the
host star induced by the presence of the exoplanets. In contrast, _direct
imaging_ performs the observation and collects photons of both the star and
its close environment. This latter observational technique offers the unique
capability of acquiring multi-spectral images of the so-called circumstellar
environment in the vicinity of young stars. Objects of high astrophysical
interest located in the circumstellar environment include exoplanets and
circumstellar disks made of gas and dust in which exoplanets can form by
accretion. In this context, direct imaging is a method of choice to study the
interactions between exoplanets and the associated circumstellar disks, as
well as the planet formation processes.

Direct imaging requires to face two observational challenges: the objects of
interest (disks and exoplanets) have a very low contrast with respect to their
host star (typically lower than $10^{-4}-10^{-3}$ in the infrared) and are
located in the immediate surrounding of the star, which requires a high
angular resolution to separate them (disks are generally seen inside an angle
of less than 1 arcsecond). The angular resolution requirement can be achieved
using large ground-based telescopes equipped with a high-performance adaptive
optics system (to compensate in real-time for atmospheric turbulence). The
contrast is improved by masking out most of the star light with a coronagraph.
However, this is not sufficient to recover interpretable images of the
circumstellar environment: residual light from the star still dominates, see
Fig. 1. In order to further reduce the impact of the star light, differential
imaging is used: several images are captured in configurations that induce a
relative motion of the objects of interest with respect to the structures due
to the star light (namely, the speckles created by the residual phase
aberrations in the pupil of the telescope). Two configurations can be used: in
angular differential imaging (ADI) [2], a sequence of images is acquired
during a few hours of observation. During the acquisition, the pupil of the
telescope keeps a constant orientation (so-called _pupil-tracking mode_) while
the field-of-view rotates due to Earth’s rotation. This leads to a rotation of
the objects of interest about the optical axis in the images, with quasi-
static speckle created by uncorrected optical aberrations. In spectral
differential imaging (SDI) [3], images in several spectral bands are captured.
Due to diffraction, the speckle pattern scales linearly with the wavelength,
in first approximation. After proper spectral rescaling of the images, the
speckle patterns are aligned while the objects of interest undergo a radial
motion due to the scaling transform. ADI and SDI can be associated to form
Angular and Spectral Differential Imaging (ASDI) sequences. The images
recorded in ADI, SDI, or ASDI should be combined, in a processing step, to
improve the contrast and to obtain interpretable images of the circumstellar
environnement. The classical processing pipeline performs first an estimation
of the component due to the star (e.g., by averaging the stack of images with
aligned speckles), a removal of this component from the data, and finally an
alignment and stacking of the residuals (so as to compensate for the rotations
and scaling of the field-of-view).

Beyond a simple average, the component due to the star can be estimated using
a median (cADI method [2]), a weighted linear combination [4, 5, 6], or a
principal component analysis [7, 8]. However, these approaches suffer from a
common pitfall: part of the signal of interest is included in the star
component, it is then lost when the star component is subtracted from the
data, a phenomenon described as the _self-subtraction problem_ [9, 10]. It is
particularly important close to the star, where the rotations and scalings
have the least effect, and the disentangling between the component of interest
and the star light is the most difficult. Other artifacts can also be
observed: partial replicas, suppression of some smooth extended-structures,
smearing or non-uniform attenuations.

To reduce the impact of self-subtraction, several approaches were considered.
Some works consider a parametric model of a disk and iteratively adjust its
parameters [9, 11]. This kind of approach is only applicable to simple disk
structures such as ellipses. Another technique, called Reference Differential
Imaging (RDI), uses additional images of an isolated star (i.e., without any
exoplanet or disk) either captured simultaneously (the so-called star-hopping
technique [12]) or from previous observations. The effectiveness of this
approach depends on the similarity between the reference and the actual
observations (similar star brightness and spectrum, close observation
conditions). Finally, the most recent and promising approaches address jointly
the problem of estimating the star light residuals and reconstructing the
spatial distribution of the disk and exoplanets. In ADI, two approaches based
on an inverse problems formulation were recently proposed: MAYONNAISE [13] and
REXPACO [14]. These two algorithms use different strategies for the separation
of the two components and the regularization of the inversion. MAYONNAISE uses
a first step to identify a low-dimension subspace that contains the residual
star light (after speckle alignment, the sequence of images is indeed highly
correlated). A truncated singular value decomposition (SVD) is applied to the
sequence with aligned speckles to estimate the star component, then the
residuals are derotated (to align the disk and exoplanet component) and
averaged to produce an estimate of the component of interest. The estimation
of each component is then alternately refined. The rank of the truncated SVD
is progressively increased from 1 to 10\. At the end of this first step, the
subspace spanned by the sequence of residual star light images is identified.
In a second step, several components are estimated by joint minimization of a
data fidelity term: the star light residuals (restricted to lie within the
subspace identified in the first step), the disk (restricted to have a sparse
representation in a wavelet basis), the exoplanets (restricted to be sparse).
Positivity constraints are also enforced during the minimization. REXPACO
follows quite a different modeling: instead of explicitly estimating the
residual star light in each image of the sequence, it builds a statistical
model of all fluctuations other than the component of interest (i.e., noise
and star light). It learns the spatial correlations of these fluctuations at
the scale of 2D image patches, following an approach first introduced for the
detection of exoplanets (the so-called PACO algorithm, based on PAtch
COvariances [15]). The component of interest is deconvolved with an edge-
preserving smoothness regularization and positivity constraint.

This paper addresses the problem of reconstructing circumstellar disks from
ASDI sequences. Compared to ADI sequences, this raises several challenges:

  * •

the modeling of temporal and spectral fluctuations of the residual star light,

  * •

the joint exploitation of both spectral and temporal information to
effectively extract the component of interest,

  * •

the tractability of estimating high-dimensional models from large datasets.

Typical ASDI datasets produced by the Integral Field Spectrograph (IFS) of the
Spectro-Polarimetry High-contrast Exoplanet Research instrument (SPHERE [16])
at the Very Large Telescope (VLT), in Chile, are $N=290\times 290$ pixels,
have $L=39$ spectral bands, and $T\approx 100$ exposures. Several hundred
million pixel measurements must then be combined to produce a multi-spectral
reconstruction of the component of interest.

_Our contributions:_ This paper extends the REXPACO algorithm [14] to ASDI
sequences. This extension involves several specific developments, in
particular:

  * •

a separable model of the spatial, temporal and spectral covariances,

  * •

a technique to estimate the components of the covariance model,

  * •

a strategy to jointly refine the model of the residual star light and to
reconstruct the disk,

  * •

a spatio-spectral regularization of the reconstructed multi-spectral images.

While the observation model is quite specific to direct imaging in astronomy,
our statistical modeling of the nuisance terms is very general and could be
applied to other imaging modalities for component unmixing problems.

Section II develops the statistical model for the residual star light and
different noise contributions. Around this model, section III builds a
reconstruction method that extracts and deconvolves the component of interest:
the multi-spectral image of the disk that surrounds the star. Results on
several ASDI sequences obtained by the SPHERE-IFS instrument at the VLT are
presented in section IV.

##  II Statistical model of all nuisance terms

In contrast to other methods in the literature [2, 10, 13], we do not
explicitly extract the residual star light component from the data but rather
develop a statistical model to capture both the residual star light (i.e., the
speckles) and the various noise contributions (photon noise, electronic noise,
sensor deficiencies). Thanks to the speckle alignment techniques (namely,
pupil tracking and image rescaling according to the wavelength), residual star
light is very similar from one image to the next (up to some chromatic
factor). There are, however, some fluctuations due to noise, chromatic
phenomena, and the evolution of the phase aberrations during observation.
These fluctuations display some spatial and spectral correlations and are
highly non-stationary: in particular, they are much stronger close to the star
(i.e., the image center).

Figure 1: Illustration of an ASDI dataset: a sequence of multi-spectral cubes
of the star SAO 206462 collected by SPHERE-IFS at the VLT, Chile. (a) images
captured at different wavelengths; (b) spatio-spectral slices along the two
lines –1– and –2– drawn in (a); (c) spatio-temporal slices along the lines –1–
and –2–. The four square areas define four regions studied in more details in
Fig.2. The component of interest, a spiral-shaped circumstellar disk, is shown
in (d) based on our reconstruction given in Fig.6. In the first channel, shown
in blue, the signal of the disk is faint compared to the stellar leakages.

Figure 1 shows an ASDI dataset of a star (SAO 206462) surrounded by a bright
disk. Slices along different dimensions of this 4D dataset are displayed. The
coronagraphic mask is aligned with the star, at the center of the field-of-
view (center of the images shown in Fig.1(a)). Residual star light dominates
the central area and extends over most of the field of view. It takes the form
of granular intensity structures (speckles). During a pre-processing step, all
images were rescaled by a wavelength-specific factor
$\lambda_{\text{ref}}/\lambda$ to compensate for diffraction and spatially
align the speckles. The solid line –1– drawn in the $x$ direction in Fig.1(a)
crosses a bright speckle. This speckle is visible in the left part of Fig.1(b)
and the first row of Fig.1(c). It remains at the same spatial location for all
wavelengths $\lambda$ and all capture times $t$. Structures of interest, such
as the disk that surrounds the star SAO 206462, undergo a rotation about the
image center throughout time and a scaling with the wavelength (due to the
rescaling applied in the preprocessing step). These spatial transformations
are visible in the slices along the dotted line –2– drawn in the images of
Fig.1(a): the line crosses the disk (as well as an area with strong residual
star light, close to the image center). The spatio-spectral slice shown at the
right of Fig.1(b) displays a scaling of the disk with respect to the
wavelength (shorter wavelengths are dilated due to the speckle-aligning
preprocessing), whereas the rotation motion can be noted in the spatio-
temporal slices shown at the bottom of Fig.1(c), in particular for a bright
structure of the disk put in a box, which is moving closer to the image center
during the sequence. Fig.1(d) shows only the component of interest: the
circumstellar disk. The images were obtained with the reconstruction method
introduced in this paper, see Fig.6 for a different visualization of the
reconstructed disk. By comparing Fig.1(a) and Fig.1(d), it appears that high-
contrast observations suffer from a strong nuisance component which has to be
numerically suppressed in order to reconstruct the component of interest.

The modeling of residual star light and noise has a strong impact on the
reconstruction of the component of interest $\boldsymbol{u}$, as discussed in
section IV. In the following of the section, we assume that the object
$\boldsymbol{u}$ has a negligible impact on the statistical distribution of
the nuisance term, i.e., the statistical distribution of the aligned data
$\text{p}_{V}(\boldsymbol{v})$ in the absence of disk or exoplanet is
identical to the distribution of the residuals
$\text{p}_{V}(\boldsymbol{v}-\mathbf{M}\boldsymbol{u})$ obtained when the
modeled contribution $\mathbf{M}\boldsymbol{u}$ of the component of interest
has been subtracted from the data (the direct model, $\mathbf{M}$, is
presented in section III-A, Eq.(23)). This assumption is reasonable in high
contrast imaging since residual light from the star is the dominant component
in the data: the photon noise level is then almost independent of the presence
of other sources in the field of view. Besides, we introduce in section III-C
a strategy to jointly estimate the statistical distribution of the nuisance
terms and reconstruct the component of interest. This strategy extends the
applicability of our modeling to data with signal-dependent perturbations.

###  II-A Patch-based statistical modeling

Image patches (i.e., neighborhoods of a few tens to a hundred pixels) offer an
interesting trade-off between locality (small enough to capture a local
behavior) and complexity (they include enough pixels to collect geometrical
and textural information). Their use has been very successful in image
restoration, from methods based on image self-similarity [17], collaborative
filtering [18], sparse coding [19, 20], mixture models [21, 22], or Gaussian
models [23]. Whereas deep neural networks have become the state-of-the-art
approach to learn rich models of natural images (either generative or
discriminative), patch-based models retain serious advantages when the number
of training samples is limited or in the case of highly non-stationary images.
As can be seen in Fig.1(a), images in an ASDI dataset are far from stationary.
Observations made during separate nights around different stars also often
display significantly different structures because of changes in the observing
conditions (which impact the residual aberrations uncorrected by adaptive
optics, and hence the spatial distribution of speckles due to star light) and
star brightness. This limits the possibility to use external observations to
learn a model to process a specific ASDI sequence and motivates the
development of a patch-based approach based solely on the ASDI sequence of
interest.

Under our patch-based model, the distribution of an ASDI sequence
$\boldsymbol{v}\in\mathbb{R}^{NLT}$ _without disk or exoplanet_ (formed by the
collection of $T$ multi-spectral images with $L$ spectral bands and $N$ pixels
in each band) is given by:

$\displaystyle\text{p}_{V}(\boldsymbol{v})\approx\prod_{n\in\mathbb{K}}\text{p}_{V_{n}}(\mathbf{E}_{n}\boldsymbol{v})\,,$
(1)

where $\text{p}_{V}$ is the joint distribution of the whole ASDI dataset,
$\mathbf{E}_{n}$ is the linear operator that extracts a $K\times L\times
T$-pixel spatio-spectro-temporal patch centered at the $n$-th spatial location
of the field-of-view (i.e., $\boldsymbol{v}_{n}=\mathbf{E}_{n}\boldsymbol{v}$
is a 4D-patch, we do not differentiate the $x$ and $y$ spatial dimensions to
simplify the notations but rather use 2D spatial indices $n$). The set of
spatial locations $\mathbb{K}$ is defined to prevent patch overlapping while
tiling the whole field-of-view (i.e., $\text{Card}(\mathbb{K})\times K=N$ and
juxtaposed square patches are used).

The model (1) assumes a statistical independence between patches, which is a
simplifying hypothesis that eases a data-driven learning of the distribution
$\text{p}_{V}$. In the sequel, each distribution $\text{p}_{V_{n}}$ is modeled
by a different multivariate Gaussian in order to capture the correlations
between observations within a spatio-spectro-temporal patch. By adapting the
parameters of these Gaussian distributions to the spatial location $n$, a non-
stationary model is obtained, with the capability to capture the variations
between areas close to the star (at the center of the image) and areas farther
away. The statistical model of a patch is thus given by:

$\displaystyle\text{p}_{V_{n}}(\boldsymbol{v}_{n})=\frac{1}{\sqrt{|2\pi\mathbf{C}_{n}|}}\exp\left(-\tfrac{1}{2}\bigl{\|}\boldsymbol{v}_{n}-\boldsymbol{\mu}_{n}\bigr{\|}_{\mathbf{C}_{n}^{-1}}^{2}\right)\,,$
(2)

with
$\|\boldsymbol{a}\|_{\mathbf{B}}^{2}=\boldsymbol{a}^{\top}\mathbf{B}\boldsymbol{a}$
and $|\mathbf{C}_{n}|$ the determinant of matrix $\mathbf{C}_{n}$. The
Gaussian distribution $\text{p}_{V_{n}}$ is parameterized by the mean patch
$\boldsymbol{\mu}_{n}\in\mathbb{R}^{KLT}$ and the covariance matrix
$\mathbf{C}_{n}\in\mathbb{R}^{KLT\times KLT}$. In order to estimate these two
quantities at each location $n$, additional hypotheses and an estimation
technique are required.

###  II-B Constraining the structure of the average vector and the covariance
matrix

Estimating and handling different Gaussian parameters for each patch location
is not doable owing to the number of parameters involved: the set of all mean
vectors $\\{\boldsymbol{\mu}_{n}\\}_{n\in\mathbb{K}}$ has as many free
parameters as the total number of measurements in $\boldsymbol{v}$ (i.e.,
$NLT$) and the set of all covariance matrices
$\\{\mathbf{C}_{n}\\}_{n\in\mathbb{K}}$ represents many times the number of
measurements in $\boldsymbol{v}$ (more precisely, $NLT(KLT+1)/2$ free
parameters, which represents more than $300\,000$ times the size of
$\boldsymbol{v}$ for typical values of $K\approx 13\times 13$, $L\approx 39$,
and $T\approx 100$).

There are two options to reduce the number of parameters in the Gaussian
models of equations (1) and (2): (i) assume some level of stationarity for the
means or covariances with respect to the spatial location $n$; (ii) enforce a
structure on the mean $\boldsymbol{\mu}_{n}$ and the covariance
$\mathbf{C}_{n}$ to limit the number of free parameters. Beyond obtaining more
tractable models, these assumptions are also indispensable, for a single ASDI
dataset $\boldsymbol{v}$, to constrain the estimator of the parameters of the
Gaussian models.

The strong spatial non-stationarity of ASDI datasets led us to favor option
(ii). We considered several ways to select a structure suitable to ASDI
observations and built on our experience of point-source detection in ASDI
datasets (i.e., exoplanet detection [15, 24, 25]). We found that it is
preferable to use a common mean vector $\boldsymbol{\mu}_{n}$ for all times
$t$ rather than a time-specific mean vector common to all wavelengths (the
spectral variations being stronger than the temporal fluctuations):

$\displaystyle\boldsymbol{\mu}_{n}=\text{vec}\\!\left(\begin{pmatrix}|\\\
\boldsymbol{\mu}_{n}^{\text{spec }}\\\
|\end{pmatrix}\overset{\text{\scriptsize$\longleftarrow
T\longrightarrow$}}{\begin{pmatrix}1&\cdots&1\end{pmatrix}}\right)\,,$ (3)

where $\text{vec}(\cdot)$ performs the vectorization of a matrix by stacking
its columns (it transforms a $KL\times T$ matrix into a vector of dimension
$KLT$) and $\boldsymbol{\mu}_{n}^{\text{spec}}$ is a $KL$-pixel multi-spectral
vector that represents the temporal average of the multi-spectral patches.

To capture the structures of both the spatial and the spectral covariances, we
model the covariance between two pixels of the patch $\boldsymbol{v}_{n}$ by:

$\text{Cov}\\!\left[\boldsymbol{v}_{n,(k_{1},\ell_{1},t_{1})},\,\boldsymbol{v}_{n,(k_{2},\ell_{2},t_{2})}\right]=0\text{
if }t_{1}\neq t_{2}\,,\\\
=\sigma_{n,t}^{2}\mathbf{C}_{n}^{\text{spat}}(k_{1},k_{2})\mathbf{C}_{n}^{\text{spec}}(\ell_{1},\ell_{2})\text{
if }t_{1}=t_{2}=t\,,$ (4)

where $\sigma_{n,t}^{2}$ is a scalar that encodes the global level of
fluctuation in the multi-spectral slice at time $t$,
$\mathbf{C}_{n}^{\text{spat}}$ is a $K\times K$ covariance matrix encoding the
spatial structure of the fluctuations (a $K$-pixel spatial patch corresponds
to a 2D square window, so this covariance matrix contains information about 2D
spatial structures), and matrix $\mathbf{C}_{n}^{\text{spec}}$ is an $L\times
L$ covariance matrix encoding spectral correlations. To prevent a degeneracy
by multiplicative factors, we normalize covariance matrices
$\mathbf{C}_{n}^{\text{spat}}$ and $\mathbf{C}_{n}^{\text{spec}}$ such that
their trace be equal to $K$ and $L$, respectively. In the covariance model of
equation (4), multi-spectral slices at different times $t_{1}$ and $t_{2}$ are
considered uncorrelated (and, thus, independent given the joint Gaussian
assumption of equation (2)). The time-varying variance parameter
$\sigma_{n,t}^{2}$ plays the role of a scale parameter in a compound-Gaussian
model [26] (also known as a Gaussian scale mixture model, see e.g. [27]). A
large value of parameter $\sigma_{n,t}^{2}$ almost discards the time frame $t$
from the $n$-th 4D patch, which limits the impact of possible outliers [28].

The covariance structure given in equation (4) corresponds to the following
separable covariance matrix:

$\displaystyle\text{Cov}\\!\left[\boldsymbol{v}_{n}\right]=\text{diag}(\boldsymbol{\sigma}_{n}^{2})\otimes\mathbf{C}_{n}^{\text{spec}}\otimes\mathbf{C}_{n}^{\text{spat}}\,,$
(5)

where $\text{diag}(\boldsymbol{\sigma}_{n}^{2})$ is a $T\times T$ diagonal
matrix whose $t$-th diagonal entry is $\sigma_{n,t}^{2}$ and $\otimes$ is
Kronecker matrix product: $\mathbf{A}\otimes\mathbf{B}$, with
$\mathbf{A}\in\mathbb{R}^{n\times n}$ and $\mathbf{B}\in\mathbb{R}^{m\times
m}$, is the $nm\times nm$ matrix with a $n\times n$ block structure such that
the $ij$-th block is the $m\times m$ matrix $A_{ij}\mathbf{B}$. Note that this
is equivalent to modeling each multi-spectral slice
$\boldsymbol{v}_{n,t}\in\mathbb{R}^{KL}$ of $\boldsymbol{v}_{n}$ as random
vectors following the compound-Gaussian model
$\mathcal{N}(\boldsymbol{\mu}_{n}^{\text{spec}},\sigma_{n,t}^{2}\mathbf{C}_{n}^{\text{spec}}\otimes\mathbf{C}_{n}^{\text{spat}})$,
or that the scaled and centered vectors
$\frac{1}{\sigma_{n,t}}(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\text{spec}})$
are independent and identically distributed for all $1\leq t\leq T$ according
to the centered Gaussian
$\mathcal{N}(\boldsymbol{0},\mathbf{C}_{n}^{\text{spec}}\otimes\mathbf{C}_{n}^{\text{spat}})$.

With the structure of the mean vector $\boldsymbol{\mu}_{n}$ given in equation
(3), corresponding to a multi-spectral patch constant through time, there are
only $NL$ free parameters to estimate all mean vectors from
$\boldsymbol{v}\in\mathbb{R}^{NLT}$. The covariance structure defined in
equations (4) and (5) leads to $T+K(K+1)/2+L(L+1)/2-2$ free parameters per 4D
patch (the -2 comes from the two normalization constraints), which leads to
approximately $NK/2$ free parameters for the whole ASDI dataset (because $K\gg
L$), and is typically one to two orders of magnitudes smaller than the total
number of measurements in $\boldsymbol{v}$ ($K/2$ is typically less than one
hundred whereas $LT$ is several thousands). Jointly with an adequate
estimation method, the structures assumed in equations (3), (4), and (5) can
thus be used to derive a _non-stationary_ model of the nuisance terms.

###  II-C Estimation of the model parameters $\boldsymbol{\mu}_{n}$ and
$\mathbf{C}_{n}$

The estimation of the parameters of a separable covariance model has been
studied by several previous works, see for example [29, 30, 31]. We build on
these works and introduce several additional elements: (i) whereas most works
consider decompositions of the covariance matrix as a Kronecker product of two
factors, we also include in equations (4) and (5) the temporal scaling factors
$\sigma_{n,t}$ for increased robustness [28]; (ii) given the limited number of
samples, we replace maximum likelihood estimates by shrinkage covariance
estimators [32, 33] to ensure that all estimated covariance matrices are
definite positive and to reduce estimation errors; (iii) to account for the
superimposition of a component of interest and nuisance terms, we develop a
joint estimation strategy in section III-C based on the estimation technique
developed in this paragraph.

Figure 2: Spatial and spectral correlation matrices computed in the four
regions depicted in Fig.1(a). The angular separation with respect to the star
(i.e., image center) increases from the region in (a), close to the star, to
the region in (d), which is farther.

We first derive the maximum likelihood estimates (MLEs) of the means
$\boldsymbol{\mu}_{n}^{\text{spec}}$, variance factors $\sigma_{n,t}^{2}$ and
matrices $\mathbf{C}_{n}^{\text{spec}}$ and $\mathbf{C}_{n}^{\text{spat}}$:

######  Proposition 1 (Maximum Likelihood Parameter Estimates).

The maximum likelihood estimates, under the Gaussian model of section II-B
are:

$\displaystyle\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}$
$\displaystyle=\left(\sum_{t=1}^{T}\widehat{\sigma}_{n,t}^{-2}\boldsymbol{v}_{n,t}\right)\biggl{/}\left(\sum_{t=1}^{T}\widehat{\sigma}_{n,t}^{-2}\right)$
(6) $\displaystyle\widehat{\sigma}_{n,t}^{2}$
$\displaystyle=\tfrac{1}{KL}\left\|\boldsymbol{v}_{n,t}-\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}\right\|_{\widehat{\mathbf{C}}_{n}^{\mathrm{spec}-1}\otimes\widehat{\mathbf{C}}_{n}^{\mathrm{spat}-1}}^{2}$
(7) $\displaystyle\widehat{\mathbf{C}}_{n}^{\mathrm{spec}}$
$\displaystyle=\tfrac{1}{TK}\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\overline{\mathbf{V}}_{t}^{\top}\widehat{\mathbf{C}}_{n}^{\mathrm{spat}-1}\overline{\mathbf{V}}_{t}$
(8) $\displaystyle\widehat{\mathbf{C}}_{n}^{\mathrm{spat}}$
$\displaystyle=\tfrac{1}{TL}\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\overline{\mathbf{V}}_{t}\widehat{\mathbf{C}}_{n}^{\mathrm{spec}-1}\overline{\mathbf{V}}_{t}^{\top}\,,$
(9)

where $\overline{\mathbf{V}}_{t}$ is a $K\times L$ matrix corresponding to the
centered multi-spectral patch at time $t$: at row $k$ and column $\ell$ it is
equal to
$[\boldsymbol{v}_{n,t}]_{k,\ell}-[\boldsymbol{\mu}_{n}^{\mathrm{spec}}]_{k,\ell}$.

###### Proof.

See appendix A. ∎

The multi-spectral mean $\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}$ is
obtained by weighted averaging, with weights inversely proportional to the
variance $\sigma_{n,t}^{2}$: this limits the impact of outliers. The variance
$\sigma_{n,t}^{2}$ corresponds to the average squared deviation to the mean,
computed after proper spatial and spectral whitening. By factoring matrices
$\widehat{\mathbf{C}}_{n}^{\mathrm{spat}-1}$ and
$\widehat{\mathbf{C}}_{n}^{\mathrm{spec}-1}$ under the form
$\mathbb{W}_{n}^{\mathrm{spat}}\mathbb{W}_{n}^{\mathrm{spat}}\,{}^{\top}$ and
$\mathbb{W}_{n}^{\mathrm{spec}}\mathbb{W}_{n}^{\mathrm{spec}}\,{}^{\top}$ (for
example by Cholesky factorization), products of the form
$\mathbb{W}_{n}^{\mathrm{spat}}\,{}^{\top}\,\overline{\mathbf{V}}_{t}$ can be
interpreted as a spatial whitening operation applied to each column of the
centered multi-spectral patch $\overline{\mathbf{V}}_{t}$, and products of the
form
$\mathbb{W}_{n}^{\mathrm{spec}}\,{}^{\top}\,\overline{\mathbf{V}}_{t}^{\top}$
can be interpreted as a spectral whitening applied to each row of
$\overline{\mathbf{V}}_{t}$. Spectral and spatial covariance matrices
$\widehat{\mathbf{C}}_{n}^{\mathrm{spec}}$ and
$\widehat{\mathbf{C}}_{n}^{\mathrm{spat}}$ are thus estimated by forming the
covariance of spatially (resp. spectrally) whitened vectors, with weights
inversely proportional to the variance $\sigma_{n,t}^{2}$.

From the expressions in Proposition 1, it is not possible to derive a closed-
form expression of each parameter that does not depend on other parameters.
Yet, these formulae can be applied alternately until convergence, a method
called _flip-flop_ in [29] where a faster convergence is reported compared to
maximizing the log-likelihood using an iterative optimization algorithm
(Newton’s method).

Figure 2 illustrates the spatial and spectral covariance matrices estimated
from an ASDI dataset captured at the VLT, in Chile. These matrices are
computed from patches extracted at the center of the four regions of the field
of view indicated by color boxes in Fig.1(a). To compare matrices with very
different variances, we normalized each covariance Cov[a,b] by
$\sqrt{\text{Cov[a,a]}\text{Cov[b,b]}}$, i.e., we show the correlation
coefficients. Due to the vectorization of 2D spatial patches, the spatial
correlations display a blocky structure. The spatial correlations within a
patch globally decrease with the 2D distance between pixels. They are stronger
in the area (a) which is the closest to the star. Spectral correlations are
also much stronger close to the star. As can be observed in Fig.1(a), after
the scaling transform applied in the pre-processing step, regions far from the
star are not seen at the longest wavelengths. The size of multi-spectral
patches extracted in these regions is reduced from $KL$ to $KL_{\text{eff}}$
pixels (with $L_{\text{eff}}<L$ the effective number of wavelengths seen at
location $n$) and the size of the spectral covariance matrix
$\widetilde{\mathbf{C}}_{n}^{\mathrm{spec}}$ is reduced accordingly, from
$L\times L$ to $L_{\text{eff}}\times L_{\text{eff}}$.

To reduce the estimation error on the covariance matrices and ensure their
definite-positiveness, shrinkage techniques combine the maximum likelihood
estimator with another estimator of smaller variance. Like in our previous
works [15, 25, 14], we consider the convex combination between the sample
covariance and a diagonal covariance, with a non-constant diagonal. This
corresponds to shrinking off-diagonal entries of the covariance matrix by a
factor $1-\tilde{\rho}$, where $\tilde{\rho}$ is estimated from the data, see
[32, 33, 15]. We recall in the following Proposition the general expression of
a shrinkage covariance estimator:

######  Proposition 2 (Shrinkage Estimator of Covariance).

Let $\\{\boldsymbol{v}_{n,t}\\}_{t=1..T}$ be a collection of $T$ independent
random vectors, such that
$(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})/\sigma_{n,t}$ be i.i.d. for all
$t$ according to a Gaussian distribution
$\mathcal{N}(\boldsymbol{0},\mathbf{C}_{n})$. The bias-variance trade-off is
improved by replacing the maximum likelihood covariance estimator:

$\displaystyle\widehat{\mathbf{C}}_{n}=\tfrac{1}{T}\sum_{t=1}^{T}\frac{1}{{\sigma}_{n,t}^{2}}(\boldsymbol{v}_{n,t}-{\boldsymbol{\mu}}_{n})(\boldsymbol{v}_{n,t}-{\boldsymbol{\mu}}_{n})^{\top}\,,$
(10)

by the shrinkage estimator:

$\displaystyle\widetilde{\mathbf{C}}_{n}=(1-\tilde{\rho}_{n})\widehat{\mathbf{C}}_{n}+\tilde{\rho}_{n}\mathrm{Diag}(\hat{\mathbf{C}}_{n})\,,$
(11)

with $\mathrm{Diag}(\mathbf{C})$ the diagonal matrix such that
$[\mathrm{Diag}(\mathbf{C})]_{i,i}=[{\mathbf{C}}]_{i,i}$ and
$\tilde{\rho}_{n}\in[0,1]$. Parameter $\tilde{\rho}_{n}$ controls the level of
shrinkage applied to off-diagonal values of $\widehat{\mathbf{C}}_{n}$. The
data-driven value for parameter
$\tilde{\rho}_{n}=\rho(\widehat{\mathbf{C}}_{n};\,\widetilde{T}_{n})$ given by
the oracle-approximating method of Chen _et al_ [33] is:

$\displaystyle\rho\left(\mathbf{C};\,T\right)$
$\displaystyle=\frac{\operatorname{tr}({\mathbf{C}}^{2})+\operatorname{tr}^{2}({\mathbf{C}})-2\operatorname{tr}[\mathrm{Diag}(\mathbf{C})^{2}]}{({T}+1)\left(\operatorname{tr}({\mathbf{C}}^{2})-\operatorname{tr}[\mathrm{Diag}(\mathbf{C})^{2}]\right)}\,,$
(12)

where $\widetilde{T}_{n}$ corresponds to the equivalent number of terms in the
weighted sum (10):

$\displaystyle\widetilde{T}_{n}$
$\displaystyle=\frac{\left(\sum_{t}1/\sigma_{n,t}^{2}\right)^{2}}{\sum_{t}1/\sigma_{n,t}^{4}}\,.$
(13)

###### Proof.

The derivation of the shrinkage factor is based on [33] and is given in [15],
appendix A. The equivalent number of terms was first introduced in [28]. ∎

The estimator $\widetilde{\mathbf{C}}_{n}$ defined in Proposition 2 shrinks
off-diagonal values of $\widehat{\mathbf{C}}_{n}$ towards 0 (by multiplication
by the factor $1-\tilde{\rho}_{n}$) and leaves diagonal values unchanged
(i.e., the variances). Beyond reducing the estimation error, the shrinkage
estimator also ensures that $\widetilde{\mathbf{C}}_{n}$ is positive definite,
which is essential in section III since the cost function that is minimized to
perform image reconstruction involves the inverse of the covariance matrix.
Compared to other covariance matrix regularization strategies such as diagonal
loading (i.e., adding a fraction of the identity matrix $\epsilon\mathbf{I}$
to $\widehat{\mathbf{C}}_{n}$), Proposition 2 is attractive because it is
data-driven: it adapts to the fluctuations observed in the data and the number
of samples. It is thus well-suited to imaging systems suffering from non-
stationary perturbations.

To apply Proposition 2, the true mean $\boldsymbol{\mu}_{n}$ and scale factors
$\sigma_{n,t}$ are replaced by their estimates, which requires that the
estimation error on these parameters be low. The expression of
$\widetilde{\mathbf{C}}_{n}$ is general, it does not include the specific
structure of our problem, namely, the spatio-spectral separability. Before
specializing the expression of the shrinkage covariance estimator to this
structure, we give in Proposition 3 an equivalent definition of the shrinkage
estimator.

######  Proposition 3 (Loss function associated to the shrinkage covariance
estimator).

The shrinkage covariance matrix estimator $\widetilde{\mathbf{C}}_{n}$ defined
in Proposition 2 is the minimizer of a modified Gaussian co-log-likelihood:

$\displaystyle\widetilde{\mathbf{C}}_{n}=\operatorname*{arg\,min}_{\mathbf{C}}\;\frac{T}{2}\log|\mathbf{C}|+\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}\|_{\mathbf{\Gamma}}^{2}\,,$
(14)

with $\mathbf{\Gamma}=\mathbf{\Psi}_{n}\odot(\mathbf{C}^{-1})$ the shrunk
precision matrix, defined as the Hadamard (element-wise) product between a
weighting matrix $\mathbf{\Psi}_{n}$ and the precision matrix
$\mathbf{C}^{-1}$. The weighting matrix $\mathbf{\Psi}_{n}$ has 1 on the
diagonal and the value $(1-\tilde{\rho}_{n})$ off the diagonal
($\tilde{\rho}_{n}$ is defined in Prop.2).

###### Proof.

This is a generalization of Eq.(30) in [15], the proof is given in appendix B.
∎

Note that the expression of the shrinkage estimator in Eq.(11) can be more
compactly written
$\widetilde{\mathbf{C}}_{n}=\mathbf{\Psi}_{n}\odot\widehat{\mathbf{C}}_{n}$.
Proposition 3 shows how to modify the co-log-likelihood so that the covariance
matrix that minimizes the loss function changes from the maximum likelihood
estimator $\widehat{\mathbf{C}}_{n}$ to the shrinkage estimator
$\widetilde{\mathbf{C}}_{n}$.

Based on Proposition 3, we introduce a loss function suitable for separable
covariance matrices:

######  Proposition 4 (Loss function for separable covariance matrix
estimation).

Let us define the loss function $\mathscr{D}_{n}$:

$\mathscr{D}_{n}(\boldsymbol{v}_{n};\,\boldsymbol{\mu}_{n}^{\mathrm{spec}},\boldsymbol{\sigma}_{n},\mathbf{C}_{n}^{\mathrm{spat}},\mathbf{C}_{n}^{\mathrm{spec}})=\\\
\tfrac{T}{2}\log|\mathrm{diag}(\boldsymbol{\sigma}_{n}^{2})\otimes\mathbf{C}_{n}^{\mathrm{spec}}\otimes\mathbf{C}_{n}^{\mathrm{spat}}|+\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\mathrm{spec}}\right\|_{\mathbf{\Gamma}}^{2},$
(15)

where
$\mathbf{\Gamma}={\left[\left(\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\mathbf{C}_{n}^{\mathrm{spec}}\right)\otimes\left(\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\mathbf{C}_{n}^{\mathrm{spat}}\right)\right]^{-1}}$
is the shrunk spatio-spectral precision matrix. We define the weighting
matrices like in Proposition 3: matrix $\mathbf{\Psi}_{n}^{\mathrm{spec}}$ is
the $L\times L$ matrix with 1 on the diagonal and the value
$(1-\tilde{\rho}_{n}^{\mathrm{spec}})$ off-diagonal, and matrix
$\mathbf{\Psi}_{n}^{\mathrm{spat}}$ is the $K\times K$ matrix with 1 on the
diagonal and the off-diagonal value $(1-\tilde{\rho}_{n}^{\mathrm{spat}})$.
Shrinkage parameters are computed with:

$\displaystyle\tilde{\rho}_{n}^{\mathrm{spec}}$
$\displaystyle=\rho(\widehat{\mathbf{C}}_{n}^{\mathrm{spec}};\,\widetilde{P}_{n}^{\mathrm{spec}})\,,$
(16) $\displaystyle\tilde{\rho}_{n}^{\mathrm{spat}}$
$\displaystyle=\rho(\widehat{\mathbf{C}}_{n}^{\mathrm{spat}};\,\widetilde{P}_{n}^{\mathrm{spat}})\,,$
(17) and the equivalent number of terms defined by:
$\displaystyle\widetilde{P}_{n}^{\mathrm{spec}}$
$\displaystyle=K\frac{\left(\sum_{t}1/\hat{\sigma}_{n,t}^{2}\right)^{2}}{\sum_{t}1/\hat{\sigma}_{n,t}^{4}}\text{
and
}\widetilde{P}_{n}^{\mathrm{spat}}=L\frac{\left(\sum_{t}1/\hat{\sigma}_{n,t}^{2}\right)^{2}}{\sum_{t}1/\hat{\sigma}_{n,t}^{4}}\,.$
(18)

The set of parameters
$\\{\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}},\,\widehat{\boldsymbol{\sigma}}_{n}^{2},\,\widetilde{\mathbf{C}}_{n}^{\mathrm{spec}},\,\widetilde{\mathbf{C}}_{n}^{\mathrm{spat}}\\}$
that minimize $\mathscr{D}_{n}$, for fixed values of the shrinkage parameters
$\tilde{\rho}_{n}^{\mathrm{spec}}$ and $\tilde{\rho}_{n}^{\mathrm{spat}}$,
verifies the following equations:

$\displaystyle\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}$
$\displaystyle=\left(\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\boldsymbol{v}_{n,t}\right)\biggl{/}\left(\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\right)$
(19) $\displaystyle\widehat{\sigma}_{n,t}^{2}$
$\displaystyle=\tfrac{1}{KL}\left\|\boldsymbol{v}_{n,t}-\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}\right\|_{\mathbf{\Gamma}}^{2}$
(20) $\displaystyle\widetilde{\mathbf{C}}_{n}^{\mathrm{spec}}$
$\displaystyle=\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\tfrac{1}{TK}\sum_{t=1}^{T}\tfrac{1}{\widehat{\sigma}_{n,t}^{2}}\overline{\mathbf{V}}_{t}^{\top}(\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\widetilde{\mathbf{C}}_{n}^{\mathrm{spat}-1})\overline{\mathbf{V}}_{t}$
(21) $\displaystyle\widetilde{\mathbf{C}}_{n}^{\mathrm{spat}}$
$\displaystyle=\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\tfrac{1}{TL}\sum_{t=1}^{T}\tfrac{1}{\widehat{\sigma}_{n,t}^{2}}\overline{\mathbf{V}}_{t}(\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\widetilde{\mathbf{C}}_{n}^{\mathrm{spec}-1})\overline{\mathbf{V}}_{t}^{\top}$
(22)

where $\overline{\mathbf{V}}_{t}$ is defined as in Proposition 1.

###### Proof.

See appendix C. ∎

The loss function $\mathscr{D}_{n}$ is a direct application of Proposition 3
to the structured covariance matrix defined in Eq.(5). As in Proposition 1,
Eqs.(19)–(22) form a set of non-linear equations that can be applied
iteratively until a fixed point is reached.

##  III Reconstruction of the component of interest

###  III-A Direct model

We extend the forward model developed for ADI in [14] by including the
spectral dimension. Since the whole ASDI sequence is acquired within a short
time (a few hours during a single night), the component that contains the
circumstellar sources such as the disk and possibly exoplanets _does not
evolve_ during the observations (rotation of the disk or exoplanets about
their host star is negligible at such short time scales). The multi-spectral
image of this component is simply described by the vector
$\boldsymbol{u}\in\mathbb{R}_{+}^{N^{\prime}L}$ of its pixel values (the
number $N^{\prime}$ of pixels in each spectral band of the reconstruction may
differ from $N$, note that there is no temporal dimension in the
reconstruction). The contribution of $\boldsymbol{u}$ to the data
$\boldsymbol{v}$ is modeled by the linear operator:

$\displaystyle\mathbf{M}\boldsymbol{u}\,,\text{ with
}\mathbf{M}=\begin{pmatrix}\mathbf{M}_{1}\\\ \vdots\\\ \mathbf{M}_{T}\\\
\end{pmatrix}\text{ and
}\mathbf{M}_{t}=\mathbf{S}\mathbf{Z}\mathbf{A}\mathbf{R}_{t}\mathbf{B}_{t}\,,$
(23)

where $\mathbf{M}_{t}$, the model for the $t$-th frame, accounts for several
instrumental effects:

  * •

a _blur_ $\mathbf{B}_{t}$ due to the instrumental blurring modeled as a 2D
discrete convolution by the off-axis point spread function (PSF),

  * •

a _rotation_ $\mathbf{R}_{t}$ applied to all off-axis sources due to the
pupil-tracking mode (the field-of-view rotates while the residual star light
remains fixed), implemented as a sparse interpolation matrix,

  * •

an _attenuation_ $\mathbf{A}$, very strong on the optical axis, then quickly
decreasing (due to the coronagraph), modeled as a diagonal matrix [14],

  * •

the absence of measurements outside the spatial extension of the sensor (a
non-square area due to the instrumental design of the integral field
spectrograph), modeled as a diagonal matrix $\mathbf{Z}$ that replaces values
outside the sensor area by _zeros_ and keeps other values unchanged.

  * •

the image _scaling_ applied during the pre-processing step produces a last
transform $\mathbf{S}$ (time-invariant), corresponding to a sparse
interpolation matrix.

With the SPHERE-IFS instrument, the off-axis point spread function (PSF) is
quite stable and its core is almost rotation invariant. The model given in
equation (23) can thus be approximated by:

$\displaystyle\mathbf{M}\boldsymbol{u}\approx\begin{pmatrix}\mathbf{F}_{1}\\\
\vdots\\\ \mathbf{F}_{T}\\\ \end{pmatrix}\\!\mathbf{B}\boldsymbol{u}\,.$ (24)

where $\mathbf{B}$ is a time-invariant blurring operator and
$\\{\mathbf{F}_{t}\\}_{t=1..T}$ are sparse matrices that perform rotations,
scalings, and attenuations according to the transmission of the coronagraph
and the sensor field-of-view. The model of equation (24) is only approximate:
it neglects possible anisotropies or temporal evolutions of the PSF. Thanks to
this approximation, a single convolution of the multi-spectral cube is
performed instead of $T$ convolutions, which leads to a dramatic acceleration
of the numerical evaluation of the forward model (by 1 to 2 orders of
magnitude) which is critical to achieving reconstructions on datasets in the
order of a few hours. We checked that the impact of these approximations on
the reconstructions was negligible in practice.

###  III-B Regularized inversion

We reconstruct the component of interest using a penalized maximum likelihood
approach, i.e., by solving the following numerical optimization problem:

$\displaystyle\widehat{\boldsymbol{u}}=\operatorname*{arg\,min}_{\boldsymbol{u}\geq\mathbf{0}}\;-\log\text{p}_{V}\bigl{(}\,\boldsymbol{v}\,|\,\mathbf{M}\boldsymbol{u};\,\\{\boldsymbol{\mu}_{n},\mathbf{C}_{n}\\}_{n\in\mathbb{K}}\,\bigr{)}+\mathscr{R}(\boldsymbol{u})\,,$
(25)

where $\mathscr{R}(\boldsymbol{u})$ is a regularization term to favor
plausible reconstructions $\boldsymbol{u}$. We selected a combination of two
regularization functions: an edge-preserving one that favors smooth images
with sharp edges co-located at all wavelengths and a sparsity-inducing
$\text{L}^{1}$ norm. The regularization writes:

$\mathscr{R}(\boldsymbol{u})=\beta_{\text{smooth}}\sum_{n=1}^{N^{\prime}}\sqrt{\tfrac{1}{L}\sum_{\ell=1}^{L}\|\mathbf{\nabla}_{\\!n}\boldsymbol{u}_{\cdot,\ell}\|_{2}^{2}+\epsilon^{2}}\\\
+\beta_{\text{sparse}}\sum_{n=1}^{N^{\prime}}\sum_{\ell=1}^{L}|u_{n,\ell}|\,,$
(26)

with $\mathbf{\nabla}_{\\!n}\boldsymbol{u}_{\cdot,\ell}$ the 2D spatial
gradient at pixel $n$ of the $\ell$-th spectral channel, evaluated by finite
differences, and $\epsilon$ a parameter chosen so as to be negligible compared
to the average norm of the spatial gradient where there is a sharp edge (the
regularization then approaches an isotropic vectorial total variation [34])
and similar to the gradient magnitude in smoothly-varying areas (this prevents
the apparition of the staircasing effect common with total variation). Hyper-
parameters $\beta_{\text{smooth}}$ and $\beta_{\text{sparse}}$ balance the
weight of each regularization term with respect to the data-fitting term. Note
that, due to the positivity constraint in equation (25), the $\text{L}^{1}$
norm $\|\boldsymbol{u}\|_{1}$ corresponds to the simple differentiable term
$\sum_{n=1}^{N^{\prime}}\sum_{\ell=1}^{L}u_{n,\ell}$ for any feasible object
$\boldsymbol{u}$, and the regularization $\mathscr{R}(\boldsymbol{u})$ is
differentiable.

###  III-C Joint estimation of the statistics of nuisance terms and
reconstruction of the component of interest

The estimation of the parameters of the statistical model for nuisance terms
should be performed on data free of the component of interest, i.e., on the
residuals $\boldsymbol{v}-\mathbf{M}\boldsymbol{u}$. Starting from an initial
guess $\widehat{\boldsymbol{\mu}}_{n}^{\text{spec }(0)}$,
$\widehat{\boldsymbol{\sigma}}_{n}^{2\;(0)}$,
$\widetilde{\mathbf{C}}_{n}^{\text{spec }(0)}$, and
$\widetilde{\mathbf{C}}_{n}^{\text{spat }(0)}$ computed under the assumption
$\boldsymbol{u}=0$ (i.e., directly on the data), reconstruction and estimation
steps can be alternated to successively refine the contribution of the objects
and the background statistics. In order to speed up the process, we consider a
joint reconstruction and estimation approach. Since in section II-C we derived
our estimators from the loss function $\mathscr{D}_{n}$ (Proposition 4), we
can perform jointly the reconstruction and the estimation of the parameters of
the nuisance model by minimizing a common cost function:

$\displaystyle\widehat{\boldsymbol{u}}=\operatorname*{arg\,min}_{\boldsymbol{u}\geq\mathbf{0}}\\!\\!\underbrace{\min_{\begin{array}[]{c}\scriptstyle\\{\widehat{\boldsymbol{\mu}}_{n}^{\text{spec}}\\}_{n\in\mathbb{K}}\\\
\scriptstyle\\{\widehat{\boldsymbol{\sigma}}_{n}^{2}\\}_{n\in\mathbb{K}}\\\
\scriptstyle\\{\widetilde{\mathbf{C}}_{n}^{\text{spec}}\\}_{n\in\mathbb{K}}\\\
\scriptstyle\\{\widetilde{\mathbf{C}}_{n}^{\text{spat}}\\}_{n\in\mathbb{K}}\end{array}}\sum_{n\in\mathbb{K}}\mathscr{D}_{n}(\boldsymbol{v}_{n}-[\mathbf{M}\boldsymbol{u}]_{n})+\mathscr{R}(\boldsymbol{u})}_{\mathscr{C}(\boldsymbol{u})}\,,$
(31)

where the dependence of $\mathscr{D}_{n}$ on parameters
$\widehat{\boldsymbol{\mu}}_{n}^{\text{spec}}$,
$\widehat{\boldsymbol{\sigma}}_{n}^{2}$,
$\widetilde{\mathbf{C}}_{n}^{\text{spec}}$, and
$\widetilde{\mathbf{C}}_{n}^{\text{spat}}$ is not explicitly written. In the
expression of $\mathscr{D}_{n}$ given in Proposition 4, the spatio-spectro-
temporal patch $\boldsymbol{v}_{n}$ is replaced by the residual patch
$\boldsymbol{v}_{n}-[\mathbf{M}\boldsymbol{u}]_{n}$ obtained when removing the
contribution of the component of interest.

To solve this hierarchical minimization problem, we suggest using a large-
scale continuous optimization technique to solve the outer minimization
problem and the fixed-point method based on Eqs.(19)–(22) to solve the inner
problem. The following proposition gives the expression of the gradient of
$\mathscr{C}(\boldsymbol{u})$.

######  Proposition 5 (Gradient of the cost function with respect to
reconstructed component $\boldsymbol{u}$).

The gradient of the cost function $\mathscr{C}(\boldsymbol{u})$ with respect
to $\boldsymbol{u}$ is:

$\nabla_{\boldsymbol{u}}\mathscr{C}(\boldsymbol{u})=\mathbf{M}^{\top}\sum_{n\in\mathbb{K}}\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}(\boldsymbol{u})}\mathbf{P}_{n,t}^{\top}\widehat{\mathbf{\Gamma}}_{n}(\boldsymbol{u})\bigl{[}\mathbf{P}_{n,t}\mathbf{M}\boldsymbol{u}\hfill\\\\[-8.61108pt]
\hfill-(\boldsymbol{v}_{n,t}-\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}(\boldsymbol{u}))\bigr{]}\\\\[-4.30554pt]
+\beta_{\mathrm{smooth}}\sum_{n=1}^{N^{\prime}}\frac{1}{2\sqrt{\sum_{\ell=1}^{L}\|\mathbf{D}_{\\!n,\ell}\boldsymbol{u}\|_{2}^{2}+\epsilon^{2}}}\mathbf{D}_{\\!n,\ell}^{\top}\mathbf{D}_{\\!n,\ell}\boldsymbol{u}+\beta_{\mathrm{sparse}}\boldsymbol{1}\\\\[-30.1388pt]
$ (32)

where $\mathbf{P}_{n,t}$ is the $KL\times KLT$ matrix that extracts a multi-
spectral patch at spatial location $n$ and time $t$,
$\widehat{\mathbf{\Gamma}}_{n}(\boldsymbol{u})$ is equal to
$\left(\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\widetilde{\mathbf{C}}_{n}^{\mathrm{spec}-1}(\boldsymbol{u})\right)\otimes\left(\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\widetilde{\mathbf{C}}_{n}^{\mathrm{spat}-1}(\boldsymbol{u})\right)$,
and $\mathbf{D}_{n,\ell}$ is the $2\times N^{\prime}L$ matrix that computes
the horizontal and vertical finite differences at pixel $n$ for the wavelength
$\ell$ of the object.

###### Proof.

The expression follows directly from differentiation rules by noting that the
differentials with respect to parameters
$\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}(\boldsymbol{u})$,
$\widehat{\boldsymbol{\sigma}}_{n}^{2}(\boldsymbol{u})$,
$\widetilde{\mathbf{C}}_{n}^{\mathrm{spec}}(\boldsymbol{u})$,
$\widetilde{\mathbf{C}}_{n}^{\mathrm{spat}}(\boldsymbol{u})$ are all null due
to Proposition 4. When applying the chain rule, the terms implying
differentials of $\widehat{\sigma}_{n,t}^{2}(\boldsymbol{u})$,
$\widehat{\mathbf{\Gamma}}(\boldsymbol{u})$, or
$\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}(\boldsymbol{u})$ thus cancel
out and only the terms where these parameters are considered fixed remain. ∎

To solve the smooth optimization problem under bound constraints of Eq.(31),
we use a limited-memory quasi-Newton method with bound constraints, VMLM-B
[35], which is a more efficient variant of L-BFGS-B [36]. This algorithm
evaluates the cost function $\mathscr{C}(\boldsymbol{u})$ and the gradient
$\nabla_{\boldsymbol{u}}\mathscr{C}(\boldsymbol{u})$ given in Proposition 5
iteratively until reaching convergence to a local optimum.

##  IV Results

Figure 3: Empirical distribution of the centered and whitened patches, for
different covariance models.

###  IV-A Validation of the statistical model of the nuisance component

Before evaluating the reconstruction method on high-contrast observations of
circumstellar disks, we show that our statistical model of the nuisances is
relevant. We consider the same ASDI dataset as in [25] (star HR 8799, with its
three known exoplanets masked out so that the data correspond only to the
nuisance term) and perform a similar analysis: we represent in Fig.3 the
empirical distribution of the patches after performing different pre-
processing. If random vectors $\boldsymbol{v}_{n}$ are accurately modeled by a
Gaussian distribution with mean $\boldsymbol{\mu}_{n}$ and covariance
$\mathbf{C}_{n}$, as in equation (2), then the centered and whitened vectors
$\mathbf{C}_{n}^{-1/2}(\boldsymbol{v}_{n}-\boldsymbol{\mu}_{n})$ should follow
$\mathcal{N}(\boldsymbol{0},\mathbf{I})$. Therefore, we compare a standard
Gaussian distribution with the empirical marginal distribution of
$\mathbf{C}_{n}^{-1/2}(\boldsymbol{v}_{n}-\boldsymbol{\mu}_{n})$ for several
covariance models: the three models considered in [25], drawn in gray dashed-
lines: (i) no covariance ($\mathbf{C}_{n}=\mathbf{I}$); (ii) only spatial
covariances; (iii) spatial covariances + temporal and spectral weighting; and
four additional models: (iv) diagonal spatial and spectral covariances (i.e.,
spatial, spectral, and temporal weighting via a separable model); (v) full
spatial covariance, diagonal spectral covariance, and temporal weighting; (vi)
full spectral covariance, diagonal spatial covariance, and temporal weighting;
and finally (vii) the full separable model introduced in this paper (see
equation (4)). The centered and whitened patches are best modeled with the
full model (the green curve is the closest to the red dashed-line
corresponding to the standard Gaussian distribution), which justifies its use
in our loss function $\mathscr{D}_{n}$. We give a more detailed study of the
empirical distribution of patches in the Supplementary Material associated
with this paper.

Figure 4: Spatial distribution of the spatial and spectral shrinkage
parameters.

The shrinkage parameters $\rho_{n}^{\text{spat}}$ and $\rho_{n}^{\text{spec}}$
may have a significant impact on the statistical model. We monitor their
values by displaying maps of the spatial and spectral shrinkage parameters in
Fig.4 for the dataset shown in Fig.1. Values of $\rho_{n}^{\text{spat}}$ and
$\rho_{n}^{\text{spec}}$ remain small (below 0.13), which indicates a moderate
bias toward 0. The spatial shrinkage is stronger at the borders of the field-
of-view where there are fewer samples available due to the spectral scaling.
Spectral shinkage is stronger for this dataset at some locations of the field-
of-view which shows that setting a global shrinkage value for the whole field-
of-view would be sub-optimal.

###  IV-B Reconstruction of disks from VLT/SPHERE-IFS data

Now, we consider two datasets resulting from the observation of two different
stars hosting a known circumstellar disk: SAO 206462 (already shown in Fig.1,
observation date: 2015-05-15, ESO program: 095.C-0298) which contains a spiral
disk and HR 4796 (observation date: 2015-02-02, ESO program: 095.C-0298) with
a thin elliptical disk.

Figure 5: Ablation study: neglecting the spatial and/or spectral correlations
of the nuisance component negatively affects the reconstruction results. The
images shown correspond to the deconvolved reconstructions
$\hat{\boldsymbol{u}}$. Fake RGB images are shown with colors selected to
cover the infrared spectrum according to the colormaps on the right.

We first perform an ablation study by showing how the reconstruction results
are impacted if simpler covariance models are considered rather than the full
model of equation (4). Fig.5 displays the reconstructed disk component in four
cases, from left to right: (i) no modeling of spatial or spectral
correlations; (ii) modeling only of spatial correlations; (iii) modeling only
of spectral correlations; and (iv) the full model proposed in this paper. The
dataset HR 4796 suffers from a strong nuisance component. Failing to account
for the spatial correlations leads to severe artifacts: a ghost circular
structure is reconstructed. By closely inspecting the center of the field-of-
view of SAO 206462, the presence of spurious structures can be noted in all
reconstructions except the reconstruction with the full model. These
qualitative assessments indicate that the modeling of correlations is
beneficial to the reconstructions.

Figure 6: Disk reconstruction results. The images shown for the method REXPACO
ASDI in the fourth column are re-blurred to be easier to compare with the
reference methods (i.e., they correspond to $\mathbf{B}\hat{\boldsymbol{u}}$).
Deblurred results $\hat{\boldsymbol{u}}$ are shown on the last column (they
correspond to the last column of Fig.5). Fake colors are used as in the
previous figure.

Fig.6 shows reconstructions obtained by several reference methods along with
reconstructions obtained with our method. Since other methods do not perform a
deconvolution, we show in the fourth column of Fig.6 our reconstruction re-
blurred at the resolution of the instrument (i.e.,
$\mathbf{B}\hat{\boldsymbol{u}}$ rather than $\hat{\boldsymbol{u}}$). Three
methods were considered as reference: cADI [2] which estimates the nuisance
component by temporally stacking the observations (using medians), PCA ADI [7]
which uses a principal component analysis to remove the nuisance component
(the best number of modes has been manually selected by trial and error), and
PACO ASDI [25] which was originally designed to extract exoplanets from ASDI
datasets but can also partially reconstruct thin disks. For cADI and PCA ADI,
we use the official implementations of the SPHERE data center [37, 38], while
we use our unsupervised pipeline of PACO ASDI [25]. The reconstructions
obtained with the reference methods all suffer from artifacts, in particular
in the center of the field-of-view. In comparison, the blurred and the
deblurred reconstructions shown in the last two columns of Fig. 6 are much
more satisfactory: the structure of the disks is in good match with standard
models of circumstellar disks [39, 40]. Besides, the REXPACO ASDI
reconstruction of HR 4796 (respectively, SAO 206462) can be compared
qualitatively with the reconstruction given in Fig.4 (respectively, Fig.1
bottom-left) of [39] (respectively, of [40]) obtained using customed routines
of cADI on the same datasets considered in this paper. The REXPACO ASDI
reconstructions exhibits significantly less artifacts such as non-physical
discontinuities in the disk structures and residuals stellar leakages near the
star. The deconvolution step of the proposed method also leads to an improved
resolution of the thin structures of the disks. Reference methods cADI and PCA
ADI subtract part of the disk component when removing the nuisance term. This
leads to strong flux biases and a high-pass filtering effect. PACO ASDI is
designed to detect point-like structures, it recovers parts of the disks in
large gradient areas: it is much more successful on the thin disk of HR 4796
than on the thicker spiral disk of SAO 206462. Finally, the multi-spectral
REXPACO ASDI reconstructions of Fig.6 can be compared with the mono-spectral
reconstructions obtained with the REXPACO algorithm [14] (see Fig.11 of [14])
on mono-spectral datasets of the same target stars. These mono-spectral
datasets were recorded by the InfraRed Dual Imaging Spectrograph (IRDIS) of
the SPHERE instrument operating simultaneously with the IFS but in a different
spectral band and resolution. A joint processing of multi-spectral datasets is
particularly relevant for disks having a circular symmetry, like SAO 206462,
since it helps to unmix the light coming from the disk and the star light
(since the two components not always superimposed due to the chromatic scaling
of speckles). This joint multi-processing leads to a better rejection of the
noise thus limiting the non-physical reconstruction artifacts like
discontinuities in the spiral arms.

##  V Conclusion

Correlations of the noise (or of a nuisance component) are generally neglected
when solving an image reconstruction problem. If not, a stationarity
assumption is typically used to model and compensate for these correlations.
In this paper, we considered an application in which the nuisance component
has a strong non-stationarity. We resorted to a local modeling of the spatial
and spectral correlations to achieve a more accurate statistical description
of the data and obtain reconstructions without noticeable artifacts. Our
modeling of these correlations is based on a separable approximation to reduce
the number of free parameters. We developed estimators based on shrinkage
covariance estimators and a compound-Gaussian model to reduce estimation
variance and improve robustness. These mean and covariance estimators of our
statistical model of nuisances are evaluated jointly to the reconstruction of
the component of interest. This leads to an algorithm that scales to the size
of typical ASDI datasets.

We have shown on several datasets obtained with the SPHERE-IFS instrument at
the VLT, in Chile, that much improved reconstructions of circumstellar disks
could be obtained with our method. This will help astrophycists better
understand the phenomena that govern the formation of planetary systems and
the interactions between exoplanets and dust disks. Beyond its use in
astronomy, we think that this general approach of modeling spatial and
spectral correlations in an image reconstruction problem could also be applied
to other fields of imaging such as microscopy.

##  Appendix A Proof of the MLEs given in Proposition 1

Under the assumptions of section II-B, the co-log-likelihood of the 4D patch
$\boldsymbol{v}_{n}$, in the absence of disk or exoplanet, writes:

$\displaystyle\mathscr{L}$
$\displaystyle=-\log\text{p}_{V_{n}}\\!(\boldsymbol{v}_{n})$
$\displaystyle=\tfrac{TK}{2}\log|\mathbf{C}_{n}^{\text{spec}}|+\tfrac{TL}{2}\log|\mathbf{C}_{n}^{\text{spat}}|+\sum_{t=1}^{T}\biggl{(}\tfrac{KL}{2}\log\sigma_{n,t}^{2}+$
$\displaystyle\hskip
50.00008pt+\frac{1}{2\sigma_{n,t}^{2}}\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\text{spec}}\right\|_{\left(\mathbf{C}_{n}^{\text{spec}}\otimes\mathbf{C}_{n}^{\text{spat}}\right)^{-1}}^{2}\biggr{)}\,,$
(33)

where we used the following properties of the Kronecker product:
$|\mathbf{A}\otimes\mathbf{B}|=|\mathbf{A}|^{m}|\mathbf{B}|^{n}$ (with
$\mathbf{A}$ an $n\times n$ matrix and $\mathbf{B}$ an $m\times m$ matrix) and
$(\mathbf{A}\otimes\mathbf{B})^{-1}=\mathbf{A}^{-1}\otimes\mathbf{B}^{-1}$.

To obtain the MLEs, we differentiate $\mathscr{L}$:

$\partial\mathscr{L}=\tfrac{TK}{2}\operatorname{tr}(\mathbf{C}_{n}^{\text{spec
}-1}\partial\mathbf{C}_{n}^{\text{spec}})+\tfrac{TL}{2}\operatorname{tr}(\mathbf{C}_{n}^{\text{spat
}-1}\partial\mathbf{C}_{n}^{\text{spat}})\\\
+\sum_{t=1}^{T}\biggl{\\{}\tfrac{KL}{2}\frac{\partial\sigma_{n,t}^{2}}{\sigma_{n,t}^{2}}-\frac{\partial\sigma_{n,t}^{2}}{2\sigma_{n,t}^{4}}\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\text{spec}}\right\|_{\mathbf{C}_{n}^{\text{spec
}-1}\otimes\mathbf{C}_{n}^{\text{spat }-1}}^{2}\\\
+\frac{1}{\sigma_{n,t}^{2}}\\!\left(\boldsymbol{\mu}_{n}^{\text{spec}}-\boldsymbol{v}_{n,t}\right)^{\top}\\!\left(\mathbf{C}_{n}^{\text{spec
}-1}\otimes\mathbf{C}_{n}^{\text{spat
}-1}\right)\\!\partial\boldsymbol{\mu}_{n}^{\text{spec}}\\\
-\frac{1}{2\sigma_{n,t}^{2}}\operatorname{tr}\\!\left(\mathbf{C}_{n}^{\text{spec
}-1}\partial\mathbf{C}_{n}^{\text{spec}}\mathbf{C}_{n}^{\text{spec
}-1}\overline{\mathbf{V}}_{t}^{\top}\mathbf{C}_{n}^{\text{spat
}-1}\overline{\mathbf{V}}_{t}\right)\\\
-\frac{1}{2\sigma_{n,t}^{2}}\operatorname{tr}\\!\left(\mathbf{C}_{n}^{\text{spat
}-1}\partial\mathbf{C}_{n}^{\text{spat}}\mathbf{C}_{n}^{\text{spat
}-1}\overline{\mathbf{V}}_{t}\mathbf{C}_{n}^{\text{spec
}-1}\overline{\mathbf{V}}_{t}^{\top}\right)\biggr{\\}}$ (34)

where we obtained the last two terms by rewriting the squared norm
$\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\text{spec}}\|_{\mathbf{C}_{n}^{\text{spec
}-1}\otimes\mathbf{C}_{n}^{\text{spat }-1}}^{2}$ under the form
$\operatorname{tr}(\overline{\mathbf{V}}_{t}^{\top}\mathbf{C}_{n}^{\text{spat
}-1}\overline{\mathbf{V}}_{t}\mathbf{C}_{n}^{\text{spec }-1})$, with
$\overline{\mathbf{V}}_{t}$ the $K\times L$ matrix whose element at row $k$
and column $\ell$ is
$[\boldsymbol{v}_{n,t}]_{k,\ell}-[\boldsymbol{\mu}_{n}^{\text{spec}}]_{k,\ell}$.

The following set of conditions is sufficient for $\partial\mathscr{L}$ to be
equal to zero:

$\displaystyle\sum_{t=1}^{T}\frac{1}{\sigma_{n,t}^{2}}\\!\left(\boldsymbol{\mu}_{n}^{\text{spec}}-\boldsymbol{v}_{n,t}\right)=0$
(35a)
$\displaystyle\tfrac{KL}{2}\frac{1}{\sigma_{n,t}^{2}}-\frac{1}{2\sigma_{n,t}^{4}}\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\text{spec}}\right\|_{\mathbf{C}_{n}^{\text{spec
}-1}\otimes\mathbf{C}_{n}^{\text{spat }-1}}^{2}=0$ (35b)
$\displaystyle\tfrac{TK}{2}\mathbf{I}-\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}\mathbf{C}_{n}^{\text{spec
}-1}\overline{\mathbf{V}}_{t}^{\top}\mathbf{C}_{n}^{\text{spat
}-1}\overline{\mathbf{V}}_{t}=\mathbf{0}$ (35c)
$\displaystyle\tfrac{TL}{2}\mathbf{I}-\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}\mathbf{C}_{n}^{\text{spat
}-1}\overline{\mathbf{V}}_{t}\mathbf{C}_{n}^{\text{spec
}-1}\overline{\mathbf{V}}_{t}^{\top}=\mathbf{0}$ (35d)

Condition (35a) holds if:

$\displaystyle\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}$
$\displaystyle=\left(\sum_{t=1}^{T}\sigma_{n,t}^{-2}\boldsymbol{v}_{n,t}\right)\biggl{/}\left(\sum_{t=1}^{T}\sigma_{n,t}^{-2}\right).$
(36) Condition (35b) holds if: $\displaystyle\widehat{\sigma}_{n,t}^{2}$
$\displaystyle=\tfrac{1}{KL}\left\|\boldsymbol{v}_{n,t}-\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}\right\|_{\widehat{\mathbf{C}}_{n}^{\mathrm{spec}-1}\otimes\widehat{\mathbf{C}}_{n}^{\mathrm{spat}-1}}^{2}\,.$
(37) Condition (35c) holds if:
$\displaystyle\widehat{\mathbf{C}}_{n}^{\mathrm{spec}}$
$\displaystyle=\tfrac{1}{TK}\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\overline{\mathbf{V}}_{t}^{\top}\widehat{\mathbf{C}}_{n}^{\mathrm{spat}-1}\overline{\mathbf{V}}_{t}\,.$
(38) Condition (35d) holds if:
$\displaystyle\widehat{\mathbf{C}}_{n}^{\mathrm{spat}}$
$\displaystyle=\tfrac{1}{TL}\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\overline{\mathbf{V}}_{t}\widehat{\mathbf{C}}_{n}^{\mathrm{spec}-1}\overline{\mathbf{V}}_{t}^{\top}\,.$
(39)

These correspond to the expressions given in Proposition 1.

##  Appendix B Proof of Proposition 3

We start by rewriting the $\text{L}^{2}$ norm under the form of a trace:

$\displaystyle\sum_{t=1}^{T}$
$\displaystyle\frac{1}{2\sigma_{n,t}^{2}}\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}\|_{\mathbf{\Gamma}}^{2}$
$\displaystyle=\operatorname{tr}\\!\left(\\!\bigl{(}\mathbf{\Psi}_{n}\odot(\mathbf{C}^{-1})\bigr{)}\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})^{\top}\\!\right)$
(40)
$\displaystyle=\operatorname{tr}\\!\left[\mathbf{C}^{-1}\\!\left(\mathbf{\Psi}_{n}\odot\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})^{\top}\right)\right]\\!.$
(41)

The differential of the cost function introduced in Proposition 3, when all
other parameters are fixed, writes:

$\partial\\!\left[\frac{T}{2}\log|\mathbf{C}|+\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}\|_{\mathbf{\Gamma}}^{2}\right]=\frac{T}{2}\operatorname{tr}(\mathbf{C}^{-1}\partial\mathbf{C})\\\
+\operatorname{tr}\\!\left[\\!\mathbf{C}^{-1}\partial\mathbf{C}\mathbf{C}^{-1}\\!\left(\mathbf{\Psi}_{n}\odot\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})^{\top}\right)\\!\right]$
(42)

A sufficient condition for the differential to be equal to zero is:

$\displaystyle\frac{T}{2}\mathbf{I}+\mathbf{C}^{-1}\\!\left(\mathbf{\Psi}_{n}\odot\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})^{\top}\right)=\mathbf{0}$
(43)

which holds if:

$\displaystyle\mathbf{C}=\mathbf{\Psi}_{n}\odot\frac{1}{T}\sum_{t=1}^{T}\frac{1}{\sigma_{n,t}^{2}}(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})(\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n})^{\top}\,.$
(44)

This is equivalent to the definition of the shrinkage covariance matrix
estimator of Proposition 2 since the weighting matrix $\mathbf{\Psi}_{n}$ is
defined as the $D\times D$ matrix with 1 on the diagonal and the value
$(1-\tilde{\rho}_{n})$ off the diagonal, and $\tilde{\rho}_{n}$ is computed
according to Eq.(12).

##  Appendix C Proof of Proposition 4

As stated in Proposition 4, we consider the shrinkage matrices
$\mathbf{\Psi}^{\text{spec}}$ and $\mathbf{\Psi}^{\text{spat}}$ to be fixed
(i.e., $\tilde{\rho}_{n}^{\text{spec}}$ and $\tilde{\rho}_{n}^{\text{spat}}$
are estimated once for all).

Before computing the differential of $\mathscr{D}_{n}$, note that its last
term,
$\frac{1}{2\sigma_{n,t}^{2}}\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\mathrm{spec}}\right\|_{\mathbf{\Gamma}}^{2}$,
can be rewritten in several different ways using the following identities:

$\|\boldsymbol{a}\|_{(\mathbf{\Psi}_{1}\otimes\mathbf{\Psi}_{2})\odot(\mathbf{\Gamma}_{1}\otimes\mathbf{\Gamma}_{2})}^{2}=\|\boldsymbol{a}\|_{(\mathbf{\Psi}_{1}\odot\mathbf{\Gamma}_{1})\otimes(\mathbf{\Psi}_{2}\odot\mathbf{\Gamma}_{2})}^{2}$
$\displaystyle=\operatorname{tr}[((\mathbf{\Psi}_{1}\odot\mathbf{\Gamma}_{1})\otimes(\mathbf{\Psi}_{2}\odot\mathbf{\Gamma}_{2}))\boldsymbol{a}\boldsymbol{a}^{\top}]$
(45)
$\displaystyle=\operatorname{tr}[((\mathbf{\Psi}_{1}\otimes\mathbf{\Psi}_{2})\odot(\mathbf{\Gamma}_{1}\otimes\mathbf{\Gamma}_{2}))\boldsymbol{a}\boldsymbol{a}^{\top}]$
(46)
$\displaystyle=\operatorname{tr}[(\mathbf{\Gamma}_{1}\otimes\mathbf{\Gamma}_{2})((\mathbf{\Psi}_{1}\otimes\mathbf{\Psi}_{2})\odot(\boldsymbol{a}\boldsymbol{a}^{\top}))]$
(47)
$\displaystyle=\operatorname{tr}[\mathbf{A}^{\top}(\mathbf{\Psi}_{2}\odot\mathbf{\Gamma}_{2})\mathbf{A}(\mathbf{\Psi}_{1}\odot\mathbf{\Gamma}_{1})]$
(48)
$\displaystyle=\operatorname{tr}\bigl{[}\mathbf{\Gamma}_{2}\bigl{(}\mathbf{\Psi}_{2}\odot(\mathbf{A}(\mathbf{\Psi}_{1}\odot\mathbf{\Gamma}_{1})\mathbf{A}^{\top})\bigr{)}\bigr{]}$
(49)
$\displaystyle=\operatorname{tr}\bigl{[}\mathbf{\Gamma}_{1}\bigl{(}\mathbf{\Psi}_{1}\odot(\mathbf{A}^{\top}(\mathbf{\Psi}_{2}\odot\mathbf{\Gamma}_{2})\mathbf{A})\bigr{)}\bigr{]}\,,$
(50)

where $\mathbf{A}$ is a $m\times n$ matrix obtained by reshaping the column
vector $\boldsymbol{a}$ according to the dimensions of matrices
$\mathbf{\Gamma}_{1}\in\mathbb{R}^{n\times n}$ and
$\mathbf{\Gamma}_{2}\in\mathbb{R}^{m\times m}$, such that
$\text{vec}(\mathbf{A})=\boldsymbol{a}$.

The total differential of the cost function $\mathscr{D}_{n}$, for fixed
matrices $\mathbf{\Psi}^{\text{spec}}$ and $\mathbf{\Psi}^{\text{spat}}$, is:

$\partial\mathscr{D}_{n}=\tfrac{TK}{2}\operatorname{tr}(\mathbf{C}_{n}^{\text{spec
}-1}\partial\mathbf{C}_{n}^{\text{spec}})+\tfrac{TL}{2}\operatorname{tr}(\mathbf{C}_{n}^{\text{spat
}-1}\partial\mathbf{C}_{n}^{\text{spat}})\\\
+\sum_{t=1}^{T}\biggl{(}\tfrac{KL}{2}\frac{\partial\sigma_{n,t}^{2}}{\sigma_{n,t}^{2}}-\frac{\partial\sigma_{n,t}^{2}}{2\sigma_{n,t}^{4}}\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\mathrm{spec}}\right\|_{\mathbf{\Gamma}}^{2}\\\
+\frac{1}{\sigma_{n,t}^{2}}\\!\left(\boldsymbol{\mu}_{n}^{\text{spec}}-\boldsymbol{v}_{n,t}\right)^{\top}\mathbf{\Gamma}\partial\boldsymbol{\mu}_{n}^{\text{spec}}\\\
-\frac{1}{2\sigma_{n,t}^{2}}\operatorname{tr}\biggl{[}\mathbf{C}_{n}^{\mathrm{spec}-1}\partial\mathbf{C}_{n}^{\mathrm{spec}}\mathbf{C}_{n}^{\mathrm{spec}-1}\bigl{\\{}\hfill\\\
\hfill\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\bigl{(}\overline{\mathbf{V}}_{t}^{\top}(\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\mathbf{C}_{n}^{\mathrm{spat}-1})\overline{\mathbf{V}}_{t}\bigr{)}\bigr{\\}}\biggr{]}\\\
-\frac{1}{2\sigma_{n,t}^{2}}\operatorname{tr}\biggl{[}\mathbf{C}_{n}^{\mathrm{spat}-1}\partial\mathbf{C}_{n}^{\mathrm{spat}}\mathbf{C}_{n}^{\mathrm{spat}-1}\bigl{\\{}\hfill\\\
\hfill\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\bigl{(}\overline{\mathbf{V}}_{t}(\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\mathbf{C}_{n}^{\mathrm{spec}-1})\overline{\mathbf{V}}_{t}^{\top}\bigr{)}\bigr{\\}}\biggr{]}\biggr{)}\,,\\\
$ (51)

where $\overline{\mathbf{V}}_{t}$ represents the spatio-spectral patch in
matrix form (see Proposition 1).

The following set of conditions is sufficient for $\partial\mathscr{D}_{n}$ to
be equal to zero:

$\displaystyle\sum_{t=1}^{T}\frac{1}{\sigma_{n,t}^{2}}\\!\left(\boldsymbol{\mu}_{n}^{\mathrm{spec}}-\boldsymbol{v}_{n,t}\right)=0$
(52a) $\displaystyle
KL\sigma_{n,t}^{2}-\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\mathrm{spec}}\right\|_{\mathbf{\Gamma}}^{2}=0$
(52b)
$\displaystyle\tfrac{TK}{2}\mathbf{I}-\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}\mathbf{C}_{n}^{\mathrm{spec}-1}\bigl{\\{}\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot$
$\displaystyle\hskip
71.13188pt\bigl{(}\overline{\mathbf{V}}_{t}^{\top}(\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\mathbf{C}_{n}^{\mathrm{spat}-1})\overline{\mathbf{V}}_{t}\bigr{)}\bigr{\\}}=\mathbf{0}$
(52c)
$\displaystyle\tfrac{TL}{2}\mathbf{I}-\sum_{t=1}^{T}\frac{1}{2\sigma_{n,t}^{2}}\mathbf{C}_{n}^{\mathrm{spat}-1}\bigl{\\{}\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot$
$\displaystyle\hskip
71.13188pt\bigl{(}\overline{\mathbf{V}}_{t}(\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\mathbf{C}_{n}^{\mathrm{spec}-1})\overline{\mathbf{V}}_{t}^{\top}\bigr{)}\bigr{\\}}=\mathbf{0}$
(52d)

Condition (52a) holds if:

$\displaystyle\widehat{\boldsymbol{\mu}}_{n}^{\mathrm{spec}}$
$\displaystyle=\left(\sum_{t=1}^{T}\frac{1}{\sigma_{n,t}^{2}}\boldsymbol{v}_{n,t}\right)\biggl{/}\left(\sum_{t=1}^{T}\frac{1}{\sigma_{n,t}^{2}}\right).$
(53) Condition (52b) holds if: $\displaystyle\widehat{\sigma}_{n,t}^{2}$
$\displaystyle=\tfrac{1}{KL}\left\|\boldsymbol{v}_{n,t}-\boldsymbol{\mu}_{n}^{\mathrm{spec}}\right\|_{\mathbf{\Gamma}}^{2}\,.$
(54) Condition (52c) holds if:
$\displaystyle\widehat{\mathbf{C}}_{n}^{\mathrm{spec}}$
$\displaystyle=\tfrac{1}{TK}\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\bigl{(}\overline{\mathbf{V}}_{t}^{\top}(\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\mathbf{C}_{n}^{\mathrm{spat}-1})\overline{\mathbf{V}}_{t}\bigr{)}\,.$
(55) Condition (52d) holds if:
$\displaystyle\widehat{\mathbf{C}}_{n}^{\mathrm{spat}}$
$\displaystyle=\tfrac{1}{TL}\sum_{t=1}^{T}\frac{1}{\widehat{\sigma}_{n,t}^{2}}\mathbf{\Psi}_{n}^{\mathrm{spat}}\odot\bigl{(}\overline{\mathbf{V}}_{t}(\mathbf{\Psi}_{n}^{\mathrm{spec}}\odot\mathbf{C}_{n}^{\mathrm{spec}-1})\overline{\mathbf{V}}_{t}^{\top}\bigr{)}\,.$
(56)

These correspond to the expressions given in Proposition 4.

## Acknowledgment

This work has been supported by the Région Auvergne-Rhône-Alpes under the
project DIAGHOLO. This work has made use of the SPHERE Data Center, jointly
operated by OSUG/IPAG (Grenoble, France), PYTHEAS/LAM/CESAM (Marseille,
France), OCA/Lagrange (Nice, France), Observatoire de Paris/LESIA (Paris,
France), and Observatoire de Lyon/CRAL (Lyon, France).

## References

  * [1] M. Perryman, _The exoplanet handbook_. Cambridge University Press, 2018.
  * [2] C. Marois, D. Lafrenière, R. Doyon, B. Macintosh, and D. Nadeau, “Angular differential imaging: A powerful high-contrast imaging technique,” _The Astrophysical Journal_ , vol. 641, no. 1, p. 556, 2006.
  * [3] R. Racine, G. A. Walker, D. Nadeau, R. Doyon, and C. Marois, “Speckle noise and the detection of faint companions,” _Publications of the Astronomical Society of the Pacific_ , vol. 111, no. 759, p. 587, 1999.
  * [4] D. Lafrenière, C. Marois, R. Doyon, D. Nadeau, and E. Artigau, “A new algorithm for point-spread function subtraction in high-contrast imaging: a demonstration with angular differential imaging,” _The Astrophysical Journal_ , vol. 660, no. 1, p. 770, 2007.
  * [5] C. Marois, C. Correia, J.-P. Véran, and T. Currie, “TLOCI: A fully loaded speckle killing machine,” _International Astronomical Union_ , vol. 8, no. S299, pp. 48–49, 2013.
  * [6] C. Marois, C. Correia, R. Galicher, P. Ingraham, B. Macintosh, T. Currie, and R. De Rosa, “GPI PSF subtraction with TLOCI: the next evolution in exoplanet/disk high-contrast imaging,” in _SPIE Astronomical Intrumentation + Telescopes_ , vol. 9148. International Society for Optics and Photonics, 2014, p. 91480U.
  * [7] R. Soummer, L. Pueyo, and J. Larkin, “Detection and characterization of exoplanets and disks using projections on Karhunen-Loève eigenimages,” _The Astrophysical Journal Letters_ , vol. 755, no. 2, p. L28, 2012.
  * [8] A. Amara and S. P. Quanz, “PYNPOINT: an image processing package for finding exoplanets,” _Monthly Notices of the Royal Astronomical Society_ , vol. 427, no. 2, pp. 948–955, 2012.
  * [9] 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_, vol. 545, p. A111, 2012.
  * [10] B. Pairet, L. Jacques, and F. Cantalloube, “Iterative low-rank and rotating sparsity promotion for circumstellar disks imaging,” in _Signal Processing with Adaptive Sparse Structured Representations_ , vol. 1, 2019, p. 1.
  * [11] T. M. Esposito, M. P. Fitzgerald, J. R. Graham, and P. Kalas, “Modeling self-subtraction in angular differential imaging: Application to the HD 32297 debris disk,” _The Astrophysical Journal_ , vol. 780, no. 1, p. 25, 2013.
  * [12] Z. Wahhaj, J. Milli, C. Romero, L. Cieza, A. Zurlo, A. Vigan, E. Peña, G. Valdes, F. Cantalloube, J. Girard _et al._ , “A search for a fifth planet around HR 8799 using the star-hopping RDI technique at VLT/SPHERE,” _Astronomy & Astrophysics_, vol. 648, p. A26, 2021.
  * [13] B. Pairet, F. Cantalloube, and L. Jacques, “MAYONNAISE: a morphological components analysis pipeline for circumstellar discs and exoplanets imaging in the near-infrared,” _Monthly Notices of the Royal Astronomical Society_ , vol. 503, no. 3, pp. 3724–3742, 2021.
  * [14] O. Flasseur, S. Thé, L. Denis, É. Thiébaut, and M. Langlois, “REXPACO: An algorithm for high contrast reconstruction of the circumstellar environment by angular differential imaging,” _Astronomy & Astrophysics_, vol. 651, p. A62, 2021.
  * [15] O. Flasseur, L. Denis, É. Thiébaut, and M. Langlois, “Exoplanet detection in angular differential imaging by statistical learning of the nonstationary patch covariances - The PACO algorithm,” _Astronomy & Astrophysics_, vol. 618, p. A138, 2018.
  * [16] J.-L. Beuzit, A. Vigan, D. Mouillet, K. Dohlen, R. Gratton, A. Boccaletti, J.-F. Sauvage, H. M. Schmid, M. Langlois, C. Petit _et al._ , “SPHERE: the exoplanet imager for the Very Large Telescope,” _Astronomy & Astrophysics_, vol. 631, p. A155, 2019.
  * [17] A. Buades, B. Coll, and J.-M. Morel, “A non-local algorithm for image denoising,” in _IEEE Computer Society Conference on Computer Vision and Pattern Recognition_ , vol. 2. IEEE, 2005, pp. 60–65.
  * [18] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, “Image denoising by sparse 3-D transform-domain collaborative filtering,” _IEEE Transactions on Image Processing_ , vol. 16, no. 8, pp. 2080–2095, 2007.
  * [19] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” _IEEE Transactions on Signal Processing_ , vol. 54, no. 11, pp. 4311–4322, 2006.
  * [20] J. Mairal, F. Bach, J. Ponce, G. Sapiro, and A. Zisserman, “Non-local sparse models for image restoration,” in _IEEE International Conference on Computer Vision_. IEEE, 2009, pp. 2272–2279.
  * [21] D. Zoran and Y. Weiss, “From learning models of natural image patches to whole image restoration,” in _IEEE International Conference on Computer Vision_. IEEE, 2011, pp. 479–486.
  * [22] G. Yu, G. Sapiro, and S. Mallat, “Solving inverse problems with piecewise linear estimators: From Gaussian mixture models to structured sparsity,” _IEEE Transactions on Image Processing_ , vol. 21, no. 5, pp. 2481–2499, 2011\.
  * [23] M. Lebrun, A. Buades, and J.-M. Morel, “A nonlocal Bayesian image denoising algorithm,” _SIAM Journal on Imaging Sciences_ , vol. 6, no. 3, pp. 1665–1688, 2013.
  * [24] O. Flasseur, L. Denis, É. Thiébaut, and M. Langlois, “An unsupervised patch-based approach for exoplanet detection by direct imaging,” in _IEEE International Conference on Image Processing_. IEEE, 2018, pp. 2735–2739.
  * [25] ——, “PACO ASDI: an algorithm for exoplanet detection and characterization in direct imaging with integral field spectrographs,” _Astronomy & Astrophysics_, vol. 637, p. A9, 2020.
  * [26] E. Conte, M. Lops, and G. Ricci, “Asymptotically optimum radar detection in compound-Gaussian clutter,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 31, no. 2, pp. 617–625, 1995.
  * [27] M. J. Wainwright and E. P. Simoncelli, “Scale mixtures of Gaussians and the statistics of natural images.” in _Neural Information Processing Systems_ , vol. 12, 1999, pp. 855–861.
  * [28] O. Flasseur, L. Denis, É. Thiébaut, and M. Langlois, “Robustness to bad frames in angular differential imaging: a local weighting approach,” _Astronomy & Astrophysics_, vol. 634, p. A2, 2020.
  * [29] N. Lu and D. L. Zimmerman, “The likelihood ratio test for a separable covariance matrix,” _Statistics & Probability Letters_, vol. 73, no. 4, pp. 449–457, 2005.
  * [30] M. G. Genton, “Separable approximations of space-time covariance matrices,” _Environmetrics: The official Journal of the International Environmetrics Society_ , vol. 18, no. 7, pp. 681–695, 2007.
  * [31] K. Werner, M. Jansson, and P. Stoica, “On estimation of covariance matrices with Kronecker product structure,” _IEEE Transactions on Signal Processing_ , vol. 56, no. 2, pp. 478–491, 2008.
  * [32] O. Ledoit and M. Wolf, “A well-conditioned estimator for large-dimensional covariance matrices,” _Journal of Multivariate Analysis_ , vol. 88, no. 2, pp. 365–411, 2004.
  * [33] Y. Chen, A. Wiesel, Y. C. Eldar, and A. O. Hero, “Shrinkage algorithms for MMSE covariance estimation,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 10, pp. 5016–5029, 2010.
  * [34] X. Bresson and T. F. Chan, “Fast dual minimization of the vectorial total variation norm and applications to color image processing,” _Inverse Problems & Imaging_, vol. 2, no. 4, p. 455, 2008.
  * [35] É. Thiébaut, “Optimization issues in blind deconvolution algorithms,” in _Astronomical Data Analysis II_ , vol. 4847. International Society for Optics and Photonics, 2002, pp. 174–183.
  * [36] C. Zhu, R. H. Byrd, P. Lu, and J. Nocedal, “Algorithm 778: L-BFGS-B: Fortran subroutines for large-scale bound-constrained optimization,” _ACM Transactions on Mathematical Software_ , vol. 23, no. 4, pp. 550–560, 1997.
  * [37] P. Delorme, N. Meunier, D. Albert, E. Lagadec, H. L. Coroller, R. Galicher, D. Mouillet, A. Boccaletti, D. Mesa, J.-C. Meunier _et al._ , “The sphere data center: a reference for high contrast imaging processing,” _arXiv preprint arXiv:1712.06948_ , 2017.
  * [38] R. Galicher, A. Boccaletti, D. Mesa, P. Delorme, R. Gratton, M. Langlois, A.-M. Lagrange, A.-L. Maire, H. Le Coroller, G. Chauvin _et al._ , “Astrometric and photometric accuracies in high contrast imaging: The SPHERE speckle calibration tool (SpeCal),” _Astronomy & Astrophysics_, vol. 615, p. A92, 2018.
  * [39] J. Milli, A. Vigan, D. Mouillet, A.-M. Lagrange, J.-C. Augereau, C. Pinte, D. Mawet, H. M. Schmid, A. Boccaletti, L. Matra _et al._ , “Near-infrared scattered light properties of the HR 4796 a dust ring - A measured scattering phase function from 13.6° to 166.6°,” _Astronomy & Astrophysics_, vol. 599, p. A108, 2017.
  * [40] A.-L. Maire, T. Stolker, S. Messina, A. Müller, B. A. Biller, T. Currie, C. Dominik, C. A. Grady, A. Boccaletti, M. Bonnefoy _et al._ , “Testing giant planet formation in the transitional disk of SAO 206462 using deep VLT/SPHERE imaging,” _Astronomy & Astrophysics_, vol. 601, p. A134, 2017.

See pages 1 of REXPACO_ASDI_SUPPLEMENTARY.pdf See pages 2 of
REXPACO_ASDI_SUPPLEMENTARY.pdf