# Using Data Imputation for Signal Separation in High contrast Imaging

Bin Ren, Laurent Pueyo, Christine Chen, Elodie Choquet, John H. Debes, Gaspard Duchêne, Francois Menard, Marshall D. Perrin

► To cite this version:

Bin Ren, Laurent Pueyo, Christine Chen, Elodie Choquet, John H. Debes, et al.. Using Data Imputation for Signal Separation in High contrast Imaging. The Astrophysical Journal, 2020, 892 (2), pp.74. 10.3847/1538-4357/ab7024 . hal-03584367

HAL Id: hal-03584367

https://amu.hal.science/hal-03584367

Submitted on 6 Dec 2022

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
# Using Data Imputation for Signal Separation in High Contrast Imaging

Bin Ren (任彬) <sup>1,2,3</sup> Laurent Pueyo, <sup>4</sup> Christine Chen <sup>4,2</sup> Élodie Choquet <sup>5</sup> John H. Debes <sup>4</sup> Gaspard Duchêne <sup>6,7</sup>
François Ménard <sup>7</sup> and Marshall D. Perrin <sup>7</sup>

<sup>1</sup>Department of Astronomy, California Institute of Technology, 1216 East California Boulevard, Pasadena, CA 91125, USA

<sup>2</sup>Department of Physics and Astronomy, The Johns Hopkins University, 3701 San Martin Drive, Baltimore, MD 21218, USA

<sup>3</sup>Department of Applied Mathematics and Statistics, The Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA

<sup>4</sup>Space Telescope Science Institute (STScI), 3700 San Martin Drive, Baltimore, MD 21218, USA

<sup>5</sup>Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

<sup>6</sup>Astronomy Department, University of California, Berkeley, CA 94720, USA

<sup>7</sup>Université Grenoble-Alpes, CNRS Institut de Planétologie et d'Astrophysique (IPAG), F-38000 Grenoble, France

(Received 2019 August 19; Revised 2020 January 20; Accepted 2020 January 25; Published 2020 March 31)

## Abstract

To characterize circumstellar systems in high contrast imaging, the fundamental step is to construct a best point spread function (PSF) template for the non-circumstellar signals (i.e., star light and speckles) and separate it from the observation. With existing PSF construction methods, the circumstellar signals (e.g., planets, circumstellar disks) are unavoidably altered by over-fitting and/or self-subtraction, making forward modeling a necessity to recover these signals. We present a forward modeling-free solution to these problems with data imputation using sequential non-negative matrix factorization (DI-sNMF). DI-sNMF first converts this signal separation problem to a “missing data” problem in statistics by flagging the regions which host circumstellar signals as missing data, then attributes PSF signals to these regions. We mathematically prove it to have negligible alteration to circumstellar signals when the imputation region is relatively small, which thus enables precise measurement for these circumstellar objects. We apply it to simulated point source and circumstellar disk observations to demonstrate its proper recovery of them. We apply it to Gemini Planet Imager (GPI) K1-band observations of the debris disk surrounding HR 4796A, finding a tentative trend that the dust is more forward scattering as the wavelength increases. We expect DI-sNMF to be applicable to other general scenarios where the separation of signals is needed.

*Unified Astronomy Thesaurus concepts:* Astronomy data analysis (1858); Coronagraphic imaging (313); Debris disks (363); Dimensionality reduction (1943); Multivariate analysis (1913); Protoplanetary disks (1300)

## 1. Introduction

High contrast imaging of circumstellar systems in visible and near infrared light offers direct spectroscopic and astrometric information on the planet or the distribution of the dust in the system. With the advancements in both instrumental design (e.g., Gemini/GPI: Macintosh et al. 2008, 2014; VLT/SPHERE: Beuzit et al. 2008, 2019; Vigan et al. 2016) and data post-processing (e.g., Lafrenière et al. 2007; Soummer et al. 2012; Amara & Quanz 2012; Milli et al. 2012; Pueyo et al. 2012; Pueyo 2016. See Pueyo 2018 for a review), direct imaging techniques have been rapidly developing in the past five years. As a result, many objects including both

exoplanets (e.g., 51 Eri: Macintosh et al. 2015; Rajan et al. 2017; PDS 70: Keppler et al. 2018; Müller et al. 2018; Christiaens et al. 2019) and circumstellar disks (e.g., HD 106906: Kalas et al. 2015; Lagrange et al. 2016; HD 61005: Esposito et al. 2016; HD 114082: Wahhaj et al. 2016; HIP 73145 [HD 131835]: Feldt et al. 2017; HD 129590: Matthews et al. 2017) have been discovered and/or characterized.

In the post-processing of direct imaging data, the fundamental step is to properly model and remove the star light and speckles (i.e., point spread function, PSF) from the observations. Existing PSF modeling methods have been optimized for unresolved exoplanets where the point spread function is known (Marois et al. 2010; Pueyo 2016), or for circumstellar disks in space-based observations where the PSF is relatively stable (Ren et al. 2018). For ground-based observations, po-
larization differential imaging (PDI) has revealed complex circumstellar structures (e.g., Avenhaus et al. 2018; Monnier et al. 2019; Garufi et al. 2020). However, PDI requires linear polarization of the signals, which is typically not present for planets. For disks, such signals can be too faint and thus beyond detection. For the disks detected through PDI, total intensity signals will help better understand the dust properties.

To characterize circumstellar disks in total intensity, post-processing strategies such as reference differential imaging (RDI: Smith & Terrile 1984), angular differential imaging (ADI: Marois et al. 2006) and spectral differential imaging (SDI: Marois et al. 2000; Sparks & Ford 2002) are adopted. These strategies are accompanied with contaminations from quasi-static and rapidly varying speckles (Pueyo 2018), and the contaminations are expected be reduced with advanced post-processing methods (e.g., Locally Optimized Combination of Images [LOCI]: Lafrenière et al. 2007; Karhunen-Loève Image Projection [KLIP]: Soummer et al. 2012; Amara & Quanz 2012). Nevertheless, existing post-processing methods are limited by over-fitting and self-subtraction issues, and forward modeling is required to recover the surface brightness distribution of the circumstellar objects. However, forward modeling is currently only optimized for unresolved point sources (Marois et al. 2010; Pueyo 2016). For resolved circumstellar disks, forward modeling needs assumptions on both morphology and flux (e.g., space-based: Choquet et al. 2016, 2017, 2018; ground-based: Follette et al. 2017; Esposito et al. 2018), which is model-dependent thus likely unable to capture the minute irregularly-shaped structures (e.g., the disks in Avenhaus et al. 2018; Monnier et al. 2019; Garufi et al. 2020).

To properly separate circumstellar disk signals from observations, Ren et al. (2018) has demonstrated the applicability of sequential non-negative matrix factorization (sNMF) to space-based observations. For these observations with stable PSFs, sNMF is able to resolve the over-fitting problem for RDI observations after adopting a scaling factor. However, sNMF is still limited by not only unstable PSFs where multiple scaling factors are needed, but also ADI or SDI post-processing where self-subtraction still persists.

In this paper, we aim at extracting thus characterizing circumstellar signals in high contrast imaging observations using a data imputation approach. In multivariate statistical analysis, data imputation is a process which fills missing data based on their relationship with known data (e.g., Johnson & Wichern 2007). For high contrast imaging data, specifically, we first flag the regions that host circumstellar signals as “missing data” based on prior knowledge of their location. We then replace the “missing data” regions with PSF-only signals (i.e., data imputation). The circumstellar signal then resides in the residual image when the data imputation model

is removed from the original image. In this paper, we study the data imputation property of sNMF, i.e., DI-sNMF.

The structure of this paper is as follows: in Section 2 we describe the DI-sNMF method. In Section 3, we apply it to a simulated planet dataset. In Section 4, we apply it to circumstellar disk simulation and observation. We discuss and summarize the paper in Section 5. We present our mathematical derivations and proofs in Appendix A and Appendix B, and the auxiliary figures in Appendix C.

## 2. Methods

The NMF method was first introduced by Paatero & Tapper (1994), then Lee & Seung (2001) provided multiplicative update rules to guarantee its convergence. Its application to astronomical data started with Blanton & Roweis (2007), where weighting terms were added to take into account of observational uncertainties. Zhu (2016) then proposed a parallel form in studying spectroscopic data in astronomy to increase computational efficiency. Ren et al. (2018) demonstrated the applicability of NMF to high contrast imaging with a sequential construction of the components (i.e., sNMF).

### 2.1. sNMF for High Contrast Imaging

In the post-processing of high contrast imaging data, we first have $N_{\text{ref}}$ reference images with each image having $N_{\text{pix}}$ pixels. For each image, we flatten it to a 1-dimensional row matrix of length $N_{\text{pix}}$ by reordering the elements from the original 2-dimensional matrix. Following the definition of symbols in Ren et al. (2018), we denote the real-valued $N_{\text{ref}}$-by-$N_{\text{pix}}$ reference matrix by $R \in \mathbb{R}^{N_{\text{ref}} \times N_{\text{pix}}}$, which contains $N_{\text{ref}}$ flattened reference images. The NMF method performs dimension reduction via approximating $R$ with the product of two matrices (i.e., $R \approx WH$): a coefficient matrix $W \in \mathbb{R}^{N_{\text{ref}} \times n}$, and a component matrix $H \in \mathbb{R}^{n \times N_{\text{pix}}}$, where $n$ is the number of components. In the component matrix, each row is an NMF component that is expected to be a representation of at least a fraction of the signals in $R$.

*Component Construction*—The decomposition of $R$ is achieved with the following vectorized weighted multiplicative update rules (Zhu 2016):

$$W^{(k+1)} = W^{(k)} \circ \frac{(V \circ R)H^{(k)T}}{[V \circ W^{(k)}H^{(k)}]H^{(k)T}}, \quad (1)$$

$$H^{(k+1)} = H^{(k)} \circ \frac{W^{(k)T}(V \circ R)}{W^{(k)T}[V \circ W^{(k)}H^{(k)}]}, \quad (2)$$

where superscript (k) and (k+1) are iteration steps k and (k+1), respectively. Superscript T denotes matrix transpose. $V \in \mathbb{R}^{N_{\text{ref}} \times N_{\text{pix}}}$ is a weighting matrix (which is the element-wise inverse of the squared uncertainties of R in Zhu 2016).
Symbol ◦ and fraction bar (· · ·) are element-wise multiplication and division of two matrices of the identical dimension, respectively. In the iteration, the initialization values are uniformly drawn from 0 to 1, and the convergence of the above iteration rules is guaranteed when the iteration number increases (Zhu 2016). Consequently, this initialization guarantees non-negativity for W and H. For an intuitive understanding of the component construction update rules, if we ignore the weighting matrix V and the coefficient matrix W (k) on the righthand side of Equation (1), we have

$$W = \frac{RH^T}{HH^T},$$

which is a vector projection of R to H. Similarly, the intuitive understanding of Equation (2) is the symbolic conjugate of Equation (1). In other words, the update rules performs least square estimation but with all the matrix elements non-negative (Ren et al. 2018). Nevertheless, we note that such an approach is not mathematically established, and it is for the purpose of intuitive understanding of NMF update rules only.

To rank the contribution from each component, Ren et al. (2018) studied the sequential construction of NMF components (i.e., sNMF). Specifically, given m constructed components ($m = 1, 2, \dots, n-1$), the construction of the ($m+1$)-th component is initialized with the converged values of the first m components in the first m rows of H, with the ($m+1$)-th row of H still uniformly drawn from 0 to 1.

*Target Modeling*—For a flattened target image $T \in \mathbb{R}^{1 \times N_{pix}}$, it can be modeled by the sNMF components H constructed above. The coefficients are obtained with the following update rule:

$$\omega^{(k+1)} = \omega^{(k)} \circ \frac{(v \circ T)H^T}{[v \circ \omega^{(k)}H]H^T}, \quad (3)$$

where $v \in \mathbb{R}^{1 \times N_{pix}}$ is the weighting matrix for T, and $\omega \in \mathbb{R}^{1 \times n}$ is the coefficient matrix for T. This modeling procedure is expected to be linear to the first order (Ren et al. 2018), i.e., the signal of the circumstellar disk or planet can be linearly separated from the PSF. For an intuitive understanding of the target modeling update rule, we can ignore the v and $\omega^{(k)}$ terms in Equation (3) and obtain a least square estimation in the form of

$$\omega = \frac{TH^T}{HH^T},$$

which is obtaining the coefficients the projection of the target T onto the component basis H, and thus it is effectively a special case of Equation (1). Nevertheless, similarly as for the intuitive understanding of the component construction update rules, we readdress such an approach is not mathematically founded and is for illustration purpose only.

sNMF with “Missing Data”—When there is “missing data” in R or T, Equations (1)–(3) can handle and ignore them during component construction or target modeling. This is achieved by first flagging the missing data with a binary mask, and then ignoring them during component construction and/or target modeling. Specifically, the mask is set to 0 when there is circumstellar signal and 1 otherwise. This binary mask concept is identical to the mask in previous methods in high contrast imaging (e.g., Pueyo et al. 2012; Milli et al. 2012). For sNMF, we neglect the contribution from circumstellar signals by treating their corresponding regions as mock “missing data” in the sequential construction of the sNMF component basis. We note that this approach ignores the matrix elements, rather than affecting the update rules by replacing these elements with 0. We prove that this element ignoring approach does not bias the sNMF components with non-PSF signals or zero values in Appendix A.1. We quantify in Appendix A.2 that the bias introduced by missing data is negligible when the masked out region is relatively small.

## 2.2. Data Imputation for “Missing Data”

The data imputation approach imputes the mock “missing data” back by modeling an observation with the sNMF components obtained in the previous subsection. Specifically, in target modeling, we obtain the coefficients for a target using Equation (3) while ignoring the “missing data” region for both the target and the component basis (i.e., the partial target and the partial component basis). We denote the corresponding coefficients by $\omega'$.

We apply $\omega'$ to the whole component basis H and obtain the data imputation model,

$$\tilde{T} = \omega' H \quad (4)$$

where $\tilde{T} \in \mathbb{R}^{1 \times N_{pix}}$ is the DI-sNMF model of T. The missing data region is now imputed with PSF-only signals.

We now establish the mathematical foundation for the data imputation procedure for DI-sNMF. On one hand, we prove in Appendix B.1 that when the non-PSF signals are neglected, they do not bias the coefficients obtained from Equation (3), i.e., $\omega'$. These coefficients are applied to the full sNMF component basis in Equation (4) to obtain an empirical model of the “missing data” region using PSF-only signals. On the other hand, when the “missing data” region is relatively small, these coefficients are expected to have negligible deviation from the coefficients in the scenario where there is only PSF signals in the entire field of view, see Appendix B.2 for the proof.

We develop the nmf_imaging package (Ren 2018)<fnref n="1" /> for the above steps, i.e., component construction and target modeling with “missing data” using binary masks.

<footnote n="1">¹ https://github.com/seawander/nmf_imaging</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Aperture photometry RDI recovery percentages of planet signals injected to HST/STIS BAR5 observations of HD 38393. The injected planets are (a) 10×, (b) 1×, and (c) 0.1× the brightness of the speckles at $r = 1''25$. In (b), DI-sNMF is able to recover the injected planet signal at $1\sigma$ levels within ~5% with decreasing bias when the number of components increases, which outperforms sNMF and KLIP. For the bright planet in (a), DI-sNMF robustly recovers it with small fractional uncertainty. For the faint planet in (c), DI-sNMF is still able to recover the signal albeit with larger fractional uncertainty. Note: (1) The photometry aperture is a 3×3 pixel square, and we conservatively estimate the uncertainties using the square root of the quadratic sum of the element-wise standard deviation across the 90 residual images inside the aperture. (2) To reduce the cluttering of the plots, we only present the uncertainties for DI-sNMF in shaded areas, and the typical uncertainties for other reduction results are presented at the top of the corresponding plots. (3) The normalized DI-sNMF images (dimension: 87×87 pixel or 4''41×4''41), obtained by dividing the reduction results by the total brightness of the corresponding injected planet, are shown in log scale from 0.1% to 10% of the total brightness of the planet models. The numbers of components used are presented at the bottom of the images. (4) At the location of the injected point source, the photon noise of the PSF wing is ~6% per pixel.</figcaption></figure>

## 2.3. Circumstellar Signal Extraction

When the DI-sNMF model in Equation (4) is removed from the original target image, the residual map,

$$ \tilde{D} = T - \tilde{T}, \qquad (5) $$

where $\tilde{D} \in \mathbb{R}^{1 \times N_{\text{pix}}}$, is expected to host the circumstellar signals (i.e., non-PSF signals).

# 3. Exoplanet Imaging

## 3.1. Data Reduction

We first demonstrate the DI-sNMF method using RDI post-processing of a point source. We use the *Hubble Space Telescope (HST)/STIS* coronagraphic imaging observations of HD 38393 behind the BAR5 occulter (Schneider et al. 2017; Debes et al. 2019) for this purpose. HD 38393 was visited with 9 orbits during HST Cycle 23 (Proposal ID: 14426,

PI: J. Debes), each with a distinct telescope orientation and a 3×3 sub-pixel dithering pattern. In each dithering position of a visit, a 0.2 s readout is repeated for 10 times. These data have been used in Ren et al. (2018), here we only focus on the observations at the central dithering location, i.e., a total of 90 images (= 9 telescope orientations × 1 dithering position × 10 readouts). The 90 images are of dimension 87×87 pixel (4''41×4''41), and we focus only on the pixels with angular separation between 1''0 and 1''5. To simulate STIS' instrumental response to an unresolved planet signal in visible light, we obtain the *TinyTim* point source model (Krist et al. 2011)<fnref n="2" /> for a G0 spectral energy distribution.

We simulate RDI post-processing as follows. We first inject a TinyTim point source at certain brightness with Pois-

<footnote n="2">² http://www.stsci.edu/software/tinytim/</footnote>
son noise to the observations of one telescope orientation (i.e., a total of 10 images). We position the point source at a physical location with radial separation *r* = 1.25 and position angle PA = 70° (North to East) to avoid occultation of the planet by the occulter or by the diffraction spikes, thus ensuring the maximum coverage of the planet and forming the target exposures. We then construct the KLIP and sNMF component bases using the rest of the 80 reference images that are planet-free. The 10 target images are then reduced with KLIP, sNMF, and DI-sNMF. For DI-sNMF, specifically, we use the basis to impute the “missing data” region—a 3×3 pixel region<fnref n="3" /> where the planet is expected to reside—with PSF-only signals. We repeat this injection subtraction procedure for 9 times for all of the telescope observations, and calculate the element-wise median and standard deviation for the 90 PSF-subtracted residual images.

We repeat the above procedure 3 times for a planet at different brightness levels—one that is as bright as the speckles (i.e., in the 3×3 pixel box, the total planet signal is equal to the total speckle signal), and two that are either 10 or 0.1 times of that. For the results obtained from all the reduction methods, we measure the brightness of the point source following Debes et al. (2019): we first integrate over the 3×3 pixel box for the element-wise median of the 90 PSF-subtracted images, then subtract that value by 9 times the median of the surrounding pixel values (that are within 2 and 4 pixels from the point source center) to remove the background. We calculate a conservative uncertainty<fnref n="4" /> by taking the square root of the quadratic sum of the element-wise standard deviation of the 90 PSF-subtracted images in the 3×3 pixel box. We present the recovery fraction results, as well as some recovery fraction images from DI-sNMF, in Figure 1.

## 3.2. Results

In Figure 1, we observe that DI-sNMF is able to robustly recover the injected planet at various brightness levels. When more DI-sNMF components are used, the recovery is better. Specifically, for a planet that is as bright as the speckles, we notice that DI-sNMF is able to properly recover its flux within 5% at $1\sigma$. We notice that the fractional noise for the point source that is 10× the brightness of the speckles is similar but slightly smaller than that of 1×. This indicates that the variation introduced by PSF subtraction overwhelms the Poisson noise from the injected planet. When the planet is 0.1× the brightness of the speckles, both PSF subtraction variation and Poisson noise contribute to the uncertainty, thus resulting to larger fractional uncertainties.

For KLIP, we confirm its expected over-subtraction behavior, i.e., the recovered planet flux decreases when more components are used (Soummer et al. 2012). In this way, KLIP requires forward modeling to recover the original flux for the planet (i.e., KLIP-FM: Pueyo 2016). Forward modeling requires prior knowledge of the 2-dimensional distribution of the unresolved planet signal (Pueyo 2016), however this information is not required in DI-sNMF. Given that DI-sNMF only requires a 2-dimensional mask to flag the “missing data” region, it is therefore another independent planet characterization method with less prior requirement. In addition, DI-sNMF is able to achieve smaller uncertainties than KLIP.

For sNMF, we observe a less aggressive over-subtraction than KLIP, and sNMF still requires a scaling factor to recover the flux as studied in Ren et al. (2018). However, the adoption of a scaling factor is expected to increase the overall brightness distribution in the field of view, thus increasing the overall background noise. For unstable PSF images such as ground-based observations, multiple scaling factors are needed to better recover the point source and its surrounding background, however such approach is currently not well-founded (Ren et al. 2018). Instead, we study the DI-sNMF method to avoid this problem in a mathematically well-established way.

In this section, we only performed DI-sNMF for the RDI target modeling process. For point source observations with ADI post-processing, we can mimic the RDI procedure by masking out the point sources during both component construction and target modeling (see similar approaches in Pueyo 2016 and references therein). For observations with SDI post-processing, we can mimic RDI by spatially scaling the observations at different wavelengths according to their wavelengths and masking out the planet sources. We do not demonstrate the ADI, SDI, or ADI + SDI capabilities of DI-sNMF for point sources since there is no mathematical difference between them and RDI. For the DI-sNMF applicability in the component construction process, we demonstrate it in a more complicated scenario with circumstellar disks in the following section.

# 4. Circumstellar Disk Imaging

## 4.1. Simulation (ADI)

### 4.1.1. Data Reduction

We now explore the capability of DI-sNMF in terms of frame diversity and mask overlap, which can be translated to a simulated ADI post-processing scenario. We first assemble from the original HST/STIS 810 images of HD 38393 in Section 3.1, but instead sample one image from each dithering position at each telescope orientation, totaling 81 images (= 9 telescope orientations × 9 dithering position × 1 readout) to obtain moderate variation of the signals. To simulate observations for ADI post-processing, we treat these 81 images

<footnote n="3">³ The size of this region is identical to region that is used to determine the point source contrast for STIS BAR5 in Debes et al. (2019).</footnote>

<footnote n="4">⁴ The uncertainties in this paper are $1\sigma$ unless otherwise specified.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Data imputation for simulated ADI post-processing of an isotropic ring (dimension of each image: 87×87 pixel or 4.''41×4.''41). Different rows correspond to the disks that are 10×, 1×, 0.1×, and 0.01× the brightness of the PSF wing, respectively. When more components are used the disks are better recovered. When the disk is 0.01× the brightness of the PSF, its recovery is marginal. See Figure 3 for the corresponding residual maps. Note: each row is at different brightness levels, and each color bar is only for the corresponding row.</figcaption></figure>

as independent PSF realizations that do not depend on tele-
scope orientations. We assign 1.°5 separations among adja-
cent images, thus obtaining a total field rotation of 120° Dur-
ing data reduction, we reduce the impact from STIS’ diffrac-
tion spikes and the BAR5 occulter using a mask created by
Debes et al. (2017). For a specific application to a ground-
based dataset, see Section 4.2.

*Simple Ring*—To simulate ADI data reduction on mock “missing data”, we first demonstrate the DI-sNMF method using an isotropic scattering ring generated by the *MCFOST* radiative transfer software (Pinte et al. 2006, 2009). We inject the ring to the 81 images, construct the components by masking out the region containing the ring, and impute the ring region with PSF signals using DI-sNMF. We calculate the element-wise median of the derotated and reduced images to obtain the final result, and calculate the element-wise standard deviation of these images to obtain the noise map. We generate the ring at four brightness levels—the disk are 10×, 1×, 0.1×, and 0.01× the brightness of the PSF wing

per pixel, respectively. We present the ring models and re-
duction results in Figure 2, and the residual maps in Figure 3.

*Complex Structures*—Being aware of different circumstellar structures, we generalize the ring study and mask out four different regions—bar, ring, clump, and spiral arms. Noticing the fact that the recovery quality of circumstellar signals depends only on the data imputation quality for the PSF signals (see Appendix A), we do not perform the injection-recovery process as for the planet or the ring model. Instead, we focus on the recovery quality of the PSF signals when the original values are known for a specific frame. For each circumstellar structure, we construct the component basis using the masked PSF signals, then compare the recovered image with the original (i.e., non-masked) image of the frame. We present the results from the ADI reduction with DI-sNMF for these structures in Figure 4.

4.1.2. Results

For the ring disk model, we notice in Figure 2 that DI-
sNMF with ADI can detect the disk from 10× down to 0.1×
the brightness of the PSF wing. When the disk is 0.01×
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Residual maps for Figure 2. The residual maps are obtained by comparing the median of the DI-sNMF results with the ring models, and divide that by the standard deviation of the DI-sNMF results for different scenarios. In all scenarios, the disk model can be better recovered when more components are used. When all the components are used, the disks are recovered within $0.5\sigma$.</figcaption></figure>

the PSF wing, the detection is marginal since it is below the photon noise of ~0.06× the PSF wing. For all brightness levels, more components are likely better recover the ring brightness. In the residual maps in Figure 3, we notice that when more components are used, the detection of the ring is reached within $0.5\sigma$.

Noticing the ring model might be a rare case of well-recoverable geometry, we present the PSF recovery fractions for a single frame under the four different scenarios in Figure 4. Even though the ideal approach is to perform inject-recovery exercises for these scenarios, we do not implement this approach since we expect the recovery of circumstellar signals to depend only on the recovery of the PSF signals (Appendix A). When we exclude the circumstellar structures during component construction, most of the field of view is covered by at least 80% of the images (Figure 4c), which thus allows for the capture of the PSF signals with DI-sNMF.

Using the DI-sNMF components, we notice that when we increase the number of components used in modeling, the recovered values for the pixels are closer to the observation. When 10 components are used (Figure 4f), most of the pixels are recovered with an error of less than ~5%, and the masked regions are recovered with an error of less than ~10%. In other words, a recovery precision of ~10% for the PSF signals thus allows for the recovery of circumstellar signals that are at least ~10% the brightness of the PSF wing or speckles.

For the point source simulation in Section 3, we observe that a few components are able to properly recover the flux. However for extended structures, we observe that using too few components can lead to insufficient recovery for the PSF signals, see Figure 4 where using too few components may cause systematic offsets in PSF recovery thus biasing circumstellar structure extraction.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Data imputation for simulated ADI post-processing in four scenarios: bar, ring, clump, and spiral. Columns (a) and (b): original and masked exposures. Column (c): coverage fraction of the field of view when the structures are excluded in the simulated ADI sequence. Most of the pixels are covered by >80% of the 81 images, thus allowing for the PSF signal capture during component construction. Columns (d), (e) and (f): fractional difference between a single frame (i.e., column a) and the DI-sNMF PSF models for the masked frame (i.e., column b) with *n* = 1, 5, and 10 components, respectively. The recovery quality increases as more components are used. When 10 components are used, the masked region in column (b) are recovered within ~10%. Note: all of the images here are of dimension 87×87 pixel (4.′′41×4.′′41).</figcaption></figure>

## 4.2. Application (ADI + SDI)

### 4.2.1. Data Reduction

We demonstrate the ADI + SDI capability of DI-sNMF with the ground-based observations of a debris disk—the HR 4796A system that hosts a narrow disk (Mouillet et al. 1997; Augereau et al. 1999; Schneider et al. 1999, 2009, 2018; Perrin et al. 2015; Debes et al. 2008; Milli et al. 2015, 2017, 2019). We obtained the Gemini/GPI K1-band coronagraphic observations of the system on 2015 April 03 (“K1-coron”; GS-2015A-Q-27, PI C. Chen). The GPI K1-coron data at 1.89 µm–2.20 µm (*R* ≈ 62–70, with 37 channels) are taken with 90s exposure, totaling 48 exposures with a total field rotation of 78°5 (a total of 48 × 37 = 1776 images).

To recover the disk with DI-sNMF, which is marginally visible at far separations but dominated by speckles in the close-in regions in Figure C1, we assume it with an inclination of 76° and a position angle for the major axis of 26° from North to East (Perrin et al. 2015; Milli et al. 2017). To create a mask where the disk is excluded for data imputation, we exclude the region that are between 0.′′99 and 1.′′23 (70 and 87 pixels. 14.166 ± 0.007 mas pixel⁻¹: De Rosa et al. 2015) from the star in the deprojected plane. For each observation, we rotate the original mask to match the corresponding time-averaged parallactic angle. To make the maximum use of available PSF signal, we spatially scale the 37 channel images based on their wavelength ratios with that of the 1st channel. For each exposure, we then collapse the 37 channels to 5 bins via pixel-wise median combination (6 channel
per bin, excluding the first 5 and last 2 channels) to increase the signal quality.

After the above procedure, we have a total of 240 images (5 bins each with 48 exposures). Identical procedures are performed on the masks, and the non-zero elements in the masks are assigned to be 1. See Figure C1 in Appendix C for the raw spatially scaled images before collapsing, and the coverage fraction map after the above disk-exclusion procedure. After the above procedure, most of the field is covered by PSF-only signals. We use the 240 scaled-and-collapsed images for sequential NMF component construction, with the 240 masks excluding the disks at different detector locations in each of the images. In Figure C2 in Appendix C, we present component examples that are constructed in this way—no disk signal is captured in the components as expected from the mathematical derivation in Appendix A.1.

Using the above components that cover most of the field of view (that includes the region where the disk resides), we model the disk-free regions with DI-sNMF. For each image, we flag the region where the disk resides as “missing data” using the corresponding binary mask created above, then obtain the DI-sNMF coefficients for the disk-excluded regions. The coefficients are then used to impute the disk region with the constructed sNMF components. We then obtain the residual image by subtracting the DI-sNMF model from the original target image. The residual images are then spatially scaled back to their original size on detector, derotated to north up and east left. We obtain the final result by taking the element-wise median of these images, and the noise map using the element-wise standard deviation of these images. The corresponding SNR map is then the ratio between the final result and the noise map.

### 4.2.2. Results

*Image*—Using the above approach, we are able to obtain the disk image at 5 different bins (i.e., 5 wavelengths). We present the reduction result<fnref n="5" /> using 10 components<fnref n="6" /> and the corresponding signal to noise ratio (SNR) map for the central wavelength (2.05 µm) in Figure 5.

We observe the ring structure with SNR > 5, and the disk image does not have self-subtraction effects as expected for DI-sNMF. For comparison, we reduce the data using KLIP ADI with 10 components, obtain and present the final image and SNR map in the identical scale as DI-sNMF. In the reduction results, there is noticeable self-subtraction effects in the KLIP result in the forms of negative halos around the disk and reduced surface brightness for the disk, however these effects

are not detectable in the DI-sNMF result. What is more, DI-sNMF is able to recover the minor axis of the ring.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Reduction results of the GPI K1-coron observations of HR 4796A from DI-sNMF and KLIP. (a) Disk image at the central wavelength (2.05 µm) obtained from DI-sNMF. (b) SNR map for the disk image, obtained by dividing the image by the element-wise standard deviation of the composing 48 individual images. (c) and (d) are the corresponding KLIP results, with the negative regions around the disk demonstrating the self-subtraction effect. In comparison, the DI-sNMF result not only has no noticeable negative regions around the disk, which demonstrates that it has no self-subtraction effect, but also has higher SNRs for the disk.</figcaption></figure>

We expect that DI-sNMF is able to properly recover the HR 4796A disk in two aspects. On one hand, there are no detectable self-subtraction effects in the DI-sNMF result in Figure 5a, indicating that self-subtraction is not impacting the result. On the other hand, the SNRs for the background regions (i.e., non-disk regions) are well-behaved with values of ~±1 in Figure 5b, i.e.,

$$ \frac{\tilde{D}_i}{\sigma_i} = \pm 1 + o(1), \qquad (6) $$

where *i* is the pixel index in the non-disk regions<fnref n="7" />, $\sigma_i$ is the corresponding standard deviation, and $o(\cdot)$ is the little o notation (i.e., $|o(x)| \ll |x|$). For independent and normally distributed noise (usually assumed for the PSF subtraction residuals), $D_i \sim N(\mu_i, \sigma_i)$, where $\mu_i$ is the corresponding expectation. In this way, $(\frac{D_i - \mu_i}{\sigma_i})^2$ follows a chi-squared

<footnote n="5">⁵ To reduce confusions from image spatial scaling, the images are all presented in detector counts per pixel in this paper unless otherwise specified.</footnote>

<footnote n="6">⁶ The final disk image does not have significant variation when more than 5 components are used.</footnote>

<footnote n="7">⁷ This definition of *i* is only applicable within Section 4.2.2. In other Sections, *i* refers to any pixel inside the field of view.</footnote>
distribution with 1 degree of freedom (i.e., $\chi_1^2$), whose corresponding expectation is 1. Assuming ideal PSF subtraction (i.e., $\mu_i = 0$), we have the expectation of $(\frac{D_i}{\sigma_i})^2$ to be 1. This value is consistent with the squared version of Equation (6) to the first order, and thus we expect that in our regime of calculating the noise map, the PSF-only signals in the background regions are consistent with ideal PSF reconstruction. Consequently, we expect that the PSF-only signals in the disk-hosting regions are also properly modeled and subtracted, since the constructed PSFs do not have significant deviations between non-masked and masked scenarios when the masked out regions are relatively small (Appendix B.2). Nevertheless, we estimate that the deviation has an upper limit of ~10%, since we masked out ~15% of the field of view for the HR 4696A observations, see the notes in Appendix B.2 for more specific discussion.

**Scattering Phase Function**—Based on the photometry recovery results in Section 3, we do not focus on the brightness of the disk to reduce redundancy. Instead, we extract the surface brightness distribution as a function of scattering angle (i.e., the scattering phase function, SPF: Hansen & Travis 1974. See Hughes et al. 2018 for a recent collection of SPFs for debris disks). Assuming the disk is circular and there is no offset between the disk center and the star, we obtain the original SPF curves using the equation in Ren et al. (2019). We then normalize the curves at 90° ± 10°, and divide them by the original SPF of an isotropic scattering disk created by MCFOST to correct for limb-brightening and convolution effects (Milli et al. 2017). We present the normalized SPFs in Figure 6.

The K1 SPFs are generally consistent with the SPF in the H2 band observed using the SPHERE instrument (1.67 µm: Milli et al. 2017). Nevertheless, for the probed scattering angles, we observe a tentative trend that the SPF is more forward scattering as the wavelength increases. Although the differences are less than 2σ, this trend is expected for the dust in the geometrical scattering regime (e.g., Figure 11 in Schröder et al. 2014).

We expect the above trend to be not biased by the symmetric scattering in the most forward and most backward scattering regions (see Milli et al. 2017 for the corresponding discussion), since it does not change when we proportionally remove the most backward scattering signals from the most forward scattering regions. In addition, we argue that the halo structure in HR 4796A (Schneider et al. 2018) is likely not biasing the measurements, since it is ~10% the brightness than the ring, and the 5 image are reduced in an identical way thus the halo creates a systematic offset for all of them. What is more, the difference of the most forward scattering part in the 5 SPFs is larger than 10%, making it unlikely for the deviation to be originated from the halo.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Normalized SPFs of the HR 4796A disk in K1 observation with DI-sNMF. The 1.67 µm data (black) is the averaged SPHERE H2 SPF in Milli et al. (2017). In the probed scattering angles, the inset image suggests that the dust are more forward scattering as wavelength increases.</figcaption></figure>

To further establish the above SPF trend as a function of wavelength, the comparison between the observations at multiple bands (e.g., H, J, K1, K2) is preferred (e.g., Chen et al. 2020). We do not perform a multi-band experiment since the flux calibration between the different bands need to be properly treated, and more sky rotation are needed to properly recover the SPFs at the other bands (Chen et al. 2020).

# 5. Discussion and Summary

In this paper, we present a new data reduction method, DI-sNMF, for high contrast imaging obervaitons. It is a data imputation approach in both the component construction and the target modeling processes in the extraction of circumstellar objects (exoplanets, disks, etc.). Using one specific property of non-negative matrix factorization, that a fraction of signals can be masked out during both component construction and target modeling, we are able to attribute PSF-only signals to the region where circumstellar objects reside.

## 5.1. Exoplanets

We demonstrate in this paper that planet signals can be properly recovered in an RDI approach. We construct the sNMF components from full field observations of PSF-only signals, and apply the data imputation approach to the target modeling procedure by masking out the planet and recover the PSF signal at the location of the planet. In our simulation in Section 3, we are able to demonstrate that DI-sNMF is able to not only properly recover the injected flux of the planet,
but also obtain smaller uncertainties than KLIP. In this way, DI-sNMF is another independent method to recover planets in addition to the KLIP forward modeling method (KLIP-FM: Pueyo 2016). Nevertheless, DI-sNMF has less assumptions than KLIP-FM: it only requires a 2-dimensional binary mask, rather than a 3-dimensional distribution of light (i.e., a PSF model for a point source).

Similarly as for KLIP-FM, a forward model matched filter (Ruffio et al. 2017) can be applied to the entire field of view to discover and characterize planets for DI-sNMF. In addition, although we do not demonstrate the ADI + SDI capability of DI-sNMF for point sources, rotating and spatially scaling the 2-dimensional binary mask will transform the ADI + SDI problem to a classical RDI problem (Pueyo 2016). Nevertheless, we demonstrate the ADI + SDI scenario for a more complicated problem—circumstellar disks—to show that the DI-sNMF approach is able to achieve the goal of properly recovering the signals, and a point source is just a subset of the more complex circumstellar disk geometry.

## 5.2. Circumstellar Disks

We first demonstrate in this paper that a simple disk model (i.e., an isotropic scattering disk) can be properly recovered in ADI post-processing using DI-sNMF. We then generalize this study, and demonstrate that PSF signals can be properly recovered in various circumstellar structures using simulated ADI approaches. In comparison with point source extraction, more components are needed to better recover the PSF signals.

As an application to disk observations, we perform a precise measurement of the SPF trend for our ADI + SDI post-processing of the debris disk surrounding HR 4796A, finding a tentative trend that the SPFs are more forward scattering as the observation wavelength increases within the *K*1 band. Although a small region on the disk that is closest to the star is not covered due to the limited field rotation (Figure C1d), the final DI-sNMF image is still able to cover the entire field of view, since this region is then fully covered when we combine the derotated reduced images. In comparison with Chen et al. (2020) where no clear trend of the SPFs is observed at different bands using KLIP ADI forward modeling, we argue that it is the limitation of the KLIP ADI forward modeling method: the minor axis (i.e., the most forward scattering region for the disk surrounding HR 4796A) region suffers from self-subtraction with KLIP ADI (e.g., Figure 5c), and forward modeling is unable to probe these values and thus reduce their uncertainties. Therefore the KLIP ADI forward modeling process extrapolates the SPFs for these self-subtraction regions, and such extrapolation depends on the input SPF models which might not be representative of the HR 4796A disk's SPFs. In addition, the uncertainties obtained from forward modeling depends on the

likelihood function, however the correlated noises have to be take into account (e.g., Wolff et al. 2017). Nevertheless, the SPFs extracted from the DI-sNMF images are empirical measurements that have least dependency on such limitations.

To fully eliminate the not covered region for data impputation in ADI + SDI post-processing, the addition of PSF-only exposures (e.g., HST/NICMOS: Choquet et al. 2014; HST/STIS: Ren et al. 2017; Keck/NIRC2: Ruane et al. 2019) is needed, and a careful selection of them (e.g., Ruane et al. 2019) is expected to increase the computational efficiency of DI-sNMF. This joint ADI + SDI + RDI approach is expected to recover the disk image with theoretically negligible alteration of the circumstellar signals when the data impputation region is relatively small (Appendix B). Although we demonstrate the DI-sNMF capability for circumstellar structures that have well-defined geometry, we expect that more complex disk structures could be properly recovered—especially for protoplanetary disks which are currently only properly recovered mainly through PDI (e.g., spiral arm systems: Benisty et al. 2015; Stolker et al. 2016; Monnier et al. 2019).

Another promising way to eliminate the not covered region is the iterative approach presented in Pairet et al. (2018), which sequentially captures the non-PSF structures since the PSF does not rotate with respect to the detector as the telescope rotates. The combination of sNMF with this approach is likely able to remove the not covered regions in the ADI + SDI approach. However such approach is still not able to handle azimuthally symmetric structures since they do not have distinguishable relative motion from the PSF. In this case, using the RDI references is the unique advantage of DI-sNMF, and the face-on structures are expected to be recovered with minor alteration of their signals.

## 5.3. Limitation and Generalization

In this paper, we only masked out ≲15% of the entire field of view to ensure that the DI-sNMF models do not have significant deviation from the true PSF-only signals. In Appendix B.2 we have derived that the area fraction of the masked region has a second order effect on the deviation. Mathematically speaking, assuming the variance contribution from all the pixels in the field of view are equal, they will have a total deviation of less than ~2% (≈ 0.15×0.15). When more areas are masked out (i.e., ≳50%), the DI-sNMF method may lead to significant deviations (i.e., ≳25%)—we do not further investigate this effect in this paper. However, under this scenario, a possible solution is to increase the field of view to reduce the corresponding impact (for example, the clump example in Figure 4).

We note that sufficient field rotation is needed for DI-sNMF in ADI post-processing. To construct the sNMF components for the PSF-only signals, the regions hosting
non-PSF signals have to be flagged by the binary masks. In addition, the PSF signals within the masked regions should not have variation timescale that is faster than the exposure time of each image, otherwise the DI-sNMF components will not be able to represent the PSF signals and fall into the regime of inadequate data imputation. If these two conditions are not met, a possible solution is the sequential capture of rotating (or non-rotating) signals presented in Pairet et al. (2018). In addition, known PSF exposures from available surveys (e.g., Choquet et al. 2014; Ren et al. 2017; Ruane et al. 2019) will help to loosen these requirements. Despite these limitations, we expect that the most promising application of DI-sNMF is to the thin edge-on disk systems where the required field rotation is smaller.

In summary, we have demonstrated that the DI-sNMF approach is expected to properly recover the signals for high contrast imaging observations. It not only is another method for planet characterization, but also opens a new gate to characterize circumstellar disks for integral field spectroscopic observations, and even obtaining their wavelength-dependent image (i.e., reflectance spectra, e.g., Debes et al. 2008; Bhowmik et al. 2019). More general applications of DI-sNMF are to the scenarios where the separation of signals is needed—for example, extracting the spectra for the substellar companions in coronagraphic long-slit spectroscopic observations (e.g., Vigan et al. 2008; Hinkley et al. 2015).

Facilities: HST (STIS), Gemini:South (Gemini Planet Imager)

Software: MCFOST (Pinte et al. 2006, 2009), nmf_imaging (Ren 2018), TinyTim (Krist et al. 2011)

## Acknowledgments

We appreciate the suggestions in both mathematical expression and main text from the anonymous referee that improved this paper. B.R. thanks Nicolas Charon and Bingxiao Xu for useful discussions, and Christophe Pinte for sharing the MCFOST software and commenting on this paper. Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive at the Space Telescope Science Institute. STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. Based on observations obtained at the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnología e Innovación Productiva (Argentina), Ministério da Ciência, Tecnologia e Inovação (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This research project (or part of this research project) was conducted using computational resources (and/or scientific computing services) at the Maryland Advanced Research Computing Center (MARCC: https://www.marcc.jhu.edu). This work was in part funded by contract 61362448-122362 from WFIRST NASA project office via a subcontract through Stanford.

## Appendix A sNMF Component Construction

In this Appendix, we convert the component construction problem for both the ADI and SDI post-processing scenarios to a “missing data”<fnref n="8" /> problem in statistics. This is achieved by assigning zero weight to the known regions containing circumstellar signals (i.e., circumstellar disks or exoplanets).

### A.1 Component Construction with Missing Data

In an ADI post-processing setup, the sky is rotated relatively to the detector frame on the telescope. In the detector frame, the non-rotating non-circumstellar signals are quasi-static; while the circumstellar signals are rotated around the center of the frame with the sky. Given the position and morphology of these circumstellar signals, we minimize their contribution in the component construction procedure for the non-circumstellar signals in Equations (1) and (2), and this

<footnote n="8">⁸ Note: Unless otherwise specified, the “missing data” in this paper refers to the regions that are expected to contain circumstellar signals.</footnote>
is achieved by assigning zero weight to these positions in the weighting matrix V.

The assignment of zero weight for positions hosting circumstellar signals thus converts the ADI problem to a missing data problem in statistics, and the data imputation procedure will be able to fill these regions with PSF-only signals. In this subsection, we prove that the zero-weighted “missing data” has no contribution to the NMF component construction procedure.

**Lemma 1 (Components):** Circumstellar signals do not influence the sNMF component basis of PSF signals when the circumstellar signals are treated as missing data.

**Proof:** We prove this Lemma by way of induction. Let component matrix $H = [H_1^T, H_2^T, \dots, H_n^T]^T$, where $H_i \in \mathbb{R}^{1 \times N_{pix}}$ ($i = 1, 2, \dots, n$) is the i-th NMF component.

a) When $n = 1$, we have row vector $H = H_1 \in \mathbb{R}^{1 \times N_{pix}}$ and column vector $W = W_1 \in \mathbb{R}^{N_{ref} \times 1}$, and the iteration rule in Equation (2), when we ignore the V matrix as in Ren et al. (2018)<fnref n="9" />, is rewritten as

$$H_1^{(k+1)} = H_1^{(k)} \circ \frac{W_1^{(k)T} R}{W_1^{(k)T} W_1^{(k)} H_1^{(k)}}, \quad (A1)$$

where the j-th element in $H_1$ ($j = 1, 2, \dots, N_{pix}$) has its iteration rule:

$$H_{1j}^{(k+1)} = H_{1j}^{(k)} \frac{\sum_{i=1}^{N_{ref}} W_{1i}^{(k)} R_{ij}}{\sum_{i=1}^{N_{ref}} W_{1i}^{(k)} (W_{1i}^{(k)} H_{1j}^{(k)})_{ij}}. \quad (A2)$$

The numerator in Equation (A2) is a real-valued number, since it is the product of a row vector $W_1^{(k)T}$ and the j-th column of R (i.e., $R_{(\cdot),j}$). To illustrate the effect of the weighting matrix, especially on its ability of excluding circumstellar signals, we introduce a binary 0-1 indicator matrix of the same dimension as R to represent the exclusion of the circumstellar signals in the weighting matrix, i.e., $1 \in \mathbb{R}^{N_{ref} \times N_{pix}}$, and for each matrix element indexed $ij$ ($i = 1, 2, \dots, N_{ref}$; $j = 1, 2, \dots, N_{pix}$), $1_{ij} = 0$ means it contains both non-circumstellar and circumstellar signals and $1_{ij} = 1$ means it contains only non-circumstellar signals (i.e., PSF signal only). With this representation, we have $R \equiv R \circ 1$. In this way, the numerator can be rewritten

as

$$\sum_{i=1}^{N_{ref}} W_{1i}^{(k)} R_{ij} = \sum_{i=1}^{N_{ref}} W_{1i}^{(k)} R_{ij} 1_{ij} \quad (A3)$$

$$= \sum_{i=1, 1_{ij}=1}^{N_{ref}} W_{1i}^{(k)} R_{ij} 1_{ij} + \sum_{i=1, 1_{ij}=0}^{N_{ref}} W_{1i}^{(k)} R_{ij} 1_{ij} \quad (A4)$$

$$= \sum_{i=1, 1_{ij}=1}^{N_{ref}} W_{1i}^{(k)} R_{ij}, \quad (A5)$$

where the second term in Equation (A4) is 0, i.e., the zero-valued elements in $R_{(\cdot),j}$ (i.e., when $1_{ij} = 0$) do not influence the numerator, and consequently they do not have influence on the iteration of $H_{1j}$. Given that the zero-valued elements in $1_{(\cdot),j}$ corresponds to the missing data in $R_{(\cdot),j}$, we have the missing data in $R_{(\cdot),j}$ do not have influence on the iteration of $H_{1j}$.

Combining all the columns in R, the zero-weighted values therefore do not have impact on $H_1$.

b) For $n > 1$, assuming the Lemma holds for the m-th component ($m = 1, 2, \dots, n-1$), then for the $(m+1)$-th component, the sequential construction of the NMF components follows (adapted from Appendix C in Ren et al. 2018):

$$H_{m+1}^{(k+1)} = H_{m+1}^{(k)} \circ \frac{W_{m+1}^{(k)T} R}{W_{m+1}^{(k)T} W_{m'}^{(k)} H_{m'}^{(k)} + W_{m+1}^{(k)T} W_{m+1}^{(k)} H_{m+1}^{(k)}}, \quad (A6)$$

where $W_{m'} = [W_1, W_2, \dots, W_m] \in \mathbb{R}^{N_{ref} \times m}$ and $H_{m'} = [H_1^T, H_2^T, \dots, H_m^T]^T \in \mathbb{R}^{m \times N_{pix}}$ are the first m constructed coefficient and component matrices, respectively.

The numerator in Equation (A6) has the same form as that in Equation (A2), therefore following the identical argument as for $H_1$ from Equation (A1) to Equation (A5), then if the Lemma holds for $H_m$, the zero-weighted elements in R do not have impact on $H_{m+1}$.

c) Combining the arguments in a) and b), we have proven that if the components are constructed sequentially, the zero-weighted elements in the reference matrix R do not have impact on the component matrix H. ■

In addition, the circumstellar signal-excluded observations in an observation sequence can be added with their references (i.e., the RDI references) to form a new set of PSF-only observations. Following the identical arguments as above, the sNMF component basis for this set contains only PSF signals.

## A.2 Component Variation Induced by Missing Data

In the previous subsection, we proved that the sNMF components are not influenced by the circumstellar signals when

<footnote n="9">⁹ The proof with V is in principle achieved with the substitutions in Appendix A3 of Blanton & Roweis (2007).</footnote>
they are treated as missing data. However, that does not guarantee that missing data does not influence the component basis. In this subsection, we estimate the sNMF basis difference between two scenarios: one that is constructed when there is no missing data, the other that is constructed with missing data. Specifically, we assume the missing data region of the latter is filled with PSF-only signals in the former scenario.

**Theorem 1 (Components from Missing Data):** When the signal contribution from the masked region is relatively small with respect to the entire field of view, the sNMF component basis does not have significant variation between the non-masked and masked scenarios.

**Proof:** Without loss of generality, we ignore the *V* matrix as in Appendix A.1, and focus on the first column in the coefficient matrix *W* and the first row in the coefficient matrix *H* (i.e., *W*₁ and *H*₁, respectively) when *n* = 1.

First, the *i*-th element of *W*₁ in Equation (1) has a form of

$$ 
\begin{aligned}
W_{1i}^{(k+1)} &= W_{1i}^{(k)} \circ \frac{\sum_{j=1}^{N_{\text{pix}}} R_{ij} H_{1j}^{(k)}}{W_{1i}^{(k)} \left(\sum_{j=1}^{N_{\text{pix}}} H_{1j}^{(k)2}\right)} \\
&= \frac{\sum_{j=1}^{N_{\text{pix}}} R_{ij} H_{1j}^{(k)}}{\sum_{j=1}^{N_{\text{pix}}} H_{1j}^{(k)2}}.
\end{aligned} 
\quad (A7) $$

When we denote the corresponding coefficient matrix generated by missing data with *W*', we have

$$ 
\begin{aligned}
W'_{1i}^{(k+1)} &= \frac{\sum_{j=1,1_{ij}=1}^{N_{\text{pix}}} R_{ij} H_{1j}^{(k)}}{\sum_{j=1}^{N_{\text{pix}}} H_{1j}^{(k)2}} \\
&= W_{1i}^{(k+1)} - \frac{\sum_{j=1,1_{ij}=0}^{N_{\text{pix}}} R_{ij} H_{1j}^{(k)}}{\sum_{j=1}^{N_{\text{pix}}} H_{1j}^{(k)2}}.
\end{aligned} 
\quad (A8) $$

Noticing $W_{1i}$ is applied to $H_1$ to approximate the *i*-th row of *R*, i.e., the *i*-th reference image, when the missing data region in the reference image $R_{i(\cdot)}$ is relatively small and does not have significant variation across different images, i.e.,

$$ \epsilon_{\text{image}} \equiv \frac{\sum_{j=1,1_{ij}=0}^{N_{\text{pix}}} 1}{N_{\text{pix}}} \ll 1, \quad (A9) $$

we can rewrite Equation (A8) as

$$ 
\begin{aligned}
W'_{1i}^{(k+1)} &= W_{1i}^{(k+1)} - o(W_{1i}^{(k+1)}) \\
&= W_{1i}^{(k+1)}[1 - o(1)],
\end{aligned} 
\quad (A10) $$

where $o(\cdot)$ is the little *o* notation, i.e., $|o(x)| \ll |x|$.

Second, the *j*-th element of $H_1$ in Equation (2), similarly as for $W_{1i}$, has a form of

$$ 
\begin{aligned}
H_{1j}^{(k+1)} &= H_{1j}^{(k)} \circ \frac{\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)} R_{ij}}{\left(\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)2}\right) H_{1j}^{(k)}} \\
&= \frac{\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)} R_{ij}}{\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)2}}.
\end{aligned} 
\quad (A11) $$

When we denote the corresponding component matrix generated by missing data with $H'$, we have

$$ H'_{1j}^{(k+1)} = \frac{\sum_{i=1,1_{ij}=1}^{N_{\text{ref}}} W_{1i}^{(k)} R_{ij}}{\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)2}} \quad (A12) $$

$$ = H_{1j}^{(k+1)} - \frac{\sum_{i=1,1_{ij}=0}^{N_{\text{ref}}} W_{1i}^{(k)} R_{ij}}{\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)2}}. \quad (A13) $$

Noticing $H_{1j}$ is applied to $W_1$ to approximate the *j*-th column of *R*, i.e., the *j*-th element in each flattened reference image, when the data missing rate of the *j*-th element in the reference images $R_{i(\cdot)j}$ is relatively small and does not have significant variation across different pixels, i.e.,

$$ \epsilon_{\text{pixel}} \equiv \frac{\sum_{i=1,1_{ij}=0}^{N_{\text{ref}}} 1}{N_{\text{ref}}} \ll 1, \quad (A14) $$

we can rewrite Equation (A13) as

$$ H'_{1j}^{(k+1)} = H_{1j}^{(k+1)}[1 - o(1)]. \quad (A15) $$

Third, we note that in the above connection between $H'_{1j}^{(k+1)}$ and $H_{1j}^{(k+1)}$ in Equation (A15), we assumed that the coefficient matrix $W_1$ was not impacted by the missing data. We now replace $W_1$ with $W'_1$ in Equation (A12) using Equation (A10) to reflect the actual procedure, then we have

$$ 
\begin{aligned}
H'_{1j}^{(k+1)} &= \frac{\sum_{i=1,1_{ij}=1}^{N_{\text{ref}}} W'_{1i}^{(k)} R_{ij}}{\sum_{i=1}^{N_{\text{ref}}} W'_{1i}^{(k)2}} \\
&= \frac{\sum_{i=1,1_{ij}=1}^{N_{\text{ref}}} W_{1i}^{(k)} R_{ij} [1 - o(1)]}{\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)2} [1 - o(1)]^2} \\
&= [1 + o(1)] \frac{\sum_{i=1,1_{ij}=1}^{N_{\text{ref}}} W_{1i}^{(k)} R_{ij}}{\sum_{i=1}^{N_{\text{ref}}} W_{1i}^{(k)2}},
\end{aligned}
$$

where the derivation is kept up to the first order, and this fraction is the identical expression with Equation (A12). When we follow the derivation of Equation (A15) for the above equation, and replace the small *o* notations with $\epsilon_{\text{image}}$ and $\epsilon_{\text{pixel}}$ in Equations (A9) and (A14), we have

$$ H'_{1j}^{(k+1)} = (1 + \epsilon_{\text{image}}) H_{1j}^{(k+1)} (1 - \epsilon_{\text{pixel}}). \quad (A16) $$

In the above equation, we expect that the deviation between $H'$ and $H$ to be first order when either Equation (A9) or Equation (A14) dominates the fraction of missing data. Nevertheless, when the two terms are equal, i.e., $\epsilon_{\text{image}} = \epsilon_{\text{pixel}} = \epsilon$, Equation (A16) can be rewritten as

$$ 
\begin{aligned}
H'_{1j}^{(k+1)} &= (1 - \epsilon^2) H_{1j}^{(k+1)} \\
&= [1 - o^2(1)] H_{1j}^{(k+1)},
\end{aligned} 
\quad (A17) $$
which is a second order deviation. As an example, this is consistent with the Zhu (2016) observation of a few percent deviation in NMF components when 20% of the data were randomly discarded in their Section 3.3, since 20% × 20% = 4%. This deviation is consequently propagated to both W and H in Equation (A8) and Equation (A13), respectively. As a result, we conclude that the deviation of the sNMF components induced by missing data is between the first order and the second order.

## Appendix B sNMF Target Modeling

In this Appendix, we convert the target modeling problem for RDI, ADI, and SDI post-processing scenarios to a “missing data” problem in statistics. Identical to the component construction procedure in Appendix A, this is achieved by assigning zero weight to the known regions containing circumstellar signals.

### B.1 Data Imputation for Missing Data

**Lemma 2 (Modeling):** Zero-weighted “missing data” region in an image can be imputed with PSF-only signals using the sNMF components.

**Proof:** In target modeling, for a given target exposure $T \in \mathbb{R}^{1 \times N_{\text{pix}}}$ with element-wise squared uncertainty $v \in \mathbb{R}^{1 \times N_{\text{pix}}}$, we use the sNMF components to model it with Equation (3). In this step, the regions containing circumstellar signals are ignored by assigning zero weights to them, and we ignore the $v$ matrix for simplicity of proof as in Appendix A. With the introduction of another indicator row vector $1_T \in \mathbb{R}^{1 \times N_{\text{pix}}}$ as in Equation (A3), the coefficient row matrix is updated as follows:

$$ 
\begin{aligned}
\omega^{(k+1)} &= \omega^{(k)} \circ \frac{T H^T}{\omega^{(k)} H H^T} \\
&\equiv \omega^{(k)} \circ \frac{(T \circ 1_T) H^T}{[\omega^{(k)} H \circ 1_T] H^T} \\
&= \omega^{(k)} \circ \frac{(T \circ 1_T) H^T}{[\omega^{(k)} H \circ 1_T] H^T}.
\end{aligned}
\quad (B1) $$

The $i$-th element of $\omega$ ($i = 1, 2, \dots, n$) is thus

$$ \omega_i^{(k+1)} = \omega_i^{(k)} \frac{\sum_{j=1}^{N_{\text{pix}}} T_j 1_{T,j} H_{ij}}{\sum_{j=1}^{N_{\text{pix}}} (\omega^{(k)} H)_j 1_{T,j} H_{ij}}, \quad (B2) $$

where the numerator,

$$ 
\begin{aligned}
\sum_{j=1}^{N_{\text{pix}}} T_j 1_{T,j} H_{ij} &= \sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} T_j H_{ij} 1_{T,j} + \sum_{j=1,1_{T,j}=0}^{N_{\text{pix}}} T_j H_{ij} 1_{T,j} \\
&= \sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} T_j H_{ij}.
\end{aligned}
\quad (B3) $$

Therefore, following the same argument on the relationship between Equation (A5) and Equation (A1), we conclude that

the circumstellar signals do not have impact on the modeling of a target exposure, and both $\omega$ and $H$ contain the contribution from only the non-circumstellar signals.

### B.2 Coefficient Variation Induced by Missing Data

**Theorem 2 (Imputation):** When the signal contribution from the masked region is relatively small, for a given target and NMF basis, the coefficients for the target do not have significant variation between the non-masked and masked scenarios.

**Proof:** When the NMF iteration converges, Equation (35) in Ren et al. (2018) states, for component $H_i$, its coefficient $\omega_i$ has a form of

$$ \omega_i = \frac{T H_i^T}{H_i H_i^T} \frac{1}{1 + \sum_{j=1, j \neq i}^n \frac{\omega_j}{\omega_i} \frac{H_j H_i^T}{H_i H_i^T}}. \quad (B4) $$

Given the sparse coefficients in NMF modeling (Ren et al. 2018), for different components, we have two categories.

First, the coefficient is zero. Due to the non-negativity constraint in Equation (B4), the proof is trivial.

Second, the coefficient is non-zero. If we denote the masked coefficient by $w'_i$, then

$$ \omega'_{i} = \frac{(T \circ 1_T) H_i^T}{(H_i \circ 1_T)(H_i \circ 1_T)^T} \frac{1}{1 + \sum_{j=1, j \neq i}^n \frac{\omega'_j (H_j \circ 1_T)(H_i \circ 1_T)^T}{\omega'_i (H_i \circ 1_T)(H_i \circ 1_T)^T}}. \quad (B5) $$

In this category, the theorem is equivalent to proving

$$ |\omega_i - \omega'_i| = o(\omega_i). \quad (B6) $$

We first focus on the absolute difference of the first term on the righthand side of Equations (B4) and (B5). When we substitute Equation (B3) into the difference, we have

$$ 
\begin{aligned}
& \left| \frac{T H_i^T}{H_i H_i^T} - \frac{(T \circ 1_T) H_i^T}{(H_i \circ 1_T)(H_i \circ 1_T)^T} \right| \\
&= \frac{T H_i^T}{H_i H_i^T} \left| 1 - \frac{(T \circ 1_T) H_i^T}{(H_i \circ 1_T)(H_i \circ 1_T)^T} \frac{H_i H_i^T}{T H_i^T} \right| \\
&= \frac{T H_i^T}{H_i H_i^T} \left| 1 - \frac{H_i H_i^T}{(H_i \circ 1_T)(H_i \circ 1_T)^T} \frac{(T \circ 1_T) H_i^T}{T H_i^T} \right| \\
&= \frac{T H_i^T}{H_i H_i^T} \left| 1 - \frac{\sum_{j=1}^{N_{\text{pix}}} H_{ij}^2}{\sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} H_{ij}^2} \frac{\sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} T_j H_{ij}}{\sum_{j=1}^{N_{\text{pix}}} T_j H_{ij}} \right| 
\end{aligned}
\quad (B7) $$

$$ 
\begin{aligned}
&= \frac{T H_i^T}{H_i H_i^T} \left| 1 - \frac{1 - \frac{\sum_{j=1,1_{T,j}=0}^{N_{\text{pix}}} H_{ij}^2}{\sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} H_{ij}^2}}{1 - \frac{\sum_{j=1,1_{T,j}=0}^{N_{\text{pix}}} T_j H_{ij}}{\sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} T_j H_{ij}}} \right|.
\end{aligned}
\quad (B8) $$

When the summation fractions in Equation (B8) are negligible, i.e., when

$$ \frac{\sum_{j=1,1_{T,j}=0}^{N_{\text{pix}}} H_{ij}^2}{\sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} H_{ij}^2} = o(1) \quad (B9) $$
and

$$ \frac{\sum_{j=1,1_{T,j}=0}^{N_{\text{pix}}} T_j H_{ij}}{\sum_{j=1,1_{T,j}=1}^{N_{\text{pix}}} T_j H_{ij}} = o(1), \qquad (B10) $$

we can rewrite Equation (B8) as

$$ \begin{aligned} & \left| \frac{T H_i^T}{H_i H_i^T} - \frac{(T \circ 1_T) H_i^T}{(H_i \circ 1_T)(H_i \circ 1_T)^T} \right| \\ &= \frac{T H_i^T}{H_i H_i^T} \left| 1 - [1 - o(1)][1 + o(1)] \right| \\ &= o^2 \left( \frac{T H_i^T}{H_i H_i^T} \right). \end{aligned} \qquad (B11) $$

Similarly, following the derivation for Equation (B8), the other terms on the righthand side of Equations (B4) and (B5) can be proven to have negligible deviation. What is more, the proof for Equation (B11) is obtained to the second order, therefore,

$$ |\omega_i - \omega'_i| = o^2(\omega_i). \qquad (B12) $$

Combining the two categories, the proof for the theorem is complete. ■

Notes—The above theorem is valid only when Equations (B9) and (B10) hold. Given the fact that both equations are masked at the same region, when the target, T, resembles the PSF signals, Equation (B10) can be deduced from Equations (B9). We therefore only need to demonstrate that Equation (B9) is valid for direct imaging data using the HR 4796A observations in Section 4. We calculate the radial profile for the DI-sNMF components (see Figure C2 for two examples), then present the squared radial profile that is divided by the sum of the squares for each component in Figure C3.

In Figure C3, we observe that the highest fraction is ~10⁻³, and most of the fraction drops logarithmically when the radial separation increases. Based on these facts, and the fact that most of our masked pixels are far from the central region, the contribution from the masked region is expected to be smaller. For the outlier (i.e., component 3), the profile increases then reaches a plateau, which means that the entries of this component are roughly equally weighted. Equation (B9) therefore holds for this high contrast imaging data. Following similar arguments, for point source modeling in Section 3, this deviation is expected to be negligible since point sources have significantly less impacts than extended structures.

In our reduction for the HR 4796A observations, we have flagged ~15% of the region as “missing data” for data imputation. Equation (B9) is still valid in this case, since the

overall influence measured from the variance is ~10%. Furthermore, the influence is expected to be significantly smaller than ~10%, which may even reach its squared value (i.e., ~1%) since Equation (B12) is a small number to the second order. Even if the influence is at ~10%, it will not change our results in Section 4.2.2, since we expect a systematic offset due to the identical treatment of the data.

## Appendix C Auxiliary Figures

We present the spatially scaled raw exposures for HR 4796A in Figure C1, the sNMF components constructed from the “missing data” approach in Figure C2, and the fractional variance radial profiles for typical components of the HR 4796A observation in Figure C3.

<figure><img src="image_8.png" /><figcaption>Figure C1. Example spatially scaled HR 4796A K1-coron images in detector counts (a, b, and c in different sequence and/or channel numbers), and the field coverage fraction after disk exclusion (d). With the spatial scaling and rotating procedures (for SDI and ADI, respectively), most of the field is covered by PSF-only signals (median: 90%, 3σ lower limit: 5%). Note: (1) the exposure images share the same color bar; (2) only the pure dark regions are always excluded.</figcaption></figure>

## References
<figure><img src="image_9.png" /><figcaption>**Figure C2.** Normalized DI-sNMF component examples. The components are constructed from the disk-excluded and channel-collapsed 240 images of HR 4796A. No disk signal is captured in the components as expected. Note: the images share the same color bar.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure C3.** Pixel-wise fractional squared radial profiles for typical DI-sNMF components for HR 4796A observations in this paper. For each component, the fractional squared radial profile is obtained by calculating the radial profile of the squared component, and divide that profile by the sum of the squared values of the component. For each component, the squared contribution is at most at $10^{-3}$ for one pixel, and most of the values drops logarithmically when the radial separation increases. The assumption in Equation (B9) is thus valid. See the notes in Appendix B.2 for a detailed discussion.</figcaption></figure>

Avenhaus, H., Quanz, S. P., Garufi, A., et al. 2018, ApJ, 863, 44

Benisty, M., Juhasz, A., Boccaletti, A., et al. 2015, A&A, 578, L6

Beuzit, J.-L., Feldt, M., Dohlen, K., et al. 2008, Proc. SPIE, 7014, 701418

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

Bhowmik, T., Boccaletti, A., Thébault, P., et al. 2019, A&A, 630, A85

Blanton, M. R., & Roweis, S. 2007, AJ, 133, 734

Chen, C., Mazoyer, J., Poteet, C., et al. 2020, in preparation

Choquet, É., Pueyo, L., Hagan, J. B., et al. 2014, Proc. SPIE, 9143, 914357

Choquet, É., Perrin, M. D., Chen, C. H., et al. 2016, ApJL, 817, L2

Choquet, É., Milli, J., Wahhaj, Z., et al. 2017, ApJL, 834, L12

Choquet, É., Bryden, G., Perrin, M. D., et al. 2018, ApJ, 854, 53

Christiaens, V., Casassus, S., Absil, O., et al. 2019, MNRAS, 486, 5819

De Rosa, R. J., Nielsen, E. L., Blunt, S. C., et al. 2015, ApJL, 814, L3

Debes, J. H., Ren, B., & Schneider, G. 2019, JATIS, 5, 035003

Debes, J. H., Weinberger, A. J., & Schneider, G. 2008, ApJ, 673, L191

Debes, J. H., Poteet, C. A., Jang-Condell, H., et al. 2017, ApJ, 835, 205

Esposito, T. M., Fitzgerald, M. P., Graham, J. R., et al. 2016, AJ, 152, 85

Esposito, T. M., Duchêne, G., Kalas, P., et al. 2018, AJ, 156, 47

Feldt, M., Olofsson, J., Boccaletti, A., et al. 2017, A&A, 601, A7

Follette, K. B., Rameau, J., Dong, R., et al. 2017, AJ, 153, 264

Garufi, A., Avenhaus, H., Pérez, S., et al. 2020, A&A, 633, A82

Hansen, J. E., & Travis, L. D. 1974, SSRv, 16, 527

Hinkley, S., Bowler, B. P., Vigan, A., et al. 2015, ApJL, 805, L10

Hughes, A. M., Duchêne, G., & Matthews, B. C. 2018, ARA&A, 56, 541

Johnson, R., & Wichern, D. 2007, Applied Multivariate Statistical Analysis (Pearson Prentice Hall)

Kalas, P. G., Rajan, A., Wang, J. J., et al. 2015, ApJ, 814, 32

Keppler, M., Benisty, M., Müller, A., et al. 2018, A&A, 617, A44

Krist, J. E., Hook, R. N., & Stoehr, F. 2011, Proc. SPIE, 8127, 81270J

Lafrenière, D., Marois, C., Doyon, R., et al. 2007, ApJ, 660, 770

Lagrange, A.-M., Langlois, M., Gratton, R., et al. 2016, A&A, 586, L8

Lee, D. D., & Seung, H. S. 2001, in Advances in Neural Information Processing Systems 13, ed. T. K. Leen, T. G. Dietterich, & V. Tresp (MIT Press), 556

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

Macintosh, B., Graham, J. R., Barman, T., et al. 2015, Science, 350, 64

Macintosh, B. A., Graham, J. R., Palmer, D. W., et al. 2008, Proc. SPIE, 7015, 701518

Marois, C., Doyon, R., Racine, R., & Nadeau, D. 2000, PASP, 112, 91

Marois, C., Lafrenière, D., Doyon, R., et al. 2006, ApJ, 641, 556

Marois, C., Macintosh, B., & Véran, J.-P. 2010, Proc. SPIE, 7736, 77361J

Matthews, E., Hinkley, S., Vigan, A., et al. 2017, ApJL, 843, L12

Milli, J., Mouillet, D., Lagrange, A.-M., et al. 2012, A&A, 545, A111

Milli, J., Mawet, D., Pinte, C., et al. 2015, A&A, 577, A57

Milli, J., Vigan, A., Mouillet, D., et al. 2017, A&A, 599, A108

Milli, J., Engler, N., Schmid, H. M., et al. 2019, A&A, 626, A54

Monnier, J. D., Harries, T. J., Bae, J., et al. 2019, ApJ, 872, 122

Mouillet, D., Lagrange, A. M., Beuzit, J. L., & Renaud, N. 1997, A&A, 324, 1083

Müller, A., Keppler, M., Henning, T., et al. 2018, A&A, 617, L2
Paatero, P., & Tapper, U. 1994, *Environmetrics*, **5**, 111

Pairet, B., Cantalloube, F., & Jacques, L. 2018, arXiv, arXiv:1812.01333

Perrin, M. D., Duchene, G., Millar-Blanchaer, M., et al. 2015, ApJ, **799**, 182

Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, **498**, 967

Pinte, C., Ménard, F., Duchêne, G., & Bastien, P. 2006, A&A, **459**, 797

Pueyo, L. 2016, ApJ, **824**, 117
—. 2018, Direct Imaging as a Detection Technique for Exoplanets, 10

Pueyo, L., Crepp, J. R., Vasisht, G., et al. 2012, ApJS, **199**, 6

Rajan, A., Rameau, J., De Rosa, R. J., et al. 2017, AJ, **154**, 10

Ren, B. 2018, nmf_imaging, doi: 10.5281/zenodo.2424378, 1.0, Zenodo

Ren, B., Pueyo, L., Perrin, M. D., et al. 2017, Proc. SPIE, **10400**, **1040021**

Ren, B., Pueyo, L., Zhu, G. B., et al. 2018, ApJ, **852**, 104

Ren, B., Choquet, É., Perrin, M. D., et al. 2019, ApJ, **882**, 64

Ruane, G., Ngo, H., Mawet, D., et al. 2019, AJ, **157**, 118

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

Schneider, G., Gaspar, A., Debes, J., et al. 2017, Enabling Narrow est) IWA Coronagraphy with STIS BAR5 and BAR10 Occulters, Tech. rep., 3

Schneider, G., Weinberger, A. J., Becklin, E. E., et al. 2009, AJ, **137**, 53

Schneider, G., Smith, B. A., Becklin, E. E., et al. 1999, ApJ, **513**, L127

Schneider, G., Debes, J. H., Grady, C. A., et al. 2018, AJ, **155**, 77

Schröder, S. E., Grynko, Y., Pommerol, A., et al. 2014, Icarus, **239**, 201

Smith, B. A., & Terrile, R. J. 1984, *Science*, **226**, 1421

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

Sparks, W. B., & Ford, H. C. 2002, ApJ, **578**, 543

Stolker, T., Dominik, C., Avenhaus, H., et al. 2016, A&A, **595**, A113

Vigan, A., Langlois, M., Moutou, C., & Dohlen, K. 2008, A&A, **489**, 1345

Vigan, A., Bonnefoy, M., Ginski, C., et al. 2016, A&A, **587**, A55

Wahhaj, Z., Milli, J., Kennedy, G., et al. 2016, A&A, **596**, L4

Wolff, S. G., Perrin, M. D., Stapelfeldt, K., et al. 2017, ApJ, **851**, 56

Zhu, G. 2016, arXiv, arXiv:1612.06037