Streaming Factor Trajectory Learning for Temporal Tensor Decomposition
Authors: Shikai Fang
Xin Yu
Shibo Li
Zheng Wang
Robert Kirby
Shandian Zhe
Venue: NeurIPS 2023 poster (2023)
Keywords: Tensor Decomposition
streaming method
Bayesian model

Abstract: Practical tensor data is often along with time information. Most existing temporal decomposition approaches estimate a set of fixed factors for the objects in each tensor mode, and hence cannot capture the temporal evolution of the objects' representation. More important, we lack an effective approach to capture such evolution from streaming data, which is common in real-world applications.  To address these issues, we propose Streaming Factor Trajectory Learning (SFTL) for temporal tensor decomposition. We use Gaussian processes (GPs) to model the trajectory of  factors so as to flexibly estimate their temporal evolution. To address the computational challenges in handling streaming data, we convert the GPs into a state-space prior by constructing an equivalent stochastic differential equation (SDE).  We develop an efficient online filtering algorithm to estimate a decoupled running posterior of the involved factor states upon receiving new data. The decoupled estimation enables us to conduct standard Rauch-Tung-Striebel smoothing to compute the full posterior of all the  trajectories in parallel, without the need for revisiting any previous data. We have shown the advantage of SFTL in both synthetic tasks and real-world applications.

Streaming Factor Trajectory Learning for Temporal Tensor Decomposition
Shikai Fang
Xin Yu
Shibo Li
Kahlert School of Computing The University of Utah shikai fangautah. edu
Kahlert School of Computing The University of Utah yuxwinde gma 1. com
Kahlert School of Computing The University Utah shiboli i cs cs@gmailllll com
Zheng Wang
Robert M. Kirby Kahlert School of Computing The University of Utah ki kirbyêcs. utah. edu
Shandian Zhe
Kahlert School Computing The University of Utah u120884708 a0t.4hh.hhh....mm... edu
Kahlert School Computing The University of Utah zhe@cs utah.edu
Abstract
Practical tensor data is often along with time information. Most existing tem- poral decomposition approaches estimate set of fixed factors for the objects in each tensor mode, and hence canmot capture the temporal evolution of the objects' representation. More important, lack effective approach to capture such evolution from streaming data, which is common in weal-world applications. To address these issues, propose Streaming Factor Trajectory Learning (SFTL) for temporal tensor decomposition. We use Gaussian processes (GPs) to model the trajectory of factors as flexibly estimate their temporal evolution. To address the computational challenges handling streaming data, we convert the GPs into state-space prior by constructing an equivalent stochastic differential equation (SDE). We develop an efficient online filtering algorithm to estimate a decoupled running posterior of the involved factor states upon receiving new data. The decou- pled estimation enables us to conduct standard Rauch-Tung Turne-ttreeell smoothing compute the full posterior of all the trajectories in parallel, without the need for visiting any previous data. We have shown the advantage of SFTL in both synthetic tasks and real-world applications. The code available at httttss////gie com/ xuangu- fang/ Stramiieinnn Facttorrir--aaoooooooooooooooo Tra Traeooorrrrreeetror.. Learning
Introduction
Tensor data is common in real-world applications. For example, (patient, drug, clinic) from medical service records and four-mode tensor (customer, commodity, seller, web-page) from the database of an online shopping platform. Tensor decomposition fundamental tool for tensor data analysis. It introduces a set of factors to represent the objects in each mode, and estimate these factors by reconstructing the observed entry values. These factors can be viewed as the underlying properties of the objects. We can use them to search for interesting structures within the objects (e.g., communities and outliers) or as discriminate features for predictive tasks, such personalized treatment or recommendation.
extract three-mode tensor
Real-world tensor data is often accompanied with time information, namely the timestamps which the objects of different modes interact to produce the entry values. Underlying the timestamps can be rich, valuable temporal patterns. While many temporal decomposition methods are available,
Corresponding author.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

most of them estimate set of static factors for each object - they either introduce a discrete time mode [Xiong et al. 2010 Zhe al. 2016a| or inject the timestamp into the decomposition model [Zhang et al. 202 Fang et al. 2022 Li et al. 2022]. Hence, these methods cannot learn the temporal variation of the factors. Accordingly, they can miss important evolution of the objects' inner properties, such as health and income. In addition, practical applications produce data streams at a rapid pace [Du et al., 2018). Due to the resource limit, it is often prohibitively expensive to decompose the entire tensor from scratch whenever we receive new data. Many trytecyproteccing applications (e.g, SnapChat) even forbid us to preserve or re-access the previous data. Therefore, not only do we need more powerful decomposition model that estimate the factor evolution, also need an effective method to capture such evolution from fast data streams.
To address these issues, we propose SFTL, a Bayesian streaming factor trajectory learning approach for temporal tensor decomposition. Our method can efficiently handle data streams, with which to estimate posterior distribution of the factor trajectory to uncover the temporal evolution of the objects' representation. Our method never needs to keep or re-visit previous data. The contribution of our work is summarized as follows.
First, we use Gaussian process (GP) prior to sample the factor trajectory of each object as a function of time. As a nonparametric function prior, GPs are flexible enough to capture variety of complex temporal dynamics. The trajectories are combined through
a CANDECOMPPPARAFAC (CP) [Harshman 1970] or Tucker decomposition [Tucker]
[966] form sample the tensor entry values at any time point.
Second, to sidestep the expensive covariance matrix computation in the GP, which further causes challenges in streaming inference, we use spectral analysis to convert the GP into a linear time invariant (LTT) stochastic differential equation (SDE). Then we convert the SDE into an equivalent state-space prior over the factor (trajectory) states at the observed timestamps. As a result, the posterior inference becomes easier and computationally efficient. Third, we take advantage of the chain structure of the state-space prior and use the recent conditional expectation propagation framework (Wang and Zhe 2019 to develop an efficient online filtering algorithm. Whenever a collection of entries at a new timestamp arrives, method can efficiently estimate decoupled running posterior of the involved factor states, with a Gaussian product form. The decoupled Gaussian estimate enables us to run standard RTS smoothing [Särkkä 2013] to compute the full posterior of each factor trajectory independently and in parallel, without revisiting any previous data. Our method at worst has linear scalability with respect to the number of observed timestamps.
We first evaluated our method in simulation study. On synthetic datasets, SFTL successfully recovered several nonlinear factor trajectories, and provided reasonable uncertainty estimation. We then tested our method on four weal-world temporal tensor datasets for missing value prediction. In both online and final predictive performance, SFTL consistently outperforms the state-of-the-art streaming CP and Tucker decomposition algorithms by large margin. In most cases, the prediction accuracy of SFTL is even higher than the recent static decomposition methods, which have to pass through the dataset many times. Finally, we investigated the learned factor trajectories from a real-world dataset. The trajectories exhibit interesting temporal evolution.
2 Preliminaries
Tensor Decomposition. We denote an M-mode tensor by y E Rdix--xdas, where each mode m has dm dimensions, corresponding to dm objects. Each tensor entry is indexed by a tuple f= (L1y...,(N)) (M), and the value is denoted by ye- For decomposition, we introduce a set of latent factors un R.... to represent each object j in mode m (1 m M). One most pop- ular tensor decomposition model is the CANDECOMPIPARAFAC (CP) decomposition (Harsh RM R. and uses the following element- form. "yg where is the element-wist product. Another com- monly used model is Tucker decomposition [Tuckrr 1966), ye R vec(W) (n), .... @ uM) =
man 1970]. which sets R1
ge R (u), o...o C...ouM) = > II m 1 m
D ... ERM Wr IM "Y"m..... where W R x--xx RM is the tensor-core parameter, vec(-) is the vectorization, @ is the Kronecker product, and r (ri....,

---

Gaussian Process (GP)s are nonparametric function priors. For function f(x), if we place GP prior, f CP(O, <(x, x)), it means f(-) is sampled a realization of the GP with covariance function which is often chosen as a kernel function. The GP prior only models the correlation between the function values, namely, cov [fx)) = s((x, x'), and does not assume any parametric form of the function. Hence, GPs are flexible enough to estimate various complex functions from data, e.g., from multilinear to highly nonlinear. The finite projection of the GP is a Gaussian distribution. That is, given an arbitrary collection of inputs (xi,.... XN }. the corresponding function values f = ( (xx)) " follow a multi-variate Gaussian prior distribution, p(f) (f K) where is the covariance matrix and each [ mn
K: x(Xm, x,,
3 Bayesian Temporal Tensor Decomposition with Factor Trajectories
In real-world applications, tensor data is often associated with time information, namely, the times- tamps which the objects of different modes interact to generate the entry values. To capture the potential evolution of the objects' inner properties, we propose a Bayesian temporal tensor decomposi- tion model that can estimate trajectory of the factor representation. Specifically. for each object j in mode m, we model the factors as function of time, u [0, oc] + R². To flexibly capture variety of temporal evolution, we assign GP prior over each element of uf"
[uğ] (t), ...,.............. ", (t)] Le., uny (t) GP(O, mlt,, (')) (1 S S R). Given the factor trajectories, we then the CP or Tucker form to sample the entry values at different time points. For the CP form, we have p(ue(t))((t)) = (t) (u), ot).... o u (t),,-")
(I)
where u(t) (t) includes all the factor trajectories, and is the inverse noise variance, for which we assign a Gamma prior, p(7) Gam(rlag- 01). For the Tucker form, we have ((y)(t)(((()), W) ((ue(t)\ve(\W)) (u) (t) & ... w (-)) where we place standard normal prior the tensor-core, p(vec(W)) = NVvoc(W))), I). In this work, focus on continuous observations. is straightforward to extend our method for other types of observations.
Suppose we have a collection of observed entry values and timestamps, D ((l1,91,4) 41,.......... (EN: yN: [N]] where 4 S ... S Ix. We denote the sequence of timestamps when a particular object j of mode m participated in the observed entries by s711 < < Se
where oj is the participation count of the object. Note that it is a sub-sequence of {tz]. From the GP prior, the values of each u (t) at these timestamps follow a variate Gaussian distribution, p(ul)) (u,; |0, KTr where -h = [u]] (o..)................. ug (e T Faadd KT is the covariance/kermel matrix computed at these timestamps. The joint probability with the CP form is
pluueebr = II
m II II,,, Nyn ((ut.,, (l.) o... u. (Ln)),
II 1=1 Mu,
10, KY) Gam(rlao,
(2)
The joint probability with the Tucker form the same except that we use the Tucker likelihood instead and multiply with the prior of tensor-core p(W).
While this formulation is straightforward, it can introduce computational challenges. There are many multi-variate Gaussian distributions in the joint distribution (2), i.e., [N(u,) O,K)) The time and space complexity to compute each N ( 10, KJT) is O ((e))) and O ((cm) respectively. Nu
With the increase of N, the appearance count c for many objects can grow as well, making the computation cost very expensive even infeasible. The issue is particularly severe when handle streaming data - the number of timestamps grows rapidly when new data keeps coming in, so does the size of each covariance matrix.
3.1 Equivalent Modeling with State-Space Priors
To sidestep expensive covariance matrix computation and ease the inference with streaming data, we follow [Hartikainen and Särkkä 2010 to convert the GP prior into an SDE via spectral analysis.
For convenience, we abuse the notation a little bit: we denote by (.)(()) trajectory function and by ujir the values of the trajectory function at the observed timestamps.
3

---

We use a Matérn kernel Ky (1,14 =
where r(-) is the Gamma function, 0 0, K2 is the modified Bessel function of the second kind, and +
K v22
A = = -', Ep E (0.1.2.... ..... the GP covariance. Via the analysis of the power spectrum of Ky. we can show that if f(t) ~ GP(O, (t,))))), it can be characterized by linear time-invariant (LTI) SDE, with where f(k) - d /dt""
state z = (1) , ,....,(l))
dz - Az + 1 3(L); dt
(3)
where B(t) a white noise process with diffusion o².
A
1 Cp
0
Both o² and are obtained from the parameters in Ky, Due to the space limit, we leave the detailed derivation in Appendix (Section(A). The LTI-SDE is particularly useful in that its finite set of states follow a Gauss-Markov chain, i.e. the state-space prior. Given arbitrary 4 <... < Lt, we have
p(z(ti), (.).)),,.....)))).....). = p(z(ti)) I p(z((x+)))z(tt))
where p(z(t1)) (N(((t)) Px),, p(z((x+1) [z{+)) (2((tt++1 (tz(te), Q+). Px is the sta- tionary covariance matrix computed by solving the matrix Riccati equation [Lancaster and Rodman 1995], F, exp(Ax A) where A + t42 and Q P A&P cc Ag Therefore, not need the full covariance matrix as in the standard GP prior, and the computation is much more efficient. The chain structure also convenient handle streaming data as we will explain later.
do
We therefore convert the GP prior over each factor trajectory J.r (t) into an LTI-SDE. We de- note the corresponding state by (t). For example, if we choose = 1, then (t)r()) = [u] du j.r (t)/dt) For each object in mode m, we concatenate all its trajectory states into one, (t) = [<,] (t); :.si.....,,. (t)))...ii,(n((() Then on all of its timestamps sy1 < Ten we obtain
state-space prior
p(h = (h, l0, Po.), 11g4+11h33 = N(h
F m
(4)
j.k j,ky
where hm 2 (s Po
diag(Poo.. .o...... Poo), F diag
Fm j.k...
Fm j.k
Fm
m el j diag @T
and O7 = Poo
Fm Po
Fm
The joint probability of our model with the CP form now becomes
p(T) II M II p( k-
p(h
ph (33k)
II II N (ul,, (In) o... u. (L))).==")
(5)
Note that in the likelihood, each u. (tn)(1 < M) is contained in corresponding state vector hy such that symm, t, (by definition, we then have h z (t,,)). The joint probability with the Tucker form is similar, which we omit to save the space.
4 Trajectory Inference from Streaming Data
In this section, we develop an efficient, scalable algorithm for factor trajectory estimation from streaming data. In general, we assume that we receive a sequence of (small) batches observed tensor entries, (Bi, B2...... generated at different timestamps, {4, (4,,4,,...) Each batch B,, is generated timestamp l,, and l, < ln+1- Denote by D.,, all the data up to timestamp Ln Lee, Dr,, = B, U... U B, Upon receiving B,+1, intend to update our model without revisiting D,,,

---

Trajectory States State Transition
Old
New data
Factor Trajectory
n++1
θ n
Data Stream
But1
Dr, = B1 U
U UBn
t2
fn+1
Time
Figure 1: A graphical representation our factor trajectory learning, from which we can see (On+l, B are independent Dr, conditioned on G, and the noise inverse variance namely, 9nt1. Bu+1 1 Dr, 10,,T
provide the trajectory posterior estimate, (p(u ) 222,15 M, S dm
where Di..1 = Di, B
To this end, first observe that the standard state-space model with Gaussian likelihood has already provided a highly efficient streaming inference framework. Denote by x, and yn the state and observation each step n, respectively. To handle streaming observations (2,--2,---).
only need to compute and track the running posterior p(x,, y1.)) upon receiving each yx. where yin denotes the total data up to step n. This called Kalman filtering [Särkkä 2013), which only depends 1, i.e.,
on the running posterior at step processed (suppose it stops at step N), we can use Rauch-Tung-Sunn Stricbel (RTS) smoother [Särkkäl 2013] to efficiently compute the full posterior of each state, p pxx,, V1:N), ), from backward, which does not need re-access any previous observations (see Section Bin in Appendix).
and is highly efficient. After all the data
However, one cannot apply the above framework outright to our model, since the tensor decomposition likelihood of each observed entry couples the states of multiple factor trajectories, see ①1) which correspond to the participated objects at different modes. That means, the factor-state chains of different objects are dynamically intertwined through the received data. The multiplicative form of these states in the likelihood render the running posterior of each trajectory intractable to compute, mention running RTS smoother. To address this challenge, structure and use the recent conditional Expectation propagation (CEP) framework [Wang and Zhe 2019) to develop an efficient online filtering algorithm, which approximates the running posterior of the involved factor states product of Gaussian. Thereby, we can decouple the involved factor state chains, and conduct standard RTS smoothing for each chain independently.
take advantage of the chain
Specifically, denote the sequence timestamps when each object j of mode m has showed up in the data stream up to l,,, by 8111 8312
nm "iee where og,, is the object's appearance count up t. Denote by I the indexes all the objects of mode m appearing in B,, Hence, for every object j E IN we have m
J.ey, tn, and j.em (t,) the factor state of the object at t, Upon receiving each B,, we intend to approximate the running posterior of all the involved factor states and noise inverse variance 7 with the following decoupled form,
P(T, (h E I 1immsm D...) RS ((T(DD)) II M
m= II jern
(h?
jiem. (D,,),
(6)
where q( (t.)) = Gam(rlan, bn), and qha f.mm
m (D..) = Nhh P Vm Jem ). To this end, let consider given the approximation at ty, how obtain the new approximation at tn+1 (i.e., upon

---

receiving =n+ in the same form of 16. To simplify the notation, let us define the preceding states of the involved factors by θ, y.eg. J. n E 1 Im. and the current states by , 1 (h def li )m First, due to the chain structure of the prior over each (ey, 1k - 0, 2, we can see that conditioned on {θ,. 7) the current states @n+1 and the new observations B,+1 are independent of D... This is because in the graphical model representation, {O,, 7] have blocked all the paths from the old observations Dr,, to the new state and observations [Bishoop 2007]: see Fig. | for an illustration. Then, can derive that
p( On 1) & p(0n+ 2n, (P..) p(On, tet T, Dr,.) = p(,, r/Dr, where p p(O,,; (D,,)) ) is the running posterior ba>
p(On+1
II m==1 JEZ
II
phyyogr, Jo J.og.
is conditional Gaussian distribution defined in each p(h 20 S.ey.
hm
and
[e N' (2/1 (u), (fx+1) o...o u (fn+1)) )
Since
p(On, 7 r() takes the form of (6. we can analytically marginalize out each h and obtain
θ
II. M
p 7(Dmn)) Gam(T|an any bn)
N(h Megg
V.
(8)
jE7m
(e N (y/1 (ul, (tx+1) o ......
M EM (tn+1))
If we view the R.H.S of (8) as joint distribution with Ba+1, then our task amounts to estimating the posterior distribution, i.e. the I....S of (8) The product in the CP likelihood (and also Tucker likelihood) renders exact posterior computation infeasible, and we henceforth approximate
M
N [y/1 (u), (ta+1) o o.... u (tn++)))
x II N (ul", (fn+1)l7? Eg Gam (Tlae, (9) m=1
where X means approximately proportional to. To optimize these approximation terms, we use the recent conditional Expectation propagation (CEP) framework [Wang and Zhe| 2019] to develop an efficient inference algorithm. It uses conditional moment matching to update each approximation in parallel and conducts fixed point iterations, and hence can converge fast. We leave the details in the Appendix (Section[C) C). Once is done, we substitute the approximation ( into ( Then the R.H.S of (8) becomes a product of Gaussian and Gamma terms over each state and We can then immediately obtain a closed-fort estimation in the form as (6) At the beginning, when estimating
r(,,), since the preceding states On- 0, have an agg bn 01, Phop, 0, and
mm Ae
Poe in which is the prior of each 5.c%% and T (see 44)
In this way, we can continuously filter the incoming batches {B1, (3,, B2....). As result, for the factor state chain of every object in every mode m, along with each timestamp sm ".kk- we can online estimate and track running posterior approximation {g(h (D.,m, ){k = 1, 2....... which distribution. Hence, factor state, with which we can compute the posterior of the trajectory at any time point t [Bishop) 2006). Our method is summarized in Algorithm
Gaussian
the standard RTS smoother, to compute the full posterior of every
Algorithm Complexity. The time complexity of our algorithm processing a batch B, is ((BB,)))) where I. i the size. The time complexity of RTS smoother for particular object in mode m is OR cj,N). where N is the total number of timestamps. The space complexity of our algorithm is O M m== Σn tm R2), which track the running posterior of the factor state each appearing timestamp for every object. Since cN S N, the complexity of our algorithm is at worst linear in N.
5 Related Work
Many tensor decomposition methods have been developed, such as Yang and Dunson 2013 Rai
al. 2014
Zhe et al. 2015 2016b Tillinghast et al. 2020 Pan al. 2020b Fang et al. 202 602Lg

---

Algorithm 1 Streaming Factor Trajectory Learning (SFTL)
1: Input: kernel hyper-parameters a, P.
(p (0,1,2....)
2: n + 0
3: while Receiving new batch of entreis Ba+1 do
4: if' n = 0 then
5: Set an = ag, by bu: P Piegnn
0, and Vm J.
P in
6: Goto 9.
7: end if
8: Retrieve the involved preceding factor states On = { jeen lá E Iuti im and their running
posterior, p(O,,, (Dt,)) = Gam(rlan: bn.) II [Im= I Whh jnn Peeg
Vm h.gg )
9: According to 8 and D use conditional Expectation Propagation to calcu- late the running posterior of the current factor states, p(On+1, T(Dr.-1) Gam(TJ@n+1,0nnt
2
IsEzm #+1
N(h 1 Pee
Vm sieg
if Needed then
11: Run RTS smoothing on any factor state chain (baalk: = 1,2,...) of interest.
12: end if
13:
n++1
14:
end while
15: Run RTS smoothing for every factor state chain (b.))l: ...
Return: (9( [D) 2, ..... w.smmmssssssa (T( where D all the data received.
Tillinghast and Zhe 2021 Tillinghast et al., 2022 Fang al. 2022 Zhe and Du 2018 Pan et al.
2020a [ang et al. 2020 Pan et 2021 Wang et al. [022]. For temporal decomposition, most
existing methods augment the tensor with discrete time mode estimate additional factors for time steps, e.g. [Xiong et al. 2010 Rogers et al. 2013 Song et al. 2017 Du et al.
2018 Ahn
al. 2021 The most recent works have conducted continuous-time decomposition. Zhang al.
2021 used polynomial splines to model a time function as the CP coefficients, Li et al. 2022 used
neuralODE Chen al. [0188] model the entry value as a function of latent factors and time point. Fang et al. 2022) performed continuous-timg Tucker decomposition, and modeled the tensor-core as time function. To our knowledge, [Wang et al.] 2022) is the first work estimate factor trajectories. It places a GP prior in the frequency domain, and samples the factor trajectories via inverse Fourier transform. It then uses another GP sample the entry values. While successful, this method cannot handle streaming data, and the black-bos box GP decomposition lacks interpretability.
Current Bayesian streaming tensor decomposition methods include (Du et al. 2018 Fang et al. 202
Pan et al. 2020b Fang et al. 20216), which are based on streaming variational Bayes [Broderick 2013 assumed density filtering (ADF) (Boyen and Koller [998). ADF can be viewed as an instance of Expectation Propagation (EP) [Minka] 2001] for streaming data. EP approximates complex terms in the probability distribution with exponential-famill members, and uses moment matching to iteratively update the approximations, which essentially is a fixed point iteration. To address the challenge of intractable moment matching. Wang and Zhe [2019] proposed conditional EP (CEP), which uses conditional moment matching and Taylor expansion to compute the moments for factorized approximations. The theoretical guarantees and error bound analysis for EP and ADF have been studied for long time, such as [Boyen and Koller 1998 Dehaene and Barthelmé] 2015 2018). The most recent work [Fang et al., 2022 also uses SDEs to represent GPs and CEP framework for inference, but their GP prior is placed on the tensor-core, not for learning factor trajectories, and their method is only for static decomposition, and cannot handle streaming data.
6 Experiment
6.1 Simulation Study
We first conducted simulation study, for which we simulated a two-mode tensor, with two nodes per mode. Each node is represented by time-varying factor: u (t) sin sin? (2itt,, us(t)
(1 sin" sin'( (t)) sin siTt() uf(() sin(2 and w²(() cos sin (3πt) sin (2ntt.. Given
7

---

medl vajectory
Ground trub
1.0
1.0
1.0
Time
Time
Time
Time
(a) ui (f)
(b) (f)
(c) (t)
(d) nz((t)
Figure 2: The learned factor trajectories from the synthetic data. The shaded region indicates the posterior standard deviation.
RMSE
FirRecord
ServerRoom
BeijingAi-22
BeijingAir-3
PTucker
0,656 = 0.147
0.458 0.039
0.401 = 0.01
0.535 +: 0.062
Tucker-ALS
0.846 = 0.005
0.985 0.014
0.559 0.021
0.838 0.026
CP-ALS
0.882 = 0.017
0.994
0.801 0.082
0.875 0.028
CT-CP
0.664 = 0.007
0.384 +.8400
0.64 0.007
0.815 + 0.018
Statie
CT-GP
0.604 + 0.004
0.223 0.035
0.759 - 0.02
0.892 1: 0.026
BCTT
0.518 + 0.007
0.185 0.013
0.396 0.022
0.801 +: 0.02
NONFAT
0.503 +: 0.002
0.117 0.006
0.395 0.007
0.882 0.014
THIS-ODE
0.526 = 0.004
0.132 0.008
0.014
0.877 0.026
POST
0.696 + 0.019
0.516 0.028
0.658 0.103
ADF-CP
0.6-18 = 0.008
0.651
0.551 0.013
Stream
BASS-Tucker
0.976 + 0.024
1.000 0.016
1.099 0.037
0.991 0.039
SFTL-CP
0.424 = 0.014
0.161 0.014
0.248 0.012
0.473 0.003
SFTL-Tucker
0.430 = 0.010
0.331 + 0.056
0.303 = 0.041
0.439
MAE
PTucker
0.369 = 0.009
0.259 0.008
0.26 0.006
0.263 + 0.02
Tucker-ALS
0,615 = 0.006
0.739 + 0.008
0.388 1: 0.008
0.631 0.017
CP-ALS
0.642 0.012
0.746
0.586 0.056
0.605
CT-CP
0.46 0.004
0.008
0.489 0.006
0.62
Static
CT-GP
0.414 0.001
0.165 0.033
0.55+00012
0.626
BCTT
0.355 - 0.005
0.141 0.011
0.254 0.007
0.578 0.009
NONFAT
0.341 1 0.001
0.071 0.004
0.256 1 0.004
0.626 1.0007
THIS-ODE
0.363 0.004
0.083 0.002
0.345 0.004
0.605 1.0013
POST
0.478 = 0.014
0.4022
0.352 = 0.022
0.486 0.095
ADF-CP
0.449 + 0.006
0.496 1 0.007
0.385 0.012
0.409 + 0.029
Stream
BASS
0.772 + 0.031
0.749 +: 0.01
0.931 + 0.037
0.731 + 0.02
SFTL-CP
0.242 0.006
0.108
+ 0.008
0.15 + 0.003
0.318 0.008
SFTL-Tucker
0.246 = 0.001
0.216 0.034
0.185 0.029
0.278 0.011
Table 1: Final prediction error with R = 5. The results were averaged from five runs.
these factors, an entry value time t is generated via (t) ~ A (t) 0.05). We randomly sampled 500 (irregular) timestamps from (0, 1). For each timestamp, we randomly picked two entries, and sampled their values accordingly. Overall, we sampled 1 1,000 observed tensor entry values.
We implemented SFTL with PyTorch Paszke et al. 2019 We used
and a p 0.3 for the Matérn kernel. We streamed the sampled entries according to their timestamps, and ran our streaming factor trajectory inference based on the CP form. The estimated trajectories are shown in Fig. 2 As we can see, SFTL recovers the ground-truth pretty accurately, showing that SFTL has successfully captured the temporal evolution of the factor representation for every node. It interesting observe that when t is around 0, 0.5 and 1, the posterior standard deviation (the shaded region) increases significantly. This reasonable: the ground-truth trajectories overlap at these time points, making more difficult to differentiiateestiinall their values at these time points. Accordingly, the uncertainty of the estimation increases. In Section of Appendix, we further provide the pot-meen-ssqare error (RMSE) in recovering the four trajectories, and sensitivity analysis of the kernel parameters.

---

SPIL-CP
1.00
1.2
1,00
RMSF
1.50
Macke POS7
RMSE
RMSE
RMSE
1:00
0.70
0.45
0.25
0.45
20000
40000
4000
8000
8000
16000
600D
12000
Number of Entries
Number Entries
Number of Entries
Number of Entries
(a) FitRecord
(b) ServerRoom
(c)
BeijingAi-22
(d) BeijingAir-3
Figure 3: Online prediction error with the number of processed entries (R 5).
6.2 Real-World Applications
Next, examined SFTL in four weal-world datasets: FitRecord, ServerRoom, BeijingAir-2, and BeijingAir-3. We tested 11 competing approaches. We compared with state-of-the-an streaming tensor decomposition methods based the CP Tucker model, including (1) POST (Du et al. 2018], (2) BASS-Tuckee [Fang et al.) 2021a] and (3) ADF-CP [Wang and Zhe 2019]. the state-of- the-art static decomposition algorithms, including (4) P-Tucker [Oh et al. 2018), (5) CP-ASS and (6)
Tucker-ALS [Bader and Kolda (2008). For those methods, we augment the tensor with a time mode, and convert the ordered, unique timestamps into increasing time steps. We also compared with the most recent continuous-time decomposition methods. (7) CT-CP (Zhang et al. [021]. (8) CT-GP, (9)
BCTT [Fang al. 2022]. (10) THIS-ODE (Li al. 2022]. and (TT) NONFAT Wang et al. 2022).
nonparametric factor trajectory learning, the only existing work that also estimates factor trajectories for temporal tensor decomposition. Note that the methods 4-11 cannot handle data streams. They have to iteratively access the data to update the model parameters and factor estimates. The details about the competing methods and datasets are provided in Appendix (Section
For all the competing methods, we used the publicly released implementations of the original authors. The hyper-parameter setting and turning follows the original papers. For SFTL, we chose from & 2), from [0.5, I] and p from [0.1, 0.5]. For our online filtering, the maximum number of CEP iterations was set to 50 and the tolerance level to 10-4 For numerical stability, we re-scaled the timestamps to [0, 1]. We examined the number of factors (or factor trajectories)
{2, 3, 5, 7}.
Final Prediction Accuracy. We first examined the final prediction accuracy with our learned factor trajectories. To this end, we followed [Xu al. [0122] Kang al., 2012), and randomly sampled 80% observed entry values and their timestamps for streaming inference and then tested the prediction error the remaining entries. We also compared with the static decomposition methods, which need to repeatedly access the training entries. We repeated the experiment five times, and computed the average root mean-square-emes their standard deviations. We ran our method based on both the CP and Tucker forms, denoted by SFTL-CP and SFTL-Tucker, respectively. We report the results for R in Table IT Due to the space limit, we leave the other results in the Appendix (Table 45 and 66. As we can see, SFTL outperforms all the streaming approaches by large margin. SFTL even obtains significantly better prediction accuracy than all the static decomposition approaches, except that on Server Room, SFTL second to THIS-ODE and NONFAT. Note that SFTL only went through the training entries for once. Although NONFAT can also estimate the factor trajectories. it uses GPs to perform black-box nonlinear decomposition and hence loses the interpretability. Note that NONFAT most case also outperforms the other static decomposition methods that only estimate time-invariant factors. The superior performance of SFTL and NONFAT shows the importance of capturing factor evolution.
(RMSE), average rrean-absoluteerroo (MAE), and
Online Predictive Performance. Next, we evaluated the online predictive performance of SFTL.. Whenever batch of entries a new timestamp has been processed, we examined the prediction accuracy on the test set, with our current estimate of the factor trajectories. We repeated the evaluation for five times, and examine how the average prediction error varies along with the number of processed entries. We show the results for R = 5 in Fig. 3 and the others in Appendix (Fig. 6 and 7). It clear that SFTL in most cases outperforms the competing streaming decomposition algorithms by a large margin throughout the course of running. Note that the online behavior of BASS- BASS-Tucker was quite unstable and so we excluded trajectory learning approach - even in the streaming scenario, incrementally capturing the time- variation of the factors can perform better than updating fixed, static factors. To confirm the advantage of SFTL in computational efficiency, report the running time in Section G of Appendix.
the figures. It confirms the advantage of our streaming

---

0.5
0.0
0.5
0.0
-0.5
0.0
-1.0
-0.5
-0.5
Time
Time
Time
(a) (t)
(b) 21,2
(c) vi,a(t)
0.4
0.8
0.7
-0.5
0.0
0.0
00.5
-0.7
-1.4
Time
Time
Time
(d) (t)
(c) (²)2(t)
(f) u²,s(t)
2-
0
0
-1-
Time
Time
Time
(g) (t)
(h) (1)
(i v).
Figure 4: learned factor trajectories of object 1, 2, in mode 1, 2, 3, respectively, from ServerRoom data. The shaded region shows one standard deviation and the dashed line indicates the biggest timestamp in the data.
Investigation of Learning Results. Finally, we investigated our learned factor trajectories. We set R = 3 and ran SFTL-CP on ServerRoom. In Fig. 4 we visualize the trajectory for the 1st air conditioning mode, the 2rd power usage level, and the 3rd location. The shaded region indicates the standard deviation, and the red dashed line the biggest timestamp in the data. First, we can see that posterior variance of all the trajectories grow quickly when moving to the right of the red line. This is reasonable because it is getting far away from the training region. Second, for each object, the first and second trajectories (Ist and 2nd column in Fig. 4 exhibit quite different time-varying patterns, e.g., the local periodicity in 2t,2(t) and uz (t), the gradual decreasing and increasing trends in ui_ (4) and už (4), respectively, which imply different inner properties of the object. Third, is particularly interesting to see that the third trajectory for all the objects appear to be close to zero, with relatively large posterior variance all the time. This might imply that two factor trajectories have been sufficient to represent each object. Requesting for third trajectory is redundant. More important, our model empirically able detect such redundancy and returns zero-valued trajectory.
7 Conclusion
We have presented SFTL, a probabilistic temporal tensor decomposition approach. SFTL can efficiently handle streaming data, and estimate time-varying factor representations. On four real- world applications, SFTL achieves superior online and final prediction accuracy.
Acknowledgments
This work has been supported by MURT AFOSR grant FA9550-20- 1-0358 and NSF CAREER Award ISS-2442295
10

---

References
Dawon Ahn, Jun-Gi Jang, and U Kang. Time-aware tensor decomposition for sparse tensors. Machine Learning, pages 1-22, 2021.
Brett W Bader and Tamara G Kolda. Efficient matlab computations with sparse and factored tensors.
SIAM Journal on Scientific Computing, 30(1):205-233 2008.
Christopher M Bishop. Pattern recognition and machine learning. springer, 2006.
Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer, 2007.
Xavier Boyen and Daphne Koller. Tractable inference for complex stochastic processes. In Proceedings of the Fourteenth conference on Uncertainty artificial intelligence, pages 33-42, 1998.
Tamara Broderick, Nicholas Boyd, Andre Wibisono, Ashia C Wilson, and Michael I Jordan. Stream- ing variational bayes. Advances in neural information processing systems, 26, 2013.
Chris Chatfield. The analysis of time series: an introduction. Chapman and hall/CRC, 2003.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances neural information processing systems, 31, 2018.
Guillaume Dehaene and Simon Barthelmé. Expectation propagation the large data limit. Journal of the Royal Statistical Society. Series B (Statistical Methodology), 80(1):199-217, 2018.
Guillaume Dehaene and Simon Barthelmé. Bounding errors of pnpectaaioonpropggattiin Advances in Neural Information Processing Systems, 28, 2015.
Yishuai Du, Yimin Zheng, Kuang-chih Lee, and Shandian Zhe. Probabilistic streaming tensor decomposition. 2018 IEEE International Conference on Data Mining (ICDM), pages 99-108.
IEEE, 2018.
Shikai Fang, Robert M Kirby, and Shandian Zhe. Bayesian streaming sparse Tucker decomposition. In Uncertainty in Artificial Intelligence, pages 558-567. PMLR, 2021a
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.
Shikai Fang, Akil Narayan, Robert Kirby, and Shandian Zhe. Bayesian continuous-time Tucker decomposition. In International Conference on Machine Learning, pages 6245-.245. PMLR, 2022.
R. A. Harshman. Foundations of the PARAFAC procedure: Model and conditions for factor analysis. UCLA Working Papers in Phonetics, 16:1-84, 1970.
an" an" "m"eellantoor"" iotillanttr"mmlli-nonoe
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.
Kang, Evangelos Papalexakis, Abhay Harpale, and Christos Faloutsos. Gigatensor: scaling tensor analysis by 100 times- algoortths and discoveries. In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 316-324 ACM, 2012.
Peter Lancaster and Leiba Rodman. Algebraic riccati equations. Clarendon press, 1995.
Shibo Li, Robert Kirby, and Shandian Zhe. Decomposing temporal high-order interactions via latent odes. In International Conference Machine Learning, pages 12797-12812, PMLR, 2022.
Thomas P Minka. Expectation propagation for approximate bayesian inference. In Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence, pages 362-369, 2001.
Gary W Ochlert. A note the delta method. The American Statistician, 46(1):27-29, 1992.

---

Sejoon Oh, Namyong Park, Sael Lee, and Uksong Kang. Scalable tucker factorization for sparse tensors tensors-algorithne algorithms and discoveries. 2018 IEEE 34th International Conference on Data Engineering (ICDE), pages 11 1120-1131. IEEE, 2018.
Zhimeng Pan, Zheng Wang, and Shandian Zhe. Scalable nonparametric factorization for high-order interaction events. In International Conference on Artificial Intelligence and Statistics, pages 4325-4335. PMLR, 2020a.
Zhimeng Pan, Zheng Wang, and Shandian Zhe. Streaming nonlinear Bayesian tensor decomposition. In Conference on Uncertainty in Artificial Intelligence, pages 490-499. PMLR, 2020b.
Zhimeng Pan, Zheng Wang, Jeff M Phillips, and Shandian Zhe. Self-adaptable point processes with nonparametric time decays. Advances in Neural Information Processing Systems, 34:459446600 2021.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances neural information processing systems, 32:
8026-8037, 2019.
Piyush Rai, Yingjian Wang, Shengbo Guo, Gary Chen, David Dunson, and Lawrence Carin, Scalable Bayesian low-rank decomposition of incomplete multiway tensors. In Proceedings of the 31th International Conference on Machine Learning (ICML), 2014.
Mark Rogers, Lei Li, and Stuart J Russell. Multilinear dynamical systems for tensor time series. Advances in Neural Information Processing Systems, 26:2634-2642, 2013.
Simo Särkkä. Bayesian filtering and smoothing. Number 3. Cambridge University Press, 2013.
Qingquan Song, Xiao Huang, Hancheng Ge, James Caverlee, and Xia Hu. Multi-aspect stream- ing tensor completion. In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 435-443, 2017.
Conor Tillinghast and Shandian Zhe. Nonparametric decomposition of sparse tensors. In International Conference on Machine Learning, pages 10301-10311, PMLR, 2021.
Conor Tillinghast, Shikai Fang, Kai Zhang, and Shandian Zhe. Probabilistic neural-kernel tensor decomposition. 2020 IEEE International Conference on Data Mining (ICDM), pages 531-540.
IEEE, 2020.
Conor Tillinghast, Zheng Wang, and Shandian Zhe. Nonparametric sparse tensor factorization with hicrarchical Gamma processes. In International Conference on Machine Learning, pages 21432-21448. PMLR, 2022.
Ledyard Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31: 279-311, 1966.
Zheng Wang and Shandian Zhe. Conditional expectation propagation. In UAI, page 6, 2019.
Zheng Wang, Xinqi Chu, and Shandian Zhe. Seff-modulating nonparametric event-tensor factoriza- tion. In International Conference on Machine Learning, pages 9857-9867. PMLR, 2020.
Zheng Wang, Yiming Xu, Conor Tillinghast, Shibo Li, Akil Narayan, and Shandian Zhe. Non- parametric embeddings of sparse high-order interaction events. In International Conference on Machine Learning, pages 23237-23253. PMLR, 2022.
Liang Xiong, Xi Chen, Tzo-Kuo Huang, Jeff Schneider, and Jaime G Carbonell. Temporal collabora- tive filtering with bayesian probabilistic tensor factorization. In Proceedings of the 2010 SIAM International Conference on Data Mining, pages 211-222. SIAM, 2010.
Zenglin Xu, Feng Yan, and Yuan Alan Qi. Infinite tucker decomposition: Nonparametric bayesian models for multiway data analysis. In ICML, 2012.
Y. Yang and D.B. Dunson. Bayesian conditional tensor factorizations for high-dimensional classifica- tion. Journal of the Royal Statistical Society B, revision submitted, 2013,
12

---

Yanqing Zhang. Xuan Bi, Niansheng Tang, and Annie Qu. Dynamic tensor recommender systems. Journal of Machine Learning Research, 22(65):1-35, 2021.
Shandian Zhe and Yishuai Du. Stochastic nonparametric event-tensor decomposition. Advances in Neural Information Processing Systems, pages 6856-6866, 2018.
Shandian Zhe, Zenglin Xu, Xinqi Chu, Yuan Qi, and Youngja Park. Scalable nonparametric multiway data analysis. Artificial Intelligence and Statistics, pages 1 122-11134 PMLR, 2015.
Shandian Zhe, Yuan Qi, Youngja Park, Zenglin Xu, Ian Molloy, and Suresh Chari. Dintucker: Scaling up Gaussian process models on large multidimensional arrays. Proceedings of the AAAI Conference on Artificial Intelligence, volume 30, 2016a.
Shandian Zhe, Kai Zhang, Pengyuan Wang, Kuang-chib Lee, Zenglin Xu, Yuan Qi, and Zoubin Ghahramani. Distributed flexible nonlinear tensor factorization, Advances in neural information processing systems, 29, 2b.
13

---

Appendix
A Spectral Analysis and LTI-SDE
We consider the Matérn kernel family,
VBV A
v
2v
Ku (1,
((v)2u-1 Ku
(10)
where A - 1', parameters, respectively, K, is the modified Bessel function of the second kind, 0 controls the smoothness. Since Ky is stationary kernel, i.e., Ky(t, t') Ky(t - t'), according to the Wiener-Khinchin theorem [Chatfield 2003]. if
the Gamma function, a > and > are the amplitude and length-scale
(()) GP(O, (t,))))
the energy spectrum density of f(t) can be obtained by the Fourier transform of kv (A),,
S(w) = 2 2 + v) a² (a² w²) -(++))
(11)
((v)
where w is the frequency, and v/2v We consider the commonly used choice = p + where P
p E {0, 1. Then we can observe that
o²
S(w) = (a² + w2jpt
(12)
(a + iw) iw)ptI'
a 2p+1 and i indicates an imaginary number. We expand the polynomial where 0² = =a ((+t5)
(2+x((p+1)
(a ww
2 E (kw)))* (iw)++1
(13)
where fcx|0= < k & p) the coefficients. From and 13 construct equivalent system to generate the signal f()). That is, the frequency domain, the system output's Fourier transform ((w) is given by
Σn_ k== en (1w) f(w) (tw) [Bw),
(14)
where 3 is the Fourier transform white noise process B(t) with spectral density (or diffusion) B(w) 02. The reason that by construction, f(w) (a+iw))t1, which gives exactly the same spectral
density as in 12 S(w) = [f We then conduct inverse Fourier transform on both sides of [14] obtain the representation the time domain.
d dp+1f Ck
Σ
P
+ di+tI 3(1),
(15)
dik
which is an SDE. Note that 3(1) has the density o². We can further construct new state
(f. f() ((.)..)))))) f(r)) SDE,
(where each (k) A f/dtt) and convert 15 into a linear time-invariant (LTI)
dz
dt = Az + n. 3(t))
(16)
where
1
0
A
n
1
cp-1
Cp
1
14

---

For concrete example, if we take = 1 (and so v 2), then A = 10,1; -2l], 77 [0: 1]. and 6² 4aa
The LTI-SDE particularly useful in that its finite states follow a Gauss- chain, namely the state-space prior. Specifically, given arbitrary ε < < tr, we have
p(z(t)).... .....,)))))))))))))..)))))))))) - p(z(tz)) II
k= p(z((x+1))z((a))
where p(z(t1) (2((1)) Px)), p(z((x+11)))))))))))))))))))1))) = N((((t++ F (RZ( Qs), P the sta- tionary covariance matrix computed by solving the matrix Riccati equation [Lancaster and Rodman 1995]. F2 exp (Áx A) where As f l42 and Q2 = Po0 - P oc A Therefore, we do not need the full covariance matrix as in the standard GP prior, and the computation is much more efficient. The chain structure is also convenient to handle streaming data as explained in the main
paper.
Note that for other type kernel functions, such the square exponential (SE) kernel, we can approximate the inverse spectral density 1/S(w) with a polynomial of w with negative roots, and follow the same way to construct an LTT-SDE (approximation) and state-space prior.
B RTS Smoother
Consider a standard state-space model with state x,, and observation at each time step n. The prior distribution is Gauss-Mart Markov chain,
Xxx++ (xn) = (x AA,x: Qn),
(xo) (N Po).
Suppose we have a Gaussian observation likelihood,
ply, (xa) Ny H,, Wn).
Then upon receiving each Yns we can use Kalman filtering to obtain the exact running posterior,
p(xalyisn) = Nx,,, my, Pk),
which is Gaussian, After all the data has been processed suppose it ends after step N we can use Rauch- -ung- Lung-Strriebl (RTS) smoother [Särkkää 2013) to efficiently compute the full posterior of each state from backward, which does not need to re-access any data: p(x,, lYNN) = -N(,, m P... where
m A,
P n+1 A PrP.A. + Qn,
G, A, [P- PP+il= 1- 1-
m m, + Gn (m m m
P; = P, + G. [P R+1
G.
(17)
= (- P.) and As the computation only needs the running posterior p(x, p(x,/))m)) (m))= the full posterior of the next state p(x,+1\y1:)) = (- Pnt1). It does not need to revisit previous observations y1.N
C Details about Online Trajectory Inference
In this section, we provide the details about how to update the running posterior according to X and p (in the main paper) with the conditional EP (CEP) framework (Wang and Zhe 2019
C.1
EP and CEP framework
We first give brief introduction the EP and CEP framework. Consider general probabilistic model with latent parameters O. Given the observed data D (y1,...,ya .NN). the joint probability distribution is
N
p(0, D) = p(0) II
ply, 10)
(18)
15

---

Our goal is to compute the posterior D). However, it is usually infeasible to compute the exact marginal distribution P(D). because of the complexity of the likelihood and/or prior. EP therefore seeks to approximate each term in the joint probability by an exponential-famiyy term,
yly 0) (0), p(0) co),
(19)
where c, and oo are constants to ensure the normalization consistency (they will get canceled in the inference, so we do not need to calculate them), and
(0) exp(>,, ( (0 S a N)
where A, is the natural parameter and ( is sufficient statistics. For example, if we choose
Gaussian term, Jn pm E.), then the sufficient statistics is <(() (0, 00 The moment the expectation of the sufficient statistics.
We therefore approximate the joint probability with
p(0, D) p()) II ply,
fo II fo (0) const.
(20)
Because the exponential family is closed under product operations, we can immediately obtain closed-form approximate posterior (()) x p(0[D) by merging the approximation terms in the R.H.S of (2).. which still a distribution in the exponential family.
Then the task amounts to optimizing those approximation terms (fa = n S N). EP repeatedly conducts four steps to optimize each fn.
Step 1. We obtain the calibrated distribution that integrates the context information of fx.
(())
I (0)
where (0)) the current posterior approximation.
Step 2. We construct tilted distribution to combine the true likelihood.
( x n (0) :p(y,
Note that if n = 0, we have ((0) x ()) p()).
Step 3. We project the tilted distribution back to the exponential family,
q" (0) = argmin KL()]14
where belongs to the exponential family. This can be done by moment matching,
Eg. (x)) Ep(g(0))-
(21)
That is, we compute the expected moment under with which obtain the parameters of q². For example, if q" (0) is a Gaussian distribution, then we need to compute Ep(0]
and E 00 ]. with which to obtain the mean and covariance for (0)) Hence we obtain q" (0) (( Ep 0 Ep 0 - e0 E,p
Step 4. We update the approximation term by
(0) ( q (
(22)
In practice, EP often updates all the fa's in parallel, and uses damping avoid divergence. iteratively runs the four steps until convergence. In essence, this is fixed point iteration to optimize a free energy function (a mini-mat max problem) Minka 2001].
The critical step in EP the moment matching (21). However, in many cases, it is analytically intractable to compute the moment under the tilted distribution due to the complexity of the likelihood. To address this problem, CEP considers the commonly used case that each fa has a factorized structure,
(0)
II S Jm (0m),
(23)
16

---

where each from is also in the exponential family, and {0m} are mutually disjoint. Then at the moment matching step. we need compute the moment of each pn under p, i.e. Ep The first key idea of CEP is to use the nested structure,
Ep((((m)) Bp m. pe. (0.te..))..))(((...))) (o(0.)))
(24)
where pmm = 0m Therefore, we can first compute the inner expectation, i.e, conditional moment,
E _ ()(m)))))))))))m))))))))))))(( g (m))
(25)
and then seek for computing the outer expectation, E (g(( {m)) The inner expectation is often easy to compute (e.g. with our CP/Tucker likelihood). When fn factorized individually over each element of 0, this can always be efficiently and accurately calculated by quadrature. However, the outer expectation is still difficult to obtain because p( is intractable. The second key idea of CEP that since the moment matching also between q(0,,,. and p p00,mm we can use the current marginal posterior to approximate the marginal titled distribution and then compute the outer expectation,
Exco...) (8(0m)) 2 Eg(o,..)(8(0,.
(26)
If it is still analytically intractable, we can use the delta method [Ochlert 1992 to approximate the expectation. That is, we use Taylor expansion of g(-) at the mean of 0 0mm- Take the first-order expansion an example,
pm) 2 gg ( 0mll + (0,m
Eggo (0
where J is the Jacobian of g at (0 10 [m) Then we take the expectation on the Taylor approxi- mation instead,
E (..) (B)) ~ (l
(27)
The above computation is very convenient to implement. Once obtain the conditional moment g( we simply replace the 0 fmmm by its expectation under current posterior approximation q. i.e. with which to construct q" in Step
Ey m.) (0em [m)) to obtain the matched moment g(
of EP (see (2T)). The remaining steps are the same.
C.2 Running Posterior Update
Now we use the CEP framework to update the running posterior p(On+1,7|Pr,,, in via the approximation ( To simplify the notation, let us define vin. (e, v) E Ba41, we approximate
A u (tn+1). and hence for each
M ~ II N(V hi Er Gam((7Ja, mr
N 1 (vl, o... v vM
(28)
If we substitute 9 into we can immediately obtain a Gaussian posterior approximation of each vž,, and Gamma posterior approximation of the noise inverse variance т. Then dividing the current posterior approximation with the R.H.S of 28 we can obtain the calibrated distribution,
(v .Vv 1 m S2
(T) Gam(7 ( W)
(29)
where 1 5 m S M. Next. we construct a tilted distribution.
M II
p(v) Var x
) (T)
(v())))))))))) N (v/IT (vh o... v )
(30)
m=1
To update each N N(vm mym
mm Em
in 28
we first look into the conditional tilted distribution,
p mm T) N N(v 1.Th . 2 N (vl (v%%.
ve Am
-2)
(31)
17

---

where Vir is [v 11 S S M.j + m. and
l \m
= v
vmmtt
M Exr
fm-1
The conditional tilted distribution is obviously Gaussian, and the conditional moment is straightfor- ward to obtain,
mm T.
Siv 1V2 7) 2 -1 +
Am
Vi
(32)
B[v VV m -] =(vV" 12
:T).
( -1 B + tyve km
(33)
where S denotes the conditional covariance. Next, according to we simply replace T, Ve Am
and ve mm (ve" by their expectation under the current posterior q in 32 and (33) to obtain the moments, í.e., the mean and covariance matrix, with which we can construct q² in Step 3 of the EP framework. The computation of E straightforward, and
Egive" \m
B, (vlv o o.... Eglvin EgIven o o... Egg (arl
B.y lve mm (v \m
Egg [vb, (v2)) 1 Jo... o Byy (v (m=-1
Ey E
vmttt Jo... Bg Vvsr (VM)'I
Similarly, update Gam we) in also a Gamma distribution,
we first observe that the conditional titled distribution is
V() Gam a) Gam(r W N ve,
(34)
where V2 = o... vol and
1
/
+
2
1 3 we ve + 2"
T vev v1
(35)
2
Since the conditional moments (the expectation of T and log T) are functions of a and w, when using the delta method to approximate the expected conditional moment, it is equivalent to approximating the expectation of ā and first, and then use the expected ā and to recover the moments. As result, we simply replace Ve and vev. in (35) by their expectation under the current posterior, and we obtain the approximation construct q* (T) Gam(7 [,,, [a), [,( Step 3 in EP. The remaining steps are straightforward. The running posterior update with the Tucker form likelihood follows similar way.
Ey [x]] and Eg W] With these approximated expectation, we then
D More Results on Simulation Study
D.1 Accuracy of Trajectory Recovery
We provide the quantitative result in recovering the factor trajectories. Note that there is only one competing method, NONFAT, which can also estimate factor trajectories. We therefore ran method and NONFAT on the synthetic dataset. We then randomly sampled 500 time points in the domain and evaluate the RMSE of the learned factor trajectories for each method. As shown in 2 the RMSE of NONFAT on recovering (t) and (t( is close to SFTL, showing NONFAT achieved the same (or very close) quality in recovering these trajectories. However, ut (t) and (t), the RMSE NONEAT is much larger, showing that NONFAT have failed to capture the other two trajectories. By contrast, SFTL consistently well recovered them.
D.2 Sensitive Analysis Kernel Parameters
To examine the sensitivity to the kernel parameters, we used the synthetic dataset, and randomly sampled 100 entries and new timestamps for evaluation. We then examined the length-scale p and
18

---

(t)
u2(()
wf(()
u²(1)
SFTI.
0.073
0.082
0.103
0.054
NONFAT
0.085
0.442
0.096
0.443
Table RMSE recovering trajectories on the simulation data.
0.1
0.3
0.5
0.7
0.9
SFTL-CP
0.091
0.064
0.059
0,056
0,057
Matém- 1/2
SFTL-Tucker
0.060
0.055
0.056
0.056
0.057
SFTL-CP
0.062
0.061
0.074
0.093
0.112
Matémn-32
SFTL-Tucker
0.061
0.059
0.078
0.101
0.129
(a) Prediction RMSE with = 0.3 and varying p.
a
0.1
0.3
0.5
0.7
0.9
SFTT-CP
0.056
0.064
0.057
0.059
0.063
Matérn- 1/2
SFTL--Tucker
0.065
0.055
0.054
0.055
0.055
SFTL-CP
0.072
0.061
0.063
0.060
0.059
Matérn-3/2
SFTL-Tucker
0.098
0.059
0.064
0.062
0.061
(b) Prediction RMSE with = 0.3 and varying a.
Table 3: Sensitive analysis of amplitude a and length-scale p on synthetic data.
amplitude a, for two commonly-use used Matérn kernels: Matérn- 1/2 and Matérn-3/ The study was performed on SFTL based on both the CP and Tucker forms. The results are reported in Table 3 Overall, the predictive performance of SFTL is less sensitive to the amplitude parameter a than to the length-scale parameter p. But when we use Matérn-1/2 the performance of both SFTL-CP and SFTI.-Tucke quite stable to the length-scale parameter p. When we use Matérn-3/2, the choice of the length-scale is critical.
E Real-Worle Dataset Information and Competing Methods
We tested all the methods in the following four real-world datasets.
FitRecord workout logs of EndoMondo users' health status in outdoor exercises. We extracted a three-mode tensor among 500 users, 20 sports types, and 50 altitudes. The entry values are heart rates. There are 50K observed entry values along with the timestamps.
ServerRoom temperature logs of Poznan Supercomputing and Networking Center. We extracted a three-mode tensor between 3 air conditioning modes (24%, 27 and 30%), 3 power usage levels (50% 75%, 100%) and 34 locations. We collected 10K entry values and their timestamps.
BeijingAir-2 air pollution measurement in Beijing from year 2014 to 2017. We extracted a two-mode tensor (monitoring site, pollutant), size 12 6, and collected 20K observed entry values (concentration) and their timestamps.
BeijingAir-3, extracted from the same data source as BeijingAir-2, a threc-mode tensor among 12 monitoring sites, 12 wind speeds and 6 wind directions. The entry value the PM2.5 concentration. There are 15K observed entry values at different timestamps.
We first compared with the following state-of thite--o-thheatt art streaming tensor decomposition methods based on the CP or Tucker model. (1) POST (Du et al. 2018), probabilistic streaming CP decomposition via mean-field streaming variational Bayes [Broderick al., 2013 (2) MASS-Tuucker [Fang et al.. 2021a] Bayesian streaming Tucker decomposition, which online estimates sparse tensor-core via
spike-and-slab prior to enhance the interpretability. We also implemented (3) ADF-CP, streaming CP
http:////tee.oooo.cooc//////ee/r home
Air Quuaatttttaaaalll
19

---

decomposition by combining the assumed density filtering and conditional moment matching (Wang and Zhe 2019
Next, we tested the state-of-the-art static decomposition algorithms, which have to go through the data many times. (4) P-Tucker (Oh et al. 2018]. , an efficient Tucker decomposition algorithm that performs parallel row-wise updates. (5) CP-ALS and (6) Tucker-ALS [Bader and Kolda 2008)] CP/Tucker decomposition via alternating least square (ALS) updates. The methods (1-6) are not specifically designed for temporal decomposition and cannot utilize the timestamps of the observed entries. order incorporate the time information for a fair comparison, we augment the tensor with time mode, and convert the ordered, unique timestamps into increasing time steps.
We then compared with the most recent continuous continuous-time temporal decomposition methods. Note that none of these methods can handle data streams. They have to iteratively access the data to update the model parameters and factor estimates. (7) CT-CP [Zhang al.) [021], continuous- time CP decomposition, which uses polynomial splines to model a time-varying coefficient X for each latent factor, (8) CT-GP, continuous-time GP decomposition, which extends (Zhe al.) 2016al to use GPs to learn the tensor entry value as function of the Jatent factors and time (e(() q( ..... _ GP(O, k)) (9) BCTT [Fang et al. 2022 Bayesian continuous- time Tucker decomposition, which estimates the tensor-core as a time-vayiin function, (10) THIS- model the entry value as NN(u), :.... u t) where NN is short for neural networks. (11) NONFAT [Wang et al. 2022]. nonparametric factor trajectory learning, the only existing work that also estimates factor trajectories for temporal tensor decomposition. uses a bi-level GP to estimate the trajectories in the frequency domain and applies inverse Fourier transform
ODE [Li et al. 2022]. which uses neural ODE [Chen et al. 2018
function the latent factors and time, (ya(t)
di
to return to the time domain.
F More Results about Prediction Accuracy
We report for R = 2, R and R 7, the final prediction error (after the data has been processed)
of all the methods in Table 4 Table 5 and Table 6 respectively. We report for R = 2, R = 3 and
R 7, the online predictive performance of the streaming decomposition approaches in Fig. Fig.
6. and Fig. 7 respectively.
1.50
SITL-PP
1.00
1.00
RMSE
SPTL-lucker
RISSE
RMSE
RMSE
0.60
9.7
0.70
0.25
20000
40000
4000
8000
0.2
8000
16000
12000
Number of Entries
Number of Entries
Number of Entries
Number Entries
(a) FitRecord
(b) ServerRoom
(c)
BeijingAi-2
(d) BeijingAir-3
Figure 5: Online prediction error with number of processed entries (R 2)
1.50
1.2
1.00
1.00
RMSE
1.00
AMSE
0.60
RMSE
RMSE
0.70
20000 Number of Entries
40000
0.25
4000
8000
0.2
0
8000
16000
0.45
6
6000
12000
Number of Entries
Number of Entries
Number of Entries
(a) FitRecord
(b) ServerRoom
(c) BeijingAi-22
(d) BeijingAir3
Figure 6: Online prediction error with the number of processed entries (R 3)
G
Running Time
As compared with static (non-streaming) methods, such BCTT, our method is faster and more efficient. That because whenever new data comes in, the static methods have to retrain the model from scratch and iteratively access the whole data accumulated so far, while our method only performs
20

---

RMSE
FirRecord
ServerRoom
BeijingAi-22
BeijingAi--3
PTucker
0.606 = 0.015
0.757 + 0.36
0.509 = 0.01
0.442 + 0.142
Tucker-ALS
0.914 0.01
0.991 0.016
0.586 0.016
0.896 0.032
CP-ALS
0.926 = 0.013
0.997 0.016
0.647 0.041
0.918 0.031
CT-CP
0.675 = 0.009
0.412 0.022
0.642 + 0.007
0.832 + 0.035
Static
CT-GP
0.611 0.009
0.218
0.723 0.01
0.88 0.026
BCTT
0.604 0.019
0
0.504
0.799 0.027
NONFAT
0.541 = 0.002
0.132 0.002
0.425 0.002
0.878 0.014
THIS-ODE
0.544 # 0.005
0.142 + 0.004
0.553 + 0.015
0.876 0.027
POST
0.705 = 0.013
0.767 0.155
0.539 1.9:0..00 0.01
0.695 0.135
ADF-CP
0.669 = 0.033
0.764 0.114
0.583
Stream
BASS-Tucker
0.016
0.016
1.03 = 0.05
0.982 0.058
SFTL-CP
0.437 + 0.014
0.18 0.019
0.323 0.019
0.462 + 0.009
SFTL.-Tucker
0.446 + 0.024
0.276 + 0.031
0.344 0.031
0.417 + 0.035
MAE
PTucker
0.416 + 0.005
0.388 + 0.152
0.336 0.004
0.271 + 0.053
Tucker-ALS
0,676 + 0.008
0.744 0.01
0,408 1: 0.008
0.669 + 0.02
CP-ALS
0.686 = 0.011
0.454 0.057
0.001
CT-CP
0,466 + 0.005
0.295
0.49 0.006
0.642
Static
CT-GP
0.424 0.006
0.0012
0.517
0.626 0.01
BCTT
0.419 0.015
0.534
0.343 0.003
0.579 0.018
NONFAT
0.373 = 0.001
0.083 0.000
0.282 0.002
0.622 0.006
THIS ODE
0.377 = 0.003
0.097
0.355 0.008
0.606
POST
0.485 = 0.008
0.509
0.368 0.008
0.517
ADF-CP
0.462 = 0.022
0.574 0.073
0.401 - 0.029
0,415 + 0.038
Stream
BASS
0.777 0.039
0.749 + 0.01
0.871 0.125
0.727 + 0.029
SFTL-CP
0.248 0.005
0.126 0.007
0.199 0.005
0.311 + 0.004
SFTL--Tucker
0.25 + 0.01
0.203 0.032
0.218 0.02
0.261 0.023
Table 4: Final prediction error with R = 2. The results were averaged from five runs.
RMSE
FitRecord
ServerRoom
BeijingAi-22
BeijingAir-3
PTucker
0,603 = 0.045
0,677 0.129
0.464 0.012
0.421 0.074
Tucker-ALS
0,885 = 0.007
0.989 = 0.014
0.559 0.017
0.863 0.032
CP-ALS
0.907 = 0.015
0.993 = 0.016
0.594 0.031
0.901 = 0.03
CT-CP
0.666 + 0.008
0.5 + 0.2
0.641 0.006
0.819 0.019
Static
CT-GP
0.606 = 0.008
0.217 0.025
0.749
0.895 + 0.054
BCTT
0.576 0.015
0.358 + 0.082
+.011
0.829 0.028
NONEAT
0.517 0.002
0.129 0.002
0.408 0.005
0.877
THIS-ODE
0.528 + 0.005
0.132 + 0.002
0.544
0.878 1 0.026
POST
0.706 = 0.034
0.741 0.161
0.518
0.622 = 0.123
ADF-CP
0.641 0.009
0.652 + 0.012
0.542 0.012
0.518 + 0.003
Stream
BASS-Tucker
1.008 0.017
0.016
0.038
0.99 1 0.034
SFTL-CP
0.434 + 0.014
0.178 1 0.006
0.288 1 0.017
SFTL-Tucker
0.418 1 0.01
0.289 1 0.096
0.314 1 0.049
0.41 1 0.013
MAE
PTucker
0.392 + 0.009
0.323 0.053
0.307 + 0.005
0.197 0.029
Tucker-ALS
0.648 = 0.012
0.743 0.008
0.39 0.008
0.651 +
CP-ALS
0,666 0.013
0.746
0.676
CT-CP
0.462 0.005
0.348 0.141
0.489 1 0.006
0.632 0.015
Static
CT-GP
0.419 = 0.005
0.158 + 0.022
0.544
0.627 0.015
BCTT
0.392 + 0.004
0.267 + 0.067
0.299 + 0.006
0.607 0.027
NONFAT
0.355 + 0.001
0.078 + 0.001
0.265 1 0.003
0.622 + 0.006
THIS-ODE
0.363 = 0.004
0.083 0.002
0.348 0.006
0.603 0.009
POST
0.482 = 0.022
0.102
0.351
0.009
U.442
ADF-CP
0.445 + 0.006
0.009
0.381 0.006
0.393 0.009
Stream
BASS
0.822 = 0.024
0.749 0.009
0.919 + 0.041
0.73 0.018
SFTL-CP
0.246 1 0.005
0.121 0.003
0.176 0.006
0.305 0.006
SFTL-Tucker
0.24 1 0.002
0.18 0.042
0.196 0.03
0.263 0.011
Table 5: Final prediction enoor with R 3. The results were averaged from five runs.
21

---

RMSE
FitRecord
ServerRoom
BeijingAi-2
BeijingAir3
PTucker
0.603 = 0.045
0.677 + 0.129
0.464 + 0.012
0.421 + 0.074
Tucker-ALS
0.826 0.003
0.983 0.016
0.586 0.018
0.825 0.026
CP-ALS
0.878 = 0.012
0.904
0.897 0.215
0.863 0.024
CT-CP
0.663 + 0.008
0.384 0.008
0.64 0.007
0.818 0.019
Static
CT-GP
0.603 + 0.006
0.381 0.303
0.766 0.016
0.904 +: 0.046
BCTT
0.498 = 0.011
0.368 0.01
0.813 + 0.028
NONEAT
0.497 0.003
0.128 0.002
0.394 0.004
0.88 = 0.013
THIS-ODE ODE
0.138 = 0.003
0.006
0
POST
0.675 = 0.012
0.707 0.14
0.519 0.017
0.738 + 0.068
ADF-CP
0,652 0.01
0.646
0.548
0.552 0.026
Stream
0,604 0.043
0.493 0.071
0.391 0.005
0.634
SFTL-CP
0.424 0.006
0.166 0.013
0
0,481 0.006
SFTL-Tucker
0,448 = 0.009
0.406 + 0.052
0.249 0.017
0.432 0.019
MAE
PTucker
0,353 + 0.005
1.00
0.248 0.004
0.32 1:2:::00 0.038
Tucker-ALS
0.6 0-002
0.737 0.009
0.392 + 0.011
0.619 0.015
CP-ALS
0.64 0.009
0.005
0.593
0.637
CT-CP
0.459 0.005
0.008
0.488 0.005
0.016
Static
CT-GP
0,412 1 0.004
0.282
0.557 1 0.009
162810.0
BCTT
0.342 + 0.005
0.157 0.015
0.234 0.005
0.581 4.002
NONFAT
0.335 + 0.002
0.077 0.002
0.256 1: 0.003
0.627 0.005
THIS-ODE ODE
0.362 0.002
0.089 0.002
0.357 = 0.007
POST
0.461 + 0.008
0.518 0.087
0.357 - 0.011
0.558 0.058
ADF-CP
0.451 + 0.006
0.489 + 0.009
0.381 0.014
0,411 0.025
Stream
BASS
0.745 + 0.026
0.749 + 0.01
0.903 + 0.011
0.721 + 0.038
SFTL-CP
0.243 0.003
0.111 + 0.008
0.159 0.004
0.323 + 0.003
SFTL-Tucker
0.253 = 0.004
0.273 + 0.033
0.144 0.008
0.273 0.016
Table 6: Final prediction error with R = 7. The results were averaged from five runs.
1.2
1.00
RMSE
AMSE
RMSE
RMSE
1.00
1.00
0.7
0.70
0.60
0.45
20000
40000
4000
8000
0
8000
15000
6000
12000
Number of Entries
Number of Entries
Number of Entries
Number of Entries
(a) FitRecord
(b) ServerRoom
(c)
BeijingAir22
(d) BeijingAir-
Figure 7: Online prediction error with the number of processed entries (R 7)
incremental updates and never needs to revisit the past data. To demonstrate this point, we compared the training time of our method with BCTT on BeijingAir2 dataset. All the methods were run on Linux workstation. From Table 7 we can see large speed-up of our method with both the CP and Tucker form. The higher the rank (R), the more significant the speed-up.
R R= 2
R =
R = 5
R =
SFTL-CP SFTL-Tucker
27.1
27.2
28.5
29.1
32.3
35.6
43.2
59.3
BCTT
49.5
56.1
72.1
136.7
Table 7: Running time in seconds on BeijingAir2 datasct.
22


Reviews:
[{"review_id": "ucr4xe4djQ", "reviewer": "Reviewer_DKpV", "summary": "The work proposes a new approach for the decomposition of a tensor with a temporal dimension. The developed method can apply to streaming data where decomposing the entire tensor from scratch would be prohibitive. \n\nThe key idea is to model each entry (or parameter) in the decomposition factor as a function of time, thus representing a factor trajectory in time, and to assume that this function is sampled from a Gaussian process. Then, given the factor trajectories, initial tensor values are recovered by sampling from a specific Normal distribution.", "questions": "- The authors use CP and Tucker decomposition. Why Tensor Train or Tensor Ring decompositions, are not compared to or not used?\nFor example, [1] propose a simple method incorporating temporal component into the decomposition. It is interesting to see the number of parameters comparison, compression time, and quality metrics.\n\n[1] T4DT: Tensorizing Time for Learning Temporal 3D Visual Data\nMikhail Usvyatsov, Rafael Ballester-Rippoll, Lina Bashaeva, Konrad Schindler, Gonzalo Ferrer, Ivan Oseledets", "limitations": "The limitations arise from the set of assumptions made by the methods. Is this method applicable to the streaming 3D visual data, where the change in scenes might be significant from scene to scene? How the chosen rank R affects the quality, and what is the trade-off for it? These questions are not discussed in the paper though critical for the method to be applicable in practice.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "The paper is well-written and easy to follow. The authors clearly state the necessary background, problem, and notations.\nThe considered problem of handling streaming data is relevant. However, many technical choices were left unjustified.", "weaknesses": "The primary concerns are about (implicit) assumptions made in the paper:\n- Parameters of the decomposition factor are sampled from a GP with the same kernel. That says that the smoothness of the factor trajectories is assumed to be similar. That might be a restrictive assumption on the data the method is applied for, but it is not discussed.\n- Relates to the previous concern. Continuity in the latent space of the factors translates to the continuity in the space of the initial tensor values that follow from (1). In practice, a slight change in the factors might hugely affect the initial tensor values and vice versa. It would be essential to see the results of the method on more diverse problems, e.g., streaming point cloud data. \n- Finally, the trade-off between the decomposition rank and sustaining several Gaussian processes (i.e., better modeling of the factors). The rank R used in the experiments is relatively small.\n\n\nSome of the technical choices for the method are not well justified:\n- why is Gamma distribution used?\n- how are the hyperparameters of the models chosen? Do the data dictate those, or is there a reason for them?\n\nA simple baseline would tackle the problem as a regression task, thus possibly learned by a neural network.  It would be great to add it as a baseline.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your detailed and insightlful comments! Here are our responses. Please see global responses R2 & R3 for your questions regarding regression baselines and Tensor Train and Tensor Ring decompositions. C: comments; R: responses.\n\n>C1: \"Parameters of the decomposition factor are sampled from a GP with the same kernel. That says that the smoothness of the factor trajectories is assumed to be similar. That might be a restrictive assumption on the data the method is applied for, but it is not discussed.”\n\nR1: Great comment and we do agree. For the convenience of development and experiments, we used the same kernel for all the GP priors. This can indeed restrict the family of trajectory functions that can be learned. Actually, our model is flexible enough to incorporate different kernels, and it is straightforward to modify our algorithm to support multiple kernels (namely, constructing a different LTI-SDE for GPs with a different kernel and then running the same inference procedure).\n\nOn the other hand, since we used the Matern kernel, which has already covered a large family, and is way more flexible than other popular kernels, like Gaussian kernel. For example, the Matern kernel with $\\nu=3/2$ covers all the functions that are one-time differentiable, while the Gaussian kernel only models functions that are infinitely differentiable. From our experiments, a single Matern kernel has already been able to show the advantage of our method. But we do agree that introducing multiple kernels can potentially improve the performance further. \n\nThanks for the great suggestion and we will supplement this discussion into our paper.\n\n>C2: \"Continuity in the latent space of the factors translates to the continuity in the space of the initial tensor values that follow from (1). In practice, a slight change in the factors might hugely affect the initial tensor values and vice versa ...  Is this method applicable to the streaming 3D visual data, where the change in scenes might be significant from scene to scene?\n\nR2: We appreciate the reviewer raising this interesting point. We believe SFTL is able to handle the case reasonably well.\n\nSpecifically, whenever the observation space exhibits a significant pattern change, SFTL will update the factor trajectories via the proposed online filtering and RTS smoother. The online filtering (via moment matching and CEP) integrates the information of the previous states and the newly observed data to construct a new state distribution (i.e, the running posterior). That is, the new state distribution is *not* created standalone to absorb the abrupt data change; it is melt with the previous state information. The RTS smoother further back-propagates the new state information to all the old states, and thereby *smooths out* the estimate of the entire trajectory. One can see the online filtering and RTS smoother fuse the observation information (which might contain a drastic pattern change) in a soft and global way, and the new information is spread into the entire state chain. In this way, the learning of the factor trajectories is robust and smooth.\n\nWe thank the reviewer for the great question and will investigate in-depth this problem in our future work.\n\n>C3: \"The rank R used in the experiments is relatively small. How the chosen rank R affects the quality, and what is the trade-off for it?\"\n\nR3: We appreciate the reviewer's comment. The relatively small rank R used in the experiments is due to the intrinsic low-rank property of our datasets (see details in Section E of Appendix) - either being sparse (e.g. FitRecord) or with a small number of objects in each mode (ServerRoom, BeijingAir). Hence, we believe the lower rank setting is sufficient to capture their temporal patterns.\n\nWe do appreciate with the reviewer's concern about \"the trade-off between the decomposition rank and sustaining several Gaussian processes\". In fact, we have examined SFTL with different ranks and compared the predictive performance and efficiency ($R \\in \\{2, 3, 5, 7\\}$) to show the trade-off; see Table 1 in the main paper and Table 3-6 in Appendix for the detailed results and discussions.  We will test SFTL with denser data, such as streaming point cloud data, and higher ranks to investigate this trade-off more comprehensively. We thank the reviewer for the constructive suggestion.\n\n>C4: \"why is Gamma distribution used?\"\n\nR4: In Bayesian modeling, the Gamma distribution is the mostly widely used prior distribution to model the noise (inverse variance) for Gaussian models/likelihoods. We followed this standard choice to build a full Bayesian model for factor trajectory learning. Thanks for the question. We will clarify this point in our paper.\n\n>C5: \"how are the hyperparameters of the models chosen? Do the data dictate those, or is there a reason for them?\"\n\nR5: We used cross-validation to tune and select the hyper-parameters and hence \"the data dictate those\" (see the details in Line 285-289). In addition, we investigated how sensitive the predictive performance is to the hyper-parameter choice (the amplitude $a$ and the length-scale $\\rho$). We refer the reviewer to Section D.2 (Sensitive Analysis) and Table 2 of Appendix for the detailed results.\n"}, {"review_id": "4zXbUIA8Em", "reviewer": "Reviewer_vq55", "summary": "This paper studies the dynamic tensor factorization under the streaming data setting. In particular, the authors propose to use a state-space Gaussian process (GP) to parameterize the CP or Tucker factors as functions of time. To efficiently learn the GP and handle streaming data, an equivalent SDE and a conditional expectation propagation algorithm are then derived. For experiments, the authors show that the model can learn factor and tensor entry trajectories. Besides, it achieves good competition performance, even compared with model non-streaming settings.", "questions": "1. In Equation 9, for the first Gaussian distribution, is the variance missing?\n", "limitations": "None", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1. This work presents the first model that attempts to learn factor trajectories under streaming data setting.\n\n2. The performance of completion is good. Even though the proposed model runs under the online setting, the error is smaller than many models that use the whole dataset.\n", "weaknesses": "1. Although the problem setting is new in the tensor decomposition field, most techniques used in this work seem standard, such as state-space Gaussian process and expectation propagation.\n\n2. The notations are too complex and hard to follow. There are too many super/sub-scripts and sometimes the meanings are unclear. For examples, for latent variables $z$, sometimes, the second subscript is the rank (e.g., in Line 128, $z^m_{j, r}(t)$, where $r$ is the rank). Then, in Line 141, it becomes the time stamp, e.g., $\\bar{z}^m_{j, k}$ with $k$ being the time index.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your valuable comments. Here are our responses. C: comments; R: response\n\n> C1: \"Although the problem setting is new in the tensor decomposition field, most techniques used in this work seem standard, such as state-space Gaussian process and expectation propagation.\"\n\nR1: We do agree that we did not invent new techniques from scratch. However, many works, including the many references in Related Work, also use existent techniques, via their novel application and combination, to address the problem of interest. we believe such works are still valuable contributions to the relevant fields. In addition to the new problem setting, we believe the novelty of our work is mainly threefold: (1) We are the first to apply the LTI-SDE representation of GPs to learn the factor trajectories for tensor decomposition. Not only is this modeling a nice fit to our problem, it also evades the computational bottleneck and various sparse approximations, which are the known headaches for GP learning. (2) We combined the state-space model and CEP framework to develop an efficient filtering algorithm, which overcomes the infeasibility of the standard Kalman filtering and allows fast fusion of the incremental data information. (3) Our novel approximation structure (Eq6) further enables the application of the classical RTS smoother, which can update the trajectory posterior highly efficiently and in parallel (this is critical to fast streaming inference). We have demonstrated the effectiveness of our method in a series of benchmark datasets.\n\n> C2: \"The notations are too complex and hard to follow.\"\n\nR2: Thanks for the comments and we agree. To introduce the state and trajectory values of each factor, we need to index the tensor mode, the object, the rank and the time, which indeed make the notation cluttered. We will introduce more concise notations, and avoid using similar symbols, like $z$ and $\\bar{z}$, to denote different variables. We will further polish our paper to make it more friendly to readers.\n\n> C3: \"In Equation 9, for the first Gaussian distribution, is the variance missing?\"\n\nR3: Yes, and thanks much for catching this typo! We will fix it in the paper."}, {"review_id": "SurxxIvY1v", "reviewer": "Reviewer_ZhxY", "summary": "The authors highlight that real-world tensor data often includes time information, representing the timestamps of interactions between objects in different modes. However, most existing decomposition methods for tensors with temporal information only estimate static factors for each object and cannot capture the temporal variation of the factors. To address this limitation, the authors introduce Streaming Factor Trajectory Learning (SFTL), which can efficiently handle data streams and estimate the temporal evolution of objects' representation without storing or re-accessing previous data. \n\nThey use a Gaussian process (GP) prior to model the factor trajectory of each object as a function of time. GPs are flexible nonparametric function priors that can capture complex temporal dynamics. The factor trajectories are combined using either the CANDECOMP/PARAFAC (CP) or Tucker decomposition to sample the tensor entry values at any given time.\n\nThe authors evaluate SFTL using both simulation studies and real-world temporal tensor datasets. In the simulation studies, SFTL recovers nonlinear factor trajectories and provides reasonable uncertainty estimation. When applied to real-world datasets for missing value prediction, SFTL consistently outperforms state-of-the-art streaming CP and Tucker decomposition algorithms in both online and final predictive performance. In some cases, SFTL even surpasses the prediction accuracy of recent static decomposition methods that require multiple passes through the dataset. The learned factor trajectories from a real-world dataset also exhibit interesting temporal evolution.", "questions": "Can you please explain the real-world applications (we only see the dataset names, but it is not sufficient to understand the application)?\nCan you also add how much gain you provide with STL compared to baselines, both in efficiency and accuracy, e.g. in percent?\nCan you please add more explanations to Figure 1? It is not so easy to understand how the overall method works from that figure.\n", "limitations": "The authors can discuss more the scalability and the applicability of STL. What are the applications that STL can be used in real-world?", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The main advantage of this method is its computational efficiency. To overcome the computational challenges of the GP's covariance matrix computation in streaming inference, they employ spectral analysis to convert the GP into a linear time-invariant (LTI) stochastic differential equation (SDE). This conversion enables them to transform the SDE into an equivalent state-space prior over the factor states at observed timestamps, simplifying posterior inference and making it computationally efficient.\n\nIn addition, they increase the scalability of tensor decomposition as in the following. They leverage the chain structure of the state-space prior and utilize the conditional expectation propagation framework to develop an efficient online filtering algorithm. This algorithm allows for the estimation of a decoupled running posterior of the involved factor states when new entries at a new timestamp arrive. The decoupled Gaussian estimate enables independent and parallel computation of the full posterior of each factor trajectory without revisiting previous data. The method scales linearly with the number of observed timestamps.\n", "weaknesses": "I think the paper is not so clear and it is hard to follow the contribution especially in Section 3 and 4. The simulation study looks fine, but not a big difference is observed between SOTA and the STL on the real-world applications.  Besides the real-world applications are not explained, so it is also not clear to understand what are the application areas of this method.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your comments and questions. Here are our clarification and responses. Please see global response R1 for our clarification of Fig. 1.  C: comments; R: responses\n\n> C1: \"the real-world applications are not explained, so it is also not clear to understand what are the application areas of this method. What are the applications that STL can be used in real-world? Can you please explain the real-world applications?\"\n\nR1: Thanks for the question. In fact, we did explain the application and dataset details  in **Section E of Appendix**. Due to the space limit, we put this information in Appendix, but we point out the place in Line 284 of the main paper. *FitRecord* is about the outdoor exercises, *ServerRoom* is about the temperature management for a supercomputing center, and *BeijingAir-2* and *BeijingAir-3* are about the air pollutant monitoring. Please see the details in Section E of Appendix. \n\n> C2: The simulation study looks fine, but not a big difference is observed between SOTA and the STL on the real-world applications...Can you also add how much gain you provide with STL compared to baselines, both in efficiency and accuracy, e.g. in percent?\n\nR2: In fact, we do view the difference between SOTA and SFTL is **big**.  First, as a streaming factorization approach, SFTL consistently outperforms its counterparts --- **the state-of-the-art streaming tensor decomposition algorithms** --- by a large margin in all the ranks. Please see Table 1 of the main paper, and Table 3-5 of Appendix. To highlight the difference, we follow your suggestion to list the performance gain of our method over the best competing approaches in RMSE as follows. As we can see, for every rank, SFTL achieves at least 16% improvement. In most cases, the improvement is over 30%.\n\n|  FitRecord   | R=2 |   R=3 |     R=5 |    R=7 |    \n| :---        |            ---: |---:  | ---:  | ---:  | \n| Best (streaming)    | 0.669   | 0.641   |  0.648  | 0.604   |\n| SFTL  | 0.437      | 0.418  | 0.424   | 0.424    |\n| Performance Gain  |   36%    | 35%   | 37%   | 30%   |\n\n|  ServerRoom   | R=2 |   R=3 |     R=5 |    R=7 |    \n| :---        |            ---: |---:  | ---:  | ---:  | \n| Best (streaming)    | 0.764   |  0.652   |  0.64  | 0.593  |\n| SFTL  | 0.18      | 0.178   | 0.161   | 0.166   |\n| Performance Gain   |   76%    | 72%   | 75%   | 72%   |\n\n|  BeijingAir2   | R=2 |   R=3 |     R=5 |    R=7 |    \n| :---        |            ---: |---:  | ---:  | ---:  | \n| Best (streaming)    | 0.539   |  0.518   |   0.516  |  0.491  |\n| SFTL  | 0.323     | 0.288   | 0.248   | 0.249   |\n| Performance Gain   |   41%    | 43%   | 52%   | 48%   |\n\n|  BeijingAir3   | R=2 |   R=3 |     R=5 |    R=7 |    \n| :---        |            ---: |---:  | ---:  | ---:  | \n| Best (streaming)    |  0.54   |  0.622   |   0.551  |  0.552  |\n| SFTL  | 0.417     | 0.41  | 0.439   |   0.432  |\n| Performance Gain   |   24%    | 34%   | 22%   | 16%   |\n\nSecond, even compared with non-streaming (i.e., static) decomposition methods, in most cases, SFTL still exhibits considerable improvement, except that on ServerRoom dataset, SFTL is slightly worse than NONFAT and THIS-ODE yet is still much better than all the other static decomposition methods. Note that `the static methods can repeatedly access the dataset for many training epochs, while the streaming method can only pass through the dataset for once`. Hence, the static methods have `much more privilege on accessing training information`. It is therefore not surprising to expect that the static methods should perform better in prediction accuracy than streaming methods. However, **SFTL can still outperforms the static methods in most cases and often by a large margin, which further shows the advantage of SFTL as a streaming approach**.  To highlight this point, we list the performance gain over the static methods in RMSE as follows.\n\n|  RMSE (R=2)  | FitRecord |   BeijingAir-2 |     BeijingAir-3 |\n| :---        |            ---: |---:  | ---:  | \n| Best (static methods)    | 0.543   | 0.425   |   0.442  |\n| SFTL  | 0.437      | 0.323   | 0.417   |\n| Performance Gain   |   20%    | 24%   | 8%   |\n\n|  RMSE (R=3)  | FitRecord |   BeijingAir-2 |     BeijingAir-3 |\n| :---        |            ---: |---:  | ---:  | \n| Best (static methods)    | 0.517  | 0.408   |  0.421  |\n| SFTL  | 0.418      | 0.288   | 0.41   |\n| Performance Gain   |   19%    | 30%   | 3%   |\n\n|  RMSE (R=5)  | FitRecord |   BeijingAir-2 |     BeijingAir-3 |\n| :---        |            ---: |---:  | ---:  | \n| Best (static methods)    | 0.503   | 0.395   |  0.535  |\n| SFTL  | 0.424      | 0.248   | 0.439   |\n| Performance Gain   |   16%    | 38%   | 19%   |\n\n|  RMSE (R=7)  | FitRecord |   BeijingAir-2 |     BeijingAir-3 |\n| :---        |            ---: |---:  | ---:  | \n| Best (static methods)    | 0.497   | 0.394   |  0.451  |\n| SFTL  | 0.424      | 0.249   | 0.432   |\n| Performance Gain   |   16%    | 37%   | 6%   |\n\nTo examine the efficiency gain, we extend our running time results in Table 7 of Appendix by adding the speed-up percentage. We compared with the running time of BCTT on BeijingAir2 dataset; please see Section G of Appendix for details. \n\n|  Training Time(Sec)   | R=2 |   R=3 |     R=5 |    R=7 |    \n| :---        |            ---: |---:  | ---:  | ---:  | \n| BCTT (benchmark)    | 49.5   | 56.5   |  72.1  | 136.7  |\n| SFTL-Tucker  | 32.3      | 35.6   | 43.2   | 59.3   |\n| SFTL-CP  | 27.1      | 27.2   | 28.5   | 29.1 |\n| Speed-up SFTL-Tucker   |   55%    | 59%   | 67%   | 131%   |\n| Speed-up SFTL-CP   |   83%    | 107%   | 152%   | 371%   |\n\nAs we can see, the speed-up is significant as well. The larger the rank, and the bigger the speed-up. We will add these tables into our paper."}, {"review_id": "1Xb3NrkqLK", "reviewer": "Reviewer_Dugj", "summary": "In traditional tensor decomposition, the factor of decomposition is static. Therefore, it was not possible to fit time-evolving data well. This paper develops SFTL using a Gaussian process to capture the time-developing trajectory of the factors. Their naive implementation would have increased computational complexity. Therefore, they reduce the computational complexity by replacing the Matern kernel in GP with an SDE.", "questions": "- I struggled to understand Fig 1. Since there are 3 trajectories in Fig 1, does this correspond to the SFTL for a tensor that has three 3 modes?\n- In my understanding, the size of vec(z) depends on the rank R, where R is a hyperparameter. In Algorithm 1, R is not included in Input. That point needs to be clarified.\n- In section 6.1, the word `node` suddenly appears, but I am unclear what this means. My understanding is that the input to the experiment in section 6.1 is a full-rank 2x2 time-evolving matrix with noise, is this correct?", "limitations": "The authors appropriately mention in the Appendix that the speedup of the proposed method is only applicable to Matern kernels.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "- Standard tensor decomposition can only extract linear patterns. SFTL, however, uses a Gaussian process to extract nonlinear features in the time direction.\n- Elegant techniques are used to reduce the computational complexity. It is unfortunate that these techniques are not applicable to non-matern kernels, i.e., RBF kernels.\n- They provide a theoretical discussion of the computational complexity of the proposed method and include additional experimental results in the appendix.\n- Notation of the symbols is clear and not misleading. The paper is well written.\n- It is good to have a discussion about the sensitivity of hyperparameters.", "weaknesses": "There are some known tricks to reduce the computational complexity of the Gaussian process from O(n^3), such as Inducing Variable Method, Deterministic Training Conditional (DTC) Approximation, and Fully Independent Training Conditional (FITC) Approximation. All of these tricks can be adapted to non-Matern kernels. With these tricks, I think it is possible to achieve fast computations without Hartikainen's elegant tricks. Although the elegant method of this paper is to be evaluated, it is a drawback that there is no comment on this point.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your insightful and constructive comments. Here are our responses. C: comments; R: responses\n\n>C1: \"There are some known tricks to reduce the computational complexity of the Gaussian process from O(n^3), ... I think it is possible to achieve fast computations without Hartikainen's elegant tricks. Although the elegant method of this paper is to be evaluated, it is a drawback that there is no comment on this point.\n\nR1: Great comments!  We do agree these methods, like DTC and FITC, can reduce the GP complexity as well. However, these methods have to use additional approximations, namely, a low-rank or sparse approximation of the kernel matrix to reduce the cost, which might bring additional errors and learning challenges.  Our motivation of using LTI-SDE is twofold: (1) we can completely avoid sparse/low-rank approximations and perform exact GP inference with linear time complexity. (2) With the state-space prior, we can fulfill efficient streaming inference with the online filtering and smoothing framework. On the contrary, conducting streaming inference for sparse latent GP models is much more challenging (e.g., how to dynamically update the inducing points and locations). Thanks for the constructive suggestion. We will add the discussion into our paper.\n\n>C2: \"I struggled to understand Fig 1. Since there are 3 trajectories in Fig 1, does this correspond to the SFTL for a tensor that has three 3 modes?\"\n\nR2: Yes, Fig. 1 indeed uses a tensor with 3 modes as an example to illustrate the conditional independence in our model. Please see global response R1 for our clarification of Fig. 1. We will add more detailed explanation of Fig.1 in our paper.\n\n>C3: \"In Algorithm 1, R is not included in Input. That point needs to be clarified.\"\n\nR3: Thanks for the great comment and suggestion! We do agree. We will add $R$ into the input of our algorithm table.\n\n>C4: \"In section 6.1, the word node suddenly appears, but I am unclear what this means. My understanding is that the input to the experiment in section 6.1 is a full-rank 2x2 time-evolving matrix with noise, is this correct?\"\n\nR4: We appreciate the reviewer catching this problem. Yes, your understanding is correct. \"node\" means the object in the mode, and we should use the word \"object\" to avoid any confusion. We will modify them to improve the clarify."}]

Author rebuttal:
We thank all for the reviewers for their careful, valuable and constructive comments! We post our global response to a few questions here. C:comments; R: responses 

>C1: more explanation/clarification for Figure 1. 

R1: Figure 1 explains the key step in our development of the online filtering algorithm. Whenever a new batch of observations $\mathcal{B}\_{n+1} $ is received, we want to compute the running posterior of the related trajectory states $\mathcal{\Theta}\_{n+1}$ (those are the factor states for objects showing up in $\mathcal{B}\_{n+1}$ ), which we denote by $p(\mathcal{\Theta}\_{n+1} | \mathcal{D}\_{t\_n}\cup \mathcal{B}\_{n+1})$. Here $\mathcal{D}\_{t_n}$ is all the data received before $\mathcal{B}\_{n+1}$, i.e., the past data. Our key observation is that --- as shown in Fig. 1 --- **the new states $\mathcal{\Theta}\_{n+1}$ and new data $\mathcal{B}\_{n+1}$ are  independent to the past data $\mathcal{D}\_{t_n}$, conditioned on the most-recent states (one-step earlier) $\mathcal{\Theta}\_n$ and noise inverse variance $\tau$** (in the Gaussian likelihood; see Eq1). In other words, (roughly speaking) if we know the most recent states $\mathcal{\Theta}\_n$ and $\tau$, we do not need to access the past data $\mathcal{D}\_{t_n}$ to learn the new states $\mathcal{\Theta}\_{n+1}$. This is due to the chain structure between the successive states in the state-space prior; see Eq4. Fig. 1 illustrates this observation from the graphical model representation. We can see that **conditioned on $\mathcal{\Theta}\_n$ and $\tau$, all the paths from $\mathcal{D}\_{t\_n}$ to $\mathcal{\Theta}\_{n+1}$ and $\mathcal{B}\_{n+1}$ are blocked.** According to the graphical model properties, it shows the conditional independence; please see D-separation in graphical model literature [1]. Leveraging the conditional independence, we can derive the rule of computing $p(\mathcal{\Theta}\_{n+1} | \mathcal{D}\_{t\_n}\cup \mathcal{B}\_{n+1})$ --- it only needs the running posterior of the most recent states $\Theta\_{n}$ and $\tau$; we do not need to revisit $\mathcal{D}\_{t\_n}$. Please see Eq7 and below and Eq8 for the details. We can thereby keep tack the running posterior for each factor state, and integrate with the likelihood of the new data to compute the new running posterior (via CEP) --- that is how our online filtering algorithm works.

[1] Bishop, Christopher M., and Nasser M. Nasrabadi. *Pattern recognition and machine learning*. Chapter 8.3, Springer, 2006.

>C2: A simple baseline would tackle the problem as a regression task. It would be great to add it as a baseline.

R2: Great suggestion. We supplement the experimental results in the below. We used the entry indices and timestamps as features to build three regression models: Bayesian Linear Regression(BLR), RBF-kernel support vector machine (RBF-SVM), and a neural-network(NN) with 3 MLP layers. With the same training and test datasets, the predictive performance of these models is given in the below (we also list SFTL with R=5 for comparison).

| RMSE      | FitRecord | ServerRoom     |  BeijingAir-2 |     BeijingAir-3 |    
| :---        |    :----:   |          ---: |---:  | ---:  | 
| BLR      | 0.967 ± 0.003      | 1.164 ± 0.013   | 1.118 ± 0.016   | 1.074 ± 0.028  |
| RBF-SVM   |0.881 ± 0.012        | 1.232 ± 0.002      | 1.052 ± 0.018   | 1.018 ± 0.016   |
| NN   | 0.725 ± 0.024        | 1.431 ± 0.035      | 0.856 ± 0.054   | 0.871 ± 0.022   |
|SFTL-CP   | **0.424 ± 0.014**        | **0.161 ± 0.014**      | **0.248 ± 0.012**   | 0.473 ± 0.013   |
| SFTL-Tucker   | 0.430 ± 0.010       | 0.331 ± 0.056   |  0.303 ± 0.041   | **0.439 ± 0.019**  |

We can see SFTL outperforms all the regression models by a large margin, which confirms the advantage of SFTL and the necessity of learning factor representations. Note that SFTL also exhibits large improvement on other rank settings (see Table 3-5 of Appendix). We will supplement these results and a discussion into our paper.


>C3: Tensor Train or Tensor Ring decomposition not used

R3: Our current work only focuses on the classical yet perhaps the most widely used tensor decomposition forms, CP and Tucker, to prove our modeling and algorithmic framework are effective and promising. We definitely want to extend our method to more modern decomposition forms, like Tensor Train and Tensor Ring. However, it is non-trivial to adjust our current CEP inference to Tensor Train and Tensor Ring. The reason is that the chain of matrix multiplication in their decompositions makes the standard conditional moment matching intractable to compute. In the next step, we plan to develop novel computational tricks to address this issue so as to enable streaming factor trajectory learning for Tensor train and Tensor Ring decomposition as well.

Meta-review: The paper introduces a family of dynamic tensor decomposition (CP or Tucker) models wherein the factors are time-indexed and and their trajectories are assumed to be drawn from a Gaussian (GP) process. Inference is challenging and the paper applies a clever trick to convert the GP prior into a state-space prior. The reviewers were impressed with the model and the convinced of the novelty and applicability of the techniques used to derive inference. The reviewers also engaged substantially with the authors during the discussion period and were all convinced by the authors' responses about their initial concerns and questions.
Decision: Accept (poster)