#  Functional Bayesian Tucker Decomposition for Continuous-indexed Tensor Data

Shikai Fang
Kahlert School of Computing
The University of Utah
<EMAIL_ADDRESS>
&Xin Yu
Kahlert School of Computing
The University of Utah
<EMAIL_ADDRESS>
&Zheng Wang
Kahlert School of Computing
The University of Utah
<EMAIL_ADDRESS>
&Shibo Li
Kahlert School of Computing
The University of Utah
<EMAIL_ADDRESS>
&Robert M. Kirby
Kahlert School of Computing
The University of Utah
<EMAIL_ADDRESS>
&Shandian Zhe
Kahlert School of Computing
The University of Utah
<EMAIL_ADDRESS>

###### Abstract

Tucker decomposition is a powerful tensor model to handle multi-aspect data.
It demonstrates the low-rank property by decomposing the grid-structured data
as interactions between a core tensor and a set of object representations
(factors). A fundamental assumption of such decomposition is that there were
finite objects in each aspect or mode, corresponding to discrete indexes of
data entries. However, many real-world data are not naturally posed in the
setting. For example, geographic data is represented as continuous indexes of
latitude and longitude coordinates, and cannot fit tensor models directly. To
generalize Tucker decomposition to such scenarios, we propose Functional
Bayesian Tucker Decomposition (FunBaT). We treat the continuous-indexed data
as the interaction between the Tucker core and a group of latent functions. We
use Gaussian processes (GP) as functional priors to model the latent
functions, and then convert the GPs into a state-space prior by constructing
an equivalent stochastic differential equation (SDE) to reduce computational
cost. An efficient inference algorithm is further developed for scalable
posterior approximation based on advanced message-passing techniques. The
advantage of our method is shown in both synthetic data and several real-world
applications.

##  1 Introduction

Tensor decomposition is widely used to analyze and predict multi-way or multi-
aspect data in real-world applications. For example, medical service records
can be summarized by a four-mode tensor (patients, doctor, clinics, time) and
the entry values can be the visit count or drug usage. Weather data can be
abstracted by tensors like (latitude, longitude, time) and the entry values
can be the temperature. Tensor decomposition introduces a set of low-rank
representations of objects in each mode, known as factors, and estimates these
factors to reconstruct the observed entry values. Among numerous proposed
tensor decomposition models, such as CANDECOMP/PARAFAC (CP) [Harshman, 1970]
decomposition and tensor train(TT) [Oseledets, 2011], Tucker decomposition
[Tucker, 1966] is famous for its fair capacity, flexibility and
interpretability.

Despite its success, the standard tensor decomposition framework has a
fundamental limitation. That is, it restricts the data format as structured
grids, and demands each mode of the tensor must be discrete and finite-
dimensional. That means, each mode includes a finite set of objects, indexed
by integers, such as user 1, user 2, …, and each data value is located by the
tuple of discrete indexes. However, in many real-world applications like
climate and geographic modeling, the data is represented as continuous
coordinates of modes, e.g., latitude, longitude and time.. In these cases,
each data entry is located by a tuple of continuous indexes, and tensor
decomposition cannot be applied directly.

To enable tensor models on such data, one widely-used trick is to discretize
the modes, for example, binning the timestamps into time steps (say by weeks
or months), or binning the latitude values into a set of ranges, and number
the ranges by 1, 2, 3… While doing this is feasible, the fine-grained
information is lost, and it is hard to decide an appropriate discretization
strategy in many cases [Pasricha et al., 2022]. A more natural solution is to
extend tensor models from grid-structure to continuous field. To the best of
our knowledge, all existing methods in that direction are limited to tensor-
train [Gorodetsky et al., 2015, Bigoni et al., 2016, Ballester-Ripoll et al.,
2019, Chertkov et al., 2023] or simpler CP format [Schmidt, 2009], and cannot
be applied to Tucker decomposition, which is with more compact and flexible
low-rank representation. What’s more, most of the proposed methods are
deterministic with polynomial approximation, and hard to handle noise data
well and provide probabilistic inference.

To bridge the gap, we propose FunBaT: Functional Bayesian Tucker
decomposition, which generalizes the standard Tucker decomposition to the
functional field under the probabilistic framework. We decompose the
continuous-indexed tensor data as the interaction between the Tucker core and
a group of latent functions, which map the continuous indexes to mode-wise
factors. We approximate such functions by Gaussian processes (GPs), and follow
[Hartikainen and Särkkä, 2010] to convert them into state-space priors by
constructing equivalent stochastic differential equations (SDE) to reduce
computational cost. We then develop an efficient inference algorithm based on
the conditional expectation propagation (CEP) framework [Wang and Zhe, 2019]
and Bayesian filter to approximate the posterior distribution of the factors.
We show this framework can also be applied to extend CP decomposition to
continuous-indexed data with the simplified setting. For evaluation, we
conduct experiments on both synthetic data and real-world applications of
climate modeling. The results show that FunBaT not only can outperform the
state-of-the-art methods with large margins in terms of reconstruction error,
and even can identify interpretable patterns that align with domain knowledge.

##  2 Preliminary

###  2.1 Tensor Decomposition and Function Factorization

Standard tensor decomposition operates under the following settings: given a
$K$-mode tensor $\mathcal{Y}\in\mathbb{R}^{d_{1}\times\cdots\times d_{K}}$,
where each mode $k$ contains $d_{k}$ objects, and each tensor entry is indexed
by a $K$-size tuple $\mathbf{i}=\left(i_{1},\ldots,i_{K}\right)$, denoted as
${y}_{\mathbf{i}}$. To decompose the tensor into a compact and low-rank form,
we introduce $K$ groups of latent factors, denoted as
$\mathcal{U}=\left\\{\mathbf{U}^{1},\ldots,\mathbf{U}^{K}\right\\}$. where
$\mathbf{U}^{k}=\left[\mathbf{u}_{1}^{k},\ldots,\mathbf{u}_{d_{k}}^{k}\right]^{\top}$.
Each factor $\mathbf{u}_{j}^{k}$ represents the latent embedding of the $j$-th
object in the $k$-th mode. The classic CANDECOMP/PARAFAC (CP) [Harshman, 1970]
assumes each factor $\mathbf{u}_{j}^{k}$ is a $r$-size vector, and tensor
entry is decomposed as the product-sum of the
factors:${y}_{\mathbf{i}}\approx\left(\mathbf{u}_{i_{1}}^{1}\circ\ldots\circ\mathbf{u}_{i_{K}}^{K}\right)^{\top}\mathbf{1}$,
where $\circ$ is the element-wise product. Tucker decomposition [Tucker, 1966]
takes one more step towards a more expressive and flexible low-rank structure.
It allows us to set different factor ranks $\\{r_{1},\ldots r_{K}\\}$ for each
mode, by introducing a core tensor
$\mathcal{W}\in\mathbb{R}^{r_{1}\times\cdots\times r_{K}}$, known as the
Tucker core, to model the interaction as the weighted-sum over all possible
cross-rank factors:

$\displaystyle{y}_{\mathbf{i}}\approx\operatorname{vec}(\mathcal{W})^{\top}\left(\mathbf{u}_{i_{1}}^{1}\otimes\ldots\otimes\mathbf{u}_{i_{K}}^{K}\right)$
(1)

where $\otimes$ is the Kronecker product and $\operatorname{vec}(\cdot)$ is
the vectorization operator.Tucker decomposition will degenerate to CP if we
set all modes’ ranks equal and $\mathcal{W}$ diagonal. Tensor Train(TT)
decomposition [Oseledets, 2011] is another classical tensor method. TT sets
the TT-rank $\\{r_{0},r_{1},\ldots r_{K}\\}$ firstly, where
$r_{0}=r_{K}=1$.Then factor $\mathbf{u}_{j}^{k}$ is defined as a $r_{k}\times
r_{k+1}$ matrix. The tensor entry $y_{\mathbf{i}}$ is then decomposed as a
series of the matrix product of the factors.

Function factorization refers to decomposing a complex multivariate function
into a set of low-rank separated functions [Rai, 2014, Nouy, 2015]. The most
important and widely used tensor method for function factorization is TT
[Gorodetsky et al., 2015, Bigoni et al., 2016, Ballester-Ripoll et al., 2019].
It converts the function approximation into a tensor decomposition problem. A
canonical form of applying TT to factorize a multivariate function
$f(\mathbf{x})$ with K variables $\mathbf{x}=(x_{1},\ldots x_{K})$is:

$\displaystyle f(\mathbf{x})\approx G_{1}(x_{1})\times
G_{2}(x_{2})\times\ldots\times G_{K}(x_{K})$ (2)

where $G_{k}(\cdot)$, a univariate matrix-valued function, takes the scaler
$x_{k}$ as input and output a matrix $G_{k}(x_{k})\in\mathbb{R}^{r_{k-1}\times
r_{k}}$. It’s clear to see that $G_{k}(\cdot)$ is a generalization of the
factor $\mathbf{u}_{j}^{k}$ in standard TT decomposition.

###  2.2 Gaussian Process as state space model

Gaussian process (GP) [Rasmussen and Williams, 2006] is a powerful non-
parametric Bayesian model to approximate functions. Formally, a GP is a prior
distribution over function $f(\cdot)$, such that $f(\cdot)$ evaluated at any
finite set of points $\mathbf{x}=\left\\{x_{1},\ldots,x_{N}\right\\}$ follows
a multivariate Gaussian distribution, denote as
$f\sim\mathcal{G}\mathcal{P}\left(0,\kappa\left({x},{x}^{\prime}\right)\right)$,
where $\kappa\left({x},{x}^{\prime}\right)$ is the covariance function, also
known as the kernel, measuring the similarity between two points. One of the
most widely used kernels is Matérn kernel:

$\displaystyle\kappa_{\text{Mat\'{e}rn}}=\sigma^{2}\frac{\left(\frac{\sqrt{2\nu}}{l}{||{x}-{x}^{\prime}||}^{2}\right)^{\nu}}{\Gamma(\nu)2^{\nu-1}}K_{\nu}\left(\frac{\sqrt{2\nu}}{l}{||{x}-{x}^{\prime}||}^{2}\right)$
(3)

where $\\{\sigma^{2},l,\nu,p\\}$ are hyperparameters determining the variance,
length-scale, smoothness, and periodicity of the function, $K_{\nu}$ is the
modified Bessel function, and $\Gamma(\cdot)$ is the Gamma function.

Despite the great capacity of GP, it suffers from cubic scaling complexity
$O(N^{3})$ for inference. To overcome this limitation, recent work
[Hartikainen and Särkkä, 2010] applied the spectral analysis tools to show an
equivalence between GPs with stationary kernels and linear time-invariant
stochastic differential equations (LTI-SDEs). Specifically, we can formulate a
vector-valued latent state $\mathbf{z}(x)$ comprising $f(x)$ and its
successive derivatives of $m$-th order, and the GP
$f({x})\sim\mathcal{G}\mathcal{P}\left(0,\kappa\right)$ is equivalent to the
solution of a LTI-SDE defined as:

$\displaystyle\mathbf{z}(x)=\left(f(x),\frac{\mathrm{d}f(x)}{\mathrm{d}x},\ldots,\frac{\mathrm{d}f^{m}(x)}{\mathrm{d}x}\right)^{\top},\
\frac{\mathrm{d}\mathbf{z}(x)}{\mathrm{d}x}=\mathbf{F}\mathbf{z}(x)+\mathbf{L}{w}(x)$
(4)

where $\mathbf{F}$ and $\mathbf{L}$ are time-invariant coefficients, and
${w}(x)$ is the white noise process with density $q_{s}$. At any finite
collections of points $\\{x_{1},\ldots,x_{N}\\}$ with increasing order, the
SDE in 4 can be further discretized as the Gaussian-Markov chain, also known
as the state-space model, defined as:

$\displaystyle
p(\mathbf{z}(x))=p(\mathbf{z}_{1})\prod_{n=1}^{N-1}p(\mathbf{z}_{n+1}|\mathbf{z}_{n})=\mathcal{N}(\mathbf{z}(x_{1})|\mathbf{0},\mathbf{P_{\infty}})\prod_{n=1}^{N-1}\mathcal{N}(\mathbf{z}(x_{n+1})|\mathbf{A}_{n}\mathbf{z}(x_{n}),\mathbf{Q}_{n})$
(5)

where $\mathbf{A}_{n}=\operatorname{exp}(\mathbf{F}\Delta_{n})$,
$\mathbf{Q}_{n}=\int_{t_{n}}^{t_{n+1}}\mathbf{A}_{n}\mathbf{L}\mathbf{L}^{\top}\mathbf{A}_{n}^{\top}q_{\mathbf{s}}\mathrm{d}t$,
$\Delta_{n}=x_{n+1}-x_{n}$, and $\mathbf{P_{\infty}}$ is the steady-state
covariance matrix which can be obtained by solving the Lyapunov equation
[Lancaster and Rodman, 1995]. All the above parameters in 4 and 5 are fully
determined by the kernel $\kappa$ and the time interval $\Delta_{n}$. For the
Matérn kernel 3 with half-integer smoothness $\nu$, $\mathbf{F}$, $\mathbf{L}$
and $\mathbf{P_{\infty}}$ are with closed-form formulas [Särkkä, 2013].
Specifically, when $\nu=1/2$, we have
$\\{m=0,\mathbf{F}=-1/l,\mathbf{L}=1,q_{\mathbf{s}}=2\sigma^{2}/l,\mathbf{P_{\infty}}=\sigma^{2}\\}$.
For $\nu=3/2$, we have
$m=1,\mathbf{F}=(0,1;-\lambda^{2},-2\lambda),\mathbf{L}=(0;1),q_{\mathbf{s}}=4\sigma^{2}\lambda^{3},\mathbf{P_{\infty}}=(\sigma^{2},0;0,\lambda^{2}\sigma^{2})$,
where $\lambda=\sqrt{3}/l$. With the state space format and proper design of
the observation, efficient $O(n)$inference can be performed using classical
Bayesian sequential inference techniques, like Kalman filtering and smoothing
[Hartikainen and Särkkä, 2010]. Then the original function $f(x)$ is simply
the first element of the inferred latent state $\mathbf{z}(x)$.

##  3 model

###  3.1 Functional Tucker Decomposition with Gaussian Process

Despite the successes of tensor models, the standard tensor models are
unsuitable for continuous-indexed data, such as the climate data with mode
e.g., latitude, longitude and time. The continuity property encoded in the
real-value index will be dropped while applying the discretization trick.
Additionally, it can be challenging to determine the optimal discretization
strategy, and the trained model cannot handle new objects with never-seen
indexes. The function factorization idea with tensor structure is a more
natural solution to this problem. However, the early work [Schmidt, 2009] with
the most simple CP model, is based on sampling-based inference, which is with
limited capacity and not scalable to large-scale data. The TT-based function
factorization methods [Gorodetsky et al., 2015, Bigoni et al., 2016,
Ballester-Ripoll et al., 2019, Chertkov et al., 2023] take deterministic
approximation like Chebyshev polynomials or splines, which is hard to handle
noise data well or provide uncertainty quantification.

Compared to CP and TT, Tucker decomposition is with more compact and flexible
low-rank representation. The global Tucker core can even capture more global
and interpretable patterns of the data. Thus, we aim to extend Tucker
decomposition to continuous-indexed data to fully utilize its advantages as a
Bayesian method, and propose FunBaT: FunBaT: Functional Bayesian Tucker
decomposition.

Aligned with the setting of the function factorization (2), we factorize a
$K$-variate function $f(\mathbf{i})$ in Tucker format (1) with preset latent
rank:$\\{r_{1},\ldots r_{K}\\}$:

$\displaystyle f({\mathbf{i}})=f(i_{1},\ldots
i_{K})\approx\operatorname{vec}(\mathcal{W})^{\top}\left(\mathbf{U}^{1}(i_{1})\otimes\ldots\otimes\mathbf{U}^{K}(i_{K})\right)$
(6)

where $\mathcal{W}\in\mathbb{R}^{r_{1}\times\cdots\times r_{K}}$ is the Tucker
core, which is the same as the standard Tucker decomposition.
$\mathbf{U}^{k}(\cdot)$ is a $r_{k}$-size vector-valued function, mapping the
continuous index $i_{k}$ of mode $k$ to a $r_{k}$-size latent factor.We assign
independent GP priors over each dimension of $\mathbf{U}^{k}(\cdot)$, and
model it $\mathbf{U}^{k}(\cdot)$ as the stack of groups of univariate scaler
functions. Specifically, we have:

$\displaystyle\mathbf{U}^{k}(i_{k})=[u^{k}_{1}(i_{k}),\ldots,u^{k}_{r_{k}}(i_{k})]^{T};\
u^{k}_{j}(i_{k})\sim\mathcal{GP}\left({0},\kappa(i_{k},i_{k}^{\prime})\right),j=1\dots
r_{k}$ (7)

where
$\kappa(i_{k},i_{k}^{\prime}):\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$
is the covariance (kernel) function of the $k$-th mode. In this work, we pick
the popular Matérn kernel (3) for GPs over all modes.

###  3.2  State-space-prior and Joint Probabilities

Given $N$ observed entries of a $K$-mode continuous-indexed tensor
$\mathcal{Y}$, denoted as $\mathcal{D}=\\{(\mathbf{i}^{n},y_{n})\\}_{n=1\dots
N}$, where $\mathbf{i}^{n}=(i^{n}_{1}\ldots i^{n}_{K})$ is the tuple of
continuous indexes, and $y_{n}$ is the entry values, we can assume all the
observation are sampled from the target Tucker-format function $f(\mathbf{i})$
and Gaussian noise $\tau^{-1}$. Specifically, we define the likelihood $l_{n}$
as:

$\displaystyle l_{n}\triangleq
p(y_{n}|\\{\mathbf{U}^{k}\\}_{k=1}^{K},\mathcal{W},\tau)=\mathcal{N}\left(y_{n}\mid\operatorname{vec}(\mathcal{W})^{\top}\left(\mathbf{U}^{1}(i_{1}^{n})\otimes\ldots\otimes\mathbf{U}^{K}(i_{K}^{n})\right),\tau^{-1}\right)$
(8)

We then handle the functional priors over the
$\\{\mathbf{U}^{k}\\}_{k=1}^{K}$. As introduced in Section 2.2, the dimension-
wise GP with Matérn kernel over $u^{k}_{j}(i_{k})$ is equivalent to a
$M$-order LTI-SDE (4), and then we can discrete it as a state space model (5).
Thus, for each mode’s function $\mathbf{U}^{k}$, we concatenate the state
space representation over all $r_{k}$ dimensions of $\mathbf{U}^{k}(i_{k})$ to
a joint state variable
$\mathbf{Z}^{k}(i_{k})=\operatorname{concat}[\mathbf{z}^{k}_{1}(i_{k})\ldots\mathbf{z}^{k}_{r_{k}}(i_{k})]$
for $\mathbf{U}^{k}$, where $\mathbf{z}^{k}_{j}(i_{k})$ is the state variable
of $u^{k}_{j}(i_{k})$. We can build an ordered index set
$\mathcal{I}_{k}=\\{i_{k}^{1}\ldots i_{k}^{N_{k}}\\}$, which includes the
$N_{k}$ unique indexes of mode-$k$’s among all observed entries in
$\mathcal{D}$. Then, the state space prior over $\mathbf{U}^{k}$ on
$\mathcal{I}_{k}$ is:

$\displaystyle
p(\mathbf{U}^{k})=p(\mathbf{Z}^{k})=p(\mathbf{Z}^{k}(i_{k}^{1})\ldots\mathbf{Z}^{k}(i_{k}^{N_{k}}))=p(\mathbf{Z}^{k}_{1})\prod_{s=1}^{N_{k}-1}p(\mathbf{Z}^{k}_{s+1}|\mathbf{Z}^{k}_{s})$
(9) where $\displaystyle
p(\mathbf{Z}^{k}_{1})=\mathcal{N}(\mathbf{Z}^{k}(i^{1}_{k})|\mathbf{0},\tilde{\mathbf{P}}_{\infty}^{k});\
p(\mathbf{Z}^{k}_{s+1}|\mathbf{Z}^{k}_{s})=\mathcal{N}(\mathbf{Z}^{k}(i^{s+1}_{k})|\tilde{\mathbf{A}}^{k}_{s}\mathbf{Z}^{k}(i_{k}^{s}),\tilde{\mathbf{Q}}_{s}^{k})$
(10)

The parameters in (10) are block-diagonal concatenates of the corresponding
univariate case.
Namely,$\tilde{\mathbf{P}}_{\infty}^{k}=\operatorname{BlockDiag}(\mathbf{P}_{\infty}^{k}\ldots\mathbf{P}_{\infty}^{k})$,
$\tilde{\mathbf{A}}^{k}_{s}=\operatorname{BlockDiag}(\mathbf{A}^{k}_{s}\ldots\mathbf{A}^{k}_{s})$,
$\tilde{\mathbf{Q}}^{k}_{s}=\operatorname{BlockDiag}(\mathbf{Q}^{k}_{s}\ldots\mathbf{Q}^{k}_{s})$,
where $\mathbf{P}_{\infty}^{k}$, $\mathbf{A}^{k}_{s}$ and $\mathbf{Q}^{k}_{s}$
are the corresponding parameters in (5). With $\mathbf{Z}^{k}(i_{k})$, we can
get the value of $\mathbf{U}^{k}$ by the projection matrix
$\mathbf{H}=\operatorname{BlockDiag}([1,0,\ldots],\ldots[1,0,\ldots])$:
$\mathbf{U}^{k}=\mathbf{H}\mathbf{Z}^{k}$.

With state space priors of the latent functions, we further assign a Gamma
prior over the noise precision $\tau$ and a Gaussian prior over the Tucker
core. For compact notation, we denote the set of all random variables as
$\mathbf{\Theta}=\\{\mathcal{W},\tau,\\{\mathbf{Z}^{k}\\}_{k=1}^{K}\\}$.
Finally, the joint probability of the proposed model is:

$\displaystyle
p(\mathcal{D},\mathbf{\Theta})=p(\mathcal{D},\\{\mathbf{Z}^{k}\\}_{k=1}^{K},\mathcal{W},\tau)$
$\displaystyle=p(\tau)p(\mathcal{W})\prod_{k=1}^{K}[p(\mathbf{Z}^{k}_{1})\prod_{s=1}^{N_{k}-1}p(\mathbf{Z}^{k}_{s+1}|\mathbf{Z}^{k}_{s})]\prod_{n=1}^{N}l_{n}$
(11)

where
$p(\mathcal{W})=\mathcal{N}(\operatorname*{vec}(\mathcal{W})\mid\mathbf{0},\mathbf{I})$
and $p(\tau)=\text{Gam}(\tau|a_{0},b_{0})$. $a_{0},b_{0}$ are the
hyperparameters of the Gamma prior. $l_{n}$ is the data likelihood defined at
(8).

##  4 Algorithm

The inference of the exact posterior $p(\mathbf{\Theta}|\mathcal{D})$ with
(11) is intractable, as there are multiple latent functions
$\\{\mathbf{U}^{k}\\}_{k=1}^{K}$ and the Tucker core interleave together in a
complex manner in the likelihood term $l_{n}$. To address this issue, we
propose a novel and efficient approximation algorithm, which first decouples
the likelihood term by a factorized approximation, and then applies the
sequential inference of Kalman filter and RTS smoother to infer the latent
factors $\\{\mathbf{Z}^{k}\\}_{k=1}^{K}$ at each observed index. Finally, we
employ the conditional moment matching technique to update the message factors
$\\{f_{n}\\}$ in parallel. We will introduce the details in the following
subsections.

###  4.1 Factorized Approximation with Gaussian and Gamma distribution

To estimate the intractable posterior $p(\mathbf{\Theta}|\mathcal{D})$ with a
tractable $q(\mathbf{\Theta})$, we first apply the mean-field assumption and
design the approximated posterior as a fully factorized format. Specifically,
we approximate the posterior as:

$\displaystyle p(\mathbf{\Theta}|\mathcal{D})\approx
q(\mathbf{\Theta})=q(\tau)q(\mathcal{W})\prod_{k=1}^{K}q(\mathbf{Z}^{k})$ (12)

where $q(\tau)=\operatorname{Gam}(\tau|a,b)$,
$\mathcal{N}(\operatorname{vec}(\mathcal{W})\mid\bm{\mu},\mathbf{S})$ are the
approximated posterior of $\tau$, $\mathcal{W}$. For $q(\mathbf{Z}^{k})$, we
further decompose it over the observed indexes set $\mathcal{I}_{k}$ as
$q(\mathbf{Z}^{k})=\prod_{s=1}^{N_{k}}q(\mathbf{Z}^{k}_{s})$, where
$q(\mathbf{Z}^{k}_{s})=q(\mathbf{Z}^{k}(i^{k}_{s}))=\mathcal{N}(\mathbf{Z}^{k}_{s}\mid\mathbf{m}^{k}_{s},\mathbf{V}^{k}_{s})$.
Our goal is to estimate the variational parameters
$\\{a,b,\bm{\mu},\mathbf{S},\\{\mathbf{m}^{k}_{s},\mathbf{V}^{k}_{s}\\}\\}$,
and make $q(\mathbf{\Theta})$ close to $p(\mathbf{\Theta}|\mathcal{D})$.

To do so, we use Expectation Propagation (EP) [Minka, 2001], to update
$q(\mathbf{\Theta})$. However, the standard EP cannot work due to the complex
Tucker form of the likelihood term $l_{n}$ makes it intractable to compute the
expectation of the likelihood term $l_{n}$ under $q(\mathbf{\Theta})$. Thus,
we further apply the advanced moment-matching technique, called Conditional
Expectation Propagation (CEP) [Wang and Zhe, 2019] to address this issue. With
the CEP, we put factorized approximation to decouple the likelihood term
$l_{n}$ into a group of message factors $\\{f_{n}\\}$ with closed-form:

$\displaystyle\mathcal{N}(y_{n}\mid\operatorname{vec}(\mathcal{W})^{\top}\left(\mathbf{U}^{1}(i_{1}^{n})\otimes\ldots\otimes\mathbf{U}^{K}(i_{K}^{n})\right),\tau^{-1})\approx{Z}_{n}f_{n}(\tau)f_{n}(\mathcal{W})\prod_{k=1}^{K}f_{n}(\mathbf{Z}^{k}(i_{k}^{n}))$
(13)

where $Z_{n}$ is the normalized constant,
$f_{n}(\tau)=\operatorname{Gam}(\tau|a_{n},b_{n})$,
$f_{n}(\mathcal{W})=\mathcal{N}(\operatorname{vec}(\mathcal{W})\mid\bm{\mu}_{n},\mathbf{S}_{n})$,
$f_{n}(\mathbf{Z}^{k}(i_{k}^{n}))=\mathcal{N}(\mathbf{Z}^{k}(i_{k}^{n})|\mathbf{m}^{k}_{n},\mathbf{V}^{k}_{n})$
are the message factors obtained from the conditional moment-matching of
$\mathbf{E}_{q}[l_{n}]$. Then we update the posterior $q(\tau)$ and
$q(\mathcal{W})$ that will minimize $\text{KL}(p||q)$ by simply merging the
message factors from likelihood (13) and priors:

$\displaystyle
q(\tau)=p(\tau)\prod_{n=1}^{N}f_{n}(\tau)=\operatorname{Gamma}(\tau|a_{0},b_{0})\prod_{n=1}^{N}\operatorname{Gam}(\tau|a_{n},b_{n})$
(14) $\displaystyle
q(\mathcal{W})=p(\mathcal{W})\prod_{n=1}^{N}f_{n}(\mathcal{W})=\mathcal{N}(\operatorname*{vec}(\mathcal{W})\mid\mathbf{0},\mathbf{I})\prod_{n=1}^{N}\mathcal{N}(\operatorname{vec}(\mathcal{W})\mid\bm{\mu}_{n},\mathbf{S}_{n})$
(15)

The message approximation in (13) and massage merging (14)(15) is based on the
conditional moment-matching technique and property of exponential family
present in [Wang and Zhe, 2019]. All involved computations are closed-form
updates and can be conducted in parallel. We omit the exact formulas and put
them along with the introduction of EP and CEP in the appendix.

###  4.2 Sequential State Inference with Bayesian Filter and Smoother

Handling $q(\mathbf{Z}^{k})$ is a bit more challenging. It is because the
prior of $\mathbf{Z}^{k}$ is with the chain-structured. (9)(10), and we need
to handle complicated integration over the whole chain to get marginal
posterior $q(\mathbf{Z}^{k}_{s})$ in the standard inference. However, with the
classical Bayesian filter and smoother method, we can infer the posterior
sequentially. Specifically, the state space structure over
$\mathbf{Z}^{k}_{s}$ is:

$\displaystyle
q(\mathbf{Z}^{k}_{s})=q(\mathbf{Z}^{k}_{s-1})p(\mathbf{Z}^{k}_{s}|\mathbf{Z}^{k}_{s-1})\prod_{n\in\mathcal{D}^{k}_{s}}f_{n}(\mathbf{Z}^{k}_{s})$
(16)

where $p(\mathbf{Z}^{k}_{s}|\mathbf{Z}^{k}_{s-1})$ is the transitions of the
state given in (10), and $\mathcal{D}^{k}_{s}$ is the subset of $\mathcal{D}$,
including the observation entries whose $k$-mode index is $i^{s}_{k}$. Namely,
$\mathcal{D}^{k}_{s}=\\{n:i_{k}^{n}=i_{k}^{s}\mid\mathbf{i}_{n}\in\mathcal{D}\\}$.
If we treat $\prod_{n\in\mathcal{D}^{k}_{s}}f_{n}(\mathbf{Z}^{k}_{s})$, a
group of Gaussian message factors, as the observation of the state space
model, (16) is the standard Kalman Filter(KF) [Kalman, 1960]. Thus, we can run
the KF algorithm and compute $q(\mathbf{Z}^{k}_{s})$ from $s=1$ to $N_{k}$
sequentially. After the forward pass of the states, we can run RTS smoothing
[Rauch et al., 1965] as a backward pass to compute the global posterior of
$q(\mathbf{Z}^{k}_{s})$. This sequential inference is widely used to infer
state space GP models [Hartikainen and Särkkä, 2010, Särkkä, 2013].

Thus, the whole inference algorithm is organized as follows: with the
observation $\mathcal{D}$, we initialize the approximated posteriors and
message factors for each likelihood. Then for each mode, we firstly
approximate the message factors by CEP in parallel, and then merge them to
update $q(\tau),q(\mathcal{W})$. We run the KF and RTS smoother to infer
$q(\mathbf{Z}^{k})$ at each observed index by treating the message factors as
observations. We repeat the inference until convergence. We summarize it in
Algorithm Table 1.

Algorithm Complexity The overall time complexity is $\mathcal{O}(NKR)$, where
$K$ is the number of modes, $N$ is the number of observations and $R$ is pre-
set mode rank. The space complexity is $\mathcal{O}(NK(R+R^{2}))$, as we need
to store all message factors. The linear costs of both time and space over
data size and tensor mode show the scalability of the proposed algorithm.

Probabilistic Interpolation at Arbitrary Index Although the state-space
functional prior of $\mathbf{U}^{k}$ or $\mathbf{Z}^{k}$ is defined over
finite observed index set $\mathcal{I}_{k}$, we highlight that we can handle
the probabilistic interpolation of $\mathbf{Z}^{k}$ at arbitrary index
$i_{k}^{*}\notin\mathcal{I}_{k}$ after model inference. Specifically, with
$i_{k}^{s-1}<i_{k}^{*}<i_{k}^{s}$ we can infer the
$q(\mathbf{Z}^{k}(i_{k}^{*}))$ by integrating the messages from the
transitions of its neighbors and will give a cloase-form solution:

$\displaystyle q(\mathbf{Z}^{k}(i_{k}^{*}))=\int
q(\mathbf{Z}^{k}_{s-1})p(\mathbf{Z}^{k}(i_{k}^{*})|\mathbf{Z}^{k}_{s-1}){\rm
d}\mathbf{Z}^{k}_{s-1}\int
q(\mathbf{Z}^{k}_{s})p(\mathbf{Z}^{k}_{s})|\mathbf{Z}^{k}(i_{k}^{*}){\rm
d}\mathbf{Z}^{k}_{s}=\mathcal{N}(\mathbf{m}^{*},\mathbf{V}^{*})$ (17)

We put the detailed derivation in the appendix. This fantastic property
enables us to build a continuous trajectory for each mode, and even predict
the $\mathcal{Y}$ value at any indexes, which standard discrete-mode tensor
decomposition cannot do.

Lightweight Alternative: FunBaT-CP As the CP decomposition is a special and
simplified case of Tucker decomposition, it’s easy to build a functional CP
decomposition with the proposed model and algorithm. Specifically, we only
need to set the Tucker core $\mathcal{W}$ as all-zero constant tensor except
diagonal elements as one, skip the inference step of $\mathcal{W}$, and remain
others. Then, we get FunBaT-CP, a simple and efficient functional Bayesian CP
model. For some data, we found FunBaT-CP is more robust than FunBaT, as the
dense Tucker core has more parameters and is easier to get overfitting. We
will show it in the experiment section.

Algorithm 1 FunBaT

Input: Observations $\mathcal{D}$ of a $K$-mode continuous-indexed tensor ,
kernel hyperparameters,, sorted unique indexes set $\\{\mathcal{I}_{k}\\}$ of
each mode.

Initialize approx. posterior $q(\tau)$, $q(\mathcal{W})$,
$\\{q(\mathbf{Z}^{k})\\}$ and message factors for each likelihood.

repeat

for $k=1$ to $K$ do

Approximate messages factors
$\\{f_{n}(\mathbf{Z}^{k}(i_{k}^{n})),f_{n}(\tau),f_{n}(\mathcal{W})\\}_{n=1}^{N}$
in parallel with CEP (13)

Update the approximated posterior $q(\tau)$, $q(\mathcal{W})$ by merging the
message (14)(15).

Update $q(\mathbf{Z}^{k})$ sequentially by KF and RTS smoother based on (16)

end for

until Convergence

Return: $q(\tau),q(\mathcal{W}),\\{q(\mathbf{Z}^{k})\\}_{k=1}^{K}$

##  5 Related Work

The early work to apply tensor format to function factorization is[Schmidt,
2009]. It takes the wrapped-GP to factorize the function with CP format and
infers with Monte Carlo sampling, which is not scalable to large-scale data.
Applying tensor-train(TT) to approximate the multivariate function with low-
rank structure is a popular topic [Gorodetsky et al., 2015, Bigoni et al.,
2016, Ballester-Ripoll et al., 2019, Chertkov et al., 2022, 2023], estimated
for applications in simultaneous localization and mapping(SLAM) [Aulinas et
al., 2008]. These methods mainly use Chebyshev polynomials and splines as
function basis, and rely on complex optimization methods like alternative
update, cross-interpolation[Gorodetsky et al., 2015, Bigoni et al., 2016] and
ANOVA (analysis of variance) representation[Ballester-Ripoll et al., 2019,
Chertkov et al., 2023]. Despite the compact form, they are purely
deterministic, sensitive to initial values, and hard to handle the
probabilistic inference well with data noise. Fang et al. [2022] performed
similar techniques like state-space GP and CEP to our work, but it is designed
to model temporal information in standard Tucker decomposition. Our model
takes one step to extend standard Tucker to the continuous-indexed and
functional field.

##  6 Experiment

###  6.1 Synthetic Data

We evaluated FunBaT on a synthetic task by simulating a rank-$1$ two-mode
tensor with each mode designed as a continuous function:

$\displaystyle\mathbf{U}^{1}(i_{1})=\exp(-2i_{1})\cdot\sin(\frac{3}{2}\pi
i_{1});\mathbf{U}^{2}(i_{2})=\sin^{2}(2\pi i_{2})\cdot\cos(2\pi i_{2})$ (18)

The ground truth of the continuous-mode tensor, represented as a surface in
Figure 1(a), was obtained by:
$y_{\mathbf{i}}=\mathbf{U}^{1}(i_{1})\mathbf{U}^{1}(i_{2})$. We randomly
sampled $650$ indexes entries from $[0,1]\times[0,1]$ and add Gaussian noise
$\epsilon\sim\mathcal{N}(0,0.02)$ as the observations. We used PyTorch to
implement FunBaT, which used the Matérn kernel with $\nu=3/2$, and set $l=0.1$
and $\sigma^{2}=1$. We trained the model with $R=1$ and compared the learned
continuous-mode functions with their ground truth. Our learned trajectories,
shown in Figure 1(c) and Figure 1(d) clearly revealed the real mode functions.
The shaded area is the estimated standard deviation. Besides, we used the
learned factors and the interpolation (17) to reconstruct the whole tensor
surface, shown in Figure 1(b). We found that the reconstructed surface matched
the ground truth closely with limited samples.

|  |  |
---|---|---|---

(a) Real tensor surface

(b) Estimated tensor surface

(c) $\mathbf{U}^{1}(i_{1})$

(d) $\mathbf{U}^{2}(i_{2})$

Figure 1: Results of Synthetic Data | RMSE | MAE
---|---|---
Datasets | PM2.5 | PM10 | SO2 | PM2.5 | PM10 | SO2
Resolution:$50\times 50\times 150$ |  |  |
P-Tucker | $0.805\pm 0.017$ | $0.787\pm 0.006$ | $0.686\pm 0.02$ | $0.586\pm 0.003$ | $0.595\pm 0.005$ | $0.436\pm 0.011$
Tucker-ALS | $1.032\pm 0.049$ | $1.005\pm 0.029$ | $0.969\pm 0.027$ | $0.729\pm 0.016$ | $0.741\pm 0.007$ | $0.654\pm 0.034$
Tucker-SVI | $0.792\pm 0.01$ | $0.8\pm 0.026$ | $0.701\pm 0.08$ | $0.593\pm 0.01$ | $0.605\pm 0.019$ | $0.423\pm 0.031$
Resolution:$100\times 100\times 300$ |  |  |
P-Tucker | $0.8\pm 0.101$ | $0.73\pm 0.021$ | $0.644\pm 0.023$ | $0.522\pm 0.011$ | $0.529\pm 0.013$ | $0.402\pm 0.008$
Tucker-ALS | $1.009\pm 0.027$ | $1.009\pm 0.026$ | $0.965\pm 0.023$ | $0.738\pm 0.01$ | $0.754\pm 0.007$ | $0.68\pm 0.011$
Tucker-SVI | $0.706\pm 0.011$ | $0.783\pm 0.067$ | $0.69\pm 0.086$ | $0.509\pm 0.008$ | $0.556\pm 0.031$ | $0.423\pm 0.031$
Resolution: $300\times 300\times 1000$ |  |  |
P-Tucker | $0.914\pm 0.126$ | $1.155\pm 0.001$ | $0.859\pm 0.096$ | $0.401\pm 0.023$ | $0.453\pm 0.002$ | $0.366\pm 0.015$
Tucker-ALS | $1.025\pm 0.044$ | $1.023\pm 0.038$ | $1.003\pm 0.019$ | $0.742\pm 0.011$ | $0.757\pm 0.011$ | $0.698\pm 0.007$
Tucker-SVI | $1.735\pm 0.25$ | $1.448\pm 0.176$ | $1.376\pm 0.107$ | $0.76\pm 0.033$ | $0.747\pm 0.028$ | $0.718\pm 0.023$
Resolution:$428\times 501\times 1461$ (original) |  |  |
P-Tucker | $1.256\pm 0.084$ | $1.397\pm 0.001$ | $0.963\pm 0.169$ | $0.451\pm 0.017$ | $0.493\pm 0.001$ | $0.377\pm 0.019$
Tucker-ALS | $1.018\pm 0.034$ | $1.012\pm 0.021$ | $0.997\pm 0.024$ | $0.738\pm 0.005$ | $0.756\pm 0.007$ | $0.698\pm 0.011$
Tucker-SVI | $1.891\pm 0.231$ | $1.527\pm 0.107$ | $1.613\pm 0.091$ | $0.834\pm 0.032$ | $0.787\pm 0.018$ | $0.756\pm 0.014$
Methods using continuous indexes |  |  |
FTT-ALS | $1.020\pm 0.013$ | $1.001\pm 0.013$ | $1.001\pm 0.026$ | $0.744\pm 0.007$ | $0.755\pm 0.007$ | $0.696\pm 0.011$
FTT-ANOVA | $2.150\pm 0.033$ | $2.007\pm 0.015$ | $1.987\pm 0.036$ | $1.788\pm 0.031$ | $1.623\pm 0.014$ | $1.499\pm 0.018$
FTT-cross | $0.942\pm 0.025$ | $0.933\pm 0.012$ | $0.844\pm 0.026$ | $0.566\pm 0.018$ | $0.561\pm 0.011$ | $0.467\pm 0.033$
RBF-SVM | $0.995\pm 0.015$ | $0.955\pm 0.02$ | $0.794\pm 0.026$ | $0.668\pm 0.008$ | $0.674\pm 0.014$ | $0.486\pm 0.026$
BLR | $0.998\pm 0.013$ | $0.977\pm 0.014$ | $0.837\pm 0.021$ | $0.736\pm 0.007$ | $0.739\pm 0.008$ | $0.573\pm 0.009$
FunBaT-CP | $0.296\pm 0.018$ | $0.343\pm 0.028$ | $\mathbb{0.386\pm 0.009}$ | $\mathbb{0.18\pm 0.002}$ | $0.233\pm 0.013$ | $0.242\pm 0.003$
FunBaT | $\mathbb{0.288\pm 0.008}$ | $\mathbb{0.328\pm 0.004}$ | $\mathbb{0.386\pm 0.01}$ | $0.183\pm 0.006$ | $\mathbb{0.226\pm 0.002}$ | $\mathbb{0.241\pm 0.004}$
Table 1: Prediction error over BeijingAir-PM2.5, BeijingAir-PM10, and BeijingAir-SO2 with $R=2$, which were averaged over five runs. The results for $R=3,5,7$ are in the supplementary. | RMSE | MAE
---|---|---
Mode-Rank | R=3 | R=5 | R=7 | R=3 | R=5 | R=7
P-Tucker | $1.306\pm 0.02$ | $1.223\pm 0.022$ | $1.172\pm 0.042$ | $0.782\pm 0.011$ | $0.675\pm 0.014$ | $0.611\pm 0.007$
Tucker-ALS | $>10$ | $>10$ | $>10$ | $>10$ | $>10$ | $>10$
Tucker-SVI | $1.438\pm 0.025$ | $1.442\pm 0.021$ | $1.39\pm 0.09$ | $0.907\pm 0.005$ | $0.908\pm 0.005$ | $0.875\pm 0.072$
FTT-ALS | $1.613\pm 0.0478$ | $1.610\pm 0.052$ | $1.609\pm 0.055$ | $0.967\pm 0.009$ | $0.953\pm 0.007$ | $0.942\pm 0.010$
FTT-ANOVA | $5.486\pm 0.031$ | $4.619\pm 0.054$ | $3.856\pm 0.059$ | $4.768\pm 0.026$ | $4.026\pm 0.100$ | $3.123\pm 0.0464$
FTT-cross | $1.415\pm 0.0287$ | $1.312\pm 0.023$ | $1.285\pm 0.052$ | $0.886\pm 0.011$ | $0.822\pm 0.006$ | $0.773\pm 0.014$
RBF-SVM | $2.374\pm 0.047$ | $2.374\pm 0.047$ | $2.374\pm 0.047$ | $1.44\pm 0.015$ | $1.44\pm 0.015$ | $1.44\pm 0.015$
BLR | $2.959\pm 0.041$ | $2.959\pm 0.041$ | $2.959\pm 0.041$ | $2.029\pm 0.011$ | $2.029\pm 0.011$ | $2.029\pm 0.011$
FunBaT-CP | $\mathbb{0.805\pm 0.06}$ | $\mathbb{0.548\pm 0.03}$ | $\mathbb{0.551\pm 0.048}$ | $\mathbb{0.448\pm 0.06}$ | $\mathbb{0.314\pm 0.005}$ | $\mathbb{0.252\pm 0.008}$
FunBaT | $1.255\pm 0.108$ | $1.182\pm 0.117$ | $1.116\pm 0.142$ | $0.736\pm 0.069$ | $0.647\pm 0.05$ | $0.572\pm 0.089$
Table 2: Prediction error of US-TEMP, which were averaged over five runs.

###  6.2 Real-World Applications

Datasets We evaluated FunBaT on four real-world datasets BeijingAir-PM2.5,
BeijingAir-PM10, BeijingAir-SO2 and US-TEMP. The first three are extracted
from BeijingAir111https://archive.ics.uci.edu/ml/datasets/Beijing+Multi-
Site+Air-Quality+Data, which contain hourly measurements of several air
pollutants with weather features in Beijing from 2014 to 2017. We selected
three continuous-indexed modes: (atmospheric-pressure, temperature, time) and
processed it into three tensors by using different pollutants (PM2.5, PM10,
SO2). Each dataset contains $17K$ observations across unique indexes of 428
atmospheric pressure measurements, 501 temperature measurements, and 1461
timestamps. We obtain US-TEMP from the
ClimateChange222https://berkeleyearth.org/data/. The dataset contains
temperatures of cities worldwide and geospatial features. We selected
temperature data from 248 cities in the United States from 1750 to 2016 with
three continuous-indexed modes: (latitude, longitude, time). The tensor
contains $56K$ observations across 15 unique latitudes, 95 longitudes, and 267
timestamps.

Baselines and Settings We set three groups of methods as baselines. The first
group includes the state-of-art standard Tucker models, which are P-Tucker [Oh
et al., 2018]: a scalable Tucker algorithm that performs parallel row-wise
updates,Tucker-ALS: Efficient Tucker decomposition algorithm using alternative
least square(ALS) update. [Bader and Kolda, 2008] and Tucker-SVI [Hoffman et
al., 2013]: Bayesian version of Tucker updating with stochastic variational
inference(SVI). The second group includes functional tensor-train(FTT) based
methods with different functional basis and optimization strategy:FTT-ALS
[Bigoni et al., 2016],FTT-ANOVA [Ballester-Ripoll et al., 2019],and FTT-cross
[Gorodetsky et al., 2015].As we can view the prediction task as a regression
problem and use the continuous indexe as features, we add RBF-SVM Support
Vector Machine [Hearst et al., 1998] with RBF kernels. and BLR Bayesian Linear
Regression [Minka, 2000] as competing methods as the third group.

We used the official open-source implementations of most baselines. We use the
TENEVA library [Chertkov et al., 2022, 2023] to test the FTT-based methods. We
re-scaled all continuous-mode indexes to $[0,1]$ to numerical robustness. For
FunBaT, we finetuning the models with Matérn kernels $\nu=\\{1/2,3/2\\}$ along
the kernel parameters for optimal performance for different datasets. We
examined all the methods with the factors rank $R\in\\{2,3,5,7\\}$. We assume
all mode’s rank is equal to $R$. Following [Fang et al., 2021a, Tillinghast et
al., 2020], we randomly sampled $80\%$ observed entry values for training and
then tested on the remaining. We repeated the experiments five times, and then
examined the average root mean-square-error (RMSE), average mean-absolute-
error (MAE), and their standard deviations.

To demonstrate the advantages of using continuous indexes rather than index
discretization, we set four different discretization granularities by binding
the raw continuous indexes into several discrete-indexed bins on BeijingAir-
PM2.5, BeijingAir-PM10, and BeijingAir-SO2. We then conduct the tensors with
four different resolutions: $428\times 501\times 1461$ (original resolution),
$300\times 300\times 1000$, $100\times 100\times 300$, and $50\times 50\times
150$. We tested the performance of the first group(standard tensor
decomposition) baselines on different resolutions. We tested FunBaT,FunBaT-CP,
and other baselines with continuous indexes at the original resolution.

Prediction Results Due to space limits, we only present the results with $R=2$
for BeijingAir-PM2.5, BeijingAir-PM10, and BeijingAir-SO2 datasets in Table 1,
while the other results are in supplementary. The prediction results for US-
TEMP are listed in Table 2. As shown in the tables, our approach FunBaT and
FunBaT-CP outperforms the competing methods by a significant margin in all
cases. We found the performance of standard Tucker methods varies a lot with
different discrete resolutions, showing the hardness to decide the optimal
discretization in real-world applications. We also found the interesting fact
that the performance of FunBaT-CP is much better than FunBaT on US-TEMP. As we
found their training loss is similar, it may be because the dense Tucker core
of FunBaT result in overfitting and is not suitable for the sparse US-TEMP
dataset.

|
---|---

(a) $\mathbf{U^{1}}(i_{1})$:Mode of latitude

(b) $\mathbf{U^{2}}(i_{2})$:Mode of longitude

Figure 2: Learned mode functions of US-TEMP Figure 3:
$\mathbf{U^{3}}(i_{3})$:Mode of time

Investigation of Learned Functions We explored whether the learned mode
functions reveal interpretable patterns that are consistent with domain
knowledge. We run FunBaT on US-TEMP dataset with $R=1$, and plotted the
learned latent functions of the three modes (latitude, longitude, time), shown
in Figure 23. As the first two modes correspond to latitude and longitude,
respectively, representing real geo-locations, we marked out the city groups
(the shaded circles) and some city names in Figure 2(a)2(b). The
$\mathbf{U}^{1}$ of the lantide mode in Figure 2(a) shows a clear decreasing
trend from southern cities like Miami (in Florida) to northern cities like
Anchorage (in Alaska), which is consistent with the fact that temperatures are
generally higher in the south and lower in the north. The learned longitude-
mode $\mathbf{U}^{2}$ in Figure 2(b) exhibits a steady trend from east to
west, but with a region (the blue circle) with lower values near the Western
longitude ${110}^{\circ}$. That is the Rocky Mountain Area including cities
like Denver and Salt Lake City. with higher altitudes, and results in lower
temperatures. The time-mode $\mathbf{U}^{3}$ in Figure 3 provides meaningful
insights of the climate change in history. It reveals a gradual increase over
the past 260 years, with a marked acceleration after 1950. This pattern aligns
with the observed trend of rising global warming following the
industrialization of the mid-20th century. The function around 1750-1770,
characterized by lower and oscillating values, corresponds to the Little Ice
Age, a well-documented period of global cooling in history.

##  7 Conclusion

We proposed FunBaT, a Bayesian method to generalize Tucker decomposition to
the functional field to model continuous-indexed tensor data. We adopt the
state-space GPs as functional prior and develop an efficient inference
algorithm based on CEP. The results on both synthetic and real-world tasks
demonstrate the effectiveness of our method.

## References

  * Aulinas et al. [2008] Josep Aulinas, Yvan Petillot, Joaquim Salvi, and Xavier Lladó.  The slam problem: a survey.  _Artificial Intelligence Research and Development_ , pages 363–371, 2008.
  * Bader and Kolda [2008] Brett W Bader and Tamara G Kolda.  Efficient matlab computations with sparse and factored tensors.  _SIAM Journal on Scientific Computing_ , 30(1):205–231, 2008.
  * Ballester-Ripoll et al. [2019] Rafael Ballester-Ripoll, Enrique G Paredes, and Renato Pajarola.  Sobol tensor trains for global sensitivity analysis.  _Reliability Engineering & System Safety_, 183:311–322, 2019.
  * Bickel and Doksum [2015] Peter J Bickel and Kjell A Doksum.  _Mathematical statistics: basic ideas and selected topics, volume I_ , volume 117.  CRC Press, 2015.
  * Bigoni et al. [2016] Daniele Bigoni, Allan P Engsig-Karup, and Youssef M Marzouk.  Spectral tensor-train decomposition.  _SIAM Journal on Scientific Computing_ , 38(4):A2405–A2439, 2016.
  * Chertkov et al. [2022] Andrei Chertkov, Gleb Ryzhakov, Georgii Novikov, and Ivan Oseledets.  Optimization of functions given in the tensor train format.  _arXiv preprint arXiv:2209.14808_ , 2022.  doi: 10.48550/ARXIV.2209.14808.  URL https://arxiv.org/pdf/2209.14808.pdf.
  * Chertkov et al. [2023] Andrei Chertkov, Gleb Ryzhakov, and Ivan Oseledets.  Black box approximation in the tensor train format initialized by anova decomposition.  _SIAM Journal on Scientific Computing_ , 45(4):A2101–A2118, 2023.
  * Fang et al. [2021a] Shikai Fang, Robert M Kirby, and Shandian Zhe.  Bayesian streaming sparse tucker decomposition.  In _Uncertainty in Artificial Intelligence_ , pages 558–567. PMLR, 2021a.
  * Fang et al. [2021b] 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, 2021b.
  * Fang et al. [2022] 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.
  * Gorodetsky et al. [2015] Alex A Gorodetsky, Sertac Karaman, and Youssef M Marzouk.  Function-train: a continuous analogue of the tensor-train decomposition.  _arXiv preprint arXiv:1510.09088_ , 2015.
  * Harshman [1970] R. A. Harshman.  Foundations of the PARAFAC procedure: Model and conditions for an”explanatory”multi-mode factor analysis.  _UCLA Working Papers in Phonetics_ , 16:1–84, 1970.
  * Hartikainen and Särkkä [2010] Jouni Hartikainen and Simo Särkkä.  Kalman filtering and smoothing solutions to temporal gaussian process regression models.  In _2010 IEEE international workshop on machine learning for signal processing_ , pages 379–384. IEEE, 2010.
  * Hearst et al. [1998] Marti A. Hearst, Susan T Dumais, Edgar Osuna, John Platt, and Bernhard Scholkopf.  Support vector machines.  _IEEE Intelligent Systems and their applications_ , 13(4):18–28, 1998.
  * Hoffman et al. [2013] Matthew D Hoffman, David M Blei, Chong Wang, and John Paisley.  Stochastic variational inference.  _Journal of Machine Learning Research_ , 2013.
  * Kalman [1960] Rudolph Emil Kalman.  A new approach to linear filtering and prediction problems.  1960\.
  * Lancaster and Rodman [1995] Peter Lancaster and Leiba Rodman.  _Algebraic riccati equations_.  Clarendon press, 1995.
  * Minka [2000] Thomas Minka.  Bayesian linear regression.  Technical report, Citeseer, 2000.
  * Minka [2001] Thomas P Minka.  Expectation propagation for approximate bayesian inference.  In _Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence_ , pages 362–369, 2001.
  * Nouy [2015] Anthony Nouy.  Low-rank tensor methods for model order reduction.  _arXiv preprint arXiv:1511.01555_ , 2015.
  * Oh et al. [2018] Sejoon Oh, Namyong Park, Sael Lee, and Uksong Kang.  Scalable Tucker factorization for sparse tensors-algorithms and discoveries.  In _2018 IEEE 34th International Conference on Data Engineering (ICDE)_ , pages 1120–1131. IEEE, 2018.
  * Oseledets [2011] Ivan V Oseledets.  Tensor-train decomposition.  _SIAM Journal on Scientific Computing_ , 33(5):2295–2317, 2011.
  * Pasricha et al. [2022] Ravdeep S Pasricha, Ekta Gujral, and Evangelos E Papalexakis.  Adaptive granularity in tensors: A quest for interpretable structure.  _Frontiers in Big Data_ , 5:929511, 2022.
  * Rai [2014] Prashant Rai.  _Sparse low rank approximation of multivariate functions–Applications in uncertainty quantification_.  PhD thesis, Ecole Centrale de Nantes (ECN), 2014.
  * Rasmussen and Williams [2006] Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  MIT Press, 2006.
  * Rauch et al. [1965] Herbert E Rauch, F Tung, and Charlotte T Striebel.  Maximum likelihood estimates of linear dynamic systems.  _AIAA journal_ , 3(8):1445–1450, 1965.
  * Särkkä [2013] Simo Särkkä.  _Bayesian filtering and smoothing_.  Number 3. Cambridge University Press, 2013.
  * Schmidt [2009] Mikkel N Schmidt.  Function factorization using warped gaussian processes.  In _Proceedings of the 26th Annual International Conference on Machine Learning_ , pages 921–928, 2009.
  * Tillinghast et al. [2020] Conor Tillinghast, Shikai Fang, Kai Zhang, and Shandian Zhe.  Probabilistic neural-kernel tensor decomposition.  In _2020 IEEE International Conference on Data Mining (ICDM)_ , pages 531–540. IEEE, 2020.
  * Tucker [1966] Ledyard Tucker.  Some mathematical notes on three-mode factor analysis.  _Psychometrika_ , 31:279–311, 1966.
  * Wang and Zhe [2019] Zheng Wang and Shandian Zhe.  Conditional expectation propagation.  In _UAI_ , page 6, 2019.

##  Appendix A Conditional Expectation Proprogation

###  A.1 Breif Introduction of expectation proprogation(EP)

As a general framework of Bayesian learning, expectation proprogation (EP)
[Minka, 2001] is a variational inference method for computing the posterior
distribution of a latent variable $\bm{\theta}$ given the observed data
$\mathcal{D}$. As the exact posterior $p(\bm{\theta}\mid\mathcal{D})$ can be
formed as following form:

$\displaystyle
p(\bm{\theta}\mid\mathcal{D})=\frac{1}{Z}\prod_{n}p_{n}(\mathcal{D}_{n}\mid\bm{\theta})$
(19)

where $p_{n}(\mathcal{D}_{n}\mid\bm{\theta})$ is the likelihood of $n$-th data
given $\bm{\theta}$, which may also link to the priors. $Z$ is the normalizer
to ensure (19) is a valid distribution. The exact computation of (19) is
always infeasible due to the complex form of
$p_{n}(\mathcal{D}_{n}\mid\bm{\theta})$. Therefore, EP naturally adopts a
factorized distribution approximation of real posterior:

$\displaystyle p(\bm{\theta}\mid\mathcal{D})\approx
q(\bm{\theta}\mid\mathcal{D})\sim\prod_{n}f_{n}(\bm{\theta})$ (20)

where $f_{n}(\bm{\theta})$ is an approximated factor for the actual data
likelihood factor $p_{n}(\mathcal{D}_{n}\mid\bm{\theta})$. If we can assign a
proper form of $f_{n}(\bm{\theta})$, and make $f_{n}(\bm{\theta})\approx
p_{n}(\mathcal{D}_{n}\mid\bm{\theta})$ as close as possible, then we can get a
good approximation of the posterior distribution. The EP algorithm sets the
basic form of factors using exponential-family distribution, which includes
many distributions such as Gaussian, Gamma, and Beta, and has the canonical
formula $f_{n}(\bm{\theta})\sim
exp(\bm{\lambda}_{n}^{\top}\phi\left(\bm{\theta}\right))$.
$\bm{\lambda}_{n}^{\top}$ and $\phi\left(\bm{\theta}\right)$ are regarded as
natural parameter and sufficient statistics of $f_{n}(\bm{\theta})$
respectively. The process of EP is to iteratively update the parameters of
each $f_{n}(\bm{\theta})$ until it converges. The update procedure to update
$f_{n}(\bm{\theta})$ is as follows:

(1) build a calibrating distribution (also known as context factors) by
removing current approx. factor :

$\displaystyle q^{\backslash n}(\bm{\theta})\propto
q(\bm{\theta})/f_{n}\left(\bm{\theta}\right)$ (21)

(2) build a tilted distribution by multiplying back the real data likelihood
with context factors:

$\displaystyle\hat{p}_{n}(\bm{\theta})\sim
p_{n}(\mathcal{D}_{n}\mid\bm{\theta})q^{\backslash n}(\bm{\theta})$ (22)

(3)build a new approximated posterior $q^{\star}$, estimate its parameters by
matching moments with the tilted distribution:

$\displaystyle\mathbb{E}_{q^{*}}(\phi(\bm{\theta}))=\mathbb{E}_{\hat{p}_{n}}(\phi(\bm{\theta}))$
(23)

(4)update the approximated posterior $f_{n}$ by removing context factors from
the new approximated posterior:

$\displaystyle f_{n}(\bm{\theta})\propto q^{\star}(\bm{\theta})/q^{\backslash
n}(\bm{\theta})$ (24)

###  A.2 Conditional Moment Match

With complex data likelihood where there are multi variables interleave, such
as the Tucker-based likelihood in the main paper, computation of
$\mathbb{E}_{\hat{p}}(\phi(\bm{\Theta}))$ in moment match step of classic EP
is intractable. To solve this problem, [Wang and Zhe, 2019] propose a
conditional moment match method to update the parameters of the approximated
posterior $f_{n}(\mathbf{Z}^{k}(i_{k}^{n}))$. Take $\mathbf{Z}^{k}(i_{k}^{n})$
for example, the required moments are
$\phi(\mathbf{Z}^{k}(i_{k}^{n}))=(\mathbf{Z}^{k}(i_{k}^{n}),\mathbf{Z}^{k}(i_{k}^{n}){\mathbf{Z}^{k}(i_{k}^{n})}^{T})$.
the key idea of CEP is to decompose the expectation into a nested structure:

$\displaystyle\mathbb{E}_{\hat{p}}(\phi(\bm{\Theta}))=\mathbb{E}_{\hat{p}\left(\bm{\Theta}_{\backslash\mathbf{Z}^{k}(i_{k}^{n})}\right)}\left[\mathbb{E}_{\hat{p}\left(\mathbf{Z}^{k}(i_{k}^{n})\mid\bm{\Theta}_{\backslash\mathbf{Z}^{k}(i_{k}^{n})}\right)}\left[\bm{\phi}(\bm{\Theta})\mid\bm{\Theta}_{\backslash\mathbf{Z}^{k}(i_{k}^{n})}\right]\right]$
(25)

where
$\bm{\Theta}_{\backslash\mathbf{Z}^{k}(i_{k}^{n})}\triangleq\bm{\Theta}\backslash\\{\mathbf{Z}^{k}(i_{k}^{n})\\}$.
Therefore, we can compute the conditional moment first, i.e., the inner
expectation, with an analytical form. However, the computation of outer-
expectation are under the marginal tilted distribution which is still
intractable. To solve this problem, we follow the BCTT [Fang et al., 2022] and
apply a multivariate delta method [Bickel and Doksum, 2015, Fang et al.,
2021b] to get:

$\displaystyle\mathbb{E}_{q\left(\bm{\Theta}_{\backslash\mathbf{Z}^{k}(i_{k}^{n})}\right)}\left[\bm{\rho}_{n}\right]\approx\bm{\rho}_{n}\left(\mathbb{E}_{q}\left[\bm{\Theta}_{\backslash\mathbf{Z}^{k}(i_{k}^{n})}\right]\right)$
(26)

where $\bm{\rho}_{n}$ denote the conditional moments of tilde distribution.

###  A.3 CEP update for FunBaT and FunBaT-CP

Following the above paradigm, we could work out all update formulas for all
parameters of the approximated message factors
$\\{f_{n}\\}=\\{\\{f_{n}(\mathbf{Z}^{k}(i_{k}^{n}))\\}_{k=1}^{K},f_{n}(\tau)\\}$
for FunBaT and FunBaT-CP.

For
$f_{n}(\mathbf{Z}^{k}(i_{k}^{n}))=\mathcal{N}({\mathbf{H}}\mathbf{Z}^{k}(i_{k}^{n})|\mathbf{m}^{k}_{n},\mathbf{S}^{k}_{n})$:

$\displaystyle\mathbf{S}^{k}_{n}=(\mathbb{E}_{q}[\tau]\mathbb{E}_{q}[\mathbf{a}^{\backslash
k}_{n}{\mathbf{a}^{\backslash k}_{n}}^{T}])^{-1};\
\mathbf{m}^{k}_{n}=\mathbf{S}^{k}_{n}(y_{n}\mathbb{E}_{q}[\tau]\mathbb{E}_{q}[\mathbf{a}^{\backslash
k}_{n}])$ (27)

where

$\displaystyle\mathbf{a}^{\backslash
k}_{n}=\mathcal{W}_{(k)}\left(\mathbf{U}^{K}({i_{K}^{n}})\ldots\mathbf{U}^{k+1}({i_{k+1}^{n}})\otimes\mathbf{U}^{k-1}({i_{k-1}^{n}})\ldots\mathbf{U}^{1}(i_{1}^{n})\right)$
(28)

$\mathcal{W}_{(k)}$ is the folded tucker core $\mathcal{W}$ at mode $k$.

For $f_{n}(\tau)=\operatorname{Gam}(\tau|\alpha_{n},\beta_{n})$, the update
formulas are:

$\displaystyle\alpha_{n}=\frac{3}{2};\
\beta_{n}=\frac{1}{2}y_{n}^{2}-y_{n}\mathbb{E}_{q}[\mathbf{a}_{n}]+\frac{1}{2}\operatorname{trace}[\mathbb{E}_{q}[\mathbf{a}_{n}\mathbf{a}_{n}^{T}]]$
(29)

where:

$\displaystyle\mathbf{a}_{n}=\operatorname{vec}(\mathcal{W})^{T}\left(\mathbf{U}^{1}({i_{1}^{n}})\otimes\ldots\mathbf{U}^{K}({i_{K}^{n}})\right)$
(30)

The update formulas of
$f_{n}(\mathcal{W})=\mathcal{N}(\operatorname{vec}(\mathcal{W})\mid\mathcal{\mu}_{n},\mathcal{S}_{n})$
is:

$\displaystyle\mathbf{S}_{n}=(\mathbb{E}_{q}[\tau]\mathbb{E}_{q}[\mathbf{b}_{n}{\mathbf{b}^{T}_{n}}])^{-1};\
\mathcal{\mu}_{n}=\mathbf{S}_{n}(y_{n}\mathbb{E}_{q}[\tau]\mathbb{E}_{q}[\mathbf{b}_{n}])$
(31)

Where

$\displaystyle\mathbf{b}_{n}=\mathbf{U}^{1}(i_{1}^{n})\otimes\ldots\otimes\mathbf{U}^{K}(i_{K}^{n})$
(32)

We can just use the above update formulas for FunBaT-CP by setting the
$\mathcal{W}$ as a constant diagonal tensor. However, We can also re-derive a
more convenient and elegant form based on the Hadmard product form of CP. It
will result in the similar formats on the update formulas of $f_{n}(\tau)$ and
$f_{n}(\mathbf{Z^{k}}(i_{k}^{n}))$, but with the different definitions on
$\mathbf{a}^{\backslash k}_{n}$ and $\mathbf{a}_{n}$. They are:

$\displaystyle\mathbf{a}^{\backslash
k}_{n}=\mathbf{U}^{1}(i_{1}^{n})\circ\cdots\mathbf{U}^{k-1}(i_{k-1}^{n})\circ\mathbf{U}^{k+1}(i_{k+1}^{n})\cdots\circ\mathbf{U}^{K}(i_{K}^{n})$
(33)

and

$\mathbf{a}_{n}=\mathbf{U}^{1}\left(i_{1}^{n}\right)\circ\cdots\circ\mathbf{U}^{K}\left(i_{K}^{n}\right)$
(34)

###  A.4  Derivation of the probabilistic imputation at any index

To derive the probabilistic imputation equation (17) in the main paper, we
consider a general state space model, which includes a sequence of states
${\bf x}_{1},\ldots,{\bf x}_{M}$ and the observed data $\mathcal{D}$. The
states are at time $t_{1},\ldots,t_{M}$ respectively. The key of the state
space model is that the prior of the states is a Markov chain. The joint
probability has the following form,

$\displaystyle p({\bf x}_{1},\ldots,{\bf x}_{M},\mathcal{D})=p({\bf
x}_{1})\prod_{j=1}^{M-1}p({\bf x}_{j+1}|{\bf x}_{j})\cdot p(\mathcal{D}|{\bf
x}_{1},\ldots,{\bf x}_{M}).$ (35)

Note that here we do not assume the data likelihood is factorized over each
state, like those typically used in Kalman filtering. In our point process
model, the likelihood often couples multiple states together.

Suppose we have run some posterior inference to obtain the posterior of these
states $q({\bf x}_{1},\ldots,{\bf x}_{M})$, and we can easily pick up the
marginal posterior of each state and each pair of the states. Now we want to
calculate the posterior distribution of the state at time $t^{*}$ such that
$t_{m}<t^{*}<t_{m+1}$. Denote the corresponding state by ${\bf x}^{*}$, our
goal is to compute $p({\bf x}^{*}|\mathcal{D})$. To do so, we consider
incorporating ${\bf x}^{*}$ in the joint probability (35),

$\displaystyle p({\bf x}_{1},\ldots,{\bf x}_{m},{\bf x}^{*},{\bf
x}_{m+1},\ldots,{\bf x}_{M},\mathcal{D})$ $\displaystyle=p({\bf
x}_{1})\prod_{j=1}^{m-1}p({\bf x}_{j+1}|{\bf x}_{j})\cdot p({\bf x}^{*}|{\bf
x}_{m})p({\bf x}_{m+1}|{\bf x}^{*})\cdot\prod_{j={m+1}}^{M}p({\bf
x}_{j+1}|{\bf x}_{j})\cdot p(\mathcal{D}|{\bf x}_{1},\ldots,{\bf x}_{M}).$
(36)

Now, we marginalize out ${\bf x}_{1:M\setminus\\{m,m+1\\}}=\\{{\bf
x}_{1},\ldots,{\bf x}_{m-1},{\bf x}_{m+2},\ldots,{\bf x}_{M}\\}$. Note that
since ${\bf x}^{*}$ does not appear in the likelihood, we can take it out from
the integral,

$\displaystyle p({\bf x}_{m},{\bf x}_{m+1},{\bf x}^{*},\mathcal{D})$
$\displaystyle=\int p({\bf x}_{1})\prod_{j=1}^{m-1}p({\bf x}_{j+1}|{\bf
x}_{j})\prod_{j={m+1}}^{M}p({\bf x}_{j+1}|{\bf x}_{j})\cdot p(\mathcal{D}|{\bf
x}_{1},\ldots,{\bf x}_{M}){\rm d}{\bf x}_{1:M\setminus\\{m,m+1\\}}$
$\displaystyle\cdot p({\bf x}^{*}|{\bf x}_{m})p({\bf x}_{m+1}|{\bf x}^{*})$
$\displaystyle=\frac{p({\bf x}_{m},{\bf x}_{m+1},\mathcal{D})p({\bf
x}^{*}|{\bf x}_{m})p({\bf x}_{m+1}|{\bf x}^{*})}{p({\bf x}_{m+1}|{\bf
x}_{m})}.$ (37)

Therefore, we have

$\displaystyle p({\bf x}_{m},{\bf x}_{m+1},{\bf x}^{*}|\mathcal{D})\propto
p({\bf x}_{m},{\bf x}_{m+1}|\mathcal{D})p({\bf x}^{*}|{\bf x}_{m})p({\bf
x}_{m+1}|{\bf x}^{*}).$ (38)

Suppose we are able to obtain $p({\bf x}_{m},{\bf x}_{m+1}|\mathcal{D})\approx
q({\bf x}_{m},{\bf x}_{m+1})$. We now need to obtain the posterior of ${\bf
x}^{*}$. In the LTI SDE model, we know that the state transition is a Gaussian
jump. Let us denote

$p({\bf x}^{*}|{\bf x}_{m})=\mathcal{N}({\bf x}^{*}|{\bf A}_{1}{\bf
x}_{m},{\bf Q}_{1}),\;\;p({\bf x}_{m+1}|{\bf x}_{*})=\mathcal{N}({\bf
x}_{m+1}|{\bf A}_{2}{\bf x}^{*},{\bf Q}_{2}).$

We can simply merge the natural parameters of the two Gaussian and obtain

$\displaystyle p({\bf x}_{m},{\bf x}_{m+1},{\bf x}^{*}|\mathcal{D})=p({\bf
x}_{m},{\bf x}_{m+1}|\mathcal{D})\mathcal{N}({\bf
x}^{*}|\mathbf{m}^{*},\mathbf{V}^{*}),$ (39)

where

$\displaystyle\left(\mathbf{V}^{*}\right)^{-1}$ $\displaystyle={\bf
Q}_{1}^{-1}+{\bf A}_{2}^{\top}{\bf Q}_{2}^{-1}{\bf A}_{2},$
$\displaystyle\left(\mathbf{V}^{*}\right)^{-1}\mathbf{m}^{*}$
$\displaystyle={\bf Q}_{1}^{-1}{\bf A}_{1}{\bf x}_{m}+{\bf A}_{2}^{\top}{\bf
Q}_{2}^{-1}{\bf x}_{m+1}.$ (40)

##  Appendix B More experiments results

The prediction results on $R=\\{3,5,7\\}$ BeijingAir-PM2.5, BeijingAir-PM10,
and BeijingAir-SO2 are list in Tables 3 Tables 4,Tables 5

| RMSE | MAE
---|---|---
Datasets | PM2.5 | PM10 | SO2 | PM2.5 | PM10 | SO2
Discrete resolution:$50\times 50\times 150$ |  |  |
P-Tucker | $0.812\pm 0.054$ | $0.779\pm 0.015$ | $0.668\pm 0.015$ | $0.566\pm 0.018$ | $0.56\pm 0.013$ | $0.423\pm 0.005$
Tucker-ALS | $1.063\pm 0.049$ | $1.036\pm 0.062$ | $0.965\pm 0.029$ | $0.727\pm 0.021$ | $0.738\pm 0.022$ | $0.628\pm 0.01$
Tucker-SVI | $0.758\pm 0.017$ | $0.8\pm 0.051$ | $0.652\pm 0.014$ | $0.554\pm 0.014$ | $0.583\pm 0.017$ | $0.421\pm 0.038$
Discrete resolution:$100\times 100\times 300$ |  |  |
P-Tucker | $0.794\pm 0.099$ | $0.767\pm 0.005$ | $0.683\pm 0.043$ | $0.486\pm 0.009$ | $0.512\pm 0.014$ | $0.402\pm 0.015$
Tucker-ALS | $1.02\pm 0.032$ | $1.011\pm 0.023$ | $0.97\pm 0.027$ | $0.738\pm 0.011$ | $0.751\pm 0.004$ | $0.681\pm 0.016$
Tucker-SVI | $0.681\pm 0.027$ | $0.741\pm 0.071$ | $0.714\pm 0.128$ | $0.468\pm 0.014$ | $0.526\pm 0.043$ | $0.421\pm 0.038$
Discrete Resolution: $300\times 300\times 1000$ |  |  |
P-Tucker | $1.493\pm 0.125$ | $1.439\pm 0.001$ | $1.264\pm 0.2$ | $0.532\pm 0.034$ | $0.575\pm 0.001$ | $0.46\pm 0.004$
Tucker-ALS | $1.027\pm 0.033$ | $1.027\pm 0.038$ | $1.007\pm 0.017$ | $0.743\pm 0.01$ | $0.758\pm 0.012$ | $0.699\pm 0.007$
Tucker-SVI | $1.657\pm 0.135$ | $1.408\pm 0.052$ | $1.451\pm 0.062$ | $0.79\pm 0.016$ | $0.768\pm 0.016$ | $0.732\pm 0.016$
Discrete Resolution:$428\times 501\times 1461$ (original) |  |  |
P-Tucker | $2.091\pm 0.122$ | $2.316\pm 0.001$ | $1.48\pm 0.1$ | $0.756\pm 0.002$ | $0.79\pm 0.001$ | $0.556\pm 0.017$
Tucker-ALS | $1.008\pm 0.013$ | $1.027\pm 0.036$ | $1\pm 0.023$ | $0.738\pm 0.005$ | $0.757\pm 0.009$ | $0.699\pm 0.01$
Tucker-SVI | $1.864\pm 0.03$ | $1.686\pm 0.061$ | $1.537\pm 0.121$ | $0.885\pm 0.016$ | $0.864\pm 0.015$ | $0.787\pm 0.032$
methods using continuous indexes |  |  |
FTT-ALS | $1.019\pm 0.013$ | $1.001\pm 0.013$ | $1.002\pm 0.026$ | $0.744\pm 0.007$ | $0.755\pm 0.007$ | $0.696\pm 0.011$
FTT-ANOVA | $2.151\pm 0.032$ | $2.006\pm 0.015$ | $1.987\pm 0.036$ | $1.788\pm 0.031$ | $1.623\pm 0.014$ | $1.499\pm 0.018$
FTT-cross | $0.943\pm 0.026$ | $0.933\pm 0.012$ | $0.845\pm 0.026$ | $0.566\pm 0.018$ | $0.561\pm 0.011$ | $0.467\pm 0.033$
RBF-SVM | $0.995\pm 0.015$ | $0.955\pm 0.02$ | $0.794\pm 0.026$ | $0.668\pm 0.008$ | $0.674\pm 0.014$ | $0.486\pm 0.026$
BLR | $0.998\pm 0.013$ | $0.977\pm 0.014$ | $0.837\pm 0.021$ | $0.736\pm 0.007$ | $0.739\pm 0.008$ | $0.573\pm 0.009$
FunBaT-CP | $0.294\pm 0.016$ | $\mathbb{0.347\pm 0.036}$ | $\mathbb{0.384\pm 0.01}$ | $\mathbb{0.183\pm 0.006}$ | $0.236\pm 0.014$ | $0.242\pm 0.003$
FunBaT | $\mathbb{0.291\pm 0.017}$ | $0.348\pm 0.036$ | $0.386\pm 0.011$ | $\mathbb{0.183\pm 0.01}$ | $\mathbb{0.233\pm 0.012}$ | $\mathbb{0.241\pm 0.004}$
Table 3: Prediction error over BeijingAir-PM2.5, BeijingAir-PM10, and BeijingAir-SO2 with $R=3$, which were averaged over five runs.  | RMSE | MAE
---|---|---
Datasets | PM2.5 | PM10 | SO2 | PM2.5 | PM10 | SO2
Discrete resolution:$50\times 50\times 150$ |  |  |
P-Tucker | $0.835\pm 0.078$ | $0.787\pm 0.077$ | $0.745\pm 0.046$ | $0.54\pm 0.007$ | $0.535\pm 0.007$ | $0.424\pm 0.004$
Tucker-ALS | $1.178\pm 0.055$ | $1.123\pm 0.033$ | $0.975\pm 0.029$ | $0.741\pm 0.014$ | $0.741\pm 0.009$ | $0.615\pm 0.006$
Tucker-SVI | $0.738\pm 0.022$ | $0.747\pm 0.009$ | $0.638\pm 0.01$ | $0.518\pm 0.011$ | $0.541\pm 0.007$ | $0.405\pm 0.022$
Discrete resolution:$100\times 100\times 300$ |  |  |
P-Tucker | $0.808\pm 0.065$ | $0.827\pm 0.012$ | $0.763\pm 0.023$ | $0.474\pm 0.014$ | $0.494\pm 0.008$ | $0.426\pm 0.005$
Tucker-ALS | $1.07\pm 0.03$ | $1.038\pm 0.022$ | $0.953\pm 0.015$ | $0.745\pm 0.01$ | $0.756\pm 0.008$ | $0.675\pm 0.007$
Tucker-SVI | $0.768\pm 0.105$ | $0.79\pm 0.085$ | $0.691\pm 0.087$ | $0.471\pm 0.038$ | $0.524\pm 0.035$ | $0.405\pm 0.022$
Discrete Resolution: $300\times 300\times 1000$ |  |  |
P-Tucker | $2.153\pm 0.271$ | $1.972\pm 0.001$ | $1.486\pm 0.054$ | $0.784\pm 0.054$ | $0.859\pm 0.001$ | $0.624\pm 0.02$
Tucker-ALS | $1.062\pm 0.031$ | $1.046\pm 0.029$ | $1.007\pm 0.02$ | $0.747\pm 0.01$ | $0.76\pm 0.01$ | $0.699\pm 0.007$
Tucker-SVI | $1.584\pm 0.092$ | $1.446\pm 0.035$ | $1.511\pm 0.065$ | $0.828\pm 0.037$ | $0.805\pm 0.012$ | $0.781\pm 0.026$
Discrete Resolution:$428\times 501\times 1461$ (original) |  |  |
P-Tucker | $2.359\pm 0.078$ | $2.426\pm 0.001$ | $1.881\pm 0.054$ | $1.011\pm 0.021$ | $1.094\pm 0.001$ | $0.775\pm 0.017$
Tucker-ALS | $1.045\pm 0.02$ | $1.056\pm 0.021$ | $0.999\pm 0.025$ | $0.74\pm 0.005$ | $0.761\pm 0.007$ | $0.698\pm 0.01$
Tucker-SVI | $1.574\pm 0.088$ | $1.603\pm 0.024$ | $1.536\pm 0.05$ | $0.842\pm 0.026$ | $0.879\pm 0.008$ | $0.812\pm 0.018$
methods using continuous indexes |  |  |
FTT-ALS | $1.019\pm 0.013$ | $1.000\pm 0.013$ | $1.001\pm 0.026$ | $0.744\pm 0.007$ | $0.754\pm 0.005$ | $0.695\pm 0.010$
FTT-ANOVA | $2.149\pm 0.033$ | $2.006\pm 0.014$ | $1.987\pm 0.036$ | $1.788\pm 0.031$ | $1.623\pm 0.014$ | $1.499\pm 0.018$
FTT-cross | $0.941\pm 0.024$ | $0.933\pm 0.012$ | $0.831\pm 0.015$ | $0.563\pm 0.018$ | $0.560\pm 0.011$ | $0.464\pm 0.033$
RBF-SVM | $0.995\pm 0.015$ | $0.955\pm 0.02$ | $0.794\pm 0.026$ | $0.668\pm 0.008$ | $0.674\pm 0.014$ | $0.486\pm 0.026$
BLR | $0.998\pm 0.013$ | $0.977\pm 0.014$ | $0.837\pm 0.021$ | $0.736\pm 0.007$ | $0.739\pm 0.008$ | $0.573\pm 0.009$
FunBaT-CP | $0.292\pm 0.013$ | $0.352\pm 0.035$ | $\mathbb{0.385\pm 0.009}$ | $\mathbb{0.183\pm 0.007}$ | $0.236\pm 0.013$ | $0.243\pm 0.003$
FunBaT | $\mathbb{0.288\pm 0.012}$ | $\mathbb{0.338\pm 0.03}$ | $0.388\pm 0.003$ | $0.191\pm 0.021$ | $\mathbb{0.231\pm 0.005}$ | $\mathbb{0.241\pm 0.003}$
Table 4: Prediction error over BeijingAir-PM2.5, BeijingAir-PM10, and BeijingAir-SO2 with $R=5$, which were averaged over five runs.  | RMSE | MAE
---|---|---
Datasets | PM2.5 | PM10 | SO2 | PM2.5 | PM10 | SO2
Discrete resolution:$50\times 50\times 150$ |  |  |
P-Tucker | $0.832\pm 0.037$ | $0.823\pm 0.112$ | $0.811\pm 0.058$ | $0.521\pm 0.005$ | $0.529\pm 0.009$ | $0.425\pm 0.005$
GPTF | $0.694\pm 0.013$ | $0.702\pm 0.005$ | $0.607\pm 0.016$ | $0.488\pm 0.012$ | $0.501\pm 0.005$ | $0.381\pm 0.009$
CP-ALS | $1.385\pm 0.209$ | $1.092\pm 0.04$ | $1.016\pm 0.087$ | $0.804\pm 0.048$ | $0.751\pm 0.013$ | $0.635\pm 0.031$
Tucker-ALS | $1.379\pm 0.046$ | $1.314\pm 0.047$ | $1.049\pm 0.025$ | $0.779\pm 0.014$ | $0.775\pm 0.015$ | $0.621\pm 0.01$
Tucker-SVI | $0.742\pm 0.023$ | $0.721\pm 0.01$ | $0.67\pm 0.027$ | $0.504\pm 0.014$ | $0.508\pm 0.004$ | $0.408\pm 0.011$
Discrete resolution:$100\times 100\times 300$ |  |  |
P-Tucker | $1.015\pm 0.063$ | $0.911\pm 0.063$ | $0.877\pm 0.039$ | $0.51\pm 0.012$ | $0.521\pm 0.018$ | $0.452\pm 0.009$
Tucker-ALS | $1.086\pm 0.041$ | $1.047\pm 0.015$ | $0.967\pm 0.029$ | $0.745\pm 0.013$ | $0.755\pm 0.007$ | $0.682\pm 0.011$
Tucker-SVI | $0.783\pm 0.054$ | $0.787\pm 0.052$ | $0.702\pm 0.054$ | $0.464\pm 0.017$ | $0.515\pm 0.021$ | $0.408\pm 0.011$
Discrete Resolution: $300\times 300\times 1000$ |  |  |
P-Tucker | $1.718\pm 0.155$ | $1.928\pm 0.001$ | $1.629\pm 0.05$ | $0.805\pm 0.043$ | $0.954\pm 0.001$ | $0.687\pm 0.006$
Tucker-ALS | $1.073\pm 0.039$ | $1.062\pm 0.02$ | $1.007\pm 0.018$ | $0.748\pm 0.01$ | $0.762\pm 0.01$ | $0.699\pm 0.006$
Tucker-SVI | $1.437\pm 0.051$ | $1.499\pm 0.027$ | $1.389\pm 0.042$ | $0.793\pm 0.032$ | $0.842\pm 0.013$ | $0.772\pm 0.018$
Discrete Resolution:$428\times 501\times 1461$ (original) |  |  |
P-Tucker | $2.134\pm 0.174$ | $2.483\pm 0.001$ | $2.001\pm 0.149$ | $0.924\pm 0.044$ | $1.055\pm 0.001$ | $0.786\pm 0.012$
Tucker-ALS | $1.051\pm 0.02$ | $1.054\pm 0.029$ | $1\pm 0.021$ | $0.741\pm 0.005$ | $0.76\pm 0.008$ | $0.698\pm 0.01$
Tucker-SVI | $1.292\pm 0.021$ | $1.521\pm 0.095$ | $1.454\pm 0.054$ | $0.737\pm 0.015$ | $0.872\pm 0.027$ | $0.817\pm 0.028$
Methods using continuous indexes |  |  |
FTT-ALS | $1.019\pm 0.013$ | $1.000\pm 0.013$ | $1.001\pm 0.026$ | $0.744\pm 0.007$ | $0.754\pm 0.005$ | $0.695\pm 0.010$
FTT-ANOVA | $2.149\pm 0.033$ | $2.006\pm 0.014$ | $1.987\pm 0.036$ | $1.788\pm 0.031$ | $1.623\pm 0.014$ | $1.499\pm 0.018$
FTT-cross | $0.941\pm 0.024$ | $0.933\pm 0.012$ | $0.831\pm 0.015$ | $0.563\pm 0.018$ | $0.560\pm 0.011$ | $0.464\pm 0.033$
RBF-SVM | $0.995\pm 0.015$ | $0.955\pm 0.02$ | $0.794\pm 0.026$ | $0.668\pm 0.008$ | $0.674\pm 0.014$ | $0.486\pm 0.026$
BL | $0.998\pm 0.013$ | $0.977\pm 0.014$ | $0.837\pm 0.021$ | $0.736\pm 0.007$ | $0.739\pm 0.008$ | $0.573\pm 0.009$
PCMT-CP | $\mathbb{0.296\pm 0.023}$ | $\mathbb{0.335\pm 0.02}$ | $\mathbb{0.385\pm 0.008}$ | $\mathbb{0.184\pm 0.009}$ | $\mathbb{0.231\pm 0.007}$ | $\mathbb{0.242\pm 0.003}$
PCMT-Tucker | $0.318\pm 0.014$ | $0.347\pm 0.017$ | $0.39\pm 0.009$ | $0.198\pm 0.008$ | $0.235\pm 0.008$ | $\mathbb{0.242\pm 0.003}$
Table 5: Prediction error over BeijingAir-PM2.5, BeijingAir-PM10, and
BeijingAir-SO2 with $R=7$, which were averaged over five runs.