PaperlayerView exact revision on arXiv

Source-linked AI summary

Tensor-on-tensor regression

Eric F. Lock

arXiv:1701.01037v2stat.MEstat.ML

TL;DR

Predicting one arbitrary-dimensional array from another lacked a general methodology. The paper introduces a contracted-product regression framework with reduced-CP-rank coefficients, and reports that simulations and a facial-image application illustrate the advantages of CP-rank and ridge regularization.

  • Problem

    General methodology was lacking for predicting an arbitrary-dimensional tensor from another, despite methods for several scalar-, vector-, and matrix-prediction settings.

    Source: §1 · p. 3

  • Method

    The framework uses a contracted tensor product with reduced-CP-rank coefficients and ridge-penalized least squares, with Gibbs sampling motivated by a Bayesian posterior mode.

    Source: §1 · p. 4, §3 · p. 6

  • Results

    Simulations and the facial-image application illustrate advantages of CP-rank and ridge regularization in the proposed framework.

    Source: §11 · p. 22

  • Takeaways & Limitations

    CP structure accounts for multi-way dependence in both arrays, while ridge regularization addresses auxiliary predictor high-dimensionality and multicollinearity.

    Source: §11 · p. 22

    The inference is semi-Bayesian: its prior is limited to facilitate inference under the penalized least-squares objective and is not intended to be subjective.

    Source: §11 · p. 23

Abstract

from arXiv · show

We propose a framework for the linear prediction of a multi-way array (i.e., a tensor) from another multi-way array of arbitrary dimension, using the contracted tensor product. This framework generalizes several existing approaches, including methods to predict a scalar outcome from a tensor, a matrix from a matrix, or a tensor from a scalar. We describe an approach that exploits the multiway structure of both the predictors and the outcomes by restricting the coefficients to have reduced CP-rank. We propose a general and efficient algorithm for penalized least-squares estimation, which allows for a ridge (L_2) penalty on the coefficients. The objective is shown to give the mode of a Bayesian posterior, which motivates a Gibbs sampling algorithm for inference. We illustrate the approach with an application to facial image data. An R package is available at https://github.com/lockEF/MultiwayRegression .

1 Introduction

The paper addresses the gap in predicting arbitrary-dimensional arrays from one another, using a contracted tensor product and reduced-CP-rank coefficients. It introduces the relevant array notation and defines the contraction underlying the framework.

1 Introduction

  • CP factorization represents a tensor as a linear combination of rank-1 tensors, helping capture multi-way dependence and reduce dimensionality.It is a special case of Tucker factorization with a diagonal core.Source: §1 · p. 2
  • Existing methods cover several scalar-, vector-, and matrix-prediction settings, but general arbitrary-dimensional tensor-on-tensor regression lacked methodology.The paper contrasts its framework with prior convex tensor-regularizer theory, which did not discuss estimation algorithms.Source: §1 · p. 3
  • The framework predicts an arbitrary-dimensional array from another through a contracted tensor product, with reduced-CP-rank coefficients that exploit structure across both arrays’ modes.The authors also propose ridge-penalized least-squares estimation and Gibbs sampling for inference.Source: §1 · p. 4

2 Notation and Preliminaries

  • The contracted tensor product contracts shared modes of two arrays and retains their uncontracted modes, extending matrix multiplication to higher-order arrays.The paper distinguishes arbitrary-dimensional arrays from vectors and matrices in its notation.Source: §2 · pp. 4–5

3 General framework

The framework predicts an outcome array of arbitrary dimension from a predictor array of arbitrary dimension by contracting the predictors with a coefficient array. Reduced CP-rank and an L2 penalty constrain estimation while retaining structure across both arrays.

3 General framework

  • The model contracts predictor array X with coefficient array B to produce outcome array Y across arbitrary dimensions.B’s first L modes contract predictor dimensions, and its last M modes expand along outcome dimensions.Source: §3 · p. 6
  • The matrix formulation supports all valid linear relations between the variates of X and Y.Source: §3 · p. 6

4 Estimation criteria

  • A rank-R CP representation reduces the coefficient model to R times the sum of predictor- and outcome-mode dimensions, borrowing information across both arrays.The reduced dimensionality can be several orders smaller than the unrestricted model.Source: §4 · p. 7
  • An L2 penalty addresses over-fitting and instability associated with high model dimension or predictor multicollinearity.The penalty parameter λ controls its degree; for a vector outcome and matrix predictor, the objective is ridge regression.Source: §4 · p. 7

5 Identifiability

The coefficient array is identifiable, but its CP factor components generally require constraints to resolve scaling and ordering ambiguities. The framework also recovers established regression methods in special cases and yields a connection between separable and non-separable penalties.

5 Identifiability

  • The coefficient array B is identifiable, but its factor components are not readily identified without restrictions.The component-identifiability conditions correspond to those for CP factorization.Source: §5 · p. 8
  • Scaling and ordering restrictions generally suffice when L+M ≥ 3; when L+M = 2, additional orthogonality restrictions are required.The restrictions can be imposed after estimation; for L+M = 2, singular value decomposition can identify components meeting them.Source: §5 · p. 8

6 Special cases

  • For vector outcomes from matrix predictors, the objective is standard ridge regression; for matrix-to-matrix prediction, it gives reduced rank regression at λ = 0 and reduced rank ridge regression otherwise.The stated equivalences connect the framework to ridge regression, reduced rank regression, and reduced rank ridge regression.Source: §6 · pp. 8–9
  • For scalar-on-tensor regression, the objective is equivalent to tensor ridge regression, and the proposed optimization gives a closed-form update for each U_l subproblem.The comparison is with Guo et al. (2012), whose described approach resorts to gradient-style methods for that subproblem.Source: §6 · p. 9
  • For L = 2, the separable L2 penalty without scale restrictions has the same solution as the non-separable nuclear-norm penalty.The nuclear norm is the sum of a matrix’s singular values.Source: §6 · p. 9

7 Optimization

The estimation algorithm alternates updates across the predictor and outcome modes to optimize the reduced-rank objective. Ridge regularization can be handled through modified least-squares updates, while tuning parameters can be selected by validation or DIC.

7 Optimization

  • The algorithm alternately updates each mode’s component vectors while holding the others fixed.Without ridge regularization, the component updates use separate OLS regressions.Source: §7 · p. 10, §7.1 · p. 10

7.1 Least-squares

  • For each predictor-mode update, the contracted predictor components form a design matrix for an OLS regression on the vectorized outcome.Outcome-mode components are updated analogously using separate OLS regressions.Source: §7.1 · p. 10

7.2 Ridge-regularization

  • Ridge-regularized optimization can use the least-squares algorithm with augmented predictor and outcome arrays, or equivalent efficient updates on the original data.Each sub-step improves the objective until convergence to a potentially local optimum; tempered regularization may improve robustness to initial values and local optima.Source: §7.2 · p. 11

7.3 Tuning parameter selection

  • Select λ and R using predictive accuracy on training and test sets or K-fold cross-validation; DIC offers a model-based alternative under the Bayesian framework.The cross-validation approach is described as free of model assumptions.Source: §7.3 · p. 12

8 Inference

The penalized least-squares objective is connected to a Bayesian posterior, motivating Gibbs sampling for inference and uncertainty quantification. The approach supports posterior inference on the coefficient array and predictions, while interpretation of individual CP factors requires care with identifiability.

8 Inference

  • The point estimator is consistent under a wide variety of distributional assumptions, including settings with correlated responses or predictors.Source: §8 · p. 12
  • The penalized objective gives the posterior mode under a rank-R Gaussian prior on B, motivating MCMC for uncertainty quantification.The unregularized objective is the maximum-likelihood estimate under the rank restriction; when λ = 0, the prior on B is flat over rank-R tensors.Source: §8 · pp. 12–13
  • Posterior draws support inference on B, conditional means, and predictive distributions for new response arrays without requiring identifiability of the CP factors.The paper also states that inference on the factors themselves is possible under the conditions of its identifiability section.Source: §8 · pp. 12, 14
  • Gibbs sampling draws from the full conditionals of the CP factors and error variance; the factor conditionals are multivariate normal.Source: §8 · p. 13
  • Interpreting individual factor components requires careful attention to their identifiability.A post-hoc transformation of components is one stated approach.Source: §8 · p. 14

9 Simulation study

The simulation study evaluates tensor-on-tensor prediction across sample sizes, signal levels, regularization choices, and assumed ranks. Regularization generally improves prediction in overfitting-prone settings, while excessive shrinkage and rank misspecification can hurt performance.

9.1 Approach

  • The study predicts a three-way array Y from another three-way array X under varied conditions, evaluating 240 simulated datasets across 25 procedures by relative out-of-sample prediction error on a new dataset of 500 observations.Source: §9.1 · pp. 15–16

9.2 Results

  • With no signal, the oracle RPE is 1; larger N and λ lower RPE, while higher assumed rank raises it, consistent with overfitting from weak regularization or excess rank.The marginal mean RPE is summarized in Table 1.Source: §9.2 · p. 16, Table 1 · p. 16
  • Higher λ generally improves prediction when sample size and signal-to-noise ratio are small, but excessive λ can over-shrink coefficients and introduce bias.With higher training sample size, credible-interval coverage is generally appropriate; at N = 30 and high signal, moderate λ can yield poor coverage because inference is biased toward smaller B.Source: §9.2 · pp. 16–17, Table 2 · p. 17
  • For N = 120, SNR = 1, and λ = 0, relative prediction error is minimized when assumed rank matches the true rank; overestimating rank is generally safer than underestimating it.The rank-misspecification comparison is presented in Table 3.Source: §9.2 · p. 17, Table 3 · p. 18

10 Application

The application predicts 72 facial attributes from color images using a training set of 1,000 images and a separate 1,000-image validation set. Rank restriction with ridge regularization achieves the lowest reported relative prediction error, while posterior predictions quantify uncertainty.

10 Application

  • The model predicts 72 standardized attributes from centered 90 × 90 × 3 facial images, using 1,000 training images and a separate 1,000-image validation set.Images are sampled from unique individuals.Source: §10 · p. 19
  • RPE reaches 0.568 at R = 15 and λ = 10^5; unregularized and full-rank models perform much worse.Figure 1 compares relative prediction error across parameter choices; its top row shows separate ridge regressions without rank restriction.Source: §10 · p. 20, Figure 1 · p. 20
  • Across 1,000 test images and 72 characteristics, actual and predicted values show substantial residual variation but a clear trend, with correlation r = 0.662.Figure 2 displays the actual-versus-predicted values.Source: §10 · p. 21, Figure 2 · p. 21
  • Posterior predictive intervals cover 0.934 of values at 95% and 0.887 at 90%; Figure 3 illustrates posterior samples for five selected characteristics of one test image.The predictive distributions use 5,000 posterior samples for each of 1,000 test images.Source: §10 · pp. 21–22, Figure 3 · p. 22

11 Discussion

The Discussion explains how CP-rank and ridge regularization address multiway dependence and predictor complexity, while identifying limits of the error model and inference approach. It also notes alternative penalties and provides an R package with documented code for the methods.

11 Discussion

  • CP-rank regularization accounts for multiway dependence in predictors and outcomes, while ridge regularization accounts for auxiliary high-dimensionality and predictor multicollinearity.The authors report that simulation studies and the facial image application illustrate the advantages of both forms of regularization.Source: §11 · p. 22
  • Separable penalties, including L1 for sparse subregions within modes, are straightforward to impose, but the described alternating scheme does not readily extend to other non-separable penalties.The Discussion identifies Tucker structure as an alternative to CP structure.Source: §11 · pp. 22–23
  • Independent, homoscedastic outcome errors may be inappropriate when the outcome has auxiliary structure; array-normal errors offer a more flexible covariance model.Tensor envelope methods are also identified as a possible direction for accounting for immaterial structure in responses or predictors.Source: §11 · p. 23
  • The inference is “semi-Bayesian”: its prior is restricted to facilitate inference under the penalized least-squares objective, not to express subjective beliefs.Fully Bayesian priors for coefficient rank and the shrinkage parameter are proposed as future work.Source: §11 · p. 23

SUPPLEMENTARY MATERIAL

  • The supplementary material provides the MultiwayRegression R package with documented code for all methods described in the article.The package is distributed as a GNU zipped tar file.Source: SUPPLEMENTARY MATERIAL · p. 23

A Consistency

Under the stated fixed-dimension conditions, the penalized least-squares estimator is consistent for the true coefficient array when the fitted rank is correct. The proof uses expected squared-error minimization and convergence arguments for the empirical loss.

A Consistency

  • With the true rank R = R0 and fixed ridge penalty λ ≥ 0, the minimizer converges in probability to B0 as N → ∞.Under the Section 5 restrictions, the factorization parameters also converge in probability to θ0.Source: A Consistency · p. 24
  • The result assumes iid observations, bounded predictors, and an identifiable true rank-R0 factorization in the interior of a compact parameter space.Within an observation, errors or predictor-array elements may be correlated or have different distributions.Source: A Consistency · p. 24
  • The proof identifies B(θ0) as minimizing expected squared error, then uses almost-sure empirical-loss convergence and a Glivenko-Cantelli argument for consistency.Compactness and continuity of the loss class support the Glivenko-Cantelli step.Source: A Consistency · p. 25
  • Consistency of the coefficient array follows from consistency of the factorization parameters by the continuous mapping theorem.The theorem treats fixed dimension as N → ∞.Source: A Consistency · pp. 23, 25

B Posterior derivations

This appendix derives full conditional distributions for the CP factorization components used in the paper’s inference procedure. Under independent Gaussian errors, the updates are obtained using Bayesian linear model results, with analogous derivations across modes.

B Posterior derivations

  • The appendix derives the full conditional distributions of the factorization components used in Section 8.The coefficient array is represented through factors U1,…,UL and V1,…,VM.Source: B Posterior derivations · p. 25
  • The factor prior conditionals follow from the spherical Gaussian prior on B and are permutation invariant across its L + M modes.The derivation for U1 uses a vectorized form of the CP factorization with that factor omitted.Source: B Posterior derivations · pp. 25–26
  • With independent N(0, σ²) errors, Bayesian linear model results yield Gaussian conditional distributions for factor updates, with analogous derivations for other modes.The appendix gives conditional distributions for U1 and VM and states that the remaining factor updates follow analogously.Source: B Posterior derivations · pp. 26–27

C Proof of Proposition 1

The appendix establishes an equivalence between separable L2 penalization and nuclear-norm penalization for predicting a vector from a three-way array. The equivalence uses orthogonality of the factor columns and balanced factor norms.

C Proof of Proposition 1

  • The proposition compares separable L2 and nuclear-norm penalization for vector prediction from a three-way array.The appendix notes analogous results for matrix-from-matrix and three-way-array-from-vector prediction.Source: C Proof of Proposition 1 · p. 27
  • For each paired factor contribution, equal factor norms minimize the sum of their squared norms subject to a fixed product.The proof invokes the scalar minimization arg min with ab = c.Source: C Proof of Proposition 1 · p. 27
  • Under orthogonal columns of U1 and U2, the factor-norm sum equals 2||B||*, yielding the stated nuclear-norm equivalence.The proof relates the nonzero singular values of B to the Frobenius norms of paired factor contributions.Source: C Proof of Proposition 1 · p. 28

D Correlated data simulation

The correlated-data simulations compare predictor and response correlation under spatially structured data. Predictor correlation sometimes improves penalized prediction error, while response correlation tends to reduce accuracy and credible-interval coverage.

D Correlated data simulation

  • The simulations used Gaussian data with exponential spatial correlation, with adjacent locations correlated at ρ = 0.6 in X or E.The response grid was 5 × 10 and the predictor grid was 15 × 20.Source: D Correlated data simulation · p. 28
  • Credible intervals tended to undercover more with correlated E than with uncorrelated data or correlated X.The authors suggest this may reflect correlation in E violating the assumed likelihood model for inference, unlike correlation in X.Source: D Correlated data simulation · p. 28
  • For penalized estimation at n = 30, correlated X had significantly better RPE than uncorrelated X, while correlated E often tended toward lower accuracy.Table 4 reports mean RPE and 95% credible-interval coverage rates for correlated X or E using correct assumed ranks.Source: D Correlated data simulation · p. 28, Table 4 · p. 29
Loading 1701.01037v2…