# Bayesian Methods in Tensor Analysis

Yiyao Shi
University of California, Irvine
Irvine, CA, US
<EMAIL_ADDRESS>
&Weining Shen
University of California, Irvine
Irvine, CA, US
<EMAIL_ADDRESS>

###### Abstract

Tensors, also known as multidimensional arrays, are useful data structures in
machine learning and statistics. In recent years, Bayesian methods have
emerged as a popular direction for analyzing tensor-valued data since they
provide a convenient way to introduce sparsity into the model and conduct
uncertainty quantification. In this article, we provide an overview of
frequentist and Bayesian methods for solving tensor completion and regression
problems, with a focus on Bayesian methods. We review common Bayesian tensor
approaches including model formulation, prior assignment, posterior
computation, and theoretical properties. We also discuss potential future
directions in this field.

_K_ eywords Imaging analysis $\cdot$ Posterior inference $\cdot$ Recommender
system $\cdot$ Tensor completion $\cdot$ Tensor decomposition $\cdot$ Tensor
regression

##  1 Introduction

Tensors, also known as multidimensional arrays, are higher dimensional
analogues of two-dimensional matrices. Tensor data analysis has gained
popularity in many scientific research and business applications, including
medical imaging [5], recommender systems [78], relational learning [93],
computer vision [83] and network analysis [53]. There is a vast literature on
studying tensor-related problems such as tensor decomposition [46, 71, 89],
tensor regression [25, 86], tensor completion [83], tensor clustering [5, 86],
tensor reinforcement learning and deep learning [86]. Among them, tensor
completion and tensor regression are two fundamental problems and we focus on
their review in this article.

Tensor completion aims at imputing missing or unobserved entries in a
partially observed tensor. Important applications of tensor completion include
providing personalized services and recommendations in context-aware
recommender systems (CARS) [78], restoring incomplete images collected from
magnetic resonance imaging (MRI) and computerized tomography (CT) [20], and
inpainting missing pixels in images and videos [58, 65]. In this review, we
divide tensor completion methods into trace norm based methods and
decomposition based methods, and introduce common approaches in each category.

Different from tensor completion, tensor regression investigates the
association between tensor-valued objects and other variables. For example,
medical imaging data such as brain MRI are naturally stored as a multi-
dimensional array, and tensor regression methods are applied to analyze their
relationship with clinical outcomes (e.g., diagnostic status, cognition and
memory score) [51, 87]. Based on the role that tensor-valued object plays in
the regression model, tensor regression methods can be categorized into tensor
predictor regression and tensor response regression.

Frequentist approaches have received a big success in tensor analysis [98, 5].
In recent years, Bayesian approaches have also gained popularity as they
provide a useful way to induce sparsity in tensor models and conduct
uncertainty quantification for estimation and predictions. In this article, we
will briefly discuss common frequentist approaches to solve tensor completion
and regression problems and focus on Bayesian approaches. We also review two
commonly used tensor decompositions, i.e., CANDECOMP/PARAFAC (CP)
decomposition [42] and the Tucker decomposition [94], since they are the
foundations for most Bayesian tensor models. For example, many Bayesian tensor
completion approaches begin with certain decomposition structure on the
tensor-valued data and then use Bayesian methods to infer the decomposition
parameters and impute the missing entries. Based on the decomposition
structures being utilized, we divide these methods into CP-based, Tucker-
based, and nonparametric methods. For tensor regression methods, we classify
the Bayesian tensor regression into Bayesian tensor predictor regression and
Bayesian tensor response regression. For each category, we review the prior
construction, model setup, posterior convergence property and sampling
strategies.

The rest of this article is organized as follows. Section 2 provides a
background introduction to tensor notations, operations and decompositions.
Section 3 and 4 review common frequentist approaches for tensor completion and
regression problems, respectively. Section 5 and 6 review Bayesian tensor
completion and regression approaches, including the prior construction,
posterior computing, and theoretical properties. Section 7 provides concluding
remarks and discusses several future directions for Bayesian tensor analysis.
Figure 1 shows an outline of our review.

Figure 1: Outline of this survey.

##  2 Background

In this section, we follow [46] to introduce the notations, definitions, and
operations of tensor. We also discuss two popular tensor decomposition
approaches and highlight some challenges in tensor analysis.

Figure 2: An example of first, second and third-order tensors.

###  2.1 Basics

#### Notation:

A tensor is a multidimensional array. The dimension of a tensor is also known
as mode, way, or order. A first-order tensor is a vector; a second-order
tensor is a matrix; and tensors of order three and higher are referred to as
higher-order tensors (see Figure 2). In this review, a tensor is denoted by
Euler script letter $\mathcal{X}\in\mathbb{R}^{n_{1}\times
n_{2}\times...\times n_{d}}$. Here $d$ is the order of tensor $\mathcal{X}$,
and $n_{k}$ is the marginal dimension of the $k$th mode ($k=1,2,...,d$). The
$(i_{1},i_{2},...,i_{d})$th element of the tensor $\mathcal{X}$ is denoted by
$x_{i_{1}i_{2}...i_{d}}$ for $i_{k}=1,2,...,n_{k}$ and $k=1,2,...,d$.
Subarrays of a tensor are formed through fixing a subset of indices in the
tensor. A fiber is a vector defined by fixing all but one indices of a tensor,
and a slice is a matrix created by fixing all the indices except for those of
two specific orders in the tensor. For instance, a third-order tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ has column, row and
tube fibers, which are respectively denoted by
$\mathcal{X}_{:i_{2}i_{3}},\mathcal{X}_{i_{1}:i_{3}}$, and
$\mathcal{X}_{i_{1}i_{2}:}$ (see Figure 3(a)(b)(c)). A third-order tensor also
has horizontal, lateral, and frontal slices, denoted by
$\mathcal{X}_{i_{1}::},\mathcal{X}_{:i_{2}:}$ and $\mathcal{X}_{::i_{3}}$,
respectively (see Figure 3(d)(e)(f)).

(a) Mode-1 (column) fibers: $\mathcal{X}_{:i_{2}i_{3}}$

(b) Mode-2 (row) fibers: $\mathcal{X}_{i_{1}:i_{3}}$

(c) Mode-3 (tube) fibers: $\mathcal{X}_{i_{1}i_{2}:}$

(d) Horizontal slices: $\mathcal{X}_{i_{1}::}$

(e) Lateral slices: $\mathcal{X}_{:i_{2}:}$

(f) Frontal slices: $\mathcal{X}_{::i_{3}}$

Figure 3: Example of fibers and slices of third-order tensor. Figure 4:
Rank-$r$ CP decomposition for a third-order tensor:
$\mathcal{X}\approx\sum_{j=1}^{r}w_{j}\boldsymbol{p}_{j}^{1}\circ\boldsymbol{p}_{j}^{2}\circ\boldsymbol{p}_{j}^{3}.$

#### Tensor Operations:

The norm of a tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times
n_{2}\times...\times n_{d}}$ is the square root of the sum of the squares of
all elements, i.e.,

$\|\mathcal{X}\|=\sqrt{\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\cdots\sum_{i_{d}=1}^{n_{d}}x_{i_{1}i_{2}...i_{d}}^{2}}.$
(1)

The inner product of two same-sized tensors
$\mathcal{X},\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ is the sum
of products of their corresponding entries, i.e.,

$\langle\mathcal{X},\mathcal{Y}\rangle=\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\cdots\sum_{i_{d}=1}^{n_{d}}x_{i_{1}i_{2}...i_{d}}y_{i_{1}i_{2}...i_{d}}.$
(2)

It immediately follows that
$\langle\mathcal{X},\mathcal{X}\rangle=\|\mathcal{X}\|^{2}.$ The tensor
Hadamard product of two tensors
$\mathcal{X}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ and
$\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ is denoted by
$\mathcal{X}*_{H}\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ and
each entry of $\mathcal{X}*_{H}\mathcal{Y}$ is the product of the
corresponding entries in tensors $\mathcal{X}$ and $\mathcal{Y}$:

$(\mathcal{X}*_{H}\mathcal{Y})_{i_{1}...i_{d}}=x_{i_{1}...i_{d}}\cdot
y_{i_{1}...i_{d}}.$ (3)

The tensor contraction product, also known as the Einstein product, of two
tensors $\mathcal{X}\in\mathbb{R}^{n_{1}\times...\times n_{d}\times
p_{1}\times...\times p_{k}}$ and
$\mathcal{Y}\in\mathbb{R}^{p_{1}\times...\times p_{k}\times
m_{1}\times...\times m_{q}}$ is denoted by
$\mathcal{X}*\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{d}\times
m_{1}\times...\times m_{q}}$ and defined as

$\begin{split}(\mathcal{X}*\mathcal{Y})_{i_{1},...,i_{d},j_{1},...,j_{q}}=\sum_{c_{1}=1}^{p_{1}}\cdots\sum_{c_{k}=1}^{p_{k}}x_{i_{1},...,i_{d},c_{1},...,c_{k}}y_{c_{1},...,c_{k},j_{1},...,j_{q}},\end{split}$
(4)

where $i_{g}=1,2,...,n_{g}$ for $g=1,2,...,d$, and $j_{s}=1,2,...,m_{s}$ for
$s=1,2,...,q$. Moreover, a $d$th-order tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times...\times n_{d}}$ is rank
one if it can be written as the outer product of $d$ vectors, i.e,

$\mathcal{X}=\boldsymbol{p}^{1}\circ\boldsymbol{p}^{2}\circ\cdots\circ\boldsymbol{p}^{d},$

where
$\boldsymbol{p}^{k}=(p_{1}^{k},p_{2}^{k},...,p_{n_{k}}^{k})\in\mathbb{R}^{n_{k}}$
$(k=1,2,...,d)$ is a vector, and the symbol “$\circ$” represents the vector
outer product. Each element of the tensor $\mathcal{X}$ is the product of
corresponding vector elements:
$x_{i_{1}i_{2}...i_{d}}=p_{i_{1}}^{1}p_{i_{2}}^{2}...p_{i_{d}}^{d}$ for
$i_{k}=1,2,...,n_{k}$ and $k=1,2,...,d$. A tensor $\mathcal{X}$ is rank $r$ if
$r$ is the smallest number such that $\mathcal{X}$ is the sum of $r$ outer
products of vectors:
$\mathcal{X}=\sum_{j=1}^{r}\boldsymbol{p}_{j}^{1}\circ\boldsymbol{p}_{j}^{2}\circ\cdots\circ\boldsymbol{p}_{j}^{d}$.

Tensor matricization, also known as tensor unfolding or flattening, is an
operation that transforms a tensor into a matrix. Given a tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times...\times n_{d}}$, the
$k$th-mode matricization arranges the mode-$k$ fibers to be columns of the
resulting matrix, which is denoted by $\boldsymbol{X}_{(k)}$ ($k=1,2,...,d$).
The element $(i_{1},i_{2},...,i_{d})$ of tensor $\mathcal{X}$ corresponds to
the entry $(i_{k},j)$ of $\boldsymbol{X}_{(k)}$, where $j=1+\sum_{t=1,t\neq
k}^{d}(i_{t}-1)J_{t}$ with $J_{t}=\prod_{m=1,m\neq k}^{t-1}n_{m}$. In
addition, a tensor can be transformed into a vector through tensor
vectorization. For a tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times...\times
n_{d}}$, the vectorization of $\mathcal{X}$ is denoted by
vec($\mathcal{X})\in\mathbb{R}^{\prod_{i=1}^{d}n_{i}}$. The element
$(i_{1},i_{2},...,i_{d})$ of tensor $\mathcal{X}$ corresponds to the element
$1+\sum_{t=1}^{d}(i_{t}-1)M_{t}$ of vec($\mathcal{X}$), where
$M_{t}=\prod_{m=1}^{t-1}n_{m}$.

The $k$-mode tensor matrix product of a tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{d}}$ with a
matrix $\boldsymbol{A}\in\mathbb{R}^{m\times n_{k}}$ is denoted by
$\mathcal{X}\times_{k}\boldsymbol{A}$, which is of size
$n_{1}\times\cdots\times n_{k-1}\times m\times n_{k+1}\times\cdots\times
n_{d}$. Elementwisely, we have
$(\mathcal{X}\times_{k}\boldsymbol{A})_{i_{1},\ldots,i_{k-1},j,i_{k+1},\ldots,i_{d}}=\sum_{i_{k}=1}^{n_{k}}\mathcal{X}_{i_{1},\ldots,i_{d}}\boldsymbol{A}_{ji_{k}}$.
The $k$-mode vector product of a tensor $\mathcal{X}\in\mathbb{R}^{n_{1}\times
n_{2}\times\cdots\times n_{d}}$ with a vector
$\boldsymbol{a}\in\mathbb{R}^{n_{k}}$ is denoted by
$\mathcal{X}\bar{\times}_{k}\boldsymbol{a}$, which is of size
$n_{1}\times\cdots\times n_{k-1}\times n_{k+1}\times\cdots\times n_{d}$.
Elementwisely, $(\mathcal{X}\bar{\times}_{k}\boldsymbol{a})_{i_{1}\ldots
i_{k-1}i_{k+1}\ldots
i_{d}}=\sum_{i_{k}=1}^{n_{k}}x_{i_{1}i_{2}...i_{d}}a_{i_{k}}.$

###  2.2 Tensor Decompositions

Tensor decompositions refer to methods that express a tensor by a combination
of simple arrays. Here we introduce two widely-used tensor decompositions and
discuss their applications.

#### CP decomposition:

The CANDECOMP/PARAFAC decomposition (CP decomposition) [42] factorizes a
tensor into a sum of rank-1 tensors. For a $d$th-mode tensor $\mathcal{X}$,
the rank-$r$ CP decomposition is written as

$\mathcal{X}\approx\sum_{j=1}^{r}w_{j}\boldsymbol{p}_{j}^{1}\circ\boldsymbol{p}_{j}^{2}\circ\cdots\circ\boldsymbol{p}_{j}^{d},$
(5)

where
$w_{j}\in\mathbb{R},\boldsymbol{p}_{j}^{k}\in\mathbb{S}^{n_{k}},j=1,...,r,k=1,2,...,d,\mathbb{S}^{n_{k}}=\\{\boldsymbol{a}\in\mathbb{R}^{n_{k}}|\|\boldsymbol{a}\|=1\\},$
and $\circ$ is the outer product. See Figure 4 for a graphical illustration of
CP decomposition. Sometimes the CP-decomposition is denoted by an
abbreviation:
$\mathcal{X}\approx[\\![\boldsymbol{W};\boldsymbol{P}^{1},\boldsymbol{P}^{2},...,\boldsymbol{P}^{d}]\\!],$
where $\boldsymbol{W}=\text{diag(}w_{1},...,w_{r})\in\mathbb{R}^{r\times r}$
is a diagonal matrix, and
$\boldsymbol{P}^{k}=[\boldsymbol{p}_{1}^{k},\boldsymbol{p}_{2}^{k}...,\boldsymbol{p}_{r}^{k}]\in\mathbb{R}^{n_{k}\times
r}$ are factor matrices. If tensor $\mathcal{X}$ admits a CP structure, then
the number of free parameters changes from $\prod_{i=1}^{d}n_{i}$ to
$r\times(\sum_{i=1}^{d}n_{i}-d+1)$.

If Equation (5) attains equality, the decomposition is called an exact CP
decomposition. Even for exact CP decomposition, there is no straightforward
algorithm to determine the rank $r$ of a specific tensor, and in fact the
problem is NP-hard [31]. In practice, most procedures numerically infer the
rank by fitting CP models with different ranks and choosing the one with the
best numerical performance.

Figure 5: Tucker decomposition of third-order tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$, where
$\mathcal{C}\in\mathbb{R}^{m_{1}\times m_{2}\times m_{3}}$ is the core tensor,
and $\boldsymbol{Q}^{k}\in\mathbb{R}^{n_{k}\times m_{k}}(k=1,2,3)$ are factor
matrices.

#### Tucker decomposition:

The Tucker decomposition factorizes a tensor into a core tensor multiplied by
a matrix along each mode. Given a $d$th-order tensor
$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times...\times n_{d}}$, the
Tucker decomposition is defined as

$\begin{split}\mathcal{X}\approx\mathcal{C}\times_{1}\boldsymbol{Q}^{1}\times_{2}\boldsymbol{Q}^{2}\times_{3}\cdots\times_{d}\boldsymbol{Q}^{d}=\sum_{j_{1}=1}^{m_{1}}\sum_{j_{2}=1}^{m_{2}}\cdots\sum_{j_{d}=1}^{m_{d}}c_{j_{1}j_{2}\dots
j_{d}}\boldsymbol{q}_{j_{1}}^{1}\circ\boldsymbol{q}_{j_{2}}^{2}\circ\cdots\circ\boldsymbol{q}_{j_{d}}^{d},\\\
\end{split}$ (6)

where $\mathcal{C}\in\mathbb{R}^{m_{1}\times m_{2}\times...\times m_{d}}$ is
the core tensor, $\boldsymbol{Q}^{k}\in\mathbb{R}^{n_{k}\times
m_{k}}(k=1,2,...,d)$ are factor matrices,
$c_{j_{1}j_{2}...j_{d}}\in\mathbb{R},\boldsymbol{q}_{j_{k}}^{k}\in\mathbb{S}^{n_{k}}(j_{k}=1,2,...,m_{k},k=1,2,...,d)$.
See Figure 5 for a graphical illustration of Tucker decomposition. The Tucker
decomposition can be denoted as
$\mathcal{X}\approx[\\![\mathcal{C};\boldsymbol{Q}^{1},\boldsymbol{Q}^{2},...,\boldsymbol{Q}^{d}]\\!].$
If $\mathcal{X}$ admits a Tucker structure, the number of free parameters in
$\mathcal{X}$ changes from $\prod_{i=1}^{d}n_{i}$ to
$\sum_{i=1}^{d}(n_{i}-1)\times m_{i}+\prod_{i=1}^{d}m_{i}$.

The $k$-rank of $\mathcal{X}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$,
denoted by rankk($\mathcal{X}$), is defined as the column rank of $k$th-mode
matricization matrix $\boldsymbol{X}_{(k)}$. Let
$R_{k}=$rank${}_{k}(\mathcal{X})$, then $\mathcal{X}$ is a
rank-$(R_{1},R_{2},...,R_{d})$ tensor. Trivially, $R_{k}\leq n_{k}$ for
$k=1,2,...,d$. When the equality in Equation (6) is attained, the
decomposition is called an exact Tucker decomposition. For a given tensor
$\mathcal{X}$, there always exists an exact Tucker decomposition with core
tensor $\mathcal{C}\in\mathbb{R}^{m_{1}\times m_{2}\times\cdots\times m_{d}}$
where $m_{k}$ is the true $k$-rank for $k=1,2,...,d$. Nevertheless, for one or
more $k$, if $m_{k}<R_{k}$, then the Tucker decomposition is not necessarily
exact; and if $m_{k}>R_{k}$, the model will contain redundant parameters.
Therefore, we usually want to identify the true tensor rank, i.e.,
$m_{k}=R_{k}$. While this job is easy for noiseless complete tensors, for
tensors obtained in real-world applications, which are usually noisy or
partially observed, the rank still needs to be determined by certain searching
procedure.

###  2.3 Challenges in tensor analysis

In tensor analysis, the ultrahigh dimensionality of the tensor-valued
coefficients and tensor data creates challenges such as heavy computational
burden and vulnerability to model overfitting. Conventional approaches usually
transform the tensors into vectors or matrices and utilize dimension reduction
and low-dimensional techniques. However, these methods are usually incapable
of accounting for the dependence structure in tensor entries. In past decades,
increasing number of studies impose decomposition structures on the tensor-
valued coefficients or data; thus naturally reducing the number of free
parameters, and avoiding the issues brought by high dimensionality.

In this paper, we focus on tensor regression and tensor completion problems,
where various decomposition structures including CP and Tucker have been
widely used. Specifically, a large proportion of tensor completion methods are
realized through inferring the decomposition structure based on the partially
observed tensor, and then impute the missing values through the inferred
decomposition structure. Also, tensor regression problems usually include
tensor-valued coefficients, and decomposition structures are imposed on the
coefficient tensor to achieve parsimony in parameters. In both situations, the
decomposition is not performed on a completely observed tensor, thus the rank
of decomposition cannot be directly inferred from the data. Most optimization-
based approaches determine the rank by various selection criteria, which may
suffer from low stability issues. Bayesian approaches perform automatic rank
inference through the introduction of sparsity-inducing priors. However,
efficient posterior computing and study of theoretical properties of the
posterior distributions are largely needed.

Low rankness and sparsity are commonly used assumptions in the literature to
help reduce the number of free parameters. For non-Bayesian methods,
oftentimes the task is formulated into an optimization problem, and the
assumptions are enforced by sparsity-inducing penalty functions. In
comparison, the Bayesian methods perform decompositions in the probabilistic
setting, and enforce sparsity assumptions through sparsity priors. We will
discuss more details about these approaches and how they resolve challenges in
the following sections.

##  3 Tensor Completion

Tensor completion methods aim at imputing missing or unobserved entries from a
partially observed tensor. It is a fundamental problem in tensor research and
has wide applications in numerous domains. For instance, tensor completion
techniques are extensively utilized in context-aware recommender systems
(CARS) to provide personalized services and recommendations. In ordinary
recommender systems, the user-item interaction data are collected and
formulated into a sparse interaction matrix, and the goal is to complete the
matrix and thus recommend individualized items to the users. In CARS, the
user-item interaction is collected with their contextual information (e.g.,
time and network), and the data is formulated as a high-order tensor where the
modes respectively represent users, items, and contexts. Therefore, the matrix
completion problem in ordinary recommender systems is transformed into a
tensor completion problem in CARS, and the purpose is to make personalized
recommendations to users based on the collected user-item interaction and
contextual information.

Apart from CARS, tensor completion is also applied in other research domains
including healthcare, computer vision and chemometrics. For example, medical
images collected from MRI and CT play important roles in the clinical
diagnosis process. Due to the high acquisition speed, oftentimes these high-
order images are incomplete, thus necessitating the application of tensor
completion algorithms. In the field of computer vision, the color videos can
be represented by a fourth-order tensor
(length$\times$width$\times$channel$\times$frame) by stacking the frames in
time order (see Figure 6). Tensor completion can be adopted to impute the
missing pixels and restore the lossy videos. As another example, chemometrics
is a discipline that employs mathematical, statistical and other methods to
improve chemical analysis. Tensor completion methods have been successfully
applied on various benchmark chemometric datasets including semi-realistic
amino acid fluorescence dataset [9] and flow injection dataset [66].

Tensor completion can be viewed as a generalization of matrix completion.
Since the matrix completion problems have been well-studied in the past few
decades, a natural way to conduct tensor completion is to unfold or slice the
tensor into a matrix (or matrices) and apply matrix completion methods to the
transformed matrix (or matrices). Nevertheless, the performance and efficiency
of such approaches are largely reduced by the loss of structural information
during the matricization process and excessive computational cost due to the
high dimensionality of the original tensor.

Under such circumstance, various methods that specifically focus on high-order
tensor completion have been developed. Among these techniques, a classical
group of approaches perform tensor completion through tensor decomposition.
Generally speaking, these methods impose a decomposition structure on a
tensor, and estimate the decomposition parameters based on the observed
entries of the tensor. After that, the estimated decomposition structure is
utilized to infer the missing entries of the tensor. Trace-norm based methods
are another popular class of tensor completion methods. These methods first
formulate tensor completion as a rank minimization problem, and then employ
tensor trace norm to further transform the task into a convex optimization
problem. Finally, various optimization techniques are applied to solve the
problem and thus complete the tensor. In this section we provide a brief
review of decomposition based and trace norm based tensor completion methods.
More details on these two methods and other variants of tensor completion
approaches can be found in Song et al. [83].

Figure 6: An illustration of color videos. Each frame of the video is
formulated as a third-order tensor, where the modes are length, width and
channels (RGB channels in this case). The frames are then stacked into a
fourth-order tensor according to time order.

###  3.1 Decomposition Based Methods

CP decomposition (5) and Tucker decomposition (6) are two most commonly used
decomposition-based methods for tensor completion. In [91], the authors
propose to perform CP decomposition on partially observed tensors by
iteratively imputing the missing values and estimating the latent vectors in
the CP structure. Specifically, in iteration $s~{}(s\geq 1)$, the partially
observed tensor $\mathcal{X}$ is completed by:

$\tilde{\mathcal{X}}^{(s)}=\mathcal{X}*_{H}\mathcal{M}+\mathcal{Y}^{(s)}*_{H}(\boldsymbol{1}-\mathcal{M}),$

where $*_{H}$ is the tensor Hadamard product defined in (3),
$\tilde{\mathcal{X}}^{(s)},\mathcal{X},\mathcal{Y}^{(s)},\mathcal{M}\in\mathbb{R}^{n_{1}\times...\times
n_{d}}$ are tensors of same size, $\tilde{\mathcal{X}}^{(s)}$ is the completed
tensor, $\mathcal{Y}^{(s)}$ is the interim low-rank approximation based on CP
decomposition, and $\mathcal{M}$ is the observation index tensor defined as

$\mathcal{M}_{i_{1}...i_{d}}=\begin{cases}1\quad\text{if
}\mathcal{X}_{i_{1}...i_{d}}\text{ is observed},\\\ 0\quad\text{if
}\mathcal{X}_{i_{1}...i_{d}}\text{ is unobserved}.\end{cases}$

After the tensor is completed, the decomposition parameters are estimated by
alternating least square optimization (ALS). The loop of tensor completion and
parameter estimation is repeated until convergence.

Similar approaches were adopted by Kiers et al. [43] and Kroonenberg [48] to
impute missing entries. These methods are referred to as EM-like methods,
because they can be viewed as a special expectation maximization (EM) method
when the residuals independently follow a Gaussian distribution. While the EM-
like methods are usually easy to implement, they may not perform well (e.g.,
slow convergence and converging to a local maximum) when there is a high
proportion of missing values.

Also based on the CP decomposition, Bro et al. [10] propose another type of
tensor completion methods called the Missing-Skipping (MS) method. It conducts
CP decomposition based only on the observed entries in the tensor, and is
typically more robust than the EM-like approaches when applied to tensors with
high proportion of missingness. In general, the MS methods seek to optimize
the following objective function

$L=\sum_{(i_{1},i_{2},...,i_{d})\in\Omega}\mathcal{D}(\mathcal{X}_{i_{1},...i_{d}},\mathcal{Y}_{i_{1},...,i_{d}}),$
(7)

where $\mathcal{X}\in\mathbb{R}^{n_{1}\times...n_{d}}$ is the observed tensor,
$\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ is the estimated
tensor with a CP structure, $\Omega$ is a set containing indices of all
observed entries in tensor $\mathcal{X}$, and $\mathcal{D}$ is an error
measure.

Under the optimization framework (7), Tomasi and Bro [91] define the error
measure $\mathcal{D}$ to be the squared difference between the observed and
estimated entry
$\mathcal{D}(\mathcal{X}_{i_{1},...i_{d}},\mathcal{Y}_{i_{1},...,i_{d}})=(\mathcal{X}_{i_{1},...i_{d}}-\mathcal{Y}_{i_{1},...,i_{d}})^{2}$,
and employ a modified Gauss-Newton iterative algorithm (i.e., Levenberg-
Marquardt method) [50, 63] to solve the optimization problem. Acar et al. [1]
utilize a weighted error and minimize the objective function based on first-
order gradient, which is shown to be more scalable to larger problem sizes
than the second-order optimization method in [91]. Moreover, the optimization
problem can be analyzed in a Bayesian setting by treating the error measure
$\mathcal{D}$ to be the log-likelihood function. We will discuss more details
about these probabilistic methods in Section 5.

Tucker decomposition is another widely utilized tool to conduct tensor
completion. While the CP-based completion approaches enjoy nice properties
including uniqueness (with the exception of elementary indeterminacies of
scaling and permutation) and nice interpretability of latent vectors, methods
that employ Tucker structure are able to accommodate more complex interaction
among latent vectors and are more effective than CP-based methods. Therefore,
in some real-world applications where the completion accuracy is prioritized
over the uniqueness and latent vector interpretation, Tucker-based approaches
are potentially more suitable than the CP-based methods.

Similar to CP-based methods, EM-like approaches and MS approaches are still
two conventional ways for Tucker-based tensor completion algorithms. Walczak
and Massart [96] and Andersson and Bro [2] discuss the idea of utilizing EM-
like Tucker decomposition to solve tensor completion in their earlier works.
This method is further combined with higher-order orthogonal iteration to
impute missing data [22]. As an example of MS Tucker decomposition,
Karatzoglou et al. [40] employ a stochastic gradient descent algorithm to
optimize the loss function based only on the observed entries. There are also
researches that develop the MS-based approaches under the Bayesian framework,
see Section 5 for more details.

In recent years, several studies utilize hierarchical tensor (HT)
representations to provide a generalization of classical Tucker models. Most
of the HT representation based methods are implemented using projected
gradient methods. For instance, Rauhut et al. [76, 77] employ Riemannian
gradient iteration method to establish an iterative hard thresholding
algorithm in their model. Besides, the Riemannian optimization is utilized to
construct the manifold for low-rank tensors in [14, 41, 47].

###  3.2 Trace Norm Based Methods

In [58] and a subsequent paper [57], the authors generalize matrix completion
to study tensors and solve the tensor completion problem by considering the
following optimization:

$\begin{split}\min_{\mathcal{Y}}&:\|\mathcal{Y}\|_{*},\\\
\text{s.t.}&:\mathcal{Y}_{\Omega}=\mathcal{X}_{\Omega},\\\ \end{split}$ (8)

where $\mathcal{X}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ is the observed
tensor, $\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ is the
estimated tensor, $\Omega$ is the set containing indices of all observed
entries in tensor $\mathcal{X}$, and $\|\cdot\|_{*}$ is the tensor trace norm.
They further define the trace norm to be a convex combination of the trace
norms of all unfolding matrices

$\|\mathcal{Y}\|_{*}:=\sum_{k=1}^{d}\alpha_{k}\|\boldsymbol{Y}_{(k)}\|_{*},$
(9)

where $\alpha_{i}$’s are non-negative weights satisfying
$\sum_{k=1}^{d}\alpha_{k}=1$. The optimization problem (8) with norm defined
in (9) is called a sum of nuclear norm (SNN) model. If the noise is Gaussian
distributed, the SNN model is equivalent with

$\min_{\mathcal{Y}}\frac{\lambda}{2}\|\mathcal{P}_{\Omega}(\mathcal{Y}-\mathcal{X})\|^{2}+\sum_{k=1}^{d}\alpha_{k}\|\boldsymbol{Y}_{(k)}\|_{*},$
(10)

where $\lambda>0$ is a tuning parameter, $\mathcal{P}_{\Omega}(\cdot)$ denotes
all the entries in the observed index set $\Omega$, $\|\cdot\|$ is the tensor
norm defined in (1), and $\|\cdot\|_{*}$ is the matrix trace norm. This
optimization problem can be solved by block coordinate decent algorithm [58]
and splitting methods (e.g., Alternating Direction Method of Multipliers,
ADMM) [20, 92, 82].

Using a similar model as (8), Mu et al. [65] propose to apply trace norm on a
balanced unfolding matrix instead of utilizing the summation of trace norms in
(9). In literature, it is also common to consider alternative norms such as
incoherent trace norm [103] and tensor nuclear norm [44, 106]. There are other
studies that impose trace norms on the factorized matrices rather than
unfolding matrices [59, 102, 62]; these approaches can be viewed as a
combination of decomposition based and trace norm based completion methods.

##  4 Tensor Regression

In this section, we review tensor regression methods, where the primary goal
is to analyze the association between tensor-valued objects and other
variables. Based on the role that the tensor plays in the regression, the
problem can be further categorized into tensor predictor regression (with
tensor-valued predictors and an univariate or multivariate response variable)
and tensor response regression (with tensor-valued response and predictors
that can be a vector, a tensor or even multiple tensors).

###  4.1 Tensor Predictor Regression

Many tensor predictor regression methods are motivated by the need of
analyzing anatomical magnetic resonance imaging (MRI) data. Usually stored in
the form of 3D images (see Figure 7 for an example), MRI presents the shape,
volume, intensity, or developmental changes in brain tissues and blood brain
barrier. These characteristics are closely related to the clinical outcomes
including diagnostic status, and cognition and memory score. It is hence
natural to formulate a tensor predictor regression to model the changes of
these scalar or vector-valued clinical outcomes with respect to the tensor-
valued MRI images.

Figure 7: An example of 3D magnetic resonance imaging (MRI). The image is
adapted with permissions from Science Photo Library. url:
https://www.sciencephoto.com/media/306963/view

In medical imaging analysis, conventional approaches are generally based on
vectorized data, either by summarizing the image data through a small number
of preidentified region of interest (ROIs), or by transforming the entire
image into a long vector. The former is highly dependent on the prior domain
knowledge and does not fully utilize the information in the raw image, and the
latter suffers from the high-dimensionality of voxels in the 3D image and
abandons important spatial information during the vectorization process. In
order to circumvent these limitations, a class of regression methods have been
developed to preserve the tensor structure. Specifically, given a univariate
response $Y$ (e.g. memory test score, disease status) and a tensor-valued
predictor $\mathcal{X}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ (e.g. 3D
image), Guo et al. [28] propose a linear regression model

$Y=\langle\mathcal{W},\mathcal{X}\rangle+b,$ (11)

where $\langle\cdot,\cdot\rangle$ is the tensor inner product defined in (2),
$\mathcal{W}$ is the coefficient tensor, and $b$ is the error. While model
(11) is a direct extension of classical linear regression model, the extension
can result in the explosion of number of unknown parameters. Specifically, the
coefficient tensor $\mathcal{W}$ includes $\prod_{i=1}^{d}n_{i}$ free
parameters, which far exceeds the typical sample size. To address this issue,
Guo et al. [28] impose a rank-$r$ CP structure (5) on $\mathcal{W}$, which
reduces the number of parameters in $\mathcal{W}$ to $r\sum_{i=1}^{d}n_{i}$.

Li et al. [55] extend model (11) to the multivariate response
$\boldsymbol{Y}=(Y_{1},Y_{2},...,Y_{q})^{\top}$ case, where each marginal
response $Y_{k}~{}(1\leq k\leq q)$ is assumed to be the summation of
$\langle\mathcal{X},\mathcal{B}_{k}\rangle$ and an error term, where
$\mathcal{X}$ is the predictor tensor, and
$\mathcal{B}_{k}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ is the coefficient
tensor. Under the assumption that the coefficients share common features, the
coefficient tensors are further formulated into a stack
$\mathcal{B}=[\mathcal{B}_{1},...,\mathcal{B}_{q}]\in\mathbb{R}^{n_{1}\times...\times
n_{d}\times q}$, on which a CP structure is imposed for parameter number
reduction.

Additionally, Zhou et al. [116] integrate model (11) with the generalized
linear regression framework, and incorporate the association between response
and other adjusting covariates into the model. Consider a scalar response $Y$,
a tensor-valued predictor $\mathcal{X}\in\mathbb{R}^{n_{1}\times...\times
n_{d}}$ and vectorized covariates $\boldsymbol{z}\in\mathbb{R}^{n_{0}}$ (e.g.,
demographic features), the generalized linear model is given by

$g\\{\mathbb{E}(Y)\\}=b+\boldsymbol{\gamma}^{\top}\boldsymbol{z}+\langle\mathcal{W},\mathcal{X}\rangle,$
(12)

where $\boldsymbol{\gamma}$ is the vector coefficient for $\boldsymbol{z}$,
$g(\cdot)$ is a link function, and $\mathcal{W}$ is the coefficient tensor
where a CP structure is assumed. Still for model (12), Li et al. [54] impose a
Tucker decomposition on $\mathcal{W}$, and demonstrate that the Tucker
structure allows for more flexibility.

In order to accommodate longitudinal correlation of the data in the imaging
analysis, Zhang et al. [105] extend model (12) under the generalized
estimating equation setting and establish asymptotic properties of the method.
Hao et al. [30] show that the linearity assumption in (11) may be violated in
some applications, and propose a nonparametric extension of (11) that
accommodates the nonlinear relationship between the response and tensor
predictor. Zhang et al. [104] use importance sketching to reduce the high
computational cost associated with the low-rank factorization in tensor
predictor regression, and establish the optimality of their method in terms of
reducing mean squared error under the Tucker structure assumption and
randomized Gaussian design. Beyond the regression framework, Wimalawarne et
al. [98] propose a binary classification method by considering a logistic loss
function and various tensor norms for regularization.

###  4.2 Tensor Response Regression

While the main focus of tensor predictor regression is analyzing the effects
of tensors on the response variables, researchers are also interested in
studying how tensor-valued outcomes change with respect to covariates. For
example, an important question in MRI studies is to compare the scans of
brains between subjects with neurological disorder (e.g., attention deficit
disorder) and normal controls, after adjusting for other covariates such as
age and sex. This problem can be formulated as a tensor response regression
problem where the MRI data, usually taking form of a three-dimensional image,
is the tensor-valued response, and other variables are predictors. Apart from
medical imaging analysis, tensor response regression is also useful in
advertisement industry. For example, the click-through rate (CTR) of digital
advertisements is often considered to be a significant indicator of the
effectiveness of an advertisement campaign. Thus an important business
question is to understand how CTR is affected by different features. Since the
CTR data can be formulated as a high-dimensional tensor (see Figure 8), we can
develop a regression model to address this problem, where the click-through
rate on target audience is the tensor-valued response, and the features of
advertisements are predictors of interest.

Figure 8: An illustration of click through rate data, which is formulated as a
third-mode tensor, where each voxel represents the click-through rate of user
$i$ reacting to advertisements from publisher $j$ at time $k$.

Given a $d$th-order tensor response
$\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{d}}$ and a vector predictor
$\boldsymbol{x}\in\mathbb{R}^{q}$, Rabusseau and Kadri [72] and Sun and Li
[87] propose a linear regression model

$\mathcal{Y}=\mathcal{B}\bar{\times}_{d+1}\boldsymbol{x}+\mathcal{E},$ (13)

where $\mathcal{B}\in\mathbb{R}^{n_{1}\times n_{2}\times...\times n_{d}\times
q}$ is an $(d+1)$th-order tensor coefficient, $\mathcal{E}$ is an error tensor
independent of $\boldsymbol{x}$, and $\bar{\times}_{d+1}$ is the $(d+1)$-mode
vector product. Without loss of generality, the intercept is set to be zero to
simplify the presentation.

Both studies [72, 87] propose to estimate the coefficients $\mathcal{B}$ by
solving an optimization problem, which consists of a squared tensor norm of
the difference between observed and estimated response
$\|\mathcal{Y}-\mathcal{B}\bar{\times}_{d+1}\boldsymbol{x}\|^{2}$ and a
sparsity structure. In Rabusseau and Kadri [72], the sparsity is obtained by a
$L_{2}$-penalty of parameters. In Sun and Li [87], the sparsity structure is
realized through a hard-thresholding constraint on the coefficients. For both
studies, decomposition structures are imposed on the tensor coefficient
$\mathcal{B}$ to facilitate parsimonious estimation of high-dimensional
parameters.

Lock [60] further extends (13) to a tensor-on-tensor regression model,
allowing a predictor of arbitrary order. Given $N$ independent samples, the
responses can be stacked into a tensor $\mathcal{Y}\in\mathbb{R}^{N\times
m_{1}\times m_{2}\times...\times m_{q}}$, and the predictors are denoted by
$\mathcal{X}\in\mathbb{R}^{N\times n_{1}\times n_{2}\times...\times n_{d}}$.
Lock [60] propose the following model:

$\mathcal{Y}=\mathcal{X}*\mathcal{B}+\mathcal{E},$ (14)

where $*$ is tensor contraction product defined in (4),
$\mathcal{B}\in\mathbb{R}^{n_{1}\times...\times n_{d}\times
m_{1}\times...\times m_{q}}$ is the coefficient tensor and $\mathcal{E}$
denotes the error. A CP structure is imposed on $\mathcal{B}$ to achieve
parsimony in parameters. The estimation of $\mathcal{B}$ is also transformed
into an optimization problem, and a $L_{2}$-penalty is included in the loss
function to prevent over-fitting. Under a similar modeling framework, Gahrooei
et al. [19] develop a multiple tenor-on-tensor regression model, where the
predictors are a set of tensors with various orders and sizes.

Based on (14), Li and Zhang [51] propose a tensor response regression that
utilizes the envelop method to remove redundant information from the response.
Raskutti et al. [75] analyze the tensor regression problem with convex and
weakly decomposable regularizers. In their regression model, either or both
the predictors and responses are tensors, and the low-rankness assumption is
realized by a nuclear norm penalty. Zhou et al. [117] focus on tensor
regression where the response is a partially observed dynamic tensor, and
impose low-rankness, sparsity and temporal smoothness constraints in the
optimization. Chen et al. [11] extend model (14) to the generalized tensor
regression setting and utilize a projected gradient descent algorithm to solve
the non-convex optimization.

##  5 Bayesian Methods in Tensor Completion

In Section 3.1, we mention that the tensor completion tasks can be realized by
performing decomposition on partially observed tensors and using the inferred
decomposition structure to impute the missing data (e.g., the Missing-Skipping
methods). The Bayesian tensor decomposition methods can be naturally applied
to study partially observed tensors. Generally, a large proportion of Bayesian
decomposition methods are based on CP (5)) or Tucker decomposition (6). A
class of nonparametric methods have also been proposed to model complex non-
linear interaction among latent factors. Recently, more decomposition
structures are analyzed under the Bayesian framework (e.g., tensor ring
decomposition [61], tensor train decomposition [38] and neural decomposition
[33]). A summary of the methods discussed in this section is given in Table 1.

###  5.1 Bayesian CP-Based Decomposition

Under the Bayesian framework, Xiong et al. [99] utilize a CP decomposition
based method to model time-evolving relational data in recommender systems. In
their study, the observed data is formed into a three-dimensional tensor
$\mathcal{R}\in\mathbb{R}^{N\times M\times K}$, where each entry
$\mathcal{R}_{ij}^{k}$ denotes user $i$’s rate on item $j$ given time $k$. A
CP structure (5) is then imposed on $\mathcal{R}$:

$\mathcal{R}\approx\sum_{d=1}^{D}\boldsymbol{U}_{d:}\circ\boldsymbol{V}_{d:}\circ\boldsymbol{T}_{d:}=[\\![\boldsymbol{U},\boldsymbol{V},\boldsymbol{T}]\\!],$
(15)

where $\boldsymbol{U},\boldsymbol{V},\boldsymbol{T}$ are latent factors
respectively corresponding to user, item, and time; and
$\boldsymbol{U}_{d:},\boldsymbol{V}_{d:},\boldsymbol{T}_{d:}$ represent the
$d$th-row of $\boldsymbol{U},\boldsymbol{V}$ and $\boldsymbol{T}$. Xiong et
al. [99] assign a Gaussian prior on the continuous entries $R_{ij}^{k}$
conditional on $\boldsymbol{U},\boldsymbol{V},\boldsymbol{T}$ as follows,

$R_{ij}^{k}|\boldsymbol{U},\boldsymbol{V},\boldsymbol{T}\sim\mathcal{N}(\langle\boldsymbol{U}_{:i},\boldsymbol{V}_{:j},\boldsymbol{T}_{:k}\rangle,\alpha^{-1}),$
(16)

where $\alpha$ is the precision, and
$\langle\boldsymbol{U}_{:i},\boldsymbol{V}_{:j},\boldsymbol{T}_{:k}\rangle$ is
the inner product of three $D$-dimensional vectors defined as

$\langle\boldsymbol{U}_{:i},\boldsymbol{V}_{:j},\boldsymbol{T}_{:k}\rangle=\sum_{d=1}^{D}U_{di}V_{dj}T_{dk}.$

A complete Bayesian setting requires full specification of the parameter
priors. In the study, multivariate Gaussian priors are put on the latent
vectors corresponding to users and items

$\displaystyle\boldsymbol{U}_{i}\sim\mathcal{N}(\boldsymbol{\mu}_{U},\boldsymbol{\Lambda}_{U}^{-1}),\quad
i=1,2,...,N,$ (17)
$\displaystyle\boldsymbol{V}_{j}\sim\mathcal{N}(\boldsymbol{\mu}_{V},\boldsymbol{\Lambda}_{V}^{-1}),\quad
j=1,2,...,M,$ (18)

and each time feature vector is assumed to depend only on its immediate
predecessor due to temporal smoothness:

$\displaystyle\boldsymbol{T}_{k}\sim\mathcal{N}(\boldsymbol{T}_{k-1},\boldsymbol{\Lambda}_{T}^{-1}),\quad
k=1,2,...,K,$ (19)
$\displaystyle\boldsymbol{T}_{0}\sim\mathcal{N}(\boldsymbol{\mu}_{T},\boldsymbol{\Lambda}_{T}^{-1}).$
(20)

Moreover, Xiong et al. [99] consider a hierarchical Bayesian structure where
the hyper-parameters
$\alpha,\boldsymbol{\Theta}_{U}\equiv\\{\boldsymbol{\mu}_{U},\boldsymbol{\Lambda}_{U}\\},\boldsymbol{\Theta}_{V}\equiv\\{\boldsymbol{\mu}_{V},\boldsymbol{\Lambda}_{V}\\},$
and
$\boldsymbol{\Theta}_{T}\equiv\\{\boldsymbol{\mu}_{T},\boldsymbol{\Lambda}_{T}\\}$
are viewed as random variables, and their prior distributions (i.e., hyper-
priors), denoted by $p(\cdot)$, are

$\begin{split}p(\alpha)=\mathcal{W}(\alpha|\tilde{W}_{0},\tilde{\nu}_{0}),\qquad\qquad\qquad\qquad\qquad\\\
p(\boldsymbol{\Theta}_{U})=p(\boldsymbol{\mu}_{U}|\boldsymbol{\Lambda}_{U})p(\boldsymbol{\Lambda}_{U})=\mathcal{N}(\boldsymbol{\mu}_{0},(\beta_{0}\boldsymbol{\Lambda}_{U})^{-1})\mathcal{W}(\boldsymbol{\Lambda}_{U}|\boldsymbol{W}_{0},\nu_{0}),\\\
p(\boldsymbol{\Theta}_{V})=p(\boldsymbol{\mu}_{V}|\boldsymbol{\Lambda}_{V})p(\boldsymbol{\Lambda}_{V})=\mathcal{N}(\boldsymbol{\mu}_{0},(\beta_{0}\boldsymbol{\Lambda}_{V})^{-1})\mathcal{W}(\boldsymbol{\Lambda}_{V}|\boldsymbol{W}_{0},\nu_{0}),\\\
p(\boldsymbol{\Theta}_{T})=p(\boldsymbol{\mu}_{T}|\boldsymbol{\Lambda}_{T})p(\boldsymbol{\Lambda}_{T})=\mathcal{N}(\boldsymbol{\mu}_{0},(\beta_{0}\boldsymbol{\Lambda}_{T})^{-1})\mathcal{W}(\boldsymbol{\Lambda}_{T}|\boldsymbol{W}_{0},\nu_{0}).\\\
\end{split}$ (21)

Here $\mathcal{W}(\boldsymbol{\Lambda}|\boldsymbol{W}_{0},\nu_{0})$ is the
Wishart distribution of a $D\times D$ random matrix $\boldsymbol{\Lambda}$
with $\nu_{0}$ degrees of freedom and a $D\times D$ scale matrix
$\boldsymbol{W}_{0}$:

$\mathcal{W}(\boldsymbol{\Lambda}|\boldsymbol{W}_{0},\nu_{0})=|\boldsymbol{\Lambda}|^{(\nu_{0}-D-1)/2}\exp\left(-\frac{\text{Tr}(\boldsymbol{W}_{0}^{-1}\boldsymbol{\Lambda})}{2}\right).$

The priors in (21) are conjugate priors for the Gaussian parameters to help
simplify the posterior computation. The parameters in the hyper-priors
$\boldsymbol{\mu}_{0},\beta_{0},\boldsymbol{W}_{0},\nu_{0},\tilde{W}_{0}$ and
$\tilde{\nu}_{0}$ can be chosen by prior knowledge or tuned by model training.

The Bayesian model in (16)–(21) is called Bayesian Probabilistic Tensor
Factorization (BPTF). The posterior distribution of the BPTF model is obtained
by Markov Chain Monte Carlo (MCMC) with Gibbs sampling [21]. While Xiong et
al. [99] use the BPTF model to perform tensor decomposition on continuous
rating data in recommender systems, similar priors have been adapted in other
applications and data types. For example, Chen et al. [12] formulate the
spatio-temporal traffic data as a third-order tensor (road
segment$\times$day$\times$time of day), where a CP structure is assumed and a
Gaussian-Wishart prior is put on the latent factors to promote conjugacy in
priors. A similar model has been used to study multi-relational network [81],
where the interaction data form a partially symmetric third-order tensor and
the tensor entries are binary indicators of whether there exists a certain
type of relationship. Correspondingly, a sigmoid function is employed in (16)
to map the outer product of latent factors onto the range $[0,1]$.

In addition, Schein et al. [79] develop a Poisson tensor factorization (PTF)
method to deal with dyadic interaction data in social networks. Specifically,
the interaction data is formulated as a fourth-order tensor $\mathcal{X}$,
where $\mathcal{X}_{ijat}$ denotes the number of interactions within a
discrete time interval $t$ involving a particular sender $i$, receiver $j$,
and action-type $a$. A Poisson distribution is employed to connect the CP
structure to the count-valued data:

$\mathcal{X}_{ijat}\sim\text{Poisson}(\sum_{k=1}^{K}\theta_{ik}^{s}\theta_{jk}^{r}\psi_{ak}\delta_{tk}).$
(22)

Gamma priors are then assigned to the latent factors,

$\begin{split}\theta_{ik}^{s}&\sim\text{Gamma}(a,b),\\\
\theta_{jk}^{r}&\sim\text{Gamma}(a,b),\\\ \psi_{ak}&\sim\text{Gamma}(c,d),\\\
\delta_{tk}&\sim\text{Gamma}(e,f).\end{split}$ (23)

Schein et al. [79] then represent the Poisson likelihood (22) as a sum of $K$
independent Poisson random variables, and derive a Variational Bayesian (VB)
algorithm to make inference on the posterior distribution.

All the aforementioned methods assume that the interactions among the latent
factors are multi-linear, which may not necessarily hold in practice. To
address this issue, Liu et al. [56] consider a neural CP decomposition that
exploits both the neural networks and probabilistic methods to capture
potential nonlinear interactions among the tensor entries. Given a tensor
$\mathcal{X}$ and the latent matrices in its CP structure
$\boldsymbol{U}^{1},...,\boldsymbol{U}^{D}$, the distribution of $\mathcal{X}$
conditional on $\boldsymbol{U}^{1},...,\boldsymbol{U}^{D}$ is given by

$p(\mathcal{X}|\\{\boldsymbol{U}^{d}\\}_{d=1}^{D})=\prod_{i_{1},...,i_{D}}\mathcal{N}(x|\mu,\sigma^{2}),$

where $x,\mu,\sigma^{2}$ are respectively short forms of
$x_{i_{1}...i_{D}},\mu_{i_{1}...i_{D}}$ and $\sigma^{2}_{i_{1}...i_{D}}$. In
order to accommodate nonlinear interaction between latent factors, $\mu$ and
$\sigma^{2}$ are defined as functions of $\boldsymbol{u}$
($\mu=\mu(\boldsymbol{u}),\sigma^{2}=\sigma^{2}(\boldsymbol{u})$), where
$\boldsymbol{u}=(U_{i_{1}:}^{1},...,U_{i_{D}:}^{D})\in\mathbb{R}^{DR}$ is a
long vector generated by concatenating one by one the $i_{d}$-th row of the
factor matrix $U^{d}$. In particular, the two functions $\mu(\cdot)$ and
$\sigma^{2}(\cdot)$ are modeled by two neural networks with the same input
$U_{i_{1}:}^{1},...,U_{i_{D}:}^{D}$

$\begin{split}\mu&=\boldsymbol{w}_{\mu}^{\top}\boldsymbol{h}(\boldsymbol{u})+b_{\mu},\\\
\log\sigma^{2}&=\boldsymbol{w}_{\sigma}^{\top}\boldsymbol{h}(\boldsymbol{u})+b_{\sigma},\end{split}$

where $\boldsymbol{h}(\boldsymbol{u})$ is a nonlinear hidden layer shared by
these two neural networks, and is defined as a tanh activation function in
[56]:

$\boldsymbol{h}(\boldsymbol{u})=tanh(\boldsymbol{W}^{\top}\boldsymbol{u}+\boldsymbol{b}).$

As discussed in Section 2.2, determining the rank of CP can be challenging in
practice. Even for a noise-free tensor, its rank specification is an NP-hard
problem [31]. In order to determine the CP rank, a common practice is to fit
models with different ranks and choose the best rank based on certain
criteria. Nevertheless, this approach may suffer from low stability issue and
high computational cost. An alternative approach is to use sparsity-inducing
priors. For example, in [74] and a subsequent work [73], the authors propose a
Bayesian low-rank CP decomposition method, which utilizes the multiplicative
gamma process (MGP) prior [4] to automatically infer the rank. Specifically,
given a CP structure

$\mathcal{X}=\sum_{r=1}^{R}\lambda_{r}\cdot\boldsymbol{u}_{r}^{(1)}\circ\boldsymbol{u}_{r}^{(2)}\circ\cdots\circ\boldsymbol{u}_{r}^{(K)},$

the following priors are put on the vector
$\boldsymbol{\lambda}=(\lambda_{1},\lambda_{2},...,\lambda_{R})$:

$\displaystyle\lambda_{r}\sim\mathcal{N}(0,\tau_{r}^{-1}),\quad 1\leq r\leq R$
(24)
$\displaystyle\tau_{r}=\prod_{l=1}^{r}\delta_{l},\quad\delta_{l}\sim\text{Gamma}(a_{c},1),\quad
a_{c}>1.$ (25)

In MGP prior, as $r$ increases, the precision $\tau_{r}$ takes large values
hence shrinks $\lambda_{r}$ towards zero. Small $\lambda_{r}$ values indicate
that the term
$\lambda_{r}\cdot\boldsymbol{u}_{r}^{(1)}\circ\boldsymbol{u}_{r}^{(2)}\circ\cdots\circ\boldsymbol{u}_{r}^{(K)}$
does not have a significant impact on the CP structure, hence could be removed
from the model. Two generalizations of MGP prior are further developed,
including truncation based variant MGP-CPt and the adaptive variant MGP-CPa,
to automatically infer the rank $R$ [74, 73].

Hu et al. [37] develop a Bayesian non-negative tensor factorization that deals
with count data and automatically infers the rank of CP decomposition. In
their work, the Poisson distribution is utilized to establish a connection
between CP structure and the count-valued data. Given a tensor
$\mathcal{Y}\in\mathbb{R}^{n_{1}\times...\times n_{K}}$ and its entry
$\boldsymbol{i}=\\{i_{1},...,i_{K}\\}$, we have

$\mathcal{Y}_{\boldsymbol{i}}\sim\text{Poisson}\left(\sum_{r=1}^{R}\lambda_{r}\prod_{k=1}^{K}u_{i_{k}r}^{(k)}\right).$

The non-negativity constraints on the factor matrices
$\boldsymbol{U}^{(1)},...,\boldsymbol{U}^{(K)}$
($\boldsymbol{U}^{(k)}=[\boldsymbol{u}_{1}^{(k)},...,\boldsymbol{u}_{R}^{(k)}],k=1,2,...,K$)
are naturally satisfied by imposing Dirichlet priors on the factors
$\boldsymbol{u}_{r}^{(k)}=[u_{1r}^{(k)},...,u_{i_{k}r}^{(k)}]^{\top}$:

$\boldsymbol{u}_{r}^{(k)}\sim\text{Dir}(a^{(k)},...,a^{(k)}),$

and a gamma-beta hierarchical prior is put on $\lambda_{r}$ to promote the
automatic rank specification:

$\displaystyle\lambda_{r}\sim\text{Gamma}(g_{r},\frac{p_{r}}{1-p_{r}}),$ (26)
$\displaystyle p_{r}\sim\text{Beta}(c\epsilon,c(1-\epsilon))~{}~{}~{}\text{for
some}~{}c>0.$ (27)

Similar to the MGP prior in (24) and (25), the gamma-beta hierarchical prior
in (26) and (27) also shrinks $\lambda_{r}$ to zero as $r$ increases, and is
thus able to select the CP rank. This model is also extended to binary data by
adding an additional layer
$b_{\boldsymbol{i}}=\boldsymbol{1}(y_{\boldsymbol{i}}\geq 1)$, which takes a
count-valued entry $y_{\boldsymbol{i}}$ in $\mathcal{Y}$ and thresholds this
latent count at one to generate binary-valued entries $b_{\boldsymbol{i}}$
[36].

Instead of imposing sparsity priors on the core elements of CP structure, Zhao
et al. [108] place a hierarchical prior over the latent factors. Let
$\mathcal{X}\in\mathbb{R}^{I_{1}\times\cdots\times I_{N}}$ have a CP structure

$\mathcal{X}=[\\![\boldsymbol{A}^{(1)},...,\boldsymbol{A}^{(N)}]\\!],$

where
$\boldsymbol{A}^{(n)}=[\boldsymbol{a}_{1}^{(n)},...,\boldsymbol{a}_{I_{n}}^{(n)}]$
$(n=1,2,...,N)$ are latent factors. Let
$\boldsymbol{\lambda}=[\lambda_{1},...,\lambda_{R}]$ and
$\boldsymbol{\Lambda}=$ diag($\boldsymbol{\lambda}$). The prior distribution
of $\boldsymbol{A}^{(n)}$ is

$p(\boldsymbol{A}^{(n)}|\boldsymbol{\lambda})=\prod_{i_{n}=1}^{I_{n}}\mathcal{N}(\boldsymbol{a}_{i_{n}}^{(n)}|\boldsymbol{0},\boldsymbol{\Lambda}^{-1}),\quad
n=1,2,\ldots,N.$

A hyperprior is further defined over $\boldsymbol{\lambda}$, which is
factorized over latent dimensions

$p(\boldsymbol{\lambda})=\prod_{r=1}^{R}\text{Gamma}(\lambda_{r}|c_{0}^{r},d_{0}^{r}).$

Here $R$ is a pre-specified maximum possible rank. The latent vectors (the
$r$-th row of all latent matrices) will shrink to a zero vector as
$\lambda_{r}^{-1}$ approaches to zero. This model can also accommodate various
types of outliers and non-Gaussian noises through the introduction of a sparse
structure, and the tradeoff between the low-rank approximation and the sparse
representation can be learned automatically by maximizing the model evidence
[111].

In real-world applications including recommender systems, image/video data
analysis and internet networks, the data are sometimes produced continuously
(streaming data). Therefore it is of interest to generalize the tensor
decomposition models to analyze such data under a real-time manner, where the
model parameters can be updated efficiently upon receiving new data without
retrieving previous entries. To this end, a class of streaming tensor
decomposition methods have been developed, and some are analyzed under the
Bayesian CP network [107, 15, 18]. In general, these algorithms start with a
prior distribution of unknown parameters and then infer a posterior that best
approximates the joint distribution of these parameters upon the arrival of
new streaming data. The estimated posterior is then used as the prior for the
next update. These methods are implemented either by streaming variational
Bayes (SVB) [107, 15], or assume-density filtering (ADF) and expectation-
propagation (EP) [18].

###  5.2 Tucker-based Bayesian Decomposition Methods

Compared to the CP decomposition, the Tucker structure (6) can model more
complex interaction between latent factors. One of the early works that employ
a probabilistic Tucker structure is proposed by Chu and Ghahramani [13], where
a probabilistic framework called pTucker is developed to perform decomposition
on partially observed tensors. Given a continuous third-order tensor
$\mathcal{Y}\in\mathbb{R}^{n\times m\times d}$, a Gaussian distribution is
assigned to each entry of tensor $\mathcal{Y}$,

$\mathcal{Y}_{ijr}|\mathcal{T}\sim\mathcal{N}(\mathcal{F}_{ijr},\sigma^{2}).$

Here $\mathcal{F}$ has a Tucker structure with a core tensor $\mathcal{T}$

$\mathcal{F}_{ijr}=\text{vec}(\mathcal{T})^{\top}(\boldsymbol{v}_{r}\otimes\boldsymbol{z}_{j}\otimes\boldsymbol{x}_{i}),$

where $\otimes$ is the Kronecker product, and
$\boldsymbol{v}_{r},\boldsymbol{z}_{j}$ and $\boldsymbol{x}_{i}$ are latent
vectors. Next, independent standard normal distributions are specified over
the entries in $\mathcal{T}$ as priors:

$\mathcal{T}_{kls}\sim\mathcal{N}(0,1),\quad\forall k,l,s.$

By integrating out the core tensor $\mathcal{T}$ from the joint
$\prod_{i,j,r}p(\mathcal{Y}_{ijr}|\mathcal{T})\prod_{k,l,s}p(\mathcal{T}_{kls})$,
the distribution over the observational array still follows a Gaussian
distribution:

$\text{vec}(\mathcal{Y})\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{U}\boldsymbol{U}^{\top}+\sigma^{2}\boldsymbol{I}),$

where $\text{vec}(\mathcal{Y})$ is the vectorized tensor, $\sigma^{2}$ is the
noise level, and
$\boldsymbol{U}=\boldsymbol{V}\otimes\boldsymbol{Z}\otimes\boldsymbol{X}$
where $\boldsymbol{V},\boldsymbol{Z}$ and $\boldsymbol{X}$ are latent
matrices. To complete the Bayesian framework, standard normal distributions
are further used as priors for latent components
$\boldsymbol{X},\boldsymbol{Z}$ and $\boldsymbol{V}$. Finally, the latent
factors are estimated by maximum a posteriori (MAP) with gradient descent.

While the MAP method provides an efficient alternative to perform point
estimation for latent factors, it also has significant disadvantages including
vulnerability to overfitting and incapability to quantify the uncertainty in
the parameters. To this end, various approaches seek to provide a fully
Bayesian treatment through inferring the posterior distribution of parameters.
For instance, Hayashi et al. [32] utilize the expectation maximization (EM)
method that combines the Laplace approximation and the Gaussian process to
perform posterior inference of latent factors. They use the exponential family
distributions to connect the Tucker structure with the observed tensor, thus
developing a decomposition method compatible with various data types. In
addition, Schein et al. [80] propose a Bayesian Poisson Tucker decomposition
(BPTD) that uses MCMC with Gibbs sampling for posterior inference. Especially
focusing on the count-valued observed tensor, that method puts Poisson priors
on the Tucker structure, and assigns Gamma priors to the latent factors.
Recently, Fang et al. [16] develop a Bayesian streaming sparse Tucker
decomposition (BASS-Tucker) to deal with streaming data. BASS-Tucker assigns a
spike-and-slab prior over entries of core tensor and employs an extended
assumed density filtering (ADF) framework for posterior inference.

Similar to CP-based methods, an important task for Tucker decomposition based
method is to choose an appropriate tensor rank. Unfortunately, the problem is
still challenging when dealing with partially observed data corrupted with
noise. Zhao et al. [109] employ hierarchical sparsity-inducing priors to
perform automatic rank determination in their Bayesian tensor decomposition
(BTD) model. Specifically, the observed tensor
$\mathcal{Y}\in\mathbb{R}^{I_{1}\times...\times I_{N}}$ is assumed to follow a
Gaussian distribution with the mean following a Tucker structure:

$\text{vec}(\mathcal{Y})|\\{\boldsymbol{U}^{(n)}\\},\mathcal{G},\tau\sim\mathcal{N}((\bigotimes_{n}\boldsymbol{U}^{(n)}))\text{vec}(\mathcal{G}),\tau^{-1}\boldsymbol{I}),$

where $\\{\boldsymbol{U}^{(n)}\\}$ are latent matrices, $\mathcal{G}$ is the
core tensor, and $\tau$ is the precision. To allow a fully Bayesian treatment,
hierarchical priors are placed over all model parameters. First, a
noninformative Gamma prior is assigned to the precision parameter $\tau$

$\tau\sim\text{Gamma}(a_{0}^{\tau},b_{0}^{\tau}).$

Next, a group sparsity prior is employed over the factor matrices, i.e., each
$\boldsymbol{U}^{(n)}=[\boldsymbol{u}_{1}^{(n)},...,\boldsymbol{u}_{I_{n}}^{(n)}]^{\top}$
($\boldsymbol{u}_{i_{n}}^{(n)}$ are latent vectors) is governed by hyper-
parameters
$\boldsymbol{\lambda}^{(n)}=(\lambda_{1}^{(n)},...,\lambda_{R_{n}}^{(n)})$,
where $\lambda_{r_{n}}^{(n)}$ controls the precision related to $r_{n}$ group
(i.e., $r_{n}$-th column of $\boldsymbol{U}^{(n)}$). Let
$\boldsymbol{\Lambda}^{(n)}=$diag($\boldsymbol{\lambda}^{(n)}$), then the
group sparsity prior is given by

$\boldsymbol{u}_{i_{n}}^{(n)}|\boldsymbol{\lambda}^{(n)}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Lambda}^{(n)^{-1}}),\quad\forall
n,\forall i_{n}.$

The sparsity assumption is also imposed on the core tensor $\mathcal{G}$.
Considering the connection between latent factors and the corresponding entry
of core tensor, the precision parameter for $\mathcal{G}_{r_{1},...,r_{N}}$
can be specified as the product of precisions over $\\{\boldsymbol{u}_{\cdot
r_{n}}^{(n)}\\}_{n=1}^{N}$, which is represented by

$\mathcal{G}_{r_{1}...r_{N}}|\\{\boldsymbol{\lambda}^{(n)}\\},\beta\sim\mathcal{N}(0,(\beta\prod_{n}\lambda_{r_{n}}^{(n)})^{-1}),$

or equivalently

$\text{vec}(\mathcal{G})|\\{\boldsymbol{\lambda}^{(n)}\\},\beta\sim
N(\boldsymbol{0},(\beta\bigotimes_{n}\boldsymbol{\Lambda}^{(n)})^{-1}),$

where $\beta$ is a scalar parameter on which a Gamma prior is placed

$\beta\sim\text{Gamma}(a_{0}^{\beta},b_{0}^{\beta}).$

The hyperprior for $\boldsymbol{\lambda}^{(n)}$ plays a key role for different
sparsity-inducing priors. Two options (student-$t$ and Laplace) are commonly
used to achieve group sparsity:

$\text{Student-}t:\lambda_{r_{n}}^{(n)}\sim\text{Gamma}(a_{0}^{\lambda},b_{0}^{\lambda}),\quad\forall
n,\forall r_{n};$ $\displaystyle\text{Laplace}:$
$\displaystyle~{}\lambda_{r_{n}}^{(n)}\sim\text{IG}(1,\frac{\gamma}{2}),\quad\forall
n,\forall r_{n},$
$\displaystyle~{}\gamma\sim\text{Gamma}(a_{0}^{\gamma},b_{0}^{\gamma}).$

Table 1: Summary of Bayesian tensor decomposition methods.

Name | Decomposition | Rank Specification | Posterior | Data Type
---|---|---|---|---
Structure | Inference
BPTF [99] |  | Pre-specify | Gibbs | Continuous
PLTF [81] |  | Pre-specify | Gibbs | Binary
BGCP [12] |  | Pre-specify | Gibbs | Continuous
PTF [79] |  | Pre-specify | VB | Count
NeuralCP [56] |  | Pre-specify | AEVB | Continuous
MGP-CP [74] |  | Automatically inferred | Gibbs | Continuous/Binary
PGCP [73] | CP | Automatically inferred | Gibbs/EM | Binary/Count
BNBCP [37] | Decomposition | Automatically inferred | Gibbs/VB | Count
ZTP-CP [36] |  | Automatically inferred | Gibbs | Binary
FBCP [108] |  | Automatically inferred | VB | Continuous
BRTF [111] |  | Automatically inferred | VB | Continuous
POST [15] |  | Pre-specify | SVB | Continuous/Binary
BRST [107] |  | Automatically inferred | SVB | Continuous
SBDT [18] |  | Pre-specify | ADF&EP | Continuous/Binary
pTucker [13] |  | Pre-specify | MAP/EM | Continuous
Hayashi et al. [32] | Tucker | Pre-specify | EM | All
BPTD [80] | Decomposition | Pre-specify | Gibbs | Count
BTD [109] | Automatically inferred | VB | Continuous
BASS-Tucker [16] |  | Pre-specify | ADF&EP | Continuous
InfTucker [100] | Nonparametric | Pre-Specify | VEM | Binary/Continuous
Zhe et al. [114] | VEM
DinTucker [113] | VEM
Zhe et al. [115] | VI
SNBTD [70] | ADF&EP
POND [90] | VB
Zhe and Du [112] | VEM
Wang et al. [97] | VI
BCTT [17] | EP
TR-VBI [61] | Tensor Ring | Automatically inferred | VB | Continuous
KFT [38] | Tensor Train | N/A | VI | Continuous
He et al. [33] | Neural | N/A | AEVB | All

ADF: Assume-density filtering [8]. AEVB: Auto-Encoding Variational Bayes [45].
EM: Expectation maximization. EP: Expectation propagation [64]. Gibbs: Markov
chain Monte Carlo (MCMC) with Gibbs sampling. MAP: Maximum a posteriori. SVB:
Steaming variational Bayes. VB: Variational Bayes. VEM: Variational
expectation maximization. VI: Variational Inference. N/A: Not applicable.
Neural: Neural tensor decomposition.

###  5.3 Nonparametric Bayesian Decomposition Methods

In addition to the aforementioned linear models, a class of nonparametric
Bayesian approaches have been developed to capture the potential nonlinear
relationship between tensor entries. One of the pioneering works is InfTucker
proposed by Xu et al. [100]. Generally, InfTucker maps the latent factors into
an infinite feature space and then performs Tucker decomposition with the core
tensor of infinite size. Let $\mathcal{M}\in\mathbb{R}^{m_{1}\times...\times
m_{K}}$ be a tensor following a Tucker structure with a core tensor
$\mathcal{W}$ and latent factors
$\boldsymbol{U}^{(1)},...,\boldsymbol{U}^{(K)}$. One can assign an element-
wise standard Gaussian prior over the core tensor $\mathcal{W}$
(vec$(\mathcal{W})\sim
N(\text{vec}(\mathcal{W});\boldsymbol{0},\boldsymbol{I})$) and marginalize out
$\mathcal{W}$. The marginal distribution for tensor $\mathcal{M}$ is then
given by

$p(\mathcal{M}|\boldsymbol{U}^{(1)},...,\boldsymbol{U}^{(K)})=\mathcal{N}(\text{vec}(\mathcal{M});\boldsymbol{0},\boldsymbol{\Sigma}^{(1)}\otimes...\otimes\boldsymbol{\Sigma}^{(K)})),$
(28)

where
$\boldsymbol{\Sigma}^{(K)}=\boldsymbol{U}^{(K)}\boldsymbol{U}^{(K)^{\top}}$.
Since the goal is to capture nonlinear relationships, each row
$\boldsymbol{u}_{t}^{k}$ of the latent factors $\boldsymbol{U}^{(k)}$ is
replaced by a nonlinear feature mapping $\phi(\boldsymbol{u}_{t}^{k})$. Then a
nonlinear covariance matrix
$\boldsymbol{\Sigma}^{(k)}=k(\boldsymbol{U}^{(k)},\boldsymbol{U}^{(k)})$ can
be obtained, where $k(\cdot,\cdot)$ is a nonlinear covariance kernel function.
In InfTucker [100], $k(\cdot,\cdot)$ is chosen as the radial basis function
kernel. After feature mapping, the core tensor $\mathcal{W}$ has the size of
the mapped feature vector $\boldsymbol{u}_{t}^{k}$ on mode $k$, which is
potentially infinity. Because the covariance of vec($\mathcal{M}$) is a
function of the latent factors
$\mathcal{U}=\\{\boldsymbol{U}^{(1)},...,\boldsymbol{U}^{(K)}\\}$, equation
(28) actually defines a Gaussian process (GP) on tensor entries, where the
input is based on the corresponding latent factors $\mathcal{U}$. To encourage
sparse estimation, element-wise Laplace priors are assigned on $\mathcal{U}$:

$\boldsymbol{u}_{i}^{(k)}\sim\mathcal{L}(\lambda)\propto\text{exp}(-\lambda\|\boldsymbol{u}_{i}^{(k)}\|_{1}).$
(29)

Finally, the observed tensor $\mathcal{Y}$ is sampled from a noisy model
$p(\mathcal{Y}|\mathcal{M})$, of which the form depends on the data type of
$\mathcal{Y}$. The joint distribution is then given by

$p(\mathcal{Y},\mathcal{M},\mathcal{U})=p(\mathcal{U})p(\mathcal{M}|\mathcal{U})p(\mathcal{Y}|\mathcal{M}),$

where $p(\mathcal{U})$ is given by (29), and $p(\mathcal{M}|\mathcal{U})$ is
given by (28) with
$\boldsymbol{\Sigma}^{(k)}=k(\boldsymbol{U}^{(k)},\boldsymbol{U}^{(k)})$.

Under a similar modeling framework, Zhe et al. [114] make two modifications to
InfTucker. One is to assign a Dirichlet process mixture (DPM) prior [3] over
the latent factors that allows an undetermined number of latent clusters. The
other is to utilize a local GP assumption instead of a global GP when
generating the observed array given the latent factors, which enables fast
computation over subarrays. Specifically, the local GP-based construction is
realized by firstly breaking the whole array $\mathcal{Y}$ into smaller
subarrays $\\{\mathcal{Y}_{1},..,\mathcal{Y}_{N}\\}$. Then for each subarray
$\mathcal{Y}_{n}$, a latent real-valued subarray $\mathcal{M}_{n}$ is
generated by a local GP based on the corresponding subset of latent factors
$\mathcal{U}_{n}=\\{\boldsymbol{U}_{n}^{(1)},...,\boldsymbol{U}_{n}^{(K)}\\}$
and the noisy observation $\mathcal{Y}_{n}$ is sampled according to
$\mathcal{M}_{n}$,

$\begin{split}p(\mathcal{Y}_{n},\mathcal{M}_{n}|\mathcal{U})&\\!=\\!p(\mathcal{M}_{n}|\mathcal{U}_{n})p(\mathcal{Y}_{n}|\mathcal{M}_{n})\\!=\\!\mathcal{N}(\text{vec}(\mathcal{M}_{n});\boldsymbol{0},\boldsymbol{\Sigma}_{n}^{(1)}\otimes...\otimes\boldsymbol{\Sigma}_{n}^{(K)})p(\mathcal{Y}_{n}|\mathcal{M}_{n}),\end{split}$

where
$\boldsymbol{\Sigma}_{n}^{(k)}=k(\boldsymbol{U}_{n}^{(k)},\boldsymbol{U}_{n}^{(k)})$
is the $k$-th mode covariance matrix over the sub-factors $\mathcal{U}_{n}$.

Likewise, DinTucker [113] consider a local GP assumption and sample each of
the subarrays $\\{\mathcal{Y}_{1},...,\mathcal{Y}_{n}\\}$ from a GP based on
the latent factors
$\tilde{\mathcal{U}}_{n}=\\{\tilde{\boldsymbol{U}}_{n}^{(1)},...,\tilde{\boldsymbol{U}}_{n}^{(K)}\\}$.
Different from Zhe et al. [114], in DinTucker these latent factors are then
tied to a set of common latent factors
$\mathcal{U}=\\{\boldsymbol{U}^{(1)},...,\boldsymbol{U}^{(K)}\\}$ via a prior
distribution

$p(\tilde{\mathcal{U}}_{n}|\mathcal{U})=\prod_{k=1}^{K}\mathcal{N}(\text{vec}(\tilde{\boldsymbol{U}}_{n}^{(k)})|\text{vec}(\boldsymbol{U}^{(k)}),\lambda\boldsymbol{I}),$

where $\lambda$ is the variance parameter that controls the similarity between
$\mathcal{U}$ and $\tilde{\mathcal{U}}_{n}$. Furthermore, DinTucker divides
each subarray $\mathcal{Y}_{n}$ into $T_{n}$ smaller subarrays
$\mathcal{Y}_{n}=\\{\mathcal{Y}_{n1},...,\mathcal{Y}_{nT_{n}}\\}$ that share
the same latent factors $\\{\tilde{\mathcal{U}}_{n}\\}$, and the joint
probability is given by

$\begin{split}p(\mathcal{U},\\{\tilde{\mathcal{U}}_{n},\mathcal{M}_{n},&\mathcal{Y}_{n}\\}_{n=1}^{N})=\prod_{n=1}^{N}p(\tilde{\mathcal{U}}_{n}|\mathcal{U})\prod_{t=1}^{T_{n}}p(\mathcal{M}_{nt}|\tilde{\mathcal{U}}_{n})p(\mathcal{Y}_{nt}|\mathcal{M}_{nt}),\end{split}$

where $\mathcal{M}_{nt}$ is a latent subarray, and
$\mathcal{M}_{n}=\\{\mathcal{M}_{nt}\\}_{t=1}^{T_{n}}$. The local terms
require less memory and has a faster processing time than the global term.
More importantly, the additive nature of these local terms in the log domain
enables distributed inference, which is then realized through the MapReduce
system.

While Zhe et al. [114] and DinTucker [113] improve the scalability of their
GP-based approaches through modeling the subtensors, their methods can still
run into challenges when the sparsity level is very high in observed tensors.
To address this issue, a class of methods that does not rely on the Kronecker-
product structure in the variance (28) are proposed based on the idea of
selecting an arbitrary subset of tensor entries for training. Assume that the
decomposition is to be performed on a sparsely observed tensor
$\mathcal{Y}\in\mathbb{R}^{d_{1}\times...\times d_{K}}$. For each tensor entry
$\boldsymbol{i}=(i_{1},...,i_{K})$, Zhe et al. [115] first construct an input
$\boldsymbol{x_{i}}$ by concatenating the corresponding latent factors from
all the modes:
$\boldsymbol{x_{i}}=[\boldsymbol{u}_{i_{1}}^{(1)},...,\boldsymbol{u}_{i_{K}}^{(K)}]$,
where $\boldsymbol{u}_{i_{k}}^{(k)}$ is the $i_{k}$-th row in the latent
factor matrix $\boldsymbol{U}^{(k)}$ for mode $k$. Then each
$\boldsymbol{x}_{\boldsymbol{i}}$ is transformed to a scalar
$m_{\boldsymbol{i}}$ through an underlying function
$f:\mathbb{R}^{\sum_{j=1}^{K}d_{j}}\to\mathbb{R}$ such that
$m_{\boldsymbol{i}}=f(\boldsymbol{x_{i}})=f([\boldsymbol{u}_{i_{1}}^{(1)},...,\boldsymbol{u}_{i_{K}}^{(K)}])$.
After that, a GP prior is assigned over $f$ to learn the unknown function: for
any set of tensor entries $S=\\{\boldsymbol{i}_{1},...,\boldsymbol{i}_{N}\\}$,
the function values
$\boldsymbol{f}_{S}=\\{f(\boldsymbol{x}_{\boldsymbol{i}_{1}}),...,f(\boldsymbol{x}_{\boldsymbol{i}_{N}})\\}$
are distributed according to a multivariate Gaussian distribution with mean
$\boldsymbol{0}$ and covariance determined by
$\boldsymbol{X}_{S}=\\{\boldsymbol{x}_{\boldsymbol{i}_{1}},...,\boldsymbol{x}_{\boldsymbol{i}_{N}}\\}$:

$p(\boldsymbol{f}_{S}|\mathcal{U})=\mathcal{N}(\boldsymbol{f}_{S}|\boldsymbol{0},k(\boldsymbol{X}_{S},\boldsymbol{X}_{S})),$
(30)

where $\mathcal{U}$ is the latent factor, and $k(\cdot,\cdot)$ is a nonlinear
covariance kernel. Note that this method is equivalent to InfTucker [100] if
all entries are selected and a Kronecker-product structure in the full
covariance is applied. A standard normal prior is assigned over the latent
factors, and the observed entries
$\boldsymbol{y}=[y_{\boldsymbol{i}_{1}},...,y_{\boldsymbol{i}_{N}}]$ are
sampled from a noise model $p(\boldsymbol{y}|\boldsymbol{m})$, where
$p(\cdot)$ is selected based on the data type.

Following the sparse GP framework (30), Pan et al. [70] propose the Streaming
Nonlinear Bayesian Tensor Decomposition (SNBTD) that performs fast posterior
updates upon receiving new tensor entries. The model is augmented with feature
weights to incorporate a linear structure, and the assumed-density-filtering
(ADF) framework is extended to perform reliable streaming inference. Also
based on (30), Tillinghast et al. [90] utilize convolutional neural networks
to construct a deep kernel $k(\cdot,\cdot)$ for GP modeling, which is more
powerful in estimating arbitrarily complicated relationships in data compared
to the methods based on shallow kernel functions (e.g., RBF kernel).

Furthermore, the tensor data are sometimes observed with temporal information;
and various approaches have been proposed to preserve the accurate timestamps
and take full advantage of the temporal information. Among these methods, Zhe
and Du [112] and Wang et al. [97] perform decomposition based on event-tensors
to capture complete temporal information, and Fang et al. [17] model the core
tensor as a time-varying function, where GP prior is placed to estimate
different types of temporal dynamics.

##  6 Bayesian Methods in Tensor Regression

Similar to frequentist tensor regression methods discussed in Section 4,
Bayesian tensor regression methods can be categorized into Bayesian tensor
predictor regression and Bayesian tensor response regression. We discuss these
two classes of methods in Section 6.1 and 6.2, and their theoretical
properties in Section 6.3. We also review posterior computing in Section 6.4.
A summary of the methods discussed in this section is given in Table 2.

Table 2: Summary of Bayesian tensor regression methods.

Name | Predictor | Response | Tensor | Algorithm
---|---|---|---|---
Type | Type | Structure
Suzuki [88] | Tensor | Scalar | CP | Gibbs
BTR [26] | Tensor+Vector | Scalar | CP | Gibbs
Zhao et al. [110] | Tensor | Scalar | Nonparametric | MAP
OLGP [35] | Tensor | Scalar | Nonparametric | OLGP
AMNR [39] | Tensor | Scalar | Nonparametric | MC
Yang and Dunsun [101] | Vector (Categorical) | Scalar (Categorical) | Tucker | Gibbs
CATCH [69] | Tensor+Vector | Scalar (Categorical) | Tucker | MLE
BTRR [27] | Vector | Tensor | CP | Gibbs
Spencer et al. [84, 85] | Vector | Tensor | CP | Gibbs
SGTM [23] | Vector | Symmetric Tensor | CP | Gibbs
BSTN [49] | Vector | Tensor | Other | Gibbs
SGPRN [52] | Matrix | Tensor | Nonparametric | VI
MLTR [34] | Tensor | Tensor | Tucker | Gibbs
ART [7] | Tensor | Tensor | CP | Gibbs

Gibbs: MCMC with Gibbs sampling. MAP: Maximum a posteriori. MC: Monte Carlo
Method. MLE: Maximum likelihood estimator. OLGP: Online local Gaussian process
[68, 95]. VI: Variational Inference.

###  6.1 Bayesian Tensor Predictor Regression

In recent years, Bayesian tensor predictor regression models have gained an
increasing attention. In 2015, Suzuki [88] developed a Bayesian framework
based on the basic tensor linear regression model

$Y_{i}=\langle\mathcal{W},\mathcal{X}_{i}\rangle+\epsilon_{i},$ (31)

where $Y_{i}\in\mathbb{R}$ is a univariate response,
$\mathcal{X}_{i}\in\mathbb{R}^{M_{1}\times\cdots\times M_{K}}$ is a tensor-
valued predictor, $\mathcal{W}\in\mathbb{R}^{M_{1}\times\cdots\times M_{K}}$
is the coefficient tensor, and $\langle\cdot,\cdot\rangle$ is the tensor inner
product (2). The error terms $\epsilon_{i}$’s are assumed i.i.d. following a
normal distribution $\mathcal{N}(0,\sigma^{2})$. To achieve parsimony in free
parameters, a rank-$r$ CP structure (5) is imposed on the coefficient tensor
$\mathcal{W}$:

$\mathcal{W}=[\\![\boldsymbol{U}^{(1)},...,\boldsymbol{U}^{(K)}]\\!],$

where $\boldsymbol{U}^{(k)}\in\mathbb{R}^{r\times M_{K}}$ ($k=1,2,...,K$) are
latent factors. To complete model specification, a Gaussian prior is placed on
the latent matrices:

$\pi(\boldsymbol{U}^{(1)},...,\boldsymbol{U}^{(K)}|r)\propto\exp\Big{\\{}-\frac{r}{2\sigma^{2}_{p}}\sum_{k=1}^{K}\text{Tr}[\boldsymbol{U}^{(k)^{\top}}\boldsymbol{U}^{(k)}]\Big{\\}},$

and a prior on the rank $r$:

$\pi(r)=\frac{1}{N_{\xi}}\xi^{r(M_{1}+\cdots+M_{K})},$

where $0<\xi<1$ is a positive real number, and $N_{\xi}$ is the normalizing
constant.

In order to adjust for other covariates in the model and accommodate various
data types of the response variable, Guhaniyogi et al. [26] propose a Bayesian
method based on the generalized tensor predictor regression model (12). Given
a scalar response $y$, vectorized predictors $\boldsymbol{z}\in\mathbb{R}^{p}$
and tensor predictor $\mathcal{X}\in\mathbb{R}^{p_{1}\times
p_{2}\times...\times p_{D}}$, the regression model is given by

$y\sim
f(\alpha+\boldsymbol{z}^{\top}\boldsymbol{\gamma}+\langle\mathcal{X},\mathcal{B}\rangle,\sigma),$
(32)

where $f(\mu,\sigma)$ is a family of distributions with location $\mu$ and
scale $\sigma$, $\boldsymbol{\gamma}\in\mathbb{R}^{p}$ are coefficients for
predictors $\boldsymbol{z}$, $\mathcal{B}\in\mathbb{R}^{p_{1}\times
p_{2}\times...\times p_{D}}$ is the coefficient tensor, and
$\langle\cdot,\cdot\rangle$ is the tensor inner product (2). A CP structure is
imposed on the tensor coefficient $\mathcal{B}$:

$\mathcal{B}=\sum_{r=1}^{R}\boldsymbol{\beta}_{1}^{(r)}\circ\cdots\circ\boldsymbol{\beta}_{D}^{(r)}.$

Under the Bayesian framework, Guhaniyogi et al. [26] propose a multiway
Dirichlet generalized double Pareto (M-DGDP) prior over the latent factors
$\boldsymbol{\beta}_{j}^{(r)}$. This prior promotes the joint shrinkage on the
global and local component parameters, as well as accommodates dimension
reduction by favoring low-rank decompositions. Specifically, M-GDGP prior
first assigns a multivariate Gaussian prior on $\boldsymbol{\beta}_{j}^{(r)}$:

$\boldsymbol{\beta}_{j}^{(r)}\sim\mathcal{N}(\boldsymbol{0},(\phi_{r}\tau)\boldsymbol{W}_{jr}),~{}j=1,\ldots,D.$
(33)

The shrinkage across components is induced in an exchangeable way, with global
scale $\tau\sim\text{Gamma}(a_{\tau},b_{\tau})$ adjusted in each component by
$\phi_{r}$ for $r=1,2,...,R$, where
$\Phi=(\phi_{1},...,\phi_{R})\sim\text{Dirichlet}(\alpha_{1},...,\alpha_{R})$
that encourages shrinkage towards lower ranks in the CP structure. In
addition, $\boldsymbol{W}_{jr}=\text{diag}(w_{jr,1},\cdots,w_{jr,p_{j}})$,
$j=1,2,...,D$ and $r=1,2,...,R$ are scale parameters for each component, where
a hierarchical prior is placed,

$w_{jr,k}\sim\text{Exp}(\lambda^{2}_{jr}/2),\quad\lambda_{jr}\sim\text{Gamma}(a_{\lambda},b_{\lambda}).$
(34)

In the M-GDGP prior, flexibility in estimating
$\mathcal{B}_{r}=\\{\boldsymbol{\beta}_{j}^{(r)};1\leq j\leq D\\}$ is achieved
by modeling within-margin heterogeneity via element-specific scaling
$w_{jr,k}$. The common rate parameter $\lambda_{jr}$ shares information
between margin elements, hence encouraging shrinkage at the local scale.

Besides linear models, a class of Gaussian process (GP) based nonparametric
approaches have been proposed to model nonlinear relationships in the tensor-
valued predictors. Given a dataset of $N$ paired observations
$\mathcal{D}=\\{(\mathcal{X}_{n},y_{n})|n=1,2,...,N\\}$, Zhao et al. [110]
aggregate all $N$ tensor inputs $\mathcal{X}_{n}~{}(n=1,2,...,N)$ into a
design tensor $\mathcal{X}\in\mathbb{R}^{N\times I_{1}\times\cdots\times
I_{M}}$, and collect the responses in the vector form
$\boldsymbol{y}=[y_{1},...,y_{N}]^{\top}$. The distribution of response vector
can be factored over the observations as

$\boldsymbol{y}\sim\prod_{n=1}^{N}\mathcal{N}(y_{n}|f(\mathcal{X}_{n}),\sigma^{2}).$
(35)

Here $f(\cdot)$ is a latent function on which a GP prior is placed

$f(\mathcal{X})\sim\text{GP}(m(\mathcal{X}),k(\mathcal{X},\mathcal{X}^{\prime})|\boldsymbol{\theta}),$
(36)

where $k(\mathcal{X},\mathcal{X}^{\prime})$ is the covariance function
(kernel), $\boldsymbol{\theta}$ is the associated hyperparameter vector, and
$m(\mathcal{X})$ is the mean function which is set to be zero in [110]. The
authors further propose to use the following product kernel in (36):

$k(\mathcal{X},\mathcal{X}^{\prime})=\alpha^{2}\prod_{d=1}^{D}\exp(\frac{KL(p(\boldsymbol{x}|\Omega_{d}^{\mathcal{X}})~{}\|~{}q(\boldsymbol{x}^{\prime}|\Omega_{d}^{\mathcal{X}^{\prime}}))}{-2\beta_{d}^{2}}),$
(37)

where $\alpha$ is a magnitude hyperparameter, and $\beta_{d}$ denotes the
$d$-mode length-scale hyper-parameter. The distributions $p$ and $q$ in the
Kullback-Leibler (KL) divergence are characterized by the hyper-parameters
$\Omega_{d}$, which can be estimated from the $d$-mode unfolding matrix
$\boldsymbol{X}_{d}$ of tensor $\mathcal{X}$ by treating each
$\boldsymbol{X}_{d}$ as a generative model with $I_{d}$ variables and
$I_{1}\times\cdots\times I_{d-1}\times I_{d+1}\times\cdots\times I_{D}$
observations. Given the complete prior construction, the hyperparameters
$\boldsymbol{\theta}=\\{\alpha,\beta_{d}|d=1,2,...,D\\}$ and $\sigma$ are then
estimated by maximum a posteriori (MAP). While the computational complexity of
GP-based methods is usually excessive, Hou et al. [35] take advantage of
online local Gaussian Process (OLGP) to present a computationally-efficient
approach for the nonparametric model in (35)-(37).

To further mitigate the burden of high-dimensionality, Imaizumi and Hayashi
[39] propose an additive-multiplicative nonparametric regression (AMNR) that
concurrently decompose the functional space and the input space, and thus is
referred to as a doubly decomposing nonparametric tensor regression method.

Denote a Sobolev space by $\mathcal{W}^{\beta}(\mathcal{X})$, which is a space
of $\beta$ times differentiable functions with support $\mathcal{X}$. Let
$\mathcal{X}=\bigotimes_{k}\boldsymbol{x}_{k}:=\boldsymbol{x}_{1}\otimes\cdots\otimes\boldsymbol{x}_{K}$
be a rank-one tensor denoted by the outer product of vectors
$\boldsymbol{x}_{k}\in\mathcal{X}^{(k)}$ ($\otimes$ is the outer product). Let
$f\in\mathcal{W}^{\beta}(\bigotimes_{k}\mathcal{X}^{(k)})$ be a function on a
rank one tensor, and for any $f$ we can construct
$\tilde{f}(\boldsymbol{x}_{1},...,\boldsymbol{x}_{K})\in\mathcal{W}^{\beta}(\mathcal{X}^{(1)}\times\cdots\times\mathcal{X}^{(k)})$
such that
$\tilde{f}(\boldsymbol{x}_{1},...,\boldsymbol{x}_{K})=f(\mathcal{X})$ using
function decomposition as $\tilde{f}=f\circ h$ with
$h:(\boldsymbol{x}_{1},...,\boldsymbol{x}_{K})\to\bigotimes_{k}\boldsymbol{x}_{k}$.
Then $f$ can be decomposed into a set of local functions
$\\{f_{m}^{k}\in\mathcal{W}^{\beta}(\mathcal{X}^{(k)})\\}_{m}$ following [29]:

$f(\mathcal{X})=\tilde{f}(\boldsymbol{x}_{1},...,\boldsymbol{x}_{K})=\sum_{m=1}^{M}\prod_{k=1}^{K}f_{m}^{(k)}(\boldsymbol{x}_{k}),$
(38)

where $M$ represents the complexity of $f$ (i.e., the “rank” of the model).

Based on (38), for a rank-$R$ tensor $\mathcal{X}$, Imaizumi and Hayashi [39]
define the AMNR function as:

$f^{AMNR}(\mathcal{X}):=\sum_{m=1}^{M}\sum_{r=1}^{R}\lambda_{r}\prod_{k=1}^{K}f_{m}^{(k)}(\boldsymbol{x}_{r}^{(k)}),$
(39)

which is achieved by first writing a rank-$R$ tensor as the sum of $R$ rank-
one tensors, then for each rank-one tensor decomposing the function into a set
of local functions. Under the Bayesian framework, a GP prior is assigned to
the local functions $f_{m}^{(k)}$, and the Gaussian distribution (35) is
utilized to associate the scalar response $Y_{i}$ with the function
$f^{AMNR}(\mathcal{X}_{i})$.

While the previous studies mainly deal with the regression problems where the
response is a continuous variable, the probabilistic methods are also apply to
categorical-response regression with tensor-valued predictors, i.e., the
tensor classification problems. For example, Pan et al. [69] propose a
covariate-adjusted tensor classification model (CATCH), which jointly models
the relationship among the covariates, tensor predictors, and categorical
responses. Given categorical response $Y\in\\{1,2,...,K\\}$, vector of
adjusting covariates $\boldsymbol{U}\in\mathbb{R}^{q}$, and tensor-variate
predictors $\mathcal{X}\in\mathbb{R}^{p_{1}\times\cdots\times p_{M}}$, the
CATCH model is proposed as

$\displaystyle\boldsymbol{U}|(Y=k)\sim
N(\boldsymbol{\Phi}_{k},\boldsymbol{\Psi})$ (40)
$\displaystyle\mathcal{X}|(\boldsymbol{U}=\boldsymbol{u},Y=k)\sim\text{TN}(\boldsymbol{\mu}_{k}+\boldsymbol{\alpha}\bar{\times}_{(M+1)}\boldsymbol{u};\boldsymbol{\Sigma}_{1},...,\boldsymbol{\Sigma}_{M}),$
(41)

where
$\boldsymbol{\Phi}_{k}\in\mathbb{R}^{q},\boldsymbol{\Psi}\in\mathbb{R}^{q\times
q}$ ($\boldsymbol{\Psi}>0$ and symmetric),
$\boldsymbol{\alpha}\in\mathbb{R}^{p_{1}\times...\times p_{M}\times
q},\boldsymbol{\mu}_{k}\in\mathbb{R}^{p_{1}\times...\times p_{M}}$, and
$\boldsymbol{\Sigma}_{m}\in\mathbb{R}^{p_{m}\times p_{m}}$
($\boldsymbol{\Sigma}_{m}>0$ and symmetric, $m=1,...,M$). Here TN$(\cdot)$ is
the tensor normal distribution, and $\bar{\times}_{(M+1)}$ is the $(M+1)$-mode
tensor vector product.

In equation (40), it is assumed that $\\{Y,\boldsymbol{U}\\}$ follows a
classical LDA model, where $\boldsymbol{\Phi}_{k}$ is the mean of
$\boldsymbol{U}$ within class $k$ and $\boldsymbol{\Psi}$ is the common within
class covariance of $\boldsymbol{U}$. Similarly, in equation (41) a common
within class covariance structure of $\mathcal{X}$ is assumed (denoted by
$\boldsymbol{\Sigma}_{m},m=1,2,...,M$), which does not depend on $Y$ after
adjusting for the covariates $\boldsymbol{U}$. The tensor coefficient
$\boldsymbol{\alpha}$ characterizes the linear dependence of tensor predictor
$\mathcal{X}$ on the covariates $\boldsymbol{U}$, and $\boldsymbol{\mu}_{k}$
is the covariate-adjusted within-class mean of $\mathcal{X}$ in class $k$.

While the goal is to predict $Y$ given $\\{\boldsymbol{U},\mathcal{X}\\}$,
based on the Bayes’ rule the optimal classifier under the CATCH model is
derived by maximizing the posterior probability

$\begin{split}\hat{Y}&=\arg\max_{k=1,2,...,K}P(Y=k|\mathcal{X}=\boldsymbol{x},\boldsymbol{U}=\boldsymbol{u})\\\
&=\arg\max_{k=1,2,...,K}\pi_{k}f_{k}(\boldsymbol{x},\boldsymbol{u}),\end{split}$
(42)

where $\pi_{k}=P(Y=k)$ and $f_{k}(\boldsymbol{x},\boldsymbol{u})$ is the joint
density function of $\mathcal{X}$ and $\boldsymbol{U}$ conditional on $Y=k$.
Combining (40) and (41), equation (42) is transformed into

$\hat{Y}=\arg\max_{k=1,2,...,K}\\{a_{k}+\boldsymbol{\gamma}_{k}^{\top}\boldsymbol{U}+\langle\mathcal{B}_{k},\mathcal{X}-\boldsymbol{\alpha}\bar{\times}_{(M+1)}\boldsymbol{U}\rangle\\},$

where
$\boldsymbol{\gamma}_{k}=\boldsymbol{\Psi}^{-1}(\boldsymbol{\Phi}_{k}-\boldsymbol{\Phi}_{1}),\mathcal{B}_{k}=[\\![\boldsymbol{\mu}_{k}-\boldsymbol{\mu}_{1};\boldsymbol{\Sigma}_{1}^{-1},...,\boldsymbol{\Sigma}_{M}^{-1}]\\!]$
following a Tucker structure with core tensor
$\boldsymbol{\mu}_{k}-\boldsymbol{\mu}_{1}$ and latent matrices
$\boldsymbol{\Sigma}_{1}^{-1},...,\boldsymbol{\Sigma}_{M}^{-1}$, and
$a_{k}=\log(\pi_{k}/\pi_{1})-\frac{1}{2}\boldsymbol{\gamma}_{k}^{\top}(\boldsymbol{\Phi}_{k}+\boldsymbol{\Phi}_{1})-\langle\mathcal{B}_{k},\frac{1}{2}(\boldsymbol{\mu}_{k}+\boldsymbol{\mu}_{1})\rangle$
is a scalar that does not involve $\mathcal{X}$ or $\boldsymbol{U}$.

Given i.i.d samples
$\\{Y^{i},\boldsymbol{U}^{i},\mathcal{X}^{i}\\}_{i=1}^{n}$, the parameters
$\\{\pi_{k},\boldsymbol{\Phi}_{k},\boldsymbol{\gamma}_{k},\boldsymbol{\mu}_{k},\mathcal{B}_{k}\\}_{k=1}^{K}$
and $\\{\boldsymbol{\Sigma}_{m}\\}_{m=1}^{M}$ can be estimated to build an
accurate classifier based on the data. It is worth noting that the estimation
of $\mathcal{B}_{k}$ are penalized in order to facilitate sparsity and achieve
parsimony.

Though not modeling tensorized predictors, Yang and Dunsun [101] employ tensor
methods to deal with classification problems with categorical predictors.
Specifically, [101] develops a framework for nonparametric Bayesian
classification through performing decomposition on the tensor transformed from
the conditional probability

$P(Y=y|X_{1}=x_{1},...,X_{p}=x_{p}),$

with a categorical response $Y\in\\{1,2,...,d_{0}\\}$ and a vector of $p$
categorical predictors $\boldsymbol{X}=(X_{1},X_{2},...,X_{p})^{\top}$. The
conditional probability can be structured as a $d_{0}\times
d_{1}\times\cdots\times d_{p}$-dimensional tensor, where $d_{j}$
$(j=1,2,...,p)$ denotes the number of levels of the $j$-th categorical
predictor $X_{j}$. The tensor is referred to as conditional probability
tensor, and the set of all conditional probability tensors is denoted by
$\mathcal{P}_{d_{1},...,d_{p}}(d_{0})$. Therefore,
$\mathcal{P}\in\mathcal{P}_{d_{1},...,d_{p}}(d_{0})$ implies

$\displaystyle\mathcal{P}_{y,x_{1},...,x_{p}}\geq 0\quad\text{for
every}~{}y,x_{1},...,x_{p};$
$\displaystyle\sum_{y=1}^{d_{0}}\mathcal{P}_{y,x_{1},...,x_{p}}=1\quad\text{for
every}~{}x_{1},...,x_{p}.$

Now all the conditional probabilities are described by an entry of the
conditional probability tensor, and thus the classification problem is
converted into a tensor decomposition problem. Additionally, Yang and Dunsun
[101] prove that every conditional probability tensor
$\mathcal{P}\in\mathcal{P}_{d_{1},...,d_{p}}(d_{0})$ can be expressed by a
Tucker structure

$\begin{split}\mathcal{P}_{y,x_{1},...,x_{p}}&=P(y|x_{1},...,x_{p})=\sum_{h_{1}=1}^{k_{1}}\cdots\sum_{h_{p}=1}^{k_{p}}\lambda_{h_{1}h_{2}...h_{p}}(y)\prod_{j=1}^{p}\pi_{h_{j}}^{(j)}(x_{j}),\end{split}$

with all positive parameters satisfying

$\begin{split}&\sum_{c=1}^{d_{0}}\lambda_{h_{1}h_{2}...h_{p}}(c)=1,\quad\text{for
every}~{}h_{1},h_{2},...,h_{p},\\\
&\sum_{h=1}^{k_{j}}\pi_{h}^{(j)}(x_{j})=1,\quad\text{for every pair of
}~{}j,x_{j}.\end{split}$

The inference on the Tucker coefficients is done under the Bayesian framework.
Specifically, independent Dirichlet priors are assigned to the parameters
$\boldsymbol{\Lambda}=\\{\lambda_{h_{1},...,h_{p}}(c),c=1,2,...,d_{0}\\}$ and
$\boldsymbol{\pi}=\\{\pi_{h_{j}}^{(j)}(x_{j}),h_{j}=1,2,...,k_{j}\\}$
($x_{j}=1,2,...,d_{j},h_{j}=1,2,...,k_{j},j=1,2,...,p$):

$\displaystyle\bigg{\\{}\lambda_{h_{1},...,h_{p}}(1),...,\lambda_{h_{1},...,h_{p}}(d_{0})\bigg{\\}}\sim\text{Dirichlet}(\frac{1}{d_{0}},...,\frac{1}{d_{0}}),$
$\displaystyle\bigg{\\{}\pi_{1}^{(j)}(x_{j}),...,\pi_{k_{j}}^{(j)}(x_{j})\bigg{\\}}\sim\text{Dirichlet}(\frac{1}{k_{j}},...,\frac{1}{k_{j}}),~{}j=1,...,p.$

These priors can impose non-negative and sum-to-one constraints naturally and
lead to conditional conjugacy in posterior computation. Besides, [101] also
assign priors to the hyper-parameters in the Dirichlet priors to promote a
fully Bayesian treatment. These priors place most of the probability on a few
elements to induce sparsity in these vectors.

###  6.2 Bayesian Tensor Response Regression

Guhaniyogi and Spencer [27] propose a Bayesian regression model with tensor
response and scalar predictors. Let $\mathcal{Y}_{t}\in\mathbb{R}^{p_{1}\times
p_{2}\times...\times p_{D}}$ be a tensor valued response, and
$\boldsymbol{x}_{t}=(x_{1,t},...,x_{m,t})\in\mathcal{X}\subset\mathbb{R}^{m}$
be an $m$-dimensional vector predictor measured at time $t$. Assuming that
both response $\mathcal{Y}_{t}$ and predictors $\boldsymbol{x}_{t}$ are
centered around their respective means, the proposed regression model for
$\mathcal{Y}_{t}$ on $\boldsymbol{x}_{t}$ is given by

$\mathcal{Y}_{t}=\boldsymbol{\Gamma}_{1}x_{1,t}+\cdots+\boldsymbol{\Gamma}_{m}x_{m,t}+\mathcal{E}_{t},\quad
i=1,2,...,n,$ (43)

where $\boldsymbol{\Gamma}_{k}\in\mathbb{R}^{p_{1}\times p_{2}\times...\times
p_{D}},k=1,2,...,m$ is the tensor coefficient corresponding to predictor
$x_{k,t}$, and $\mathcal{E}_{t}\in\mathbb{R}^{p_{1}\times p_{2}\times...\times
p_{D}}$ represents the error tensor. To account for the temporal correlation
of the response tensor, the error tensor $\mathcal{E}_{t}$ is assumed to
follow a component-wise AR(1) structure across $t$:
vec$(\mathcal{E}_{t})=\kappa\text{vec}(\mathcal{E}_{t-1})+\text{vec}(\boldsymbol{\eta}_{t})$,
where $\kappa\in(-1,1)$ is the correlation coefficient, and
$\boldsymbol{\eta}_{t}\in\mathbb{R}^{p_{1}\times p_{2}\times...\times p_{D}}$
is a random tensor, with each entry following a Gaussian distribution
$\mathcal{N}(0,\sigma^{2}/(1-\kappa^{2}))$.

Next, a CP structure is imposed on each $\boldsymbol{\Gamma}_{k}$ to reduce
the dimensionality of coefficient tensors, i.e.,
$\boldsymbol{\Gamma}_{k}=\sum_{r=1}^{R}\boldsymbol{\gamma}_{1,k}^{(r)}\circ\cdots\circ\boldsymbol{\gamma}_{D,k}^{(r)}$.
Although Guhaniyogi’s previously proposed M-DGDP prior (33)(34) over the
latent factors $\boldsymbol{\gamma}_{j,k}^{(r)}$ can promote global and local
sparsity, Guhaniyogi and Spencer [27] claim that the straightforward
application of M-DGDP prior leads to inaccurate estimation due to less
desirable tail behavior of the coefficient distributions. Instead, a multiway
stick breaking shrinkage prior (M-SB) is assigned to
$\boldsymbol{\gamma}_{j,k}^{(r)}$, where the main difference compared to
M-DGDP prior is how shrinkage is achieved across ranks. The construction of
the M-SB prior is given as follows. Let
$\boldsymbol{W}_{jr,k}=\text{diag}(w_{jr,k,1},...,w_{jr,k,p_{d}})$. Then

$\boldsymbol{\gamma}_{j,k}^{(r)}\sim\mathcal{N}(0,\tau_{r,k}\boldsymbol{W}_{jr,k}).$

Further set $\tau_{r,k}=\phi_{r,k}\tau_{k}$ to be scaling specific to rank $r$
($r=1,...,R$). Then effective shrinkage across ranks is achieved by adopting a
stick breaking construction for the rank-specific parameter $\phi_{r,k}$:

$\displaystyle\phi_{r,k}=\xi_{r,k}\prod_{l=1}^{r-1}(1-\xi_{l,k}),\quad
r=1,...,R-1,$ $\displaystyle\phi_{R,k}=\prod_{l=1}^{R-1}(1-\xi_{l,k}),$

where $\xi_{r,k}\sim_{iid}\text{Beta}(1,\alpha_{k}).$ The Bayesian setting is
then completed by specifying

$\displaystyle\tau_{k}\sim\text{InvGamma}(a_{\tau},b_{\tau}),~{}~{}w_{jr,k,i}\sim\text{Exp}(\lambda_{jr,k}^{2}/2),~{}~{}\lambda_{jr,k}\sim\text{Gamma}(a_{\lambda},b_{\lambda}),$

where the hierarchical prior of $w_{jr,k,i}$ allows the local scale parameters
$\boldsymbol{W}_{jr,k}$ to achieve margin level shrinkage.

Based on the regression function (43), Spencer et al. [84, 85] develop an
additive mixed effect model that simultaneously measures the activation due to
stimulus at voxels in the $g$th brain region and connectivity among $G$ brain
regions. Let $\mathcal{Y}_{i,g,t}\in\mathbb{R}^{p_{1,g}\times\cdots\times
p_{D,g}}$ be the tensor of observed fMRI data in brain region $g$ for the
$i$-th subject at the $t$-th time point, and
$x_{1,i,t},...,x_{m,i,t}\in\mathbb{R}$ be activation-related predictors, the
regression function is given by

$\mathcal{Y}_{i,g,t}=\boldsymbol{\Gamma}_{1,g}x_{1,i,t}+\cdots\boldsymbol{\Gamma}_{m,g}x_{m,i,t}+d_{i,g}+\mathcal{E}_{i,g,t}$

for subject $i=1,2,...,n$ in region $g=1,2,...,G$ and time $t=1,2,...,T$. Here
$\mathcal{E}_{i,g,t}\in\mathbb{R}^{p_{1,g}\times\cdots\times p_{D,g}}$ is the
error tensor, of which the elements are assumed to follow a normal
distribution with zero mean and shared variance $\sigma_{y}^{2}$.
$\boldsymbol{\Gamma}_{k,g}\in\mathbb{R}^{p_{1,g}\times\cdots\times p_{D,g}}$
represents activation due to the $k$-th stimulus at $g$-th brain region. Each
$\boldsymbol{\Gamma}_{k,g}$ is assumed to follow a CP structure, and M-SB
prior is assigned to the latent factors of the CP decomposition to determine
the nature of activation. Also, $d_{i,g}\in\mathbb{R}$ are region- and
subject-specific random effects which are jointly modeled to borrow
information across regions of interest. Specifically, a Gaussian graphical
LASSO prior is imposed on these random effects:

$\displaystyle\boldsymbol{d}_{i}=(d_{i,1},...,d_{i,G})^{\top}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{\Sigma}^{-1}),\quad
i=1,2,...,n,$ $\displaystyle
p(\boldsymbol{\sigma}|\zeta)=C^{-1}\prod_{g<g_{1}}[DE(\sigma_{gg_{1}}|\zeta)]\prod_{g=1}^{G}[\text{Exp}(\sigma_{gg}|\frac{\zeta}{2})]\boldsymbol{1}_{\boldsymbol{\Sigma}\in\mathcal{P}^{+}},$

where $\mathcal{P}^{+}$ is the class of all symmetric positive definite
matrices and $C$ is a normalization constant. The covariance
$\boldsymbol{\sigma}=(\sigma_{gg_{1}}:g\leq g_{1})$ is a vector of upper
triangle and diagonal entries of the precision matrix $\boldsymbol{\Sigma}$.
By the properties of multivariate Gaussian distribution, a small value of
$\sigma_{gg_{1}}$ stands for weak connectivity between regions of interest
(ROIs) $g$ and $g_{1}$, given other ROIs. In practice, a double exponential
prior distributions is employed on the off-diagonal entries of the precision
matrix $\boldsymbol{\Sigma}$ to favor shrinkage among these entries. A full
Bayesian prior construction is completed by assigning Gamma prior on $\zeta$
and an inverse Gamma prior on the variance parameter $\sigma_{y}^{2}$.

To study brain connectome datasets acquired using diffusion weighted magnetic
resonance imaging (DWI), Guha and Guhaniyogi [23] propose a generalized
Bayesian linear modeling framework with a symmetric tensor response and scalar
predictors. Let
$\mathcal{Y}_{i}\in\mathcal{Y}\subset\mathbb{R}^{p\times...\times p}$ be a
symmetric tensor response with diagonal entries zero,
$\boldsymbol{x}_{i}=(x_{i1},...,x_{im})^{\top}$ be $m$ predictors of interest,
and $\boldsymbol{z}_{i}=(z_{i1},...,z_{il})^{\top}$ be $l$ auxiliary
predictors corresponding to the $i$th individual. Let
$\mathcal{J}=\\{\boldsymbol{j}=(j_{1},...,j_{D}):1\leq j_{1}<\cdots<j_{D}\leq
p\\}$ be a set of indices. Given that $\mathcal{Y}_{i}$ is symmetric with
dummy diagonal entries, it suffices to build a probabilistic generative
mechanism for $y_{i,\boldsymbol{j}}~{}(\boldsymbol{j}\in\mathcal{J})$. In
practice, a set of conditionally independent generalized linear models is
utilized. Let $E(y_{i,\boldsymbol{j}})=\omega_{i,\boldsymbol{j}}$, for
$\boldsymbol{j}\in\mathcal{J}$ we have

$\omega_{i,\boldsymbol{j}}=H^{-1}(\beta_{0}+B_{1,\boldsymbol{j}}x_{i1}+\cdots+B_{m,\boldsymbol{j}}x_{im}+\beta_{1}z_{i1}+\cdots+\beta_{l}z_{il}),$

where $B_{1,\boldsymbol{j}},...,B_{m,\boldsymbol{j}}$ respectively represents
the entry $\boldsymbol{j}=(j_{1},...,j_{D})$ of the $p\times\cdots\times p$
symmetric coefficient tensors $\mathcal{B}_{1},...,\mathcal{B}_{m}$ with
diagonal entries zero, $\beta_{0},\beta_{1},...,\beta_{l}\in\mathbb{R}$ are
the intercept and coefficients corresponding to variables $z_{i1},...,z_{il}$
respectively, and $H(\cdot)$ is the link function. The model formulation
implies a similar effect of any of the auxiliary variables
$(z_{i1},...,z_{il})$ on all entries of the response tensor but varying
effects of $h$-th predictor on different entries
$\boldsymbol{j}\in\mathcal{J}$ of the response tensor. To account for
associations between tensor nodes and predictors and to achieve parsimony in
tensor coefficients, a CP-like structure is imposed on symmetric coefficient
tensors $\mathcal{B}_{1},...,\mathcal{B}_{m}$, i.e.,

$B_{h,\boldsymbol{j}}=\sum_{r=1}^{R}\lambda_{h,r}u_{h,j_{1}}^{(r)}\cdots
u_{h,j_{D}}^{(r)},\quad h=1,2,...,m;~{}\boldsymbol{j}\in\mathcal{J},$ (44)

where
$\boldsymbol{u}_{h}^{(r)}=(u_{h,1}^{(r)},...,u_{h,p}^{(r)})^{\top}\in\mathbb{R}^{p}$
are latent factors and $\lambda_{h,r}\in\\{0,1\\}$ is the binary inclusion
variable determining if the $r$-th summand in (44) is relevant in model
setting. Further let
$\tilde{\boldsymbol{u}}_{h,k}=(u_{h,k}^{(1)},...,u_{h,k}^{(R)})$, then the
$h$-th predictor of interest is considered to have no impact on the $k$-th
tensor if $\tilde{\boldsymbol{u}}_{h,k}=0$. In order to directly study the
effect of tensor nodes related to the $h$-th predictor of interest, a spike-
and-slab mixture distribution prior is assigned on
$\tilde{\boldsymbol{u}}_{h,k}$:

$\displaystyle\tilde{\boldsymbol{u}}_{h,k}\sim\begin{cases}\mathcal{N}(\boldsymbol{0},\boldsymbol{M}_{h}),&\text{if
}\eta_{h,k}=1\\\ \delta_{\boldsymbol{0}},&\text{if
}\eta_{h,k}=0\end{cases},~{}~{}\eta_{h,k}\sim\text{Bern}(\xi_{h}),$
$\displaystyle\boldsymbol{M}_{h}\sim IW(\boldsymbol{S},\nu),~{}~{}\xi_{h}\sim
U(0,1),$

where $\delta_{\boldsymbol{0}}$ is the Dirac function at $\boldsymbol{0}$ and
$\boldsymbol{M}_{h}$ is a covariance matrix of order $R\times R$.
$IW(\boldsymbol{S},\nu)$ denotes an Inverse-Wishart distribution with an
$R\times R$ positive definite scale matrix $\boldsymbol{S}$ and degrees of
freedom $\nu$. The parameter $\xi_{h}$ corresponds to the probability of the
nonzero mixture component and $\eta_{h,k}$ is a binary indicator that equals
$0$ if $\tilde{\boldsymbol{u}}_{h,k}=\delta_{\boldsymbol{0}}$. Thus, the
posterior distributions of $\eta_{h,k}$’s can help identify nodes related to a
predictor.

To impart increasing shrinkage on $\lambda_{h,r}$ as $r$ grows, a hierarchical
prior is imposed on $\lambda_{h,r}$:

$\lambda_{h,r}\sim\text{Bern}(\nu_{h,r}),~{}\nu_{h,r}\sim\text{Beta}(1,r^{\zeta}),\zeta>1.$

Additionally, a Gaussian prior is placed $N(a_{\beta},b_{\beta})$ on
$\beta_{0},\beta_{1},...,\beta_{l}$.

Recently, Lee et al. [49] develop a Bayesian skewed tensor normal (BSTN)
regression, which addresses the problem of considerable skewness in the
tensorized response in a study of periodontal disease (PD). For an order-$K$
tensor response $\mathcal{Y}_{i}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$
with a vector of covariates $\boldsymbol{x}_{i}\in\mathbb{R}^{p}$, the
regression model is given by

$\mathcal{Y}_{i}=\mathcal{B}\bar{\times}_{(K+1)}\boldsymbol{x}_{i}+\mathcal{E}_{i},\quad\text{for
}i=1,2,...,n,$

where $\mathcal{B}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}\times p}$ is an
order-$(K+1)$ coefficient tensor, $\bar{\times}_{(K+1)}$ is the $(K+1)$-th
mode vector product, and
$\mathcal{E}_{i}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$ is the error
tensor. The skewness in the distribution of $\mathcal{Y}$ is modeled by

$\mathcal{E}_{i}=|\mathcal{Z}_{2i}|\times_{K}\boldsymbol{\Lambda}+\mathcal{Z}_{1i},$

where
$\boldsymbol{\Lambda}=\text{diag}(\lambda_{1},...,\lambda_{d_{K}})\in\mathbb{R}^{d_{K}\times
d_{K}}$ is a digonal matrix with skewness parameters
$\boldsymbol{\lambda}=(\lambda_{1},...,\lambda_{d_{K}})$, $|\boldsymbol{M}|$
denotes a matrix whose elements are absolute values of the corresponding
elements in matrix $\boldsymbol{M}$, and $\times_{K}$ is the mode-$K$ tensor
matrix product. The tensor
$\mathcal{Z}_{2i}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$ follows a
tensor normal distribution
$\mathcal{Z}_{2i}\sim\text{TN}(\boldsymbol{0};\boldsymbol{I}_{d_{1}},...,\boldsymbol{I}_{d_{K-1}},\boldsymbol{D}_{\boldsymbol{\sigma}}^{2})$,
and is assumed to be independent of
$\mathcal{Z}_{1i}\sim\text{TN}(\boldsymbol{0};\boldsymbol{R}_{1},...,\boldsymbol{R}_{K-1},\boldsymbol{D}_{\boldsymbol{\sigma}}\boldsymbol{R}_{K}\boldsymbol{D}_{\boldsymbol{\sigma}})$,
where $\boldsymbol{R}_{1},...,\boldsymbol{R}_{K}$ are positive-definite
correlation matrices, and
$\boldsymbol{D}_{\boldsymbol{\sigma}}=\text{diag}(\sigma_{1},...,\sigma_{d_{K}})$
is a diagonal matrix of positive scale parameters
$\sigma_{1},...,\sigma_{d_{K}}$. The parameterization for the tensor normal
$\mathcal{Z}_{1i}$ via correlation matrices
$\boldsymbol{R}_{1},...,\boldsymbol{R}_{K}$ avoids the common identifiability
issue. Only the $K$th mode of $\mathcal{Z}_{2i}$ is multiplied by a skewness
matrix $\boldsymbol{\Lambda}=\text{diag}(\lambda_{1},...,\lambda_{d_{K}})$
because it is assumed that there is a same level of skewness in all
combinations of the first $(K-1)$ modes in the PD dataset. When $\lambda_{j}$
is positive (or negative), the corresponding marginal density of
$y_{i_{1},...,i_{K-1},j}$ of tensor response $\mathcal{Y}$ is skewed to the
right (left).

Various prior distributions are put on the parameters. Specifically,
independent zero-mean normal density with pre-specified variance is utilized
as the common prior for
$\boldsymbol{\lambda}=(\lambda_{1},...,\lambda_{d_{K}})$, and common
independent inverse-gamma distribution $IG(g_{1},g_{2})$ with pre-specified
shape $g_{1}>0$ and scale $g_{2}>0$ is imposed on
$\boldsymbol{\sigma}=(\sigma_{1},...,\sigma_{d_{K}})$. Besides, the parametric
correlation matrices $\boldsymbol{R}_{1},...,\boldsymbol{R}_{K}$ are assumed
to be equicorrelation matrices with independent uniform priors $Unif(-1,1)$
for unknown off-diagonal elements. The tensor normal distribution
$\text{TN}(\boldsymbol{0};\boldsymbol{C}_{1},...,\boldsymbol{C}_{K+1})$ with
zero mean and known covariance matrices
$\boldsymbol{C}_{1},...,\boldsymbol{C}_{K+1}$ are put on tensor coefficient
$\mathcal{B}$. Lee et al. [49] also propose an alternative prior distribution
for $\mathcal{B}$, where a spike-and-slab prior is employed to introduce
sparsity.

Similar to the tensor predictor regression, Gaussian Process (GP) based
nonparametric models are also studied for regression problems with tensorized
responses. Li et al. [52] proposes a method based on the Gaussian process
regression networks (GPRN), where no special kernel structure is pre-assumed,
and tensor/matrix-normal variational posteriors are introduced to improve the
inference performance.

The previously discussed methods assume a low-dimensional structure of the
predictors (either in the form of vector or matrix), and are generally
incapable of modeling high-dimensional tensorized predictors. Under such
circumstance, various tensor-on-tensor methods are proposed to deal with
regression problems with both tensor-valued responses and predictors, and some
are analyzed under the Bayesian framework. Given a tensor response
$\mathcal{Y}_{i}\in\mathbb{R}^{p_{1}\times...\times p_{K}}$ and tensor
predictor $\mathcal{X}_{i}\in\mathbb{R}^{m_{1}\times...\times m_{K}}$, Hoff
[34] propose to associate $\mathcal{Y}_{i}$ and $\mathcal{X}_{i}$ through a
Tucker structure (6)

$\mathcal{Y}_{i}=\mathcal{X}_{i}\times_{1}\boldsymbol{B}_{1}\times_{2}\boldsymbol{B}_{2}\times_{3}\cdots\times_{K}\boldsymbol{B}_{K}+\mathcal{E}_{i},$
(45)

where $\boldsymbol{B}_{1},...,\boldsymbol{B}_{K}$ are matrices of dimension
$p_{1}\times m_{1},...,p_{K}\times m_{K}$ respectively. The error tensor
$\mathcal{E}_{i}$ are i.i.d with dimension $p_{1}\times\cdots\times p_{D}$,
and are assumed to follow a tensor normal distribution

$\mathcal{E}_{i}\sim\text{TN}(\boldsymbol{0};\boldsymbol{\Sigma}_{1},...,\boldsymbol{\Sigma}_{K}).$

Under the Bayesian framework, the matrix normal priors are assigned to
$\boldsymbol{B}_{k}|\boldsymbol{\Sigma}_{k}$, and inverse Wishart priors are
imposed on $\boldsymbol{\Sigma}_{k}$ ($k=1,2,...,K$) to deliver efficient
posterior computation.

Hoff [34] requires that the responses and predictors have same number of
modes. Lock [60] circumvent this restriction by employing the regression
structure based on the tensor contraction product in (14). Utilizing the same
structure, Billio et al. [7] develop a Bayesian dynamic regression model that
allows tensor-valued predictors and responses to be of arbitrary dimension.
Specifically, denote the tensor response by
$\mathcal{Y}_{t}\in\mathbb{R}^{p_{1}\times...\times p_{D_{1}}}$ and the tensor
predictor measured at time by
$t$$\mathcal{X}_{t}\in\mathbb{R}^{q_{1}\times...\times q_{D_{2}}}$. Billio et
al. [7] propose the following dynamic regression model:

$\mathcal{Y}_{t}=\sum_{j=1}^{q}\mathcal{B}_{j}*\mathcal{Y}_{t-j}+\mathcal{A}*\mathcal{X}_{t}+\mathcal{E}_{t},$

where $\mathcal{B}_{j}$ and $\mathcal{A}$ are coefficient tensors of dimension
$p_{1}\times\cdots\times p_{D_{1}}\times p_{1}\times\cdots\times p_{D_{1}}$
and $p_{1}\times\cdots\times p_{D_{1}}\times q_{1}\times\cdots\times
q_{D_{2}}$ respectively, and $*$ is the tensor contraction product (4). The
random error tensor $\mathcal{E}_{t}$ follows a tensor normal distribution,
$\mathcal{E}_{t}\sim\text{TN}(\boldsymbol{0};\boldsymbol{\Sigma}_{1},...,\boldsymbol{\Sigma}_{D_{1}}).$
The parsimony of coefficients is achieved by CP structures on the tensor
coefficients, and an M-DGDP prior is assigned to the latent factors to promote
shrinkage across tensor coefficients and improve computational scalability in
high-dimensional settings.

###  6.3 Theoretical Properties of Bayesian Tensor Regression

In this section, we discuss the theoretical properties for several Bayesian
tensor regression methods.

In [88], the in-sample predictive accuracy of an estimator coefficient tensor
$\hat{\mathcal{W}}$ in (31) is defined by

$\|\hat{\mathcal{W}}-\mathcal{W}^{*}\|_{n}^{2}:=\frac{1}{n}\sum_{i=1}^{n}\langle
X_{i},\hat{\mathcal{W}}-\mathcal{W}^{*}\rangle^{2},$

where $\mathcal{W}^{*}$ is the true coefficient tensor,
$\\{X_{i}\\}_{i=1}^{n}$ are the observed input samples. The out-of-sample
predictive accuracy is defined by

$\|\hat{\mathcal{W}}-\mathcal{W}^{*}\|_{L_{2}(P(X))}^{2}:=E_{X\sim
P(X)}[\langle X,\hat{\mathcal{W}}-\mathcal{W}^{*}\rangle^{2}],$

where $P(X)$ is the distribution of $X$ that generates the observed samples
$\\{X_{i}\\}_{i=1}^{n}$ and the expectation is taken over independent
realization $X$ from the observed ones.

Assume that the $l_{1}$-norm of $X_{i}$ is bounded by $1$, the convergence
rate of the expected in-sample predictive accuracy of the posterior mean
estimator $\int\mathcal{W}d\Pi(\mathcal{W}|Y_{1:n})$,

$E\bigg{[}\bigg{\|}\int\mathcal{W}d\Pi(\mathcal{W}|Y_{1:n})-\mathcal{W}^{*}\bigg{\|}_{n}^{2}\bigg{]},$

is characterized by the actual degree of freedom up to a log term.
Specifically, let $d^{*}$ be the CP-rank of the true tensor $\mathcal{W}^{*}$,
and $M_{1},...,M_{K}$ be the dimensions for each order of $\mathcal{W}^{*}$,
the rate is essentially

$O\left(\frac{\text{degree of
freedom}}{n}\right)=O(\frac{d^{*}(M_{1}+\cdots+M_{K})}{n})$

(up to a log term) and is optimal. Besides, although the true rank $d^{*}$ is
unknown, by placing a prior distribution on the rank, the Bayes estimator can
appropriately estimate the rank and give an almost optimal rate depending on
the rank. In this sense, the Bayes estimator has adaptivity to the true rank.
Additionally, the frequentist methods often assume a variant of strong
convexity (e.g., a restricted eigenvalue condition [6] and restricted strong
convexity [67]) to derive a fast convergence rate of sparse estimators such as
Lasso and the trace norm regularization estimator. In contrast, the
convergence rate in [88] does not require any strong convexity assumption in
the model.

In terms of the out-of-sample predictive accuracy, the convergence rate
achieved is also optimal up to the log-term under the infinity norm
thresholding assumption ($\|\mathcal{W}^{*}\|_{\infty}<R$, where $R>0$).
Specifically, the rate is

$O(\frac{d^{*}(M_{1}+\cdots+M_{K})}{n}(R^{2}\vee 1))$

up to a log order.

Based on equation (32), Guhaniyogi et al. [26] prove the posterior consistency
of the estimated coefficient tensor $\mathcal{B}$. Define a Kulback-Leibler
(KL) neighborhood around the true tensor $\mathcal{B}_{n}^{0}$ as

$\mathbb{B}_{n}=\bigg{\\{}\mathcal{B}_{n}:\frac{1}{n}\sum_{i=1}^{n}\text{KL}(f(y_{i}|\mathcal{B}_{n}^{0}),f(y_{i}|\mathcal{B}_{n}))<\epsilon\bigg{\\}},$

where $f(\cdot)$ is the glm density in (32). Let $\Pi_{n}$ denote posterior
probability given $n$ observations, Guhnaiyogi et al. [26] establish the
posterior consistency by showing that

$\Pi_{n}(\mathbb{B}_{n}^{c})\to 0~{}~{}\text{under
}\mathcal{B}_{n}^{0}~{}~{}a.s.~{}\text{as }n\to\infty$

when the prior $\pi_{n}(\mathcal{B}_{n})$ satisfies a concentration condition.
Based on this conclusion, Guhaniyogi et al. further establish the posterior
consistency for the M-DGDP prior employed in their study.

In a subsequent work [24], the authors relax the key assumption in [26] which
requires that both the true and fitted tensor coefficients have the same rank
in CP decomposition. Instead, the theoretical properties are obtained based on
a more realistic assumption that the rank of the fitted tensor coefficient is
merely greater than the rank of the true tensor coefficients. Under additional
assumptions, the authors prove that the in-sample predictive accuracy is upper
bounded by a quantity given below:

$E_{\mathcal{B}_{n}^{0}}\int\|\mathcal{B}_{n}-\mathcal{B}_{n}^{0}\|_{n}^{2}\Pi(\mathcal{B}_{n}|y_{1:n},X_{1:n})\leq
AH_{n}/n,$

where $H_{n}=o\\{\log(n)^{d}\\}$ and $A$ are positive constants depending on
other parameters. By applying the Jensen’s inequality

$\begin{split}E_{\mathcal{B}_{n}^{0}}&[\|E(\mathcal{B}_{n}|Y_{1:n},\mathcal{X}_{1:n})-\mathcal{B}_{n}^{0}\|_{n}^{2}]\leq
E_{\mathcal{B}_{n}^{0}}\int\|\mathcal{B}_{n}-\mathcal{B}_{n}^{0}\|_{n}^{2}\Pi(\mathcal{B}_{n}|Y_{1:n},X_{1:n}),\end{split}$

the posterior mean of the tensor coefficient,
$E(\mathcal{B}_{n}|Y_{1:n},X_{1:n})$, converges to the truth with a rate of
order $n^{-1/2}$ up to a $\log(n)$ factor, which is near-optimal. Similar to
Suzuki [88], this conclusion on convergence rate does not require the strong
convexity assumption on the model.

For the AMNR function defined in equation (39), Imaizumi and Hayashi [39]
establish the asymptotic property of the distance between the true function
and its estimator. Let $f^{*}\in\mathcal{W}^{\beta}(\mathcal{X})$
($\mathcal{W}^{\beta}(\mathcal{X})$ is the Sobolev space) be the true function
and $\hat{f}_{n}$ be the estimator for $f^{*}$. Let $M^{*}$ be the rank of the
true function. Then the behavior of the distance $\|f^{*}-\hat{f}_{n}\|$
strongly depends on $M^{*}$. Let $\|f\|_{n}$ be the empirical norm satisfying

$\|f\|_{n}^{2}:=\frac{1}{n}\sum_{i=1}^{n}f(x_{i})^{2}.$

When $M^{*}$ is finite, under certain assumptions and for some finite constant
$C>0$, by [39], it follows that

$E\|\hat{f}_{n}-f^{*}\|_{n}^{2}\leq Cn^{-2\beta/(2\beta+\max_{k}I_{k})},$

where $\max_{k}I_{k}$ is the maximum dimension of the tensor predictor
$\mathcal{X}$. This property indicates that the convergence rate of the
estimator corresponds to the minimax optimal rate of estimating a function in
$\mathcal{W}^{\beta}$ on a compact support in $\mathbb{R}^{I_{k}}$. The
convergence rate of AMNR depends only on the largest dimensionality of
$\mathcal{X}$.

When $M^{*}$ is infinite, by truncating $M^{*}$ at a finite value $M$, the
convergence rate is nearly the same as the case of finite $M^{*}$, which is
slightly worsened by a factor $\gamma/(1+\gamma)$ [39]:

$E\|\hat{f}_{n}-f^{*}\|_{n}^{2}\leq
C(n^{-2\beta/(2\beta+\max_{k}I_{k})})^{\gamma/(1+\gamma)}.$

For the CATCH model in (40)-(42), Pan et al. [69] establish the asymptotic
properties for a simplified model, where only tensor predictor $\mathcal{X}$
is collected (the covariates $\boldsymbol{U}$ are not included). They define
the classification error rate of the CATCH estimator and that of the Bayes
rule as

$\displaystyle
R_{n}=\text{Pr}(\hat{Y}(\mathcal{X}^{\text{new}}|\hat{\mathcal{B}}_{k},\hat{\pi}_{k},\hat{\boldsymbol{\mu}}_{k})\neq
Y^{\text{new}}),$ $\displaystyle
R=\text{Pr}(\hat{Y}(\mathcal{X}^{\text{new}}|\mathcal{B}_{k},\pi_{k},\boldsymbol{\mu}_{k})\neq
Y^{\text{new}}),$

where $\hat{\mathcal{B}}_{k},\hat{\pi}_{k}$ and $\hat{\boldsymbol{\mu}}_{k}$
are the estimated coefficients, and $\mathcal{B}_{k},\pi_{k}$ and
$\boldsymbol{\mu}_{k}$ are true coefficients. Under certain conditions,
$R_{n}\to R$ with probability tending to 1. In other words, CATCH can
asymptotically achieve the optimal classification accuracy.

In [101], Yang and Dunsun establish the posterior contraction rate of their
proposed classification model. Suppose that data are obtained for $n$
observations $y^{n}=(y_{1},...,y_{n})^{\top}$ ($y_{i}\in\\{1,2,...,d_{0}\\}$),
which are conditionally independent given
$\boldsymbol{X}^{n}=(\boldsymbol{x}_{1},...,\boldsymbol{x}_{n})^{\top}$ with
$\boldsymbol{x}_{i}=(x_{i1},...,x_{ip_{n}})^{\top}$, $x_{ij}\in\\{1,...,d\\}$
and $p_{n}\gg n$. Assume that the design points
$\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}$ are independent observations from
an unknown probability distribution $G_{n}$ on $\\{1,2,...,d\\}^{p_{n}}$,
denote

$\begin{split}d(P,P_{0})=\int\sum_{y=1}^{d_{0}}|P&(y|x_{1},...,x_{p})-P_{0}(y|x_{1},...,x_{p})|G_{n}(dx_{1},...,dx_{p}),\end{split}$

where $P_{0}$ is the true distribution, and $P$ is the estimated distribution.
Then under the given prior and other assumptions, it follows that

$\Pi_{n}\\{P:d(P,P_{0})\geq M\epsilon_{n}|y^{n},\boldsymbol{X}^{n}\\}\to
0~{}~{}a.s.,$

where $\epsilon_{n}\to
0~{}(n\epsilon_{n}^{2}\to\infty,\sum_{n}\exp(-n\epsilon_{n}^{2})<\infty)$, $M$
is a constant, and $\Pi_{n}(A|y^{n},\boldsymbol{X}^{n})$ is the posterior
distribution of $A$ given the observations. Based on this result, Yang and
Dunson [101] further prove that the posterior convergence of the model can be
very close to $n^{-1/2}$ under appropriate near low rank conditions.

Among tensor response regression problems, Guha and Guhaniyogi [23] establish
the convergence rate for predictive densities of their proposed SGTM model.
Specifically, let $f^{*}(\mathcal{Y}|\boldsymbol{x})$ be the true conditional
density of $\mathcal{Y}$ given $\boldsymbol{x}$ and
$f(\mathcal{Y}|\boldsymbol{x})$ be the random predictive density for which a
posterior is obtained. Define an integrated Hellinger distance between $f^{*}$
and $f$ as

$\mathcal{D}_{H}(f,f^{*})=\sqrt{\int\int(\sqrt{f(\mathcal{Y}|\boldsymbol{x})}-\sqrt{f^{*}(\mathcal{Y}|\boldsymbol{x})})^{2}\nu_{\mathcal{Y}}(d\mathcal{Y})\nu_{\boldsymbol{x}}(d\boldsymbol{x})},$

where $\nu_{\boldsymbol{x}}$ is the unknown probability measure for
$\boldsymbol{x}$ and $\nu_{\mathcal{Y}}$ is the dominating measure for $f$ and
$f^{*}$. For a sequence $\epsilon_{n}$ satisfying
$0<\epsilon_{n}<1,\epsilon_{n}\to 0$, and $n\epsilon_{n}^{2}\to\infty$, under
certain conditions it satisfies

$E_{f^{*}}\Pi_{n}\\{\mathcal{D}_{H}(f,f^{*})>4\epsilon_{n}|\\{\mathcal{Y}_{i},\boldsymbol{x}_{i}\\}_{i=1}^{n}\\}<4e^{-n\epsilon_{n}^{2}}$

for all large $n$, where $\Pi_{n}$ is the posterior density. This result
implies that the posterior probability outside a shrinking neighborhood around
the true predictive density $f^{*}$ converges to $0$ as $n\to\infty$. Under
further assumptions, the convergence rate $\epsilon_{n}$ can have an order
close to the parametric optimal rate of $n^{-1/2}$ up to a $\log(n)$ factor.

###  6.4 Posterior computation

In terms of posterior inference methods, sampling methods such as MCMC and
variational methods (e.g., Variational Expectation Maximization, Variational
Inference, and Variational Bayes) are the two popular choices for Bayesian
tensor analysis. MCMC is utilized in a majority of Bayesian tensor regression
and some Bayesian tensor completion (decomposition) problems. The ergodic
theory of MCMC guarantees that the sampled chain converges to the desired
posterior distribution, and sometimes the MAP result is utilized to initialize
the MCMC sampling for accelerating the convergence rate [99, 81]. In order to
reduce the computational cost and adapt to different situations, batch MCMC
and online MCMC are also used for posterior sampling [37, 36].

As an alternative strategy to approximate posterior densities for Bayesian
models, variational inference is very frequently employed in Bayesian tensor
completion methods. These methods do not guarantee producing samples from the
exact target density, but they are in general faster and more scalable to
large data than MCMC. In this category, Variational Expectation Maximization
(VEM) [100, 114, 113, 112], Variational Inference (VI) [115, 97, 38, 52], and
Variational Bayes (VB) [79, 37, 108, 109, 111, 90, 61] are the classical
choices, and the recently developed auto-encoding VB algorithm is employed to
deal with intractable distributions [56, 33]. Various studies also adopt
specific frameworks to reduce computational complexity (e.g., batch VB [37],
variational sparse Gaussian Processes [90, 112, 115, 97]) and accommodate
online or streaming data (e.g., online VB-EM [114], streaming VB [15, 107],
and Assumed Density Filtering / Expectation Propagation [16, 18, 70, 17]).
Additionally, Bayesian tensor completion (regression) methods also utilize
other methods including MLE [69], MAP [110] and EM [73, 32].

##  7 Conclusion

In Bayesian tensor analysis, the unique data structure and its high
dimensionality create challenges in both computation and theory. Bayesian
methods impose different decomposition structures on the tensor-valued data or
coefficients to reduce the number of free parameters. While CP, Tucker and
non-parametric decomposition are the most commonly used decomposition
structures, other decompositions have received some attention under the
Bayesian framework in recent years (e.g., tensor ring [61], tensor train [38],
neural [33]).

A full Bayesian model requires complete specification of a probabilistic model
and priors over model parameters, both of which depends on the data type. For
example, in tensor completion, when the tensor is continuous, the elements are
usually assumed to follow a Gaussian distribution with the tensor mean
following a decomposition structure [99, 56, 100]. The Gaussian distribution
can be extended to model the binary data through a link function [81]. In
terms of count data, element-wise Poisson distribution is often utilized to
relate the decomposition structure to the tensor-valued data, and a Dirichlet
or Gamma prior can be applied to latent factors or core tensor to enforce the
non-negativity in coefficients [79, 37, 80]. For tensor regression problems,
multivariate normal priors are placed over latent factors of the CP
decomposition, with Gaussian-Wishart prior on the hyper-parameters of the
normal distribution to achieve conjugacy [99, 12, 81]. Besides, specific
priors on core tensor (e.g., MGP prior [74, 73], Gamma-Beta hierarchical prior
[37]) or latent factors [109] in CP/Tucker structure can promote automatic
rank inference by letting the posterior decide the optimal rank. Sparsity
priors such as M-DGDP prior [26, 7] and M-SB prior [27] are also popular
choices for latent factors in the CP structure to promote low rankness, and
local/global sparsity. Integrating robust, interpretable and computationally
scalable Bayesian tensor methods with complex models (e.g., nonlinear machine
learning, reinforcement learning, causal inference, and dynamic models) is a
still interesting future direction.

Bayesian tensor regression has been widely used in applications, especially in
medical imaging analysis (e.g., MRI and EGG), where high resolution spatially
correlated data are produced. For both tensor-predictor and tensor-response
regressions, there is a need to model tensor-valued coefficients, which is
achieved by using CP/Tucker decomposition or nonparametric models that utilize
Gaussian process to model the non-linear relationship in the coefficient
tensor. Posterior inference is conducted by Markov Chain Monte Carlo (MCMC)
with Gibbs sampling, optimization based methods (e.g., variational Bayes) or
streaming methods (e.g., expectation propagation). It is still of interest to
develop scalable algorithms that accommodate challenging settings such as
streaming data analysis.

In terms of theoretical studies, most of the existing work is about
(near-)optimal convergence rate for posterior distributions of the tensor
coefficients in regression-related problems [88, 24, 39, 69, 101, 23]. There
are still many open problems such as theoretical analysis for Bayesian tensor
completion (and other tensor problems that we did not cover in this review)
and convergence analysis of computational algorithms.

## References

  * [1] Evrim Acar, Daniel M Dunlavy, Tamara G Kolda, and Morten Mørup.  Scalable tensor factorizations for incomplete data.  Chemometrics and Intelligent Laboratory Systems, 106(1):41–56, 2011\.
  * [2] Claus A Andersson and Rasmus Bro.  Improving the speed of multi-way algorithms:: Part i. tucker3.  Chemometrics and intelligent laboratory systems, 42(1-2):93–103, 1998.
  * [3] Charles E Antoniak.  Mixtures of dirichlet processes with applications to bayesian nonparametric problems.  The annals of statistics, pages 1152–1174, 1974.
  * [4] Anirban Bhattacharya and David B Dunson.  Sparse bayesian infinite factor models.  Biometrika, pages 291–306, 2011.
  * [5] Xuan Bi, Xiwei Tang, Yubai Yuan, Yanqing Zhang, and Annie Qu.  Tensors in statistics.  Annual review of statistics and its application, 8:345–368, 2021\.
  * [6] Peter J Bickel, Ya’acov Ritov, and Alexandre B Tsybakov.  Simultaneous analysis of lasso and dantzig selector.  2009\.
  * [7] Monica Billio, Roberto Casarin, Matteo Iacopini, and Sylvia Kaufmann.  Bayesian dynamic tensor regression.  Journal of Business & Economic Statistics, pages 1–11, 2022.
  * [8] Xavier Boyen and Daphne Koller.  Tractable inference for complex stochastic processes.  In Proceedings of the Fourteenth conference on Uncertainty in artificial intelligence, pages 33–42, 1998.
  * [9] Rasmus Bro.  Parafac. tutorial and applications.  Chemometrics and intelligent laboratory systems, 38(2):149–171, 1997.
  * [10] Rasmus Bro.  Multi-way analysis in the food industry-models, algorithms, and applications.  In MRI, EPG and EMA,” Proc ICSLP 2000. Citeseer, 1998.
  * [11] Han Chen, Garvesh Raskutti, and Ming Yuan.  Non-convex projected gradient descent for generalized low-rank tensor regression.  The Journal of Machine Learning Research, 20(1):172–208, 2019.
  * [12] Xinyu Chen, Zhaocheng He, and Lijun Sun.  A bayesian tensor decomposition approach for spatiotemporal traffic data imputation.  Transportation research part C: emerging technologies, 98:73–84, 2019.
  * [13] Wei Chu and Zoubin Ghahramani.  Probabilistic models for incomplete multi-dimensional arrays.  In Artificial Intelligence and Statistics, pages 89–96. PMLR, 2009\.
  * [14] Curt Da Silva and FJ Herrmann.  Hierarchical tucker tensor optimization-applications to tensor completion. sampta 2013.  In 10th International Conference on Sampling Theory and Application, Jacobs University Bremen, 2013.
  * [15] Yishuai Du, Yimin Zheng, Kuang-chih Lee, and Shandian Zhe.  Probabilistic streaming tensor decomposition.  In 2018 IEEE International Conference on Data Mining (ICDM), pages 99–108. IEEE, 2018.
  * [16] Shikai Fang, Robert M Kirby, and Shandian Zhe.  Bayesian streaming sparse tucker decomposition.  In Uncertainty in Artificial Intelligence, pages 558–567. PMLR, 2021.
  * [17] Shikai Fang, Akil Narayan, Robert Kirby, and Shandian Zhe.  Bayesian continuous-time tucker decomposition.  In International Conference on Machine Learning, pages 6235–6245. PMLR, 2022.
  * [18] Shikai Fang, Zheng Wang, Zhimeng Pan, Ji Liu, and Shandian Zhe.  Streaming bayesian deep tensor factorization.  In International Conference on Machine Learning, pages 3133–3142. PMLR, 2021.
  * [19] Mostafa Reisi Gahrooei, Hao Yan, Kamran Paynabar, and Jianjun Shi.  Multiple tensor-on-tensor regression: An approach for modeling processes with heterogeneous sources of data.  Technometrics, 63(2):147–159, 2021.
  * [20] Silvia Gandy, Benjamin Recht, and Isao Yamada.  Tensor completion and low-n-rank tensor recovery via convex optimization.  Inverse problems, 27(2):025010, 2011.
  * [21] Alan E Gelfand and Adrian FM Smith.  Sampling-based approaches to calculating marginal densities.  Journal of the American statistical association, 85(410):398–409, 1990.