Bayesian Tensor Regression

Rajarshi Guhaniyogi; Shaan Qamar; David B. Dunson

Source-linked AI summary

Bayesian Tensor Regression

Rajarshi Guhaniyogi, Shaan Qamar, David B. Dunson

arXiv:1509.06490v1Sep 22, 2015stat.ME

TL;DR

Tensor regression commonly vectorizes multiway predictors, losing structure and producing difficult high-dimensional estimation. This paper introduces Bayesian multiway shrinkage priors for structured tensor coefficients, with posterior consistency and efficient MCMC computation. Simulations report improved coefficient estimation over FTR and vectorized approaches, while the framework supports sparsity and uncertainty quantification.

Problem
Vectorizing tensor covariates fails to exploit their structure, while independent voxel screening does not account for their joint impact.
Method
The paper uses a rank-R parafac coefficient representation with Bayesian multiway shrinkage priors for generalized linear tensor regression.
Results
M-DGDP consistently outperforms FTR in voxel-level RMSE, including lower RMSE on true zero and non-zero voxels in low-rank settings.
Takeaways & Limitations
The framework adapts model complexity and sparsity while providing posterior uncertainty for tensor coefficient and predictive inference.
Takeaways & Limitations
Rank-1 parafac modeling severely limits flexibility by ruling out interactions among dimensions.

Abstract

from arXiv · show

This article proposes a Bayesian approach to regression with a scalar response against vector and tensor covariates. Tensor covariates are commonly vectorized prior to analysis, failing to exploit the structure of the tensor, and resulting in poor estimation and predictive performance. We develop a novel class of multiway shrinkage priors for the coefficients in tensor regression models. Properties are described, including posterior consistency under mild conditions, and an efficient Markov chain Monte Carlo algorithm is developed for posterior computation. Simulation studies illustrate substantial gains over vectorizing or using existing tensor regression methods in terms of estimation and parameter inference. The approach is further illustrated in a neuroimaging application. Keywords: Dimension reduction; multiway shrinkage prior; magnetic resonance imaging (MRI); parafac decomposition; posterior consistency; tensor regression

1 Introduction

Tensor predictors are often analyzed by independent voxel screening or vectorization, but these approaches do not adequately exploit joint tensor structure. The paper proposes Bayesian multiway shrinkage to adapt model complexity, estimate important coefficients accurately, and quantify uncertainty.

Independent voxel screening overlooks the joint impact of the overall tensor, while simultaneous-analysis methods remain sparse in the literature.
Vectorizing tensors discards spatial structure and makes it harder to learn low-dimensional relationships when sample sizes are much smaller than voxel counts.
The proposed self-calibrating procedure shrinks unimportant voxel coefficients toward zero while preserving accuracy for important coefficients.
Multiway shrinkage priors induce sparsity within and across ranks, supporting model-based rank selection and optimal region selection.
A Bayesian formulation is motivated by the need for valid uncertainty measures in settings with low or moderate sample sizes.

2 Tensor regression

The tensor regression framework models scalar responses using vector and tensor predictors while representing the tensor coefficient through a structured parafac decomposition. This reduces dimensionality, but rank restrictions and tensor-specific parameter indeterminacies shape the model’s flexibility and interpretation.

2.1 Basic notation: The parafac decomposition is a simpler special case of Tucker decomposition, using equal rank across modes and diagonal core coefficients.
2.2 Model framework: The generalized linear model adds the tensor contribution ⟨X, B⟩ to scalar predictors, with ⟨X, B⟩ = vec(X)'vec(B).
2.2 Model framework: A rank-R parafac decomposition reduces the tensor coefficient’s parameter count from the unstructured tensor size to p + R∑_j p_j parameters.
2.2 Model framework: Rank-1 parafac modeling severely limits flexibility because a single margin vector per dimension cannot represent interactions among dimensions.
2.2 Model framework: Tensor margins are not individually identifiable because scale, permutation, and, for D = 2, orthogonal-transformation indeterminacies can leave B unchanged.
2.2 Model framework: The Bayesian method targets estimation and inference for B and predictive performance, so its goals do not require identifying the individual tensor margins.

3 Multiway shrinkage priors

The section develops multiway shrinkage priors for tensor regression that combine rank adaptation with sparsity across tensor components and margins. The M-DGDP construction is designed to preserve adequate tails while shrinking unimportant voxel coefficients.

3.2 Multiway priors: The proposed multiway shrinkage priors operate on rank-decomposed tensor margins, producing simultaneous shrinkage across the tensor coefficients.Under a rank-R parafac decomposition, voxel-level coefficients are nonlinear functions of margin parameters, and the prior shrinks across components and margins.
3.2 Multiway priors: The prior is designed to adapt toward lower-rank decompositions, favor contiguous predictive tensor regions, and preserve exchangeability when no elements are known to be important.These are stated as desirable properties for a multiway prior in the absence of prior information about important tensor elements.
3.3 The multiway Dirichlet GDP prior: The M-DGDP prior uses component-specific global scales from a Gamma–Dirichlet hierarchy and element-specific local scales with shared margin-level rates.The construction sets τ_r = φ_rτ, with τ Gamma-distributed and Φ Dirichlet-distributed, while local parameters model heterogeneity within margins.
3.4 Prior hyper-parameter elicitation: Theoretical variance bounds and prior summaries guide hyperparameter selection, with default choices illustrated through voxel-level percentiles and induced prior distributions.The paper reports lower and upper variance bounds and recommends aλ = 3 with bλ = 2D√aλ to avoid overly narrow induced variance.
3.4 Prior hyper-parameter elicitation: Dirichlet concentration controls effective rank: as α decreases, realizations concentrate near simplex vertices and become increasingly sparse.The paper illustrates this behavior with realizations from a Dirichlet distribution for R = 3.

4 Posterior consistency for tensor regression

The paper establishes posterior consistency for tensor regression with a multiway shrinkage prior in a growing-dimension setting. The main theorem gives sufficient prior conditions under which posterior mass concentrates around the true tensor parameter.

4.1 Notation and framework: The consistency argument assumes a centered response, known zero fixed effects, and data generated within the specified tensor regression model class.The paper states that these simplifying assumptions ease notation and calculations and that the results generalize straightforwardly.
4.1 Notation and framework: The asymptotic framework allows tensor dimensions to grow with sample size, even though the vectorized parameter dimension can substantially exceed n.This creates theoretical challenges related to, but distinct from, high-dimensional regression and multiway contingency tables.
4.2 Main result: The main theorem is obtained from a simple sufficient condition on the prior, supported by exponentially consistent tests.Lemma 7.1 establishes the tests used in the proof, while the theorem applies the prior condition to the proposed shrinkage prior.
4.2 Main result: The M-DGDP prior yields posterior consistency under the theorem’s stated regularity conditions.The result follows because the proposed prior satisfies the sufficient prior condition used for consistency.
4.2 Main result: The theorem requires bounded true margin coefficients and sub-linear growth of the product of tensor dimensions relative to sample size.The stated conditions also include growth and regularity constraints on the tensor dimensions and parameters.

5 Posterior computation and model fitting

Posterior computation uses a blocked MCMC algorithm tailored to the multiway shrinkage prior. Marginalization, compositional sampling, and back-fitting are used to improve mixing and update tensor parameters efficiently.

5 Posterior computation and model fitting: The implementation assigns conjugate inverse-gamma and normal priors to the noise variance and fixed effects after standardizing the response and voxel predictors.The response is centered and scaled, while tensor predictor voxels are standardized to mean zero and variance one.
5 Posterior computation and model fitting: The M-DGDP prior enables efficient posterior computation through marginalization and blocking.The sampler cycles over prior hyperparameters, tensor coefficients and latent variables, and regression parameters.
5 Posterior computation and model fitting: The sampler updates [α, Φ, τ|B, W], [B, W|Φ, τ, γ, σ, y], and [γ, σ|B, y] in sequence.This block structure separates rank-specific shrinkage parameters, tensor coefficients, and regression quantities.
5 Posterior computation and model fitting: The non-trivial rank-scale block is sampled compositionally, which the authors identify as essential for good mixing under the M-DGDP prior.The Dirichlet concentration parameter is updated by griddy-Gibbs, followed by rank-specific scale updates.
5 Posterior computation and model fitting: Tensor coefficients are updated using full conditional distributions with a back-fitting procedure across margin-level conditional distributions.The procedure cycles over margin-specific coefficient and latent-variable updates across components.

6 Simulation studies

Simulations compare M-DGDP with frequentist tensor regression and vectorized Lasso across tensor dimensions, ranks, sparsity patterns, and image-like settings. M-DGDP consistently improves voxel-level estimation, provides good interval coverage, and remains comparatively robust as predictor dimensions increase.

6 Simulation studies: The simulations vary tensor dimension, parafac rank, sparsity, and signal complexity using generated tensors and three ready-made two-dimensional images.Five replicated datasets with n = 1000 are generated for the simulated settings.
6 Simulation studies: M-DGDP consistently outperforms FTR on voxel-level RMSE, with lower error for both true zero and true non-zero voxels in low-rank settings.The comparison uses generated and ready-made tensor examples and evaluates estimation across zero, non-zero, and overall voxels.
6 Simulation studies: M-DGDP adapts to varying sparsity by shrinking many coefficients near zero while accurately estimating non-zero voxels.The paper contrasts this behavior with FTR’s sensitivity to cross-validation tuning and tensor dimension.
6 Simulation studies: M-DGDP yields credible intervals with good frequentist coverage across simulated settings, both overall and for true non-zero coefficients.The study evaluates coverage and interval length for 95% credible intervals.
6 Simulation studies: As margin dimensions increase, FTR’s RMSE worsens considerably more on true zero coefficients, while M-DGDP remains the clear winner on the absolute scale for true non-zero voxels.RMSE for true non-zero voxels increases for both methods, but their relative deterioration differs substantially.

7 Simulated response with a real 3D brain image

The simulations evaluate M-DGDP on sparse 3D tensor coefficients and compare its estimation and uncertainty performance with FTR and vectorized Lasso. Across the examined cases, M-DGDP improves coefficient estimation and provides well-calibrated credible intervals.

MRI application: The 550-adolescent MRI application treats age and sex as scalar covariates and 3D MRI images as tensor covariates, with a simulated rank-2 tensor coefficient.The images have dimensions 30×30×30.
3D simulation results: 10–15% better performance is reported for M-DGDP across the three sparse 3D simulation cases.The cases contain 12%, 18%, and 30% nonzero elements, respectively.
3D simulation results: M-DGDP outperforms L1-optimized methods more strongly in sparser cases, while auto-tuning avoids sensitivity to penalty selection.L1 penalties can substantially over-shrink coefficients that differ significantly from zero.
Uncertainty quantification: M-DGDP consistently achieves coverage over 95% with reasonably short credible intervals, whereas L1-based methods generally perform worse.Coverage and interval length are summarized for 95% credible intervals.
MRI application: The analysis reports point estimates and credible intervals for age and sex and concludes that M-DGDP provides superior performance with uncertainty characterization.These scalar-covariate estimates are reported in Table 7.

MCMC algorithm

The MCMC derivations specify conditional sampling steps for the M-DGDP prior by exploiting gamma, Dirichlet, and generalized inverse Gaussian distributions. The normalized component weights are obtained by independently sampling positive variables and renormalizing them.

MCMC algorithm: The M-DGDP sampling derivation uses gamma-distributed global scale and Dirichlet-distributed component weights.The prior specifies τ ∼ Ga(aτ, bτ) and Φ ∼ Dirichlet(α1, …, αR).
MCMC algorithm: Conditional draws for each component weight use independent generalized inverse Gaussian variables with parameters determined by αr, bτ, and Cr.The sampled variables have distribution giG(αr − p0/2, 2bτ, 2Cr).
MCMC algorithm: The component weights are obtained by normalizing the independently sampled positive variables.This produces weights supported on the simplex.

Proof of lemma 3.1

The proof bounds the variance induced on voxel-level coefficients by the M-DGDP prior. It combines moment calculations for local scales, gamma-function bounds, and norm inequalities.

Variance bound: The induced voxel-level variance is bounded using the inverse moments of the local shrinkage parameter and bounds on the Dirichlet weights.The proof assumes aλ > 2 and uses ||Φ||_D^D ≤ 1 together with Hölder-type inequalities.
Variance bound: Gamma-function ratios are bounded exponentially in the Dirichlet concentration and tensor dimension to control the variance expression.The bound uses α1 = c/R and log(x+1) ≤ x.

Consistency proofs

The consistency proofs establish the testing and prior-mass ingredients needed for posterior consistency under growth and regularity conditions. Supplementary lemmas provide the probability, norm, convexity, and density results used in those arguments.

Testing argument: Posterior consistency relies in part on exponentially consistent tests for distinguishing B_n = B_0 from alternatives.The supplement defines the test sequence and states its existence under a dimension-growth condition.
Testing argument: The testing construction uses concentration bounds for Gaussian and chi-square quantities, including Laurent–Massart inequalities.These bounds control probabilities of the test events and associated remainder terms.
Prior-mass argument: The proof controls prior mass by bounding neighborhoods of the tensor factor loadings and deriving lower bounds for their probability.The argument introduces κ_n and bounds factor norms and associated prior probabilities.
Supporting lemmas: The supplementary results establish algebraic identities for rank-one tensors and norm bounds used in the consistency proof.These include a telescoping decomposition of T − F and bounds involving factor norms.