Karhunen–Loève Image Processing (KLIP)

• Karhunen–Loève Image Processing (KLIP) is a PCA-based algorithm that models and subtracts stellar PSFs to enhance high-contrast astronomical imaging.
• It projects science images onto an orthonormal basis of eigenimages derived from reference PSFs, optimizing signal recovery for exoplanets and circumstellar disks.
• Efficient implementations leverage GPU acceleration and downdate techniques, ensuring analytical throughput estimation and unbiased astrophysical measurements.

Karhunen–Loève Image Processing (KLIP) is a data-driven algorithm for high-contrast astronomical imaging, leveraging the Karhunen–Loève transform (or Principal Component Analysis, PCA) to optimally model and subtract the stellar point spread function (PSF). KLIP has emerged as a cornerstone technique for exoplanet and circumstellar disk detection and characterization, providing a unified formalism that offers both high-sensitivity detection and analytical guarantees for unbiased astrometry and photometry. KLIP projects science target images onto an orthonormal basis of eigenimages constructed from a library of reference PSFs, yielding robust PSF subtraction and enabling forward modeling, crucial for precision astrophysical inference (Soummer et al., 2012).

1. Mathematical Framework and KLIP Workflow

KLIP operates by organizing a set of $N$ reference PSF images, each vectorized into length $P$, into a data matrix $X\in\mathbb{R}^{P\times N}$, with columns mean-subtracted over the target region $\mathcal{S}$ to create a zero-mean library. The essential steps are:

1. Covariance Construction: Calculate the pixel–pixel covariance matrix

$C = \frac{1}{N-1} X X^T \in \mathbb{R}^{P\times P}$

2. Eigen-Decomposition: Diagonalize $C$,

$C = V \Lambda V^T$

where $V$ gives the orthonormal eigenimages $v_j$ and $\Lambda$ the sorted eigenvalues.

1. Basis Truncation: Select the dominant $k$ eigenimages, $V_k = [v_1, ..., v_k]$.
2. Science Image Projection: For a zero-mean target $y$, compute projection coefficients $a = V_k^T y$, reconstruct the synthetic PSF $\hat y = V_k a$, and subtract to yield the residual $r = y - \hat y$.

The choice of $k$ (number of modes) controls the suppression of stochastic PSF noise (speckles) versus the retention of astrophysical signal (planet, disk). The throughput of point-like signals remains high even for moderate $k$ due to near-orthogonality with dominant speckle modes. The signal-to-noise ratio (SNR) can be assessed analytically as a function of $k$, frequently optimized via a grid search or residual-variance plateauing (Soummer et al., 2012, Ko et al., 2024, Ren, 2023).

2. Implementation Techniques and Computational Advances

For large datasets, the covariance decomposition is performed in the lower-dimensional reference-frame space (size $N\times N$ versus pixel space $P\times P$), with KL eigenimages in pixel space constructed as scaled linear combinations of the reference images:

$$u_j = \frac{1}{\sqrt{N \lambda_j}} X_c v_j, \ \ (u_j \in \mathbb{R}^P)$$

(Ko et al., 2024). Efficient implementations use singular value decomposition (SVD) and exploit GPU acceleration, as in the torchKLIP package built on PyTorch, enabling rapid eigendecomposition, projection, and model subtraction for datasets up to $10^5$ images. PyTorch or other ML frameworks allow batched processing, seamless CPU/GPU switching, memory-efficient covariance computation, and extension towards ML-based denoisers (Ko et al., 2024).

KLIP lends itself to zone-based and annular application, in which $\mathcal{S}$ is subdivided to allow spatially varying optimization—directly analogous to the LOCI algorithm's segmentation scheme—but remains fully linear and admits analytic throughput estimation (Soummer et al., 2012).

3. Algorithmic Enhancements: Downdate and Data Imputation Variants

KLIP's computational bottleneck under Angular Differential Imaging (ADI) or Spectral Differential Imaging (SDI) arises when the reference library for each frame must exclude the target itself to prevent astrophysical self-subtraction. Naively recomputing the eigendecomposition for every frame is prohibitive for large $N$. Efficient downdate algorithms circumvent this by modifying the SVD of the full reference set to exclude a column (“downdate”), reducing per-frame cost from $O(N^2)$ to $O(1)$ in $N$ (Long et al., 2021). The downdate algorithm uses Brand's matrix identities to update the SVD with only $O(k^3+p k^2)$ operations per frame, enabling full-data-rate processing for $N\sim10^4$ observations.

KLIP has also been extended with a Data-Imputation mechanism (DIKL) addressing self-subtraction and over-subtraction of disk signals. DIKL partitions the data into “anchor” (speckle-only) and “boat” (regions of astrophysical interest) regions, performs KL projection on the anchor, and imputes the resulting PSF model onto the boat. This analytic approach achieves comparable residual variance to iterative methods (e.g., DI-sNMF) at $\sim10^3\times$ less computational cost and can be integrated into RDI pipelines with a simple modification of existing KLIP infrastructure (e.g., pyKLIP, VIP, PynPoint) (Ren, 2023).

4. Forward Modeling and Unbiased Characterization

KLIP's linear framework permits analytic forward modeling for characterization of detected sources. After identifying a candidate in the residual $F(n)$, one propagates an astrophysical model $A_\xi(n)$ parameterized by $\xi$ (e.g., flux, position) through the truncated KLIP basis:

$$A_\xi^{\text{out}}(n) = A_\xi(n) - \sum_{j=1}^k e_j(n)\langle A_\xi, e_j \rangle$$

Source parameters $\xi$ are then estimated by minimizing

$$\chi^2(\xi) = \sum_{n\in \mathcal{S}} \left[ F(n) - A_\xi^{\text{out}}(n) \right]^2$$

Yielding unbiased photometry and astrometry limited by the residual noise in $F(n)$. This provides a direct analytic alternative to forward-injection and full pipeline reprocessing, especially critical for disks and low-SNR planetary signals (Soummer et al., 2012).

5. Algorithmic Performance, Limiting Factors, and Benchmarks

KLIP matches the speckle-suppression performance of the LOCI algorithm in detection, with several operational advantages:

• Parameter robustness: The optimality of the KL basis under truncation provides stability with respect to $k$ and segmentation choices.
• Computational efficiency: Once the basis is constructed, PSF subtraction reduces to inner products; for annular KLIP, region-specific $k$ can be selected. Downdate and batched eigendecomposition implementations further accelerate large-scale reductions (Long et al., 2021).
• Signal recovery: KLIP preserves sharp disk edges and point source throughput over a broad range of $k$; LOCI can induce significant self-subtraction of disks unless highly regularized (Soummer et al., 2012, Ren, 2023).
• GPU acceleration: torchKLIP demonstrates a 2–3× speedup over pyKLIP on full-frame reductions and scalability to $N\sim 10^2$–$10^3$ frames, with GPU acceleration providing up to an additional order of magnitude in speed (Ko et al., 2024).

Empirical benchmarks using VLT/NACO β Pictoris and Keck/NIRC2 HR 8799 data sets confirm that the SNR curves, residuals, and detection limits from torchKLIP, pyKLIP, and analytic KLIP are consistent to within a few percent. Common $k$ values in the literature are $k\sim 5$ for DIKL or $k=9$–20 for planetary datasets.

6. Practical Considerations and Pipeline Integration

Best practices for high-contrast imaging with KLIP include:

• Reference library selection: Exclude frames spatially or temporally proximate to the candidate planet to avoid self-subtraction.
• Masking and zone definition: Careful masking of stellar cores and bad pixels, and annular segmentation, enhances both computational tractability and SNR.
• K optimization: Scan the SNR over $k$ to determine optimal truncation.
• Residual estimation: For detection sensitivity, false positive rate, and contrast curve determination, use injection–recovery of synthetic planets.
• Integration: Existing KLIP implementations (pyKLIP, VIP, PynPoint, ADI.jl) facilitate rapid adoption of DIKL and downdate techniques; mask and region definition are typically the only required extension for DIKL (Ren, 2023).
• Machine learning integration: KLIP implementations in PyTorch (torchKLIP) allow natural extension with ML methods (e.g., kernel PCA, autoencoders, deep-learning-based PSF predictors).

7. Limitations and Extensions

KLIP assumes stationarity of the PSF statistics across the reference library and negligible change in the mean image when frames are excluded from the set. For high S/N sources or rapidly varying PSFs, residual biases (throughput loss, modeling error) can appear; analytically estimating or empirically calibrating the throughput remains necessary for precise contrast analyses. Extensions of KLIP include:

• Batch/online SVD updating: For streaming data and real-time applications such as AO control, downdate and update formulations enable continuous incorporation of new frames (Long et al., 2021).
• Multi-column downdate: Excluding multiple frames simultaneously from the reference library, at increased SVD rank.
• Kernel and non-linear generalizations: Replacement of linear PCA with kernel PCA or deep architectures.
• Hybrid iterative-analytic methods: Combining DIKL with advanced regularization strategies (e.g., NMF) offers theoretical paths to improve disk and extended structure recovery while preserving computational efficiency (Ren, 2023).
• Machine learning-based PSF models: Initial torchKLIP prototypes in PyTorch demonstrate the feasibility of incorporating neural architectures for speckle prediction and subtraction (Ko et al., 2024).

KLIP's principled linear algebraic formulation, direct throughput/variance analysis, and operational scalability have established it as a reference PSF subtraction method for data from the Hubble Space Telescope, ground-based extreme AO instruments, and recent applications to James Webb Space Telescope imaging.