Heteroskedastic Canonical Polyadic Tensor Decomposition

Kyle Ritscher; Carlos Llosa-Vite

stat.MESep 30, 2026

Heteroskedastic Canonical Polyadic Tensor Decomposition

Authors: Kyle Ritscher, Carlos Llosa-Vite

Organizations: Computer Science Research Institute, Sandia National Laboratories, Albuquerque, NM · Department of Statistics & Data Science, University of California, Los Angeles, CA

Abstract

When minimizing the squared-error loss, the popular CP decomposition can be interpreted as parameter inference in a Gaussian model with a low-rank mean tensor and constant variance across the tensor entries. We introduce heteroskedastic-CP (HCP), which models entrywise variability with a non-constant, low-rank precision tensor, and develop an alternating block-coordinate ascent method to recover both the low-rank mean and precision tensors from noisy observations. Our procedure is computationally competitive, with the same leading-order factor-update complexity as CP-ALS. We demonstrate HCP on synthetic experiments and an EEG application.

Figures & tables

Figure 1

Figure 2 : Comparison of relative errors in estimates of $\mathcal{B}$ at varying heteroskedasticity levels.

Figure 3 : Visualization of $9$ th slice of the noisy observed tensor $\mathcal{X}$ and the CP, PGCP, and HCP estimates $\widehat{\mathcal{B}}$ , for varying heteroskedasticity levels.

Figure 4 : Panels correspond to a different true observation model, and boxplots show the $\widehat{\mathcal{B}}$ recovery in relative Frobenius error over replicates for HCP, CP, and GCP methods using Poisson, Gamma, and inverse Gaussian losses.

Panel of Figure 5 — Figure 5 : Simulation results comparing recovery performance across tensor sizes. The left panel shows relative error in recovery of $\mathcal{B}$ , while the right panel shows relative error in recovery of $\mathcal{S}$ .

Figure 6 : Estimated log-variance from HCP for the EEG tensor, averaged over channels and runs. The recovered variance is largest at lower frequencies.

Figure 7 : Visualization of the factor matrices recovered from the EEG tensor by CP and HCP with rows corresponding to components. Within each method, columns show the channel/scalp loadings, time-frequency loadings, and measurement loadings. Dashed vertical lines in measurement columns indicate recording-run boundaries.

Appendix figures & tables9 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 8 : Results for a toy 3 mode tensor ( $10\times 8\times 12$ ) with $R=3$ and known $\mathcal{S}$ . The first row displays the first horizontal slice, and the second, the first lateral slice, for the observed data, the true mean, and the estimated mean.

Figure 9 : Results for a 3 mode tensor ( $10\times 8\times 12$ ) with $R=3$ and known $\mathcal{B}=0$ . The first row displays the first horizontal slice, and the second, the first lateral slice, for the squared residuals, the true variance, and the estimated variance.

Figure 10 : Results for a 3 mode tensor ( $10\times 8\times 12$ ) with $R=3$ . For the first frontal slice we show the data, true mean tensor, and estimated mean tensor in the first row, and true vs learned variance tensors in the second row.

Panel of Figure 11 — Figure 11 : Simulation results comparing recovery performance across tensor sizes and heteroskedasticity structures. The left panel shows relative error in recovery of $\mathcal{B}$ , while the right panel shows relative error in recovery of $\mathcal{S}$ .

Figure 12 : Mean recovery when varying the fitted mean rank $R_{\mathcal{B}}$ , with the precision rank fixed at its true value $R_{\mathcal{S}}=3$ . The true mean rank is $R_{\mathcal{B}}=3$ . HCP consistently outperforms CP decomposition, with the worst performance for both methods occuring when the fitted mean rank is too small.

Figure 13 : Precision recovery when varying the fitted mean rank $R_{\mathcal{B}}$ , with the precision rank fixed at its true value $R_{\mathcal{S}}=3$ . This plot shows how misspecification of the mean rank affects estimation of the precision tensor.

Figure 14 : Mean recovery when varying the fitted precision rank $R_{\mathcal{S}}$ , with the mean rank fixed at its true value $R_{\mathcal{B}}=3$ . Importantly, even for the worst fitted precision ranks, HCP will outperform CP in mean recovery.

Figure 15 : Precision recovery when varying the fitted precision rank $R_{\mathcal{S}}$ , with the mean rank fixed at its true value $R_{\mathcal{B}}=3$ . The true precision rank is $R_{\mathcal{S}}=3$ . Unlike for the mean, if the tensor is small, the fitted rank being too large can be more detrimental than being too small.

Figure 16 : Visualization of the factor matrices recovered from the EEG tensor by HCP and CP with rows corresponding to components. Within each method, columns show the channel/scalp loadings, time-frequency loadings, and measurement loadings. Dashed vertical lines in measurement columns indicate recording-run boundaries.

Explore similar work

Bayesian Tensor Decomposition with Diffusion Model Prior
Jun 2, 2026Zerui Tao, Qibin ZhaoTensor Decomposition Diffusion Priors
RIKEN Center for Advanced Intelligence Project (AIP), Tokyo, Japan.
Spectra-Guided Neural Tucker Factorization
May 30, 2026Fusheng Wang, Yikai HouTensor Completion
School of Automation, Chongqing University of Posts and Telecommunications, Chongqing, China · College of Computer and Information Science, School of Software, Southwest University, Chongqing, China
Robust Neural Tucker Factorization with Bias Correction and Adaptive Initialization
Jun 15, 2026Yuchao Su, Yixin RanTensor Completion Neural Tangent Kernel
School of Computer Science and Engineering, Chongqing University of Science and Technology, Chongqing, China · College of Computer and Information Science, School of Software, Southwest University, Chongqing, China

Jun 2, 2026, cs.LGJ/K move · Enter open · S save

Bayesian Tensor Decomposition with Diffusion Model Prior

Zerui Tao, Qibin Zhao

RIKEN Center for Advanced Intelligence Project (AIP), Tokyo, Japan.

Tensor Decomposition Diffusion Priors Low-Rank Structure Bayesian Posterior Decomposition+1

Open paper AskPDF arXiv

Low-rank tensor decomposition (TD) is usually effective on clean, fully observed data, but it often degrades under severe missingness or noise. Low-rankness is itself a useful but limited structural prior, and additional handcrafted priors (e.g., sparsity or smoothness) still fall short of capturing the rich statistics of real-world data. To compensate for this weak inductive bias under heavy corruption, one would like to inject a learned, data-driven prior; however, the state-of-the-art diffusion models are not readily compatible with current TD and tractable posterior inference. To address these challenges, we introduce DiffBCP, a hybrid-prior Bayesian CP decomposition framework that couples a cumulative shrinkage process prior over the CP factors for automatic rank selection with an off-the-shelf pre-trained diffusion model as an implicit data prior on the reconstructed tensor. To make posterior inference tractable despite the coupling among the likelihood, low-rank constraint, and diffusion prior, we develop a split Gibbs sampler: CP factors admit conjugate updates, while the diffusion block is sampled via low-rank-guided denoising. A noise-adaptive coupling schedule further reduces sensitivity to hand-tuned annealing. Experiments on image inpainting and denoising, including high-resolution out-of-distribution images, show consistent gains over Bayesian, nonlinear, and plug-and-play TD baselines.

Abstract

Explore similar work

Bayesian Tensor Decomposition with Diffusion Model Prior

Spectra-Guided Neural Tucker Factorization

Robust Neural Tucker Factorization with Bias Correction and Adaptive Initialization