© The Authors 2023

1 Introduction

High-contrast imaging aims to detect and characterize faint signals surrounding bright central sources (e.g., Oppenheimer & Hinkley 2009; Benisty et al. 2023; Currie et al. 2023). To reach this goal, coronagraphs are used to suppress central light, and post-processing – combined with observational techniques – is implemented to remove residual light and speckles (e.g., Pueyo 2018; Follette 2023). With these setups, existing observational techniques and methods are efficiently detecting and characterizing point sources such as planets and brown dwarfs (e.g., Nielsen et al. 2019; Vigan et al. 2021), as well as circumstellar disks and quasar host galaxies in polarized light (e.g., Benisty et al. 2023; Gratadour et al. 2015). Total intensity detection for extended structures, however, is still prone to data reduction artifacts with most existing methods for ground-based observations (e.g., Ruane et al. 2019; Xie et al. 2022). Their characterization is limited by compromised data quality (e.g., Milli et al. 2012; Mazoyer et al. 2020; Olofsson et al. 2023; Xie et al. 2023).

To detect and characterize extended structures in total intensity, several statistics-based methods and their derivative methods are being used (e.g., Lafrenière et al. 2007; Soummer et al. 2012; Amara & Quanz 2012; Ren et al. 2018, 2020; Pairet et al. 2018, 2021; Flasseur et al. 2021; Samland et al. 2021; Juillard et al. 2022; Berdeu et al. 2022). Nevertheless, these methods are either prone to severe overfitting that can alter the morphology and surface brightness of such structures (e.g., Lafrenière et al. 2007; Soummer et al. 2012), or being highly selective on reference regions and computationally intensive (e.g., Ren et al. 2018, 2020), or subject to regularization terms that may provide uncertain results on the morphology and surface brightness of extended structures (e.g., Flasseur et al. 2021; Pairet et al. 2021; Juillard et al. 2022). Meanwhile, a careful selection of data reduction regions (e.g., Galicher & Marois 2011; Milli et al. 2012, 2017; Perrin et al. 2015) could provide better results than nonstatistical methods, yet these regions need to be adjusted for different systems, and such a selection is challenging when the disk morphology is face-on, complex, or even unknown. Using non-negative matrix factorization (NMF; Lee & Seung 2001), Ren et al. (2020) showed that data imputation with sequential NMF (DI-sNMF) could be a promising statistical method for simple disk morphology, and can provide high-quality results with theoretically minimized signal alteration. However, the data imputation concept has only been applied in Ren et al. (2020) for NMF in high-contrast imaging, and it can provide high-quality results only when some strict requirements are satisfied (e.g., Ren et al. 2020; Olofsson et al. 2023; Xie et al. 2023). What is more, DI-sNMF is computationally inefficient due to its iterative nature (e.g., Ren et al. 2018), calling for other analytical methods that can implement the data imputation concept.

In high-contrast imaging, the principal-component-analysis-based Karhunen–Loève image projection (KLIP; Soummer et al. 2012; Amara & Quanz 2012) method is one of the analytical method standards for data reduction. KLIP decomposes reference images into an orthogonal basis through the Karhunen–Loève (KL) transform and has been included in post-processing pipelines (e.g., Wang et al. 2015; Gomez Gonzalez et al. 2016; Stolker et al. 2019; Lucas & Bottom 2020). It is widely used in surveys of both exoplanets (e.g., Nielsen et al. 2019; Vigan et al. 2021) and circumstellar disks (e.g., Esposito et al. 2020; Xie et al. 2022; Cugno et al. 2023; Wallack et al. 2023). Due to its overfitting nature (e.g., Soummer et al. 2012; Pueyo 2016), KLIP encounters difficulty in detecting and characterizing complex and extended structures. While iterative KLIP derivatives (e.g. Pairet et al. 2018; Ginski et al. 2021; Stapper & Ginski 2022) might help to alleviate the over-subtraction and self-subtraction in certain observational setups, there are still persistent residual signals and self-subtraction artifacts that pose challenges to data interpretation. As a data replacement (i.e., imputation) attempt with KLIP, Hunziker et al. (2018) and Xuan et al. (2018) used KLIP-based background subtraction in an interpolation approach: first perform KL decomposition, then conduct KLIP while masking out the central regions. However, due to the masking of central regions, the KLIP procedure within was performed on a non-orthogonal basis. Such a KLIP treatment can introduce overfitting, which is manifested as a higher empirical model for the background, and yields negative surrounding regions even before post-processing. Noticing this overfitting, Xie et al. (2023) thus used DI-sNMF for background removal and postprocessing to extract a double-spiraled system to enable precise spiral motion measurement. After all, a proper data imputation with KLIP would offer a promising way to produce high-quality results with a high computational efficiency.

In the upcoming era of extremely large telescopes (ELTs), the reference differential imaging (RDI) observational setup, including star-hopping (Wahhaj et al. 2021), could likely dominate the observation modes to produce high-quality images of exoplanets, circumstellar disks, and quasar host galaxies (hereafter “astro-physical signals”) in total intensity. In an RDI observation, a reference star (or multiple reference stars) that does not host circumstellar signals is used to empirically capture the point spread function and speckles of the target star or quasar to reveal the target’s astrophysical surroundings. With large collecting areas (e.g., Gilmozzi & Spyromilio 2007), ELTs can reach high signal-to-noise ratios in a significantly smaller amount of time (e.g., Bowens et al. 2021) than existing telescopes. In comparison, it might be suboptimal for the classical angular differential imaging (ADI; Marois et al. 2006) observation setup to be deployed on ELTs, since ADI requires a sufficient amount of field rotation of the sky for data processing (e.g., Milli et al. 2012; Xuan et al. 2018), yet the sky rotation rate is a natural phenomenon that cannot be changed. To reduce and process the ELT data with high fidelity in RDI observations, iterative methods could potentially offer the best results. However, these iterative methods might be limited due to their significantly high computational cost compared to that of analytical methods, especially when there are high-resolution observations with a large number of datasets. Therefore, an analytical RDI data reduction method is needed to provide high-quality results in an efficient way.

In this study, by engineering the mathematical basics for KLIP, we enable it with the data imputation concept (DIKL). By modifying two steps in the standard KLIP procedure, and illustrating the DIKL results, we expect DIKL to be a promising method in RDI data reduction. In comparison with the established iterative DI-sNMF method that offers high-quality RDI results (e.g., Olofsson et al. 2023; Xie et al. 2023), such modifications allow an analytical data imputation by DIKL. It is computationally more efficient than iterative methods such as DI-sNMF by ~3 orders of magnitude, and it can provide results with a comparable quality to DI-sNMF in RDI datasets.

2 Method

In high-contrast imaging observations that require the usage of coronagraphs and adaptive optics systems, light signals from a central object (e.g., star, quasar) received on the detectors are superimposed on the astrophysical signals (e.g., Pueyo 2018; Follette 2023). The goal of post-processing is to extract the faint astrophysical signals (e.g., exoplanets, circumstellar disk, host galaxy) that are hidden behind these bright signals of non-astrophysical origin (e.g., adaptive optics halo, speckles, thermal background). For simplicity, we refer to all these nonastrophysical signals as “speckles” hereafter. We refer to an image that contains astrophysical signals as the target, t. We refer to a set of images that only contain speckle signals as reference images, R¹. In post-processing, we use the reference images to obtain the representative features of the speckles, then remove these features from a target image to reveal astrophysical signals in the residuals.

2.1 Overview of KLIP

For a reference image array containing n_ref ∊ ℕ reference images, each with n_pix ∊ ℕ pixels, we denote the array as $R\in {\text{ℝ}}^{n_{\text{pix}}\times n_{\text{ref}}}$. A column in R contains a reference image that is converted to a one-dimensional column vector (e.g., concatenating all columns in one image). The notations here follow those of statisticians and computer scientists, instead of existing publications by astronomers (e.g., Soummer et al. 2012; Pueyo 2016), in order to enable an efficient delivery of messages to non-astronomy readers in coding. Specifically, we have the reference array $R=\left[R_1,\cdots ,R_{n_{\text{ref}}}\right]$, with $R_i={\left(R_{i1},\cdots ,R_{i{n}_{\text{pix}}}\right)}^{\text{⊤}}$ for i ∊ {1, …, n_ref} and ⊤ for matrix transpose, where R_ij ∊ ℝ for j ∊{1,…, n_pix}.

For a target image, we can convert it to a column vector, $t\in {\text{ℝ}}^{n_{\text{pix}}\times 1}$. To remove the speckles from the target image, a statistics-based post-processing method first transforms the reference array, R, to obtain a basis. This basis contains the representative features of the reference array (i.e., the components). We then remove the representative features from the target vector. After post-processing, astrophysical signals are expected to reside in the residual image from feature removal.

The KLIP algorithm in Soummer et al. (2012) consists of two steps: the Karhunen-Loève (KL) transform of the reference matrix, R, to obtain the KL component basis, and the projection of the target vector, t, on the KL basis. For KLIP, the columns in both R and t have zero spatial mean. To reveal the astrophysical signals, the KLIP projection is then removed from the target image. To distinguish the difference between KLIP and DIKL later in this work, here we denote the reference array for KLIP as matrix $R^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(a\right)}\times n_{\text{ref}}}$, with $n_{\text{pix}}^{\left(a\right)}\le n_{\text{pix}}$ allowing for the selection of smaller regions in reference images. Similarly, we denote the target vector as $t^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(a\right)}\times 1}$ in KLIP processing. With these notations, we review the KLIP procedure as follows.

In the KL transform step, for a reference array, R^(a), its covariance matrix ${\text{Σ}}^{\left(a\right)}=R^{\left(a\right)\text{⊤}}{R}^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(a\right)}\times n_{\text{pix}}^{\left(a\right)}}$ is real symmetric. The spectral decomposition of Σ^(a) is $${\text{Σ}}^{\left(a\right)}=R^{\left(a\right)\top }{R}^{\left(a\right)}=Q^{\left(a\right)}{D}^{\left(a\right)}{Q}^{\left(a\right)\top },$$(1)

where $D^{\left(a\right)}=\text{diag}\left\{{\lambda }_1^{\left(a\right)},\cdots {\lambda }_{n_{\text{ref}}}^{\left(a\right)}\right\}$, a diagonal matrix, the diagonal entries of which are the eigenvalues of Σ^(a), with ${\lambda }_k^{\left(a\right)}\ge {\lambda }_{k+1}^{\left(a\right)}>0$, and $Q^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{ref}}\times n_{\text{ref}}}$ is an orthonormal matrix with its columns being the corresponding eigenvectors with Q^(a)⊤ Q^(a) = Q^(a) Q^(a)⊤ = I. Left multiply Eq. (1) by Q^(a)⊤ and right multiply it by Q^(a) and we have $${\left[R^{\left(a\right)}{Q}^{\left(a\right)}\right]}^{\top }{R}^{\left(a\right)}{Q}^{\left(a\right)}=Q^{\left(a\right)\top }{R}^{\left(a\right)\top }{R}^{\left(a\right)}{Q}^{\left(a\right)}=D^{\left(a\right)},$$

where we now have diagonalized the covariance matrix. By defining $$Z^{\left(a\right)}=R^{\left(a\right)}{Q}^{\left(a\right)},$$(2)

we have an orthogonal KL component basis, $Z^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(a\right)}\times n_{\text{ref}}}$. The columns of the KL basis can be divided by the corresponding square root of the eigenvalues (e.g., Soummer et al. 2012), and thus create an orthonormal basis, Z^(a). While we do not adopt this normalization since it does not impact the following target modeling procedure in the datasets used later in this study, such a normalization (e.g., VIP; Gomez Gonzalez et al. 2016) might be necessary for other datasets.

In the KL projection step, for a given target image vector with zero spatial mean $t^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(a\right)}\times 1}$, we can project and remove the KL component to obtain the residual $r^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(a\right)}\times 1}$, $$r^{\left(a\right)}=t^{\left(a\right)}-{\displaystyle \sum _{k=1}^{K_{\text{klip}}}{c}_k^{\left(a\right)}{Z}_k^{\left(a\right)}},$$(3)

where K_kıip ∊ {1, …, n_ref} is the cutoff count of the KL components (Soummer et al. 2012), and $$c_k^{\left(a\right)}=t^{\left(a\right)\top }{Z}_k^{\left(a\right)}$$(4)

is a scalar with $c_k^{\left(a\right)}\in \text{\hspace{0.17em}ℝ}$, which denotes the projection from the target onto the k-th KL component of R^(a).

In the next step of this study, we modify the KLIP procedure to introduce DIKL post-processing. To perform DIKL data reduction, we need to focus on different regions of the reference array, R. For the full reference array, R, we can reorder or duplicate certain rows of it to obtain two submatrices, $$\{\begin{array}{l}{R}^{\left(a\right)}=S^{\left(a\right)}R\hfill \\ R^{\left(b\right)}=S^{\left(b\right)}R\hfill \end{array},$$(5)

where $R^{\left(a\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(a\right)}\times n_{\text{ref}}}$ and $R^{\left(b\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(b\right)}\times n_{\text{ref}}}$, with $n_{\text{pix}}^{\left(a\right)}\le n_{\text{pix}}$ and $n_{\text{pix}}^{\left(b\right)}\le n_{\text{pix}}$. Selection matrices $S^{\left(a\right)}\in {\text{𝔹}}^{n_{\text{pix}}^{\left(a\right)}\times n_{\text{pix}}}$ and $S^{\left(b\right)}\in {\text{𝔹}}^{n_{\text{pix}}^{\left(b\right)}\times n_{\text{pix}}}$ are boolean, and they are used to select specific rows in R to form R^(a) and R^(b), respectively. For example, to select the jth row from R and store it in the i-th row of R^(a), we have $S_{ij}^{\left(a\right)}=1$ (and $S_{ij}^{\left(a\right)}=0$ otherwise). To be more informative on the naming conventions of the superscripts in Eq. (5), we refer the two selected matrices as an anchor matrix R^(a) and a boat matrix R^(b). The anchor matrix, R^(a), covers the regions that only host speckle signals, and the boat matrix, R^(b), can host astrophysical signals (and can also contain the anchor matrix that is adopted in this study); see Fig. 1 for an illustration of the selected regions.

Fig. 1 Partitioning of an image for DIKL. (a) An anchor image used to construct the KL basis. (b) A boat image used to construct the DIKL basis. (c) A complete image without partitioning. We note that the images here are vectorized to constitute the columns of the R(a), R(b), and R matrices, respectively.

2.2 DIKL for reference differential imaging

The DIKL method modifies KLIP in two respects. Specifically, the DIKL components are constructed by applying the eigenvectors from the anchor matrix, R^(a), to the boat matrix, R^(b), and the DIKL projection adopts the projection coefficient between the anchor target, t^(a), and the anchor KL basis, Z^(a).

On the one hand, DIKL modifies the KL transform in Eq. (2). The DIKL basis follows $${Z^{\prime }}^{\left(b\right)}=R^{\left(b\right)}{Q}^{\left(a\right)},$$(6)

where we now instead use the eigenvectors from R^(a) to construct the DIKL basis for the R^(b) regions. Here we use the prime symbol for the DIKL basis to distinguish it from a normal KL basis for the reference matrix, R^(b). By comparison, the KL basis for R^(b) is Z^(b), which is neither calculated nor used in this study. On the other hand, we modify the KL projection in Eq. (3).

For a target vector $t^{\left(b\right)}\in {\text{ℝ}}^{n_{\text{pix}}^{\left(b\right)}\times 1}$ with zero spatial mean, DIKL uses the coefficients from Eq. (4) between t^(a) and Z^(a) then applies the coefficients to the DIKL basis in Eq. (6). Specifically, the residual image after DIKL projection is $$\begin{array}{l}{r}^{\left(b\right)}=t^{\left(b\right)}-{\displaystyle \sum _{k=1}^{K_{\text{klip}}}{c}_k^{\left(a\right)}{Z^{\prime }}_k^{\left(b\right)}}\hfill \\ =t^{\left(b\right)}-{\displaystyle \sum _{k=1}^{K_{\text{klip}}}{t}^{\left(a\right)\top }{Z}_k^{\left(a\right)}{Z^{\prime }}_k^{\left(b\right)}}\hfill \\ =|S^{\left(b\right)}t\rangle -{\displaystyle \sum _{k=1}^{K_{\text{klip}}}\langle S^{\left(a\right)}t|Z_k^{\left(a\right)}\rangle |{Z^{\prime }}_k^{\left(b\right)}\rangle },\hfill \end{array}$$(7)

where the last equation adopts the bra-ket convention for the purpose of demonstration and connection with existing studies (e.g., Soummer et al. 2012).

With the above procedure, Eqs. (6) and (7) establish the DIKL post-processing procedure; see Appendix A for a corresponding pseudocode and implementation instructions. On the one hand, there is no direct projection of t^(b) onto the KL basis of R^(b), which ensures the non-overfitting of astrophysical signals for DIKL. On the other hand, the projection coefficient is obtained from where speckle-only signals are expected (in both the reference images and the target image). Both ensure that astrophysical signals in the boat regions of t^(b) are not captured in the coefficients in Eq. (7), and thus it serves as the data imputation step for DIKL. In fact, DIKL only performs the KL transform to the speckles in the anchor matrix in this work, which avoids the high variation of signals in the boat matrix (e.g., regions that are close to the coronagraph) to dominate the classical KL transform (that is performed on the entire image); see Fig. 1.

We note that Eq. (6), however, is limited in the reduction of non-RDI data. For example, in ADI datasets, we need to mask out different pixels in an observation sequence (e.g., Milli et al. 2017; Ren et al. 2020) to focus on speckle-only signals, since angular diversity is used to capture the speckles in an ADI observation sequence. We can indeed calculate an element in the covariance matrix by focusing only on the pixels that are not masked out in a pair of reference images, and perform spectral decomposition in Eq. (1) to obtain the corresponding eigenvectors. However, we cannot yet construct a basis following Eq. (2) for KLIP, nor Eq. (6) for DIKL, since weights in the eigenvector matrix, Q, are assigned to data that are masked out (i.e., artificially “missing”) in the reference array, R. To apply DIKL to ADI datasets, we can interpolate the speckle-only signals for the masked out data (e.g., Perrin et al. 2015) or iteratively fill the missing data (e.g., Bailey 2012) in R for Eq. (6), yet such treatments are beyond the current scope of this analytical study.

3 Application

To demonstrate the DIKL method using on-sky RDI data, we retrieved datasets from the SPHERE instrument at the Very Large Telescope (VLT). We adopted data from IRDIS (Dohlen et al. 2008) of SPHERE in the star-hopping mode (Wahhaj et al. 2021). The star-hopping mode hops between a target star and its reference star during an observation sequence, and thus achieves a quasi-simultaneous capture of speckle change. We retrieved the UT 2021-09-06, UT 2022-02-07, and UT 2022-04-01 observations of circumstellar disk hosts - HD 169142, PDS 201, and HD 129590 – as well as their reference stars from programs 105.209E and 105.20GP. We preprocessed the data using the IRDAP pipeline (van Holstein et al. 2017, 2020) and the stars were then located at the center of the images. In the IRDIS data, with each IRDIS pixel having a spatial scale of 12.25 mas (Maire et al. 2016), the control ring - the interior to which the adaptive optics system can perform optimal dark hole correction to reveal faint objects – spans between 85 and 115 pixels from the matrix centers. For further processing, we masked out the matrix elements that are interior to a radius of 8 pixels from the center. We display the partitioning of an image in Fig. 1 for illustration of the anchor region and boat region in this study.

For the three groups of IRDIS datasets, the images for postprocessing are the central 350 × 350 pixel from the IRDAP preprocessed data. For one exposure, there are two channels from the preprocessed data, and we added them to produce one image. There are 128, 32, and 32 target images, with 36, 12, and 32 corresponding reference images, for HD 169142, PDS 201, and HD 129590, respectively. To convert a 2D image to a column in a matrix for this study, in practice, we created a 2D binary mask, the entries of which are 1 when the corresponding pixels are selected in Figs. 1a,b. For one image, we used the where function in numpy (Harris et al. 2020) to select the pixels in the images, and created a column for the selected matrix. By performing this for all reference images, we created the selected anchor matrix, R^(a), or the boat matrix, R^(b), for further postprocessing. By performing this once for a target image, we created the selected anchor target vector, t^(a), or the boat target vector, t^(b).

In post-processing, for the three IRDIS datasets that are in matrix form, we first followed Eq. (1) to perform spectral decomposition of the matrix elements in the anchor matrix (i.e., the IRDIS control ring). We then followed Eq. (6) to conduct DIKL transform on the boat matrix – the entire image (or the matrix elements interior to the outer edge of the control ring for the HD 129590 data) – using the eigenvectors from the covariance matrix of the anchor references. At last, to remove the speckles for a target image, we followed Eq. (7) to perform DIKL reduction, whereby we adopted the KLIP projection coefficients from Eq. (4) between the anchor target and anchor references. We reshaped the 1D vectors into a 2D image, then derotated the reduction results for all target images to north-up and east-left according to their corresponding parallactic angles calculated using Pynpoint (Stolker et al. 2019). We calculated the element-wise median of the derotated reduction results as the combined image. We adopted the median-subtracted combined image as the final result.

3.1 Residual variance

The variance of residual images can be informative of the existence of astrophysical signals. By performing DIKL on the reference images of the HD 169142 dataset, we generated the corresponding fractional residual variance (FRV) curves (e.g., Soummer et al. 2012; Ren et al. 2018). Specifically, for a reduced reference image, we divided its variance by that of the corresponding original image to generate the FRV. The FRV curves are then the FRV dependence as a function of the KL components. Similarly, we generated the FRV curves for KLIP for comparison. For KLIP, FRV curves are expected to follow the fractional residual eigenvalues (Soummer et al. 2012), as can be seen in Fig. 2.

The FRV curves of DIKL converge to higher plateaus than those from KLIP. This convergence illustrates the less aggressive overfitting (or potentially non-overfitting of data; e.g., NMF; Ren et al. 2018) of DIKL. As with NMF, the FRV plateaus with DIKL are the indirect evidence that DIKL can retain more information on non-speckle signals than KLIP. In addition, the FRVs for DIKL reach stable values with ≈ 5 DIKL components, suggesting that the rest of the components have a negligible contribution in the DIKL process.

Due to the non-orthogonality of a DIKL basis, we witnessed a minor level of overfitting of the data, and thus we subtracted the median of the DIKL residuals to manually offset this effect. We present in Fig. 3 the correlation matrix for DIKL components: there exist crosstalks in the form of non-zero off-diagonal elements. However, given that the contribution of higher-order DIKL components (≳5) is negligible in the FRVs in Fig. 2, the DIKL crosstalk only impacts the first few components, and the crosstalk does not contribute to data reduction beyond the first ≈5 DIKL components due to the reaching of FRV plateaus in Fig. 2.

Fig. 2 Fractional residual variance as a function of KL components for RDI. DIKL reaches higher plateaus than KLIP, illustrating the information retention capability of DIKL.

Fig. 3 Correlation matrix of DIKL components. While KL components are mutually orthogonal, most DIKL components are not. However, with the FRVs in Fig. 2 suggesting that only the first ~5 DIKL components contribute to the removal of speckles, the crosstalk from higher order components therefore does not impact DIKL data reduction.

3.2 DIKL imaging

We used DIKL to reduce the star-hopping RDI observations of HD 169142, PDS 201, and HD 129590, which are known circumstellar disk hosts. The three disks have inclinations from nearly face-on to roughly edge-on (e.g., Pohl et al. 2017; Wagner et al. 2020; Matthews et al. 2017; Olofsson et al. 2023). We also reduced the datasets with KLIP and DI-sNMF for comparison.

Between the DIKL and KLIP results in Fig. 4 for the nearly face-on HD 169142, DIKL retrieves a two-ringed system. KLIP can nevertheless only recover the inner ring that is close to the coronagraph, with compromised data quality. The DIKL result resembles the disk image obtained in polarized light, which thus demonstrates its superiority over KLIP in conserving face-on structures. Similarly, the PDS 201 and HD 129590 results have fine extended structures only seen in DIKL and not in KLIP. By comparison, KLIP removes signals from the disks, altering disk morphology, which poses challenges for data interpretation (e.g., Wagner et al. 2020; Olofsson et al. 2023).

In comparison with the DI-sNMF method, which is mathematically well founded to deliver high-quality results (e.g., Ren et al. 2020; Olofsson et al. 2023; Xie et al. 2023), DIKL can provide satisfactory results in Fig. 5. The DIKL results not only resemble the morphology of disks from DI-sNMF (Olofsson et al. 2023; Ren et al. 2023) but also agree with the DI-sNMF surface brightness values within ~10% for the disk-hosting regions. This further shows that DIKL can be a promising method in reaching high-quality results, and with a computational efficiency that is ~3 orders of magnitude better than DI-sNMF.

With the demonstrated superiority of DIKL over the classical KLIP method in Fig. 4, as well as the consistency of its results with the established DI-sNMF method in Fig. 5, DIKL is a promising method that can produce high-quality results using a similar amount of calculation time as KLIP. However, with the DIKL basis being a non-orthogonal basis in practice, we should not yet ascribe full credibility in the fine details of DIKL results for interpretation until its application to datasets from other instruments is validated. Nevertheless, DIKL possesses a distinct advantage in producing high-quality preliminary outcomes with high computational efficiency as the first step for data analysis. Such datasets with promising results can then be reduced with more advanced methods including DI-sNMF to ensure signal quality.

Fig. 4 Comparison of reduction results using KLIP (left) and DIKL (right) for nearly face-on to roughly edge-on systems. For HD 169142 (left), PDS 201 (middle), and HD 129590 (right), DIKL reaches a data quality that might be comparable with polarized light observations.

Fig. 5 Validation of DIKL results (top) using DI-sNMF results (middle), with the former method being more computationally efficient than the latter by ~3 orders of magnitude. The difference between the two methods (bottom), obtained by first subtracting the DI-sNMF result from that of DIKL then dividing it by the DI-sNMF result, is ≲10% for the disk signals.

4 Summary

For RDI data reduction in high-contrast imaging, we demonstrated a high-efficiency and analytical approach to perform data imputation with a modified KLIP algorithm: DIKL. Specifically, we modified two steps in KLIP to reach the purpose of data imputation. On the one hand, in component construction, we used the eigenvectors of the covariance matrix – from the regions that only host speckle signals (i.e., the anchor matrix, “a”) – to generate the DIKL basis for the regions hosting non-speckle signals (i.e., the boat matrix, “b”); see Eq. (6). On the other hand, in speckle removal, we adopted the coefficients from capturing the speckles in the anchor matrix, and applied them to the DIKL basis. We thus captured the speckles for the regions that host astrophysical signals in the boat matrix; see Eq. (7). With the two modifications, we only needed to perform spectral decomposition once, in Eq. (1). What is more, the corresponding covariance matrix for the anchor matrix is faster to compute than KLIP, due to a reduced number of matrix elements. As a result, the computational complexity of DIKL is similar to or less than that of KLIP.

By avoiding the projection of astrophysical signals onto speckle features, DIKL can recover face-on structures that are normally overfit in RDI observations (see Fig. 4). On the one hand, in comparison with the mathematically well-founded iterative DI-sNMF method in Ren et al. (2020), DIKL might potentially provide analytical results of low quality due to the non-orthogonal DIKL basis in Eq. (6). However, DIKL can approach a data quality similar to that of the latter in Fig. 5, since the crosstalk of the DIKL basis is negligible for higher-order terms in Figs. 2 and 3. On the other hand, in comparison with the Hunziker et al. (2018) modification of KLIP, DIKL does not apply the KL transform to the regions interior to the control ring. Given that the KL transform captures the variance of the signals (e.g., Soummer et al. 2012), while high-contrast imaging observations have the highest variance next to the coronagraph (e.g., Pueyo 2018), the Hunziker et al. (2018) and Xuan et al. (2018) approach should not be adopted to realize the data imputation concept for KLIP. In comparison with that approach, DIKL is not sensitive to such variances, since it uses the control ring signals that are less prone to random noise than the KL transform, and thus DIKL is theoretically more plausible for data imputation. As a result, DIKL can be a promising analytical method that provides initial results, and with a computational efficiency ~3 orders of magnitude higher than the existing iterative methods (e.g., Ren et al. 2020). However, due to the crosstalk of the DIKL basis when the reference images are not stable, we recommend using DIKL for initial signal detection followed by other iterative data imputation methods for detailed characterization.

Given that DIKL is a natural extension of the KLIP algorithm, it can be implemented in the existing high-contrast imaging packages (e.g., Wang et al. 2015; Gomez Gonzalez et al. 2016; Stolker et al. 2019; Lucas & Bottom 2020) for reference differential imaging data reduction. What is more, DIKL can perform data reduction for images containing negative values (from background removal), and thus it might potentially extract fainter disks than DI-sNMF, which only takes in non-negative values. In the upcoming ELT era, we expect DIKL and its future derivative methods to provide high-quality results with high computational efficiency. Moving forward, on the one hand, we can assign different weights to the pixels in KL transform (e.g., Bailey 2012) to make it less prone to random noise or shot noise. On the other hand, for non-RDI data (e.g., angular differential imaging and spectral differential imaging), more careful derivations that can handle missing data for different regions from different images, including modifying the covariance matrix for the KL transform and the KLIP procedure, are needed in the future.

Acknowledgements

We thank the anonymous referee for their suggestions that increased the readability and clarity of this manuscript, Mamadou N’Diaye for discussions, and Chen Xie for discussions on reference differential imaging in the ELT era. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 105.209E and 105.20GP. We thank the PIs for the two programs (M. Benisty and J. Olofsson) for the datasets that validated this study. This research has received funding from the European Union’s Horizon Europe research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 101103114.

Appendix A Pseudocode for DIKL implementation

We present a pseudocode to implement DIKL in Algorithm 1. Specifically, we used the KL projection coefficient from Eq. (4), and applied it to the DIKL basis in Eq. (6), to obtain the DIKL projection. We then removed the DIKL projection from the target image to reveal the astrophysical signals in Eq. (7).

To implement a standalone DIKL function, we need two matrix operations: matrix multiplication and eigendecomposition; both are available in modern scientific programming languages. Alternatively, to implement DIKL in the field of high-contrast imaging, we can use existing KLIP frameworks (e.g., pyKLIP: Wang et al. 2015, VIP: Gomez Gonzalez et al. 2016, Pynpoint: Stolker et al. 2019, ADI.jl: Lucas & Bottom 2020), since there are only two additional commands (at the end of Algorithm 1) in addition to KLIP for RDI observations. There is one extra command to select specific regions to obtain the anchor and boat matrices, and this selection is available using the mask commands in existing frameworks.

In the actual implementation, we recommend splitting Algorithm 1 into two functions for computational efficiency. One function performs KL and DIKL transforms, and returns both the KL basis from Eq. (2) and the DIKL basis from Eq. (6). The other uses the KL basis to generate the KL projection coefficients from Eq. (4) for target t, then applies the coefficients to the DIKL basis to obtain the residuals in Eq. (7). In this way, we can use the same reference array, R, for the RDI reduction of different targets using DIKL.

To implement DIKL on existing high-contrast imaging pipelines, we detail the key modifications here. Specifically, in the KL transform we need the eigenvector matrix Q^(a) in Eq. (1), and thus such an output is needed when eigendecomposition is conducted (e.g., pyKLIP², Pynpoint³, and VIP⁴); for adi.jl⁵, where singular value decomposition is conducted, we can use the output directly. We can then multiply the boat reference matrix, R^(b), with the eigenvector matrix to obtain the DIKL basis, Z′^(b), in Eq. (6). At last, we multiply the KL projection coefficients, c^(a), from the pipelines for Eq. (4) with the DIKL basis Z′^(b), then obtain the residuals in Eq. (7) after DIKL projection.

References

All Figures

	Fig. 1 Partitioning of an image for DIKL. (a) An anchor image used to construct the KL basis. (b) A boat image used to construct the DIKL basis. (c) A complete image without partitioning. We note that the images here are vectorized to constitute the columns of the R(a), R(b), and R matrices, respectively.
In the text

	Fig. 2 Fractional residual variance as a function of KL components for RDI. DIKL reaches higher plateaus than KLIP, illustrating the information retention capability of DIKL.
In the text

	Fig. 3 Correlation matrix of DIKL components. While KL components are mutually orthogonal, most DIKL components are not. However, with the FRVs in Fig. 2 suggesting that only the first ~5 DIKL components contribute to the removal of speckles, the crosstalk from higher order components therefore does not impact DIKL data reduction.
In the text

	Fig. 4 Comparison of reduction results using KLIP (left) and DIKL (right) for nearly face-on to roughly edge-on systems. For HD 169142 (left), PDS 201 (middle), and HD 129590 (right), DIKL reaches a data quality that might be comparable with polarized light observations.
In the text

	Fig. 5 Validation of DIKL results (top) using DI-sNMF results (middle), with the former method being more computationally efficient than the latter by ~3 orders of magnitude. The difference between the two methods (bottom), obtained by first subtracting the DI-sNMF result from that of DIKL then dividing it by the DI-sNMF result, is ≲10% for the disk signals.
In the text