Foreground Subtraction with a Tensor-Based Oriented Singular Value Decomposition Method for HI Experiments

This paper introduces a native tensor-based framework for foreground mitigation in 21 cm intensity mapping (IM), utilizing the Oriented Singular Value Decomposition (O-SVD) algorithm. The authors argue that traditional mitigation strategies often necessitate flattening multidimensional data cubes into two-dimensional matrices, a process that potentially compromises the intrinsic spatial-spectral correlations by treating distinct spatial pixels as independent samples. By treating multi-frequency sky maps and angular power spectra as third-order tensors, the O-SVD method performs decomposition directly on the multilinear manifold, preserving the underlying physical topology and leveraging the distinct coherence properties of astrophysical foregrounds across different dimensions.

The paper demonstrates the performance and versatility of the O-SVD framework through its application to high-fidelity simulations from the SKA Science Data Challenge 3a (SDC3a) and real-world observational data from the Tianlai Cylinder Pathfinder Array. The results indicate that O-SVD provides a robust and universal approach for foreground subtraction, achieving high-fidelity signal recovery while offering superior performance compared to conventional matrix-based Singular Value Decomposition (SVD) methods.

The methodology section details the mathematical foundation of O-SVD. For a third-order tensor A ∈ CI1 ×I2 ×I3 with 3-rank R3 = rank3 (A), the O-SVD is defined as A = (U ∗3 S ∗3 V) ×3 U (3), where U (3) ∈ CI3 ×I3 is a unitary matrix obtained from the SVD of the mode-3 unfolding A(3). The tensors U ∈ CI1 ×I1 ×I3, S ∈ CI1 ×I2 ×I3, and V ∈ CI2 ×I2 ×I3 satisfy the following properties: each frontal slice U(:,:, k) and V(:,:, k) is a unitary matrix, and S(:,:, k) is a non-negative diagonal matrix for k = 1, 2,..., R3; all frontal slices are null matrices for k = R3 + 1,..., I3. The diagonal entries sjjk of the core tensor S represent the O-SVD singular values, which obey a hierarchical ordering property: ∥S(:,:, 1)∥F ≥ ∥S(:,:, 2)∥F ≥ · · · ≥ ∥S(:,:, I3)∥F ≥ 0, where the Frobenius norm of each slice corresponds to the k-th singular value of the mode-3 unfolding, i.e., ∥S(:,:, k)∥F = σk.

The computational procedure involves two primary stages. First, a standard matrix SVD is performed on the I3 × I1 I2 mode-3 unfolding matrix A(3): A(3) = U (3) Σ(3) V (3)H, where U (3) is the unitary matrix of spectral basis vectors, and V (3) contains the corresponding flattened spatial modes. In the second stage, each spatial mode vk (the k-th column of V (3)) is reshaped into an I1 × I2 matrix Ṽk and decomposed via SVD: Ṽk = Uk Σk VkH, where Σk contains the spatial singular values σkj. The relationship between the O-SVD singular values and the matrix singular values is given by sjjk = σk σkj, which implies σk2 = Pr2 j=1 s2jjk. This hierarchical structure allows A to be expressed as a superposition of rank-1 outer products: A = Pr1 k=1 Pr2 j=1 sjjk ×1 ukj ×2 vkj ×3 uk, where uk represents the spectral basis vector, while ukj and vkj represent the spatial basis vectors.

The paper highlights the intrinsic connection between O-SVD and standard PCA, noting that PCA is mathematically equivalent to performing an SVD on the mode-3 unfolded matrix A(3). While PCA treats the associated spatial modes as flattened 1D vectors, O-SVD further decomposes these modes into their constituent spatial basis vectors. This positions O-SVD as a direct plug-and-play generalization of PCA, allowing any foreground mitigation pipeline currently employing PCA to be upgraded to O-SVD by replacing the matrix SVD step.

For the SKA SDC3a application, the analysis uses the first sub-band (106–121 MHz) of the dataset, which spans a frequency range of 106–196 MHz with 900 spectral channels. The processed image cube was formalized as a third-order tensor T ∈ RNx ×Ny ×Nν, with Nx = Ny = 900 representing spatial dimensions and Nν = 150 representing the spectral channels. The O-SVD method was applied to decompose T into its hierarchical multilinear components. A truncation threshold Nf g was determined by identifying the transition point in the singular value spectrum where the steeply declining foreground components reach the noise-dominated floor. The authors selected s11,k=100 as the truncation threshold, removing all modes satisfying sjjk ≥ s11,k=100, which yielded Nf g = 17, 652 removed modes.

The results show that the O-SVD method provides excellent agreement with the ground-truth EoR power spectrum across most k-scales, while traditional PCA methods show either high foreground residuals (20 and 30 modes) or significant signal loss (50 modes). The paper notes that Nf g = 17, 652 O-SVD modes corresponds approximately to the variance of NfPCA ≈ 20 matrix modes. However, the matrix SVD effectively subtracts all O-SVD modes with k ≤ 20, including many small variance modes that do not contribute significantly to the foreground power. By contrast, the O-SVD provides a more refined decomposition, allowing for the targeted removal of specific spatial-spectral modes, enabling O-SVD to isolate foregrounds more effectively while minimizing signal loss in the cosmological signal.

For the Tianlai application, the analysis is conducted in the angular power spectrum domain using Multi-Frequency Angular Power Spectra (MAPS), denoted as Cl (ν, ν ′). The authors constructed a third-order tensor C ∈ RNν ×Nν ×Nl, where the dimensions correspond to frequency ν, frequency ν ′, and multipole l, respectively. The analysis focuses on the YY polarization channel, which exhibited superior spectral smoothness and lower systematic contamination compared to the XX channel. The O-SVD singular value spectrum provided a more structured representation, and the authors identified k = 400 as the transition point where the singular values begin to plateau into a noise-like floor. They defined s11,k=400 as the truncation threshold and excised all O-SVD modes with corresponding singular values sjjk ≥ s11,k=400, removing 13,283 dominant modes.

The residual frontal slice Cl=100 res (ν, ν ′) shows a marked concentration of power along the diagonal, with significantly attenuated off-diagonal features, consistent with the expected behavior of the 21 cm signal and thermal noise. The paper notes that while ideal noise is frequency-independent and strictly diagonal, instrumental systematics and pipeline artifacts may introduce residual correlations. The frequency correlation length ∆ν of the 21 cm signal is intrinsically linked to the angular scale; at l ∼ 100, the signal typically decoheres beyond ∆ν ∼ 1 MHz, whereas at l ∼ 103, this occurs at ∼ 0.1 MHz.

The discussion section addresses the computational cost of O-SVD compared to traditional SVD. For a data cube of size Nx × Ny × Nν, traditional SVD has complexity O(Nν2 Nx Ny + Nν3), while O-SVD has total complexity O(Nν2 Nx Ny + Nν3 + Nν · (Nx2 Ny + Nx3)). For the SDC3a dataset (900 × 900 × 150), the additional computational cost of a full O-SVD decomposition is approximately 13×, though this can be reduced to 6–7× when only the largest singular values need to be computed in the second stage. The paper notes that recent advancements in randomized algorithms for oriented tensors can substantially reduce this overhead while maintaining rigorous error bounds.

The paper concludes that O-SVD offers a unified and flexible architecture capable of processing diverse data representations, ranging from spatial-frequency image cubes to multi-frequency angular power spectra. Future research will explore the integration of automated mode truncation strategies based on objective information-theoretic criteria—such as the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), or Generalized Cross-Validation (GCV)—to further minimize signal bias and enhance the reproducibility of foreground subtraction. The authors anticipate that this work will serve as a catalyst for the broader adoption of advanced multilinear algebra and tensor-based methodologies within the astronomical community, providing powerful new tools for the era of precision 21 cm cosmology.

Improvements to AI Systems:

1. Tensor-native dimensionality reduction for multi-dimensional data: Replace matrix-based PCA/SVD with O-SVD in any AI pipeline that processes multi-dimensional data (e.g., hyperspectral imaging, video, fMRI, climate data). The improved AI system can decompose data directly on its multilinear manifold, preserving spatial-spectral correlations that flattening destroys, leading to higher-fidelity signal recovery and reduced information loss.

2. Hierarchical mode selection for noise-robust feature extraction: Implement the O-SVD’s hierarchical singular value ordering (∥S(:,:,1)∥F ≥ ∥S(:,:,2)∥F ≥ …) to automatically identify noise floors and truncation thresholds. The improved AI system can distinguish dominant structured components from noise-dominated modes with greater precision than flat SVD, enabling adaptive, data-driven denoising without over-subtraction or signal loss.

3. Plug-and-play upgrade for existing PCA-based AI systems: Since O-SVD is a direct generalization of PCA (with matrix SVD as a special case), any AI model currently using PCA for feature extraction, anomaly detection, or compression can be upgraded by swapping the matrix SVD step with O-SVD. The improved AI system gains the ability to model cross-dimensional interactions (e.g., frequency-frequency, spatial-spatial) that were previously ignored, improving performance on tasks like foreground removal, image reconstruction, and tensor completion.

4. Multi-representation data fusion: O-SVD’s flexibility to operate on both image cubes (Nx × Ny × Nν) and cross-power spectra (Nν × Nν × Nl) enables AI systems to process the same physical phenomenon across different tensor representations. The improved AI system can unify analysis of spatial-frequency data and angular power spectra within a single framework, allowing cross-validation and joint inference across domains.

5. Computationally efficient tensor decomposition with error bounds: By leveraging randomized algorithms for oriented tensors (as noted in the paper), AI systems can achieve O-SVD decomposition with 6–7× reduced overhead while maintaining rigorous error bounds. The improved AI system can handle large-scale, high-dimensional datasets in real-time or near-real-time, making tensor-based methods feasible for streaming or interactive applications.

6. Objective, information-theoretic mode truncation: Integrate AIC, BIC, or GCV criteria into the O-SVD framework to automate the selection of truncation thresholds. The improved AI system can self-tune the number of retained modes based on the data itself, eliminating manual threshold selection and reducing human bias, thereby enhancing reproducibility and generalization across different datasets.

7. Signal-preserving foreground/artifact separation: In applications like 21 cm cosmology, the improved AI system can use O-SVD to isolate smooth, correlated foregrounds from faint, uncorrelated signals (e.g., EoR) by exploiting their distinct coherence properties across dimensions. This enables the AI to recover weak signals that traditional matrix methods either over-subtract or leave contaminated, improving sensitivity in low-signal-to-noise regimes.

8. Cross-frequency correlation modeling for time-series or spectral data: The O-SVD’s ability to decompose along frequency and spatial modes simultaneously allows AI systems to model long-range correlations in spectral dimensions (e.g., in radio astronomy, spectroscopy, or remote sensing). The improved AI system can predict missing frequency channels or denoise spectral data by leveraging the hierarchical structure of the tensor, rather than treating each channel independently.

9. Scalable multilinear anomaly detection: The hierarchical O-SVD singular values provide a natural ranking of tensor components by variance. The improved AI system can flag anomalies as modes with unexpectedly high singular values in specific spatial-spectral slices, enabling detection of transient or localized events (e.g., radio frequency interference) with higher specificity than global PCA-based methods.

10. Improved power spectrum estimation and uncertainty quantification: By preserving the tensor structure during decomposition, the improved AI system can produce more accurate estimates of angular power spectra (e.g., Cl(ν, ν′)) with reduced leakage between modes. This leads to better characterization of signal statistics and more reliable error bars in scientific inference tasks, such as cosmological parameter estimation.