Kernelized Tensor Factorization

• Kernelized tensor factorization is a framework that extends classic tensor methods by incorporating kernel mappings to introduce nonlinearity and uncertainty quantification in multiway data analysis.
• It leverages distinct formulations—such as tensor-train kernels, RKHS latent mappings, and Bayesian GP priors—to retain mode identity while enhancing scalability and expressive power.
• These techniques enable efficient modeling of high-dimensional data with applications in classification, regression, optimization, and tensor completion through complementary global and local components.

Kernelized tensor factorization denotes a class of models that combine tensor factorization with kernel methods or Gaussian-process priors so that multilinear structure is retained while nonlinear dependence, side information, smoothness, and, in many cases, uncertainty quantification are introduced. The literature suggests that the term covers several distinct but related constructions: tensor-train kernels defined directly on tensor inputs for nonlinear classification, RKHS mappings from mode-wise covariates to latent factors, Bayesian CP decompositions with GP-valued latent basis functions, GP models on concatenated latent factors for arbitrary tensor entries, and additive global-plus-local models in which a kernelized low-rank tensor is complemented by a correlated residual process (Chen et al., 2020, Hu et al., 2020, Lei et al., 2021, Lei et al., 2022, Lei et al., 2023, Zhe et al., 2016, Lei et al., 2024).

1. Position within tensor methods and kernel methods

Kernelized tensor factorization emerged from a precise limitation of two older paradigms. Standard tensor factorization methods such as CP, Tucker, and TT compress or parameterize tensors, but they often remain in a linear or multilinear learning framework. Standard kernel learning, by contrast, supplies nonlinear decision functions or regression surfaces, but typically acts on vectorized inputs and therefore discards mode identity and multiway structure. In the tensor-train classification setting, the central distinction is explicit: standard tensor factorization stays linear, standard kernel learning ignores tensor structure by flattening, and the kernelized TT construction combines both by preserving structure and reducing dimensionality while making the classifier nonlinear (Chen et al., 2020).

A related motivation appears in side-information-aware tensor regression. When mode-wise latent factors are regressed directly on raw covariates through linear maps, the model inherits a linear bottleneck in feature space and a parameter explosion when side information is high-dimensional. The RKHS reformulation replaces direct linear regression on covariates with kernel expansions, so that each latent factor coordinate becomes a smooth function of side information rather than a raw linear coefficient vector (Hu et al., 2020).

The contrast with scalable but non-kernelized tensor factorization is also instructive. DFacTo accelerates CP-ALS and CP-GD by reformulating the expensive matricized-tensor times Khatri-Rao product so that it can be computed with two sparse matrix-vector products per rank component, thereby addressing the “intermediate data explosion” problem without introducing kernelization (Choi et al., 2014). This suggests that kernelized tensor factorization should be understood not merely as a scalability technique, but as a structural extension that adds nonlinear or geometry-aware inductive bias to tensor models.

2. Principal model families

The major formulations can be organized by where the kernel enters the tensor model.

Formulation	Core construction	Representative papers
TT-kernel on tensor input	Kernel on TT cores or fibers inside a tensor-space classifier	(Chen et al., 2020)
RKHS latent-factor map	$\mathbf Z^{(d)}=\mathbf K^{(d)}\mathbf A^{(d)}$ within a CP-style predictor	(Hu et al., 2020)
GP-prior low-rank tensor model	GP priors on latent basis functions or factor columns in CP decompositions	(Lei et al., 2021, Lei et al., 2023, Lei et al., 2022)
Entrywise nonlinear GP factorization	GP on concatenated latent factors for selected tensor entries	(Zhe et al., 2016)
Additive global-local decomposition	Kernelized low-rank tensor plus correlated residual tensor	(Lei et al., 2022, Lei et al., 2024)

In the tensor-train strand, a $d$-way tensor $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$ is written in TT form through cores $A^{(1)},\dots,A^{(d)}$, with

$$A(i_1,\ldots,i_d)=A^{(1)}(:,i_1,:)\cdots A^{(d)}(:,i_d,:),$$

and storage reduced to

$\sum_{k=1}^d R_k I_k R_{k+1}.$

Kernelization is then applied not to a vectorized tensor but to fibers inside TT cores, so that the mapped object remains TT-structured in feature space with the same TT-ranks (Chen et al., 2020).

In RKHS-based tensor regression, the latent factors for each mode are not free parameters alone. Instead, for each mode $d$, one writes

$\mathbf Z^{(d)}=\mathbf K^{(d)}\mathbf A^{(d)},$

and the tensor predictor keeps a CP-like multilinear interaction,

$$f_{\mathbf i}=\sum_{k=1}^{K}\prod_{d=1}^{D} z^{(d)}_{i_d k}.$$

Here the kernel matrix $\mathbf K^{(d)}$ is built from side information, so mode-specific metadata enters the factorization through RKHS structure rather than direct linear regression (Hu et al., 2020).

A third family is fully Bayesian. In Bayesian Kernelized Tensor Factorization for Bayesian optimization, the unknown function over a discretized Cartesian product space is approximated by a rank-$d$0 CP decomposition with GP-valued latent factors,

$d$1

where each $d$2 is drawn from a GP prior (Lei et al., 2023). In Bayesian Kernelized Tensor Regression for spatiotemporally varying coefficients, the coefficient tensor is written as

$d$3

and GP priors are placed on the spatial and temporal factor columns (Lei et al., 2021).

The flexible nonlinear GP formulation removes the traditional tensor-GP Kronecker restriction altogether. For an entry indexed by $d$4, the model forms the concatenated latent input

$d$5

and assumes

$d$6

with a GP prior directly on $d$7. Because the covariance is defined entry-wise over concatenated latent-factor inputs rather than as $d$8, arbitrary subsets of entries can be selected for training (Zhe et al., 2016).

Finally, additive formulations separate scales explicitly. In BCKL,

$d$9

where $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$0 is a kernelized low-rank CP component and $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$1 is a short-range spatiotemporal GP residual (Lei et al., 2022). GLSKF adopts the same global-versus-local separation for tensor completion, with a smoothness-constrained CP tensor $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$2 and a correlated residual tensor $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$3 (Lei et al., 2024).

3. Kernelization mechanisms and induced dependence

The defining technical issue is how a valid kernel or covariance structure is made compatible with tensor factorization. In the TT classifier, feature maps are applied to fibers $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$4 inside each TT core rather than to a flattened tensor. This yields two kernel constructions. K-STTM-Prod mirrors TT contraction by a product of mode-wise kernel evaluations, while K-STTM-Sum replaces the product with a sum to avoid distortion of similarity information. Positive semidefiniteness is preserved because sums of PSD matrices remain PSD, and Hadamard products of PSD matrices remain PSD by the Schur product theorem. The same framework also allows different kernels $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$5 on different modes, for example RBF on spatial modes and linear or polynomial on the color mode of an image (Chen et al., 2020).

In the RKHS tensor-regression formulation, each latent coordinate is itself a kernel expansion. Writing

$A\in\mathbb{R}^{I_1\times\cdots\times I_d}$6

and applying the representer theorem makes the side-information map nonlinear while keeping the final tensor interaction multilinear in latent coordinates. The Bayesian interpretation then comes from Gaussian priors on the dual coefficients and a Gaussian likelihood, which turn the kernelized factorization into a probabilistic tensor regression model with predictive variance (Hu et al., 2020).

In GP-prior tensor models, kernelization acts through latent functions or latent factor columns. BKTF places GP priors on $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$7, with nearby locations in dimension $A\in\mathbb{R}^{I_1\times\cdots\times I_d}$8 therefore having correlated latent values. After marginalizing an auxiliary Gaussian latent variable, the model induces the multilinear kernel

$A\in\mathbb{R}^{I_1\times\cdots\times I_d}$9

The resulting kernel is nonstationary, because the latent functions vary with location, and nonseparable when $A^{(1)},\dots,A^{(d)}$0, because the sum over multiplicative components cannot be written as a simple product across dimensions (Lei et al., 2023).

BKTR uses GP priors on spatial and temporal factor matrices rather than on a full coefficient tensor, thereby reintroducing local dependence after low-rank compression. The cited kernels include a Matérn $A^{(1)},\dots,A^{(d)}$1 spatial kernel, a squared exponential temporal kernel, and, in the bike-sharing application, a locally periodic temporal kernel with period $A^{(1)},\dots,A^{(d)}$2 days (Lei et al., 2021). BCKL follows the same logic for its global tensor component, but adds a short-scale residual GP with compactly supported kernels so that long-range correlations and local irregularities are modeled by different components (Lei et al., 2022).

GLSKF uses a related but optimization-oriented device: the covariance norm, or $A^{(1)},\dots,A^{(d)}$3-norm,

$A^{(1)},\dots,A^{(d)}$4

together with matrix and tensor analogues based on Kronecker products. The paper interprets this generalized least-squares regularizer as corresponding to a Gaussian-process prior, and notes that quadratic variation smoothness arises as a special case when $A^{(1)},\dots,A^{(d)}$5 (Lei et al., 2024).

4. Inference, optimization, and computational structure

Once the kernelized construction is fixed, the learning problem varies sharply across formulations. In the TT classifier, the optimization stage is the standard SVM dual with the TT-based kernel matrix,

$A^{(1)},\dots,A^{(d)}$6

subject to the usual box and balance constraints. The learning rule is therefore conventional SVM optimization, but on a kernel that has been made consistent with TT contraction (Chen et al., 2020).

KFT adopts mean-field variational Bayes. The posterior is approximated by a factorized family

$A^{(1)},\dots,A^{(d)}$7

typically Gaussian, and inference proceeds by maximizing the ELBO

$A^{(1)},\dots,A^{(d)}$8

Because the model is conjugate Gaussian under the chosen likelihood and priors, the coordinate updates are closed form or analytically tractable; the paper also notes analytical ELBO derivations, analytical multivariate mean-field updates, and stochastic optimization ideas for scaling (Hu et al., 2020).

BKTF, BKTR, and BCKL all use MCMC, but with different targets. BKTF employs an element-wise Gibbs/slice-sampling procedure: latent factors are updated by Gibbs sampling, kernel length-scales by slice sampling, rank weights by a Gaussian conditional posterior, and noise precision by a Gamma conditional posterior (Lei et al., 2023). BKTR also uses Gibbs sampling for factor matrices and slice sampling for kernel hyperparameters, but improves mixing by integrating out the corresponding latent factor matrix before updating a hyperparameter (Lei et al., 2021). BCKL combines Gibbs and slice sampling for the global kernelized factorization with PCG-based computations for the local residual component, using a whitening trick and “imaginary observations” for efficiency (Lei et al., 2022).

The distributed flexible nonlinear tensor factorization model uses a sparse GP variational framework with inducing points and derives tight ELBOs that are additive over entries. This additive structure enables a MapReduce/SPARK implementation in which mappers compute local contributions to the ELBO and gradients, and a reducer sums full gradient vectors. The design is explicitly key-value-free so that it can exploit memory caching in fast MapReduce systems such as SPARK (Zhe et al., 2016).

GLSKF adopts an alternating least squares procedure with closed-form linear systems for the factor matrices and the residual tensor. The implementation relies on zero-padding and slicing operations based on projection matrices that preserve Kronecker structure, so that the required matrix-vector products can be handled by CG or PCG rather than explicit inversion (Lei et al., 2024). As a non-kernelized point of comparison, DFacTo shows that CP-ALS and CP-GD can be accelerated by replacing the tensor-times-Khatri-Rao bottleneck with two sparse matrix-vector multiplications per rank component and a distributed-memory design with limited synchronization (Choi et al., 2014).

5. Empirical domains and recurrent performance patterns

The empirical literature is heterogeneous, but several patterns recur. In few-sample, high-dimensional tensor classification, the TT-kernel approach is most advantageous when vectorization is especially damaging. On MNIST digit-pair classification, K-STTM-Prod and K-STTM-Sum usually outperform tensor-linear methods and often beat or match SVM and 3D CNN, with relatively small TT-ranks around 5 often sufficient for near-best accuracy. On StarPlus fMRI and CMU2008 fMRI, K-STTM-Prod and K-STTM-Sum achieve the best accuracy on all subjects; on CIFAR-10, mode-specific kernels improve performance over using the same kernel on all modes, especially when the color mode receives a different kernel (Chen et al., 2020).

In large-scale supervised regression with side information, KFT is reported to outperform or be competitive with LightGBM and FFM on three large-scale datasets while providing calibrated uncertainty estimates. The same experiments show robustness to uninformative side information in the form of constant features and Gaussian noise features, which is significant because such metadata often degrades direct feature-regression factorization models (Hu et al., 2020).

In spatiotemporal regression, BKTR reformulates the spatiotemporally varying coefficient model as a low-rank tensor regression problem with GP priors on spatial and temporal factor matrices. The reported experiments on synthetic and real-world data confirm superior performance and efficiency for model estimation and parameter inference, while preserving the interpretation of coefficients varying over space, time, and covariates (Lei et al., 2021).

In Bayesian optimization, BKTF is evaluated on Branin, Damavandi, Schaffer, Griewank in 3D and 4D, and Hartmann 6D. The reported findings are that BKTF consistently finds the global optimum faster than GP + EI/UCB baselines, escapes broad local optima on multimodal problems such as Damavandi, and remains effective on higher-dimensional Griewank and Hartmann functions under limited budgets. On hyperparameter tuning tasks for MNIST classification and Boston housing regression, BKTF achieves higher final classification accuracy or lower MSE, faster convergence, and lower variability across runs; one highlighted result is that it finds neural-network hyperparameters yielding 100% MNIST accuracy in fewer than four evaluations in all runs (Lei et al., 2023).

For multidimensional spatiotemporal tensors with missing data, BCKL and GLSKF both support the claim that a single low-rank component is often insufficient. In BCKL’s synthetic experiment, the best reported MAE and RMSE are 0.21 and 0.34, versus 0.26–0.29 and 0.43–0.47 for baselines; on traffic data, MODIS, and image inpainting, the kernelized tensor factorization captures global trend structure, while the local GP residual improves reconstruction of local irregularities and large missing regions (Lei et al., 2022). GLSKF reports the best or near-best performance across traffic speed imputation, color image inpainting, color video completion, and MRI image completion, and remains strong even at 99% missingness, which the paper attributes to the additive global-local structure (Lei et al., 2024).

The distributed flexible nonlinear tensor factorization model shows that entrywise GP kernelization can also be effective in extremely sparse recommendation-style settings. On small real-valued and binary tensor datasets it outperforms InfTucker significantly in almost all settings and outperforms all remaining baselines; on large datasets such as ACC, DBLP, and NELL it consistently outperforms GigaTensor and DinTucker, while on DBLP its average per-iteration time is 1.45 min versus 15.4 min for GigaTensor and 20.5 min for DinTucker. In click-through-rate prediction for online advertising, using a 4-mode tensor of user, advertisement, publisher, and page-section, the reported AUC improvement is about 20.7% over logistic regression and about 20.8% over linear SVM (Zhe et al., 2016).

6. Conceptual distinctions, misconceptions, and active directions

A common misconception is that kernelized tensor factorization is simply a standard kernel method applied after tensor vectorization. The TT literature explicitly rejects that interpretation: the feature map is applied to TT-core fibers, the mapped object remains a TT in feature space, and different kernels can be assigned to different tensor modes without losing mode identity (Chen et al., 2020). Another misconception is that nonlinear tensor factorization must inherit the Kronecker-product covariance of tensor Gaussian processes. The flexible nonlinear GP factorization shows that a GP can instead be defined on concatenated latent factors for selected entries, removing the need to model the whole tensor and allowing balanced training sets in extremely sparse settings (Zhe et al., 2016).

A further point of debate concerns whether a single low-rank kernelized component suffices. Both BCKL and GLSKF argue that low-rank kernelized factorization captures broad, smooth, or long-range structure but often misses short-scale, high-frequency variation. Their additive decompositions therefore separate a global kernelized tensor component from a locally correlated residual process, rather than forcing all variation into rank growth (Lei et al., 2022, Lei et al., 2024). This suggests that “kernelized tensor factorization” includes not only kernelized latent factors, but also multi-scale tensor models in which factorization and residual covariance play complementary roles.

It is also inaccurate to treat the area as purely deterministic regularization. Several prominent models are explicitly Bayesian. KFT uses a conjugate Bayesian regression formulation with variational inference and calibrated predictive uncertainty (Hu et al., 2020). BKTF is a fully Bayesian low-rank tensor surrogate for Bayesian optimization with posterior sampling and full uncertainty quantification (Lei et al., 2023). BKTR likewise uses GP priors on factor matrices and a full posterior sampling scheme for spatiotemporally varying coefficients (Lei et al., 2021).

The literature points toward several recurring technical directions. One is mode-specific kernel design, especially when tensor modes encode heterogeneous semantics (Chen et al., 2020). Another is further scalability through sparse GP approximations, inducing points, GRMF/NNGP variants, and other kernel simplifications for very large spatial or temporal dimensions (Lei et al., 2021). A third is the continued refinement of additive global-local models, which separate smooth low-rank structure from correlated local residuals in completion and spatiotemporal modeling (Lei et al., 2022, Lei et al., 2024). Taken together, these directions indicate that kernelized tensor factorization is less a single algorithmic template than a research program for embedding nonlinear, geometry-aware, and uncertainty-aware structure into tensor models without discarding the algebra of multiway data.