Dynamic Tensor Decomposition via Neural Diffusion-Reaction Processes
Authors: Zheng Wang
Shikai Fang
Shibo Li
Shandian Zhe
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Tensor Decomposition
Representation Learning

Abstract: Tensor decomposition is an important tool for multiway data analysis. In practice, the data is often sparse yet associated with rich temporal information. Existing methods, however, often under-use the time information and ignore the structural knowledge within the sparsely observed tensor entries. To overcome these limitations and to better capture the underlying temporal structure, we propose Dynamic EMbedIngs fOr dynamic Tensor dEcomposition (DEMOTE). We develop a neural diffusion-reaction process to estimate dynamic embeddings for the entities in each tensor mode. Specifically, based on the observed tensor entries, we build a multi-partite graph to encode the correlation between the entities. We construct a graph diffusion process to co-evolve the embedding trajectories of the correlated entities and use a neural network to construct a reaction process for each individual entity. In this way, our model can capture both the commonalities and personalities during the evolution of the embeddings for different entities. We then use a neural network to model the entry value as a nonlinear function of the embedding trajectories. For model estimation, we combine ODE solvers to develop a stochastic mini-batch learning algorithm. We propose a stratified sampling method to balance the cost of processing each mini-batch so as to improve the overall efficiency. We show the advantage of our approach in both simulation studies and real-world applications. The code is available at https://github.com/wzhut/Dynamic-Tensor-Decomposition-via-Neural-Diffusion-Reaction-Processes.

Dynamic Tensor Decomposition via Neural Diffusion-Reaction Processes
Zheng Wang"
Shikai Fang"
Shibo Li
Kahlert School of Computing University of Utah Salt Lake City, UT 84112 u1 u120884 (ett88477tthcc..m
Kahlert School of Computing University of Utah Salt Lake City, UT 84112 shikai fang@utah.edu
Kahlert School of Computing University of Utah Salt Lake City, UT 84112 shibol cs utah.
Shandian Zhe
Kahlert School of Computing
University of Utah Salt Lake City, UT 84112 zhe@cs, utah. edu
Abstract
Tensor decomposition is an important tool for multiway data analysis. In practice, the data is often sparse yet associated with rich temporal information. Existing methods, however, often under-use the time information and ignore the struc- tural knowledge within the sparsely observed tensor entries. To overcome these limitations and to better capture the underlying temporal structure, we propose Dynamic EMbedIngs fOr dynamic Tensor dEcomposition (DEMOTE). We de- velop neural diffusion-reaction process to estimate dynamic embeddings for the entities in each tensor mode. Specifically, based on the observed tensor entries, we build a multi-partite graph to encode the correlation between the entities. We construct graph diffusion process co-evolve the embedding trajectories of the correlated entities and use neural network to construct reaction process for each individual entity. In this way, model can capture both the com- monalities and personalities during the evolution the embeddings for different neural network to model the entry value nonlin-
entities. We then ear function of the embedding trajectories. For model estimation, we combine ODE solvers to develop stochastic mini-batch learning algorithm. We propose a stratified sampling method to balance the cost of processing each mini-batch improve the overall efficiency. We show the advantage of our approach in both simulation study and real-world applications. The code available at attps://githus. www.://..ooooaaaaaaaaaaeeeeeeeee com wwt.m//whuttttt http://yyyynnneeeeeeeerreeeeeeee
via Www.al-Diffusseeee Beceee-eacctoo Processees
Introduction
Multiway data is common in real-world applications and naturally represented by tensors. For example, online shopping and promotion activities can be expressed as a three-mode tensor (customer, commodity, online merchant). Tensor decomposition is important tool for multiway data analysis. estimates embeddings for the entities each tensor mode, with which to recover the observed
"Equal contribution
t "Correspondng author.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

entry values. The embeddings can reflect the underlying structures within the entities and can be used as predictive features, such as for recommendation and ads auction.
In practice tensor data is often very sparse. That is, the observed entries only take a tiny portion of all possible entries, say, 0.01 %. In addition, the data often includes timestamps for the observed entry values, which imply rich, complex temporal variation patterns. Current tensor decomposition approaches often ignore the structure knowledge within the sparsely observed entries and under-use the temporal information, e.g., simply binning the timestamps into crude time steps (Xiong et al. 2010; Rogers al. 2013 Zhe et al., 201 2015 Du et al. 2018). More important, standard tensor decomposition estimates a static embedding for each entity. However, as the representation of entities, these embeddings summarize the underlying properties of the entities, which can naturally evolve with time, such as customer interests, user income, product popularity, and fashion. Learning static embeddings can miss capturing these interesting, important temporal knowledge. While the most recent work (Wang et al. 2022) has proposed the first decomposition method to estimate embedding trajectories, it never considers the structural knowledge within the data.
To overcome these limitations, we propose DEMOTE, a dynamic embedding approach for dynamic tensor decomposition. We construct nonlinear diffusion-reaction process in an Ordinary Differ- ential Equation (ODE) framework to estimate embedding trajectories for tensor entities. The ODE framework is known to be flexible and convenient to handle irregularly sampled timestamps and sparsely observed data (Rubanova et al. 2019). In addition, since ODE models focus on learning the dynamics (i.e., time derivatives) of the target function, they have promising potential for providing robust, accurate long-term predictions (via integration with the dynamics). Specifically, to leverage the structural knowledge within the data, we first build a multi-partite graph based on the observed entries. The graph encodes the correlations between entities at different modes in terms of their interactions. We then construct a graph diffusion process in the ODE to co-evolve the embedding trajectories of correlated entities. Next, we use neural network to construct reaction process to model the undividual-speciffe evolution for each entity. this way, our neural diffusion-reaction process captures both the commonalities and personalities of the entities in learning their dynamic embeddings. Given the embedding trajectories, we model the entry value as a latent function of the associated entities' trajectories. We use another neural network to flexibly estimate the function and to capture the complex relationships of the entities. For efficient training, we base on ODE solvers to develop a stochastic mini-batch learning algorithm. We develop stratified sampling scheme, which can balance the cost of executing the ODE solvers in each mini-batch so as to improve the efficiency.
We evaluated our method in both simulation and real-world applications. The simulation experiments show that DEMOTE can successfully capture the underlying dynamics of the entities from their temporal interactions and recover the hidden clustering structures within the trajectories. Then in three real-world applications, we tested the accuracy in predicting the tensor entry values at different time points. DEMOTE consistently outperforms the state-of tftte-of-th decomposition methods that incorporate temporal information, often by a large margin. We also demonstrated that both the diffusion and reaction processes contribute to the learning performance. Finally, we investigated the learned embedding trajectories and found interesting evolution paths and hidden structures.
2 Notations and Background
Suppose we have collected data for a K-mode tensor. Each mode ls includes de entities, which we index by 1,..., ....,dd. dd. We then index each tensor entry by tuple E = (4,,-.,,4k) where for each k. have 1 S 4e de. Suppose we observed N tensor entry values and timestamps. The dataset is denoted by D = {(Lg, (141,1,91),.... v,),...... (LN, ]N where (t,,) and { are the timestamps and entry values, respectively. Our goal is for each entity j of mode to estimate a dynamic embedding u& (t): 14, + " That is, the embedding is a time function (trajectory) with an R-dimensional output. Standard tensor decomposition introduces a static embedding representation for each entity, namely, uz is considered as time invariant. Tensor decomposition aims to estimate the embeddings (or factors) to reconstruct the tensor. For example, the classical Tucker decomposition (Tucker 1966) employs multilinear factorization model, M = W x1 U' X2... ... X K U*, where M E 1 X...xx; the entire tensor, W E IR Rg XX. is the tensor-core parameter, U comprises all the embeddings of the entities in mode k. and X is the rensor-matrix product at mode & (Kolda 2006). The popular CANDECOMP/PARAFIS
(CP) decomposition (Harshman) 1970) can be viewed as simplified version of Tucker decomposition, where we set R1
Rx R. and the tensor-core W to be
2

---

Embeddings of entities of each mode
Mode
Edges defined by observed entries (i.
u
(1)
Diffusion process along edges
Reaction process on individual entities
Time trajectory
Table of observed data
Mode
Mode
Mode
Time-
Value
k
(t)
(1
uf (t)
n²(t)
(t)
Mode
Mode 3
Figure 1: The illustration of the embedding model in DEMOTE.
diagonal. Hence, each entry value is factorized as me = =uu,, o....cufx ) where o the Hadamard (element-wise) product, and À corresponds to diag(W). While CP and Tucker decomposition are popular, their multilinear modeling nonlinear relationships of the entities, Xu et al. (2012); Zhe et al. (2015 2016a) used a Gaussian process (GP) (Rasmussen and Williams) 2006) model the entry value as random function of the embeddings, mi = _( "xx where ~ GP (0, <((e, xx:)) )), x = [ u and [u i ..... u are the embeddings of the entities in entry l and l', respectively, and k(.--) the covariance (kernel) function. Given the GP prior, any finite multi-variate Gaussian distribution, m ~ N(0, K), where m = [m....... Mey K is the N x N kernel matrix, and each [K] n(xe, xe Suppose we have collected continuous observations for the N entries
be oversimplistic for complex applications. To estimate
N entry values follow
[11 :....... Mv]. We can use Gaussian noise model: Yn = me, + Gn where Gn (O, d The marginal likelihood p(y) (y (0, I)), which we can maximize to estimate the embeddings and the model parameters.
Practical data often includes temporal information, i.e., the timestamp when each observed entry value is generated. To leverage this information, existing methods often bin the timestamps into series of steps, say, by weeks months (Xiong al. [000 Rogers et al. 2013] Zhe et al.) 60166] Song et al. 2017). The tensor is then expanded with an additional time-step mode, and one can apply any decomposition algorithm to estimate embeddings for both the entities and time steps. To capture the temporal dependency, a conditional model is often used (Xiong et al. 2010), say, p(tj+i|t,) (t [t, 71) where t, the embedding of j-th step. To leverage the continuous time information, Zhhang et al. (2021) recently developed continuous CP decomposition, where the coefficients À are modeled as time function with polynomial splines.
Model
Standard tensor decomposition assumes the embeddings static and time-invariant. However, the embeddings summarize and extract the properties of entities, which can evolve with time, such as customer interests, health status, and product popularity. Therefore, only estimating static embeddings can miss important temporal variations of the entities' properties, resulting in poor representations and predictive performance. In addition, practical tensor data are typical sparse, and only small portion of entries actually have data. Within these entries can be valuable structural knowledge. Current methods, however, are rarely aware of such knowledge. To overcome these limitations, propose DEMOTE, a novel dynamic embedding approach.
Specifically, propose ODE model to Iearn the embedding trajectories (u)(t) l s K, 1 d. The ODE framework is known to be amenable for irregularly sampled, sparsely observed data. More important, ODE models concentrate on learning the time derivative duž: /dt (i.e., dynamics), rather than trajectory function itself. Therefore, they have a promising potential to give reliable, long-term trajectory prediction (via numerical integration) even at time points far away from the training timestamps, provided the time derivative is well captured. We construct joint ODE model for all the embedding trajectories. The ODE consists of a diffusion process and

---

reaction process. The diffusion process leverages the structural knowledge in data to co-evolve the embeddings of correlated entities, so as to better overcome the data sparsity. The reaction process models the entity-specific evolution so that it can capture the individual differences in the embedding evolution. The ODE model synergizes the two processes to capture both the commonalities and personalities of these embedding trajectories.
Diffusion Process on Multi-Partite Graphs. First, we construct graph-based diffusion process to exploit the entity correlations reflected in data. Intuitively, if an observed entry involves entity (e.g., customer A) and (e.g.. commodity B), the two entities are likely correlated. Thus, we can draw an edge between A and B to express the correlation. We then generalize this intuition to create K --partite undirected graph GEE, V). to encode such correlations across all the entities the K tensor modes. Each vertex represents particular entity, and the entire collection of the entities partitioned into groups, V = V1 U... VK, where group V {uf...., o/, represents the entities of mode k. Two entities (at different modes) are connected if they were observed to interact, namely, (ufi " vf E EE 3em E D such that Ln (...... ........................ j',,... where and j' are indices at mode k and k', respectively. See Fig. for illustration. This graph naturally implies underlying information diffusion across the entities within their interactions. For example, if customer A connects to products B and C, it might mean that A distributes their interests/ interests/willlngese willingness/budgets to purchase and C. The edges between one merchant A and a list of products (B, C.... } might indicate the diffusion of willingness to increase the inventory of these products.
&.k' To flexibly estimate the diffusion rate, we introduce a weight w "Jy' for each edge (uf uf)) E. We then arrange these weights into K(K 1) adjacent matrices, W ( [1 S k. K' K #
where Wh. = (W Each Wkkk' is a sparse de x du matrix that represents the edges and
edge weights between V and VK', i.e., [Wk. = j.j if (ut uy)) E and 0 otherwise. We now construct diffusion process based on the K --partite graph. We view the embedding trajectory as a kind of concentration. For each entity j at mode k; the change rate of its concentration (embedding) uf (t) determined by the difference from the concentrations of its neighbors. Since the neighbors come from entities of all the other K - 1 modes, we have
duf dt BE(1..... K
Σ
Whe lass (uz ut())
E (w)"'
aj o..
=E(L..... K|VE
where whe his is the j-th row of Wh.s, U+(t) (t).........)..)..)....)) ..... uZ [t)] is the embeddings of all the entities mode s, of size d, R, and 0 k. = Σド ( D3.3" is the degree of vertex in Wke,, We can see that the evolution of the embeddings for different modes are coupled. Hence, it is natural to formulate the diffusion process jointly for all the embeddings, dl((t) (U' (t).))..... U (t)) /dt
(t(()) (L) (W- A))(() where
0
W1.
WLK
W2,1
0
W
WK-1K
WK11
WK.K-
A= = diag (Ee( A the degree matrix of Wkis.
k. is E.E( A K.s) and each A h.* diag(a kis ...., ad
Reaction Process of Individual Entities. Next, to capture the individual difference of each entity in evolving their embeddings, we model a local reaction process for each entity, fo, (1), L), where f(-) a neural network (NN), and 0g are the NN (reaction) parameters for mode-k entities. The metaphor from the chemical physics is as follows. While substances are being diffused across different sites, at each site a chemical reaction process happens concurrently, which varies the concentration locally. We extend the model joint diffusion-reaction process,
Ou()) = - A)U(() + F(U,t,, u(0) = Uu.
(1)
Øt

---

where F( = fee (u), fol.., (u), t).................. fex (u (uf,t) .t), .t.................... fex (u t)]T.
Entry Value Generation. Given the embedding trajectories, obtain the tensor entry value me arbitrary time I, we model me(() as a function of the relevant embeddings at time t,
me(() 9 (u), (4).
K 1x (1))
While one can follow (Xu et al. model needs to compute giant kernel matrix over all the observed entry values (see Sec. (2. It is computationally too expensive or infcasible when the number of observations is large. Hence one has to seek for complex low-rank approximations. To avoid this problem. model g with another neural network, which is not only as ficxible as GP, but is more scalable and convenient for computation. Since now, the input to g(-) consists of the trajectory values, which vary with time, our NN model for can flexibly capture the complex temporal relationship between the entities. We finally sample the observed entry values with Gaussian noise model, p(y [m] NN(y m, o²I) where y Iy/...... yN] T and m [me, (t1)..... mex (tv) ] We focus on real-valued data in this paper. However, it is straightforward to extend our approach to other types data.
2012 Zhe et al. 2016b) to assign a GP prior over g(-), the GP
4 Model Estimation
Given data D = -1: t1). (LN, tN yN),, the joint probability of our model is
p(3, () = (()) II p(0)) II,
N N (yulg (u/., (ta)..... u (lm)) o I2)
(3)
where 3 is the NN parameters for g. each 0g is the NN reaction parameters for mode-k entities, 7((3) and p(0)) are element-wise standard Gaussian, and y = (y/...... ........)"""
MN)
To obtain the trajectory
values in the Gaussian likelihood of each yn, we need to solve the ODE in ( time t,.
((t,)) ODESolve(26, 0, t,, θ)
(4)
where θ (W, .1....,00r
0x. consists of the ODE parameters. Our goal is estimate θ the initial state Lu the NN parameters B, and the noise variance o².
Stratified Mini-Bateh Sampling, We use stochastic mini-batch optimization to maximize the log joint probability so as estimate all the required parameters,
C = r(B, (0x). = log(Prior) -
n=1 log.N V (yn (xn), o01)
where log(Prior) log r(3) + Σ p(k), and x,, = (ui., (ln), u (tn) Each time, sample mini-batch of observations B, and obtain unbiased stochastic estimate of the log probability, _ = (Pgiorior) 38 Enes ((vll((x)) o²)]. We compute VL as the stochastic
gradient to update all the parameters.
For each data point n in the mini-batch, need to run ODE solving (ul., (tw).......... u (tm)) To back-propagate the gradient so as to compute the gradient w.r.t the ODE parameters and initial state 24, we can either construct a computational graph during the running of the solver (e.g., the Runge-Kutta method (Dormand and Prince| 1980)), or use the adjoint state method (Pontryagin) 1987, Chen et al. (2018) that solves an adjoint backward ODE to compute the gradient. In whichever case, we need to sort the time points in the mini-batch and solve the ODE sequentially for these time points. As a result, the number of unique time points in the mini-batch greatly influences the speed of processing the mini-batch The standard mini-batch sampling (based on the training example indices) can result in an uneven allocation of the computational cost across the mini-batches some mini-batch is fast and some including more unique time points is much slower. To address this issue, we use a simple stratified sampling approach.
obtain X,,
We collect the unique time points in the whole dataset, T = {TI, .2....} at the beginning. To conduct each stochastic update, we first sample B unique time points from T, then for each time point T₃ E C, we look at all the observed entry values produced at Tj. namely
Dr, [(Ln, lnx Un) E D)t,, = T;)-
. We randomly sample example from each Dry to collect the mini-batch B.
5

---

In this way, we ensure the cost of running ODE solvers and related gradient computation in each mini-batch is identical. There are no fluctuations in cost/running time when processing different min-batches. Empirically, we found that the overall speed of our method is much faster than vanilla stochastic mini-batch optimization (see Table 2 in Appendix).
5 Related Work
There have been many works tensor decomposition, such YYang and Dunson
2013 Rai
iet al.) 2014; [Zhe et al.] 2015 2016b a Tillinghast et al., [020]] Pan et al. 2020b Fang et al.
2021a bj (Tillinghat and Zhe, 2021 Tillinghast al. [022]] Fang et al. 2022). To integrate temporal
information. most approaches augment the tensor with a time mode (Xiong et al. 2010 Rogers et al.
2013 Zhe et al., 2016b Ahn et al.) 2021 Zhe al., 2015 Du et al.) 2018). which includes a list
of time steps. To estimate the temporal dependencies, existing methods often employ dynamic model over the time steps. such conditional Gaussian prior (Xiong et al. 2010), recurrent neural networks Wu et al., 2019), and kernel
(Ahn et al. 202T). To leverage continuous timestamps, Zhang et al. (202T) modeled the CP coefficients as a time function, and Fang et al. (2022) modeled the tensor-core as time function in the Tucker decomposition. The most recent work WWan et al., 2022) places a GP prior in the frequency domain, and construct a bi-level GP model to learn factor or embedding trajectories as a combination of Fourier bases.
Another set of works factorize the interaction events (Schein et al.) 2015 2016; Zhe and Du 2018 Pan et al. 2020a] Wang et al.] 2020; Pan et al., [021]] Wang et al. (022), and they cannol predict the interaction results (e.g., tensor entry values like purchase amount and product ratings). These works mainly leverage Poisson processes, Hawkes processes, or more general point processes to estimate the event rate. Like the standard tensor factorization, these methods also estimate static embeddings for the event participants.
Our model can be viewed only employ the reaction process for each entity, our model is a latent neural ODE (we have an additional NN that combines the latent trajectories predict the tensor entry values). However, we further leverage the structural knowledge in data to construct a multi-partite graph so as to encode the correlations of the entities. Based on the multi-partite graph, we construct diffusion process to co-evolve the embeddings of the entities. In doing so, we can better overcome the data sparsity issue. The most recent work (Li et al., 2022)) uses neural ODE to model the tensor entry value as a function of the involved embeddings and time. It differs from our method in that (1) it still learns static ttime-invarian structure knowledge within the tensor data like our method. Many other works have developed graph diffusion processes based on graph data. For example,Chanbeed (Chmmberlain et al. (202T) proposed graph neural network (GNN) by using multi-head attention to construct the adjacent matrix for graph diffusion equation over the graph nodes. Atwood and Towsley (2016) introduced a diffusion operator to develop di wwfffusionn-oovoluuiinn
extension of the neural ODE model (Chen al. 2018). If
embedding for each entity, and (2) its modeling does not use the
neural networks. |Huang et al. (2021) developed a GNN-based ODE to model both the nodes and edges in dynamic graphs.
Experiment
6.1 Simulation Study
We first examined DEMOTE on synthetic task. Specifically, considered two-mode tensor, where each mode includes 20 entities. Each entity has one underlying embedding trajectory. In the first mode, the trajectory of each entity is an exponential function, u/. (t) = of tx (1 J < 20), while in the second mode a linear function, (t) = o + 2neet.. We generated two clusters of trajectories for each mode, where those of the first ten entities form the first cluster and the remaining the second cluster. To this end, for mode we sampled the coefficients of the first ten entities' trajectories, namely, [c],.... cio] from A (i-5..... -5] 0.1I). and the remaining ten coefficients (ch,.... c from /' ((0.5.... 0.5] 0.1I). Then for mode 2, we sampled each coefficient of conditioned on its counter-pan for mode cgic) Ne c/ 0.1). That means, the coefficients of cluster- trajectories across the two modes are close, and so are those of cluster-2 trajectories. The value for particular entry& l = (2, is generated by
mt(t) = (ut, (t))
2(12 tla mod 2-0) (uf, (t)) 1(1 t1 mod 2=1)
(5)

---

(a) Ground-truth: Mode
(b) NONFAT: Mode
(c) DEMOTE: Mode 1
(d) Ground-truth: Mode
(e) NONFAT: Mode
(i) DEMOTE: Mode 2
Figure 2: The estimated embedding trajectories for each mode. The color indicates the ground-truth cluster membership.
where (.) is the indicator function. When 4 + la is even, the entry value is the trajectory value of the first entity: otherwise, it is the trajectory value of the second entity. To generate the training data, we randomly sampled entries from {41, 42] 41 I 10} U {(1,, (2)] S 4, I 203 (namely, interactions between cluster-1 entities of the two modes, and between cluster-2 entities). We then sampled Unifrom [0, 5). to obtain the corresponding entry values. We randomly generated 6,400 entry values and the timestamps for training, and another 1,600 data points for testing.
We implemented (https: / /tg/ggttttggg tuub. com wwtpe//ggiiitneeeeeeeeeeee chen /tteenttrchhii f feq) to solve ODEs and to compute the gradient w.r.t ODE parameters and initial states via automatic differentiation. For the NN of the reaction process, we used one hidden layer, with 10 neurons and tanh activation, and for the NN to predict the interaction result, used two hidden layers, 50 neurons per layer and tanh activation.
method with Pytorch (Paszke et al.] 2019). We used torchdi ffeq library
We compared with NONFAT (NONparametric Factor Trajectory learning) (Wang et al. 2022), bi-level latent GP model that uses Fourier bases to estimate factor trajectories for dynamic tensor decomposition. To our knowledge, this work is the only method (and also the most recent) that estimates trajectories. We used the original implementation (https: / / q thub. com. wwzut /NONFAT) and the default settings. We set the mini-batch size to 50, and used ADAM (Kingma and Ba 2014) algorithm for stochastic optimization. The learning rate at eau scheduler (Al-Kababji et al. 2022). The maximum
automatically adjusted
in 10 10 by the Reduce number of epochs is 2K, which is enough for convergence. The estimated trajectories are shown in Fig. f As we can see, our estimation (Fig. 2c and well matches the ground-truth and accurately recovers the cluster structure of the trajectories. The root-mean-squar error (RMSE) on the test set is 0.032. By contrast, although the test error of NONFAT is close to DEMOTE (0.034), its learned trajectories (Fig. 2b and 2ep are far from the ground-truth, and fail to reflect the cluster structure. These have shown the advantage of DEMOTE in capturing complex relationships within data to recover the underlying trajectories and their structure.
6.2 Prediction Accuracy
Datasets. We next evaluated the predictive performance of DEMOTE in three real-world applications. (1) CA Weather (Moosavi al.]]20. 2019) (https: smoosavi org/ dataasts/ssssssso sesssssssssssssssseeeeeeesssssst weather conditions in California from August 2016 to December 2020. We extracted a four-mode tensor for 7 different weather types, severity levels, 30 latitudes and 30 longitudes GPS coordinates. The entry value is the count of the particular weather condition. We collected 15K observed tensor entry values and the timestamps. (2) CA Traffic (Moosavi et al., 2019) (https: htttps://smoosawv smoosavvi corg/dataset s/1stw), traffic accidents in California from January 2018 to December 2020. We extracted four
7

---

Reaction
Diffusion
DEMOTE
0.08
0.5
nAMSE
0.7
NRMSE
NRMSE
0.4
0.06
0.65
0.3
0.04
0.6
3
5
R
(a) CA Weather
(b) CA Traffic
(c) Server Room
Figure 3: Predictive performance of the diffusion and reaction processes.
mode tensor (traffic rype, severity level, latitudes, longitude). There are traffic types, 6 severity levels, 20 latitudes and 20 longitudes. We collected 30K entry values (accident counts) at different time points. (3) Server Room (https: zenodo. http://...rrrreeee////eee./
..e.INN XINpA azaM8) temporal temperature records of Pozman Supercomputing and Networking Center. The temperatures were measured at 34 locations, under different air-condition modes (24", 27" and 30") and power usage settings (50%. 75% and 100%). ). Hence, we extracted a three-mode tensor (location, air- condition mode, power level). In total, 10K observed entry values and their timestamps were collected.
Competing Methods. The following popular and/or state-of-ther temporal decomposition ap- proaches were compared. (1) CP-DTLD, discrete-time CP decomposition with linear dynamics, where a conditional prior placed over successive time steps, p(t)+1/t,) = ( At, b. eI): A, b and u were jointly estimated during the CP decomposition. Note that (Xiong et al. 2010) is an instance this model where A and b O. (2) GP-DTLD and (3) NN-DTLD, similar to CP-DTLD, except using GP (Zhe et al. 601 and NN decomposition models (similar to (Liu et al. 2019h, respectively. (4) CP-DTND, (5) GP-DTND and (6) NN-DTND - CP, GP and NN decompo- sition with nonlinear dynamics, where the conditional prior is p( (ittttltt).. (t |or (tt,) +,,,)bbu, vl) where (.) is a nonlinear activation. The dynamics can therefore be viewed as an RNN transition. (7) CP-CT Zhang al., 2021), continuous time CP factorization, which models the CP coefficients vine-varying function, with polynomial splines. (8) GP-CT, continuous-time GP decomposition that extends (Xu et al., 2012 Zhe al. (0166) by plugging the time in the GP kernel as estimate the entry value as a function of the embeddings and time, m = g(u! , u t). (9) NN-CT, continuous- time NN decomposition, where the input consists of both the embeddings and time t. (10) THIS-ODE (Li et al. 2002)), a continuous-time decomposition, where a neural ODE is used to estimate the tensor entry values given the static embeddings and time. (11) NONFAT (Wang et al.)
2022]. a bi-level latent GP model that uses Fourier bases to estimate factor trajectories for dynamic tensor decomposition.
Settings and Results. All the approaches were implemented with PoTorch. The Square Exponential kernel was used for all the GP-related methods, including GP-(DTLD, DTND, CT). We used the same variational sparse approximation (Hensman et al.) 2013) to fulfill scalable posterior inference. Following (Zhe et al., 20166), the number of inducing point was set 100. For the NN decomposition methods, we employed a three-layer network with tnhh activation, and for THIS-ODE, we used one-layer network. The layer width chosen from (10, 25, 50, 75, 100). We used tanh as the activation function in the nonlinear dynamic baselines, including (CP, GP, NNJ-DTND. For our method, used the same NN architecture for both the reaction process and entry value prediction, which includes two hidden layers with 50 neurons per layer. For CP-CT, we employed 100 knots to fulfill the polynomial splines. For each discrete-time method, the number of time steps was chosen from (25, 50, 75, 100] via the gross-validation on the training set. We trained all the models with stochastic mini-batch optimization. We used the ADAM algorithm, and the mini-batch size was set to 100. We ran every method with 10K epochs to ensure convergence. The learning rate was automatically adjusted in [10-* 10011 by the ReduceLRon coomplateaal eau scheduler. We varied the dimension of the embeddings R from (2. 3, 5, 7). For DEMOTE, R is the number of embedding trajectories; we used computational graphs to obtain the gradient. We followed (Kang et al., 2012 Zhe et al. 201 to randomly draw 80% observed entries and their time stamps for training, with the remaining for test. We computed the normalized noot-mean-square enror (nRMSE). We repeated the evaluation five times and computed the average nRMSE and standard deviation.

---

CA Weather
R
R 3
R =
R =
CP-DTLD
0.7440 0.0035
0.7372 0.0040
0.7290 + 0.0042
0.7270
0.0044
GP-DTLD
0.7417 + 0.0031
0.7414 + 0.0036
0.7444 + 0.0036
0.7449 + 0.0039
NN-DTLD
0.7228 1: 0.0054
0.7116 + 0.0033
0.7070 1 0.0041
0.7065 1: 0.0038
CP-DTND
0.7448 1 0.0031
0.7360 0.0035
0.7273 1 0.0037
0.7280 1 0.0044
GP-DTND
0.7399 0.0034
0.7346 1 0.0032
0.7448 1: 0.0037
0.7467 0.0031
NN-DTND
0.7113 0.0045
0.6979 = 0.0126
0.6659 = 0.0122
0.6543 0.0155
CP-CT
1.0000 + 0.0096
0.9959 = 0.0067
1.0010 = 0.0017
1.0060 = 0.0034
GP-CT
0.7433 + 0.0038
0.7354 = 0.0027
0.7359 = 0.0034
0.7377 = 0.0033
NN-CT
0.8697 0.0014
0.8679 0.0022
0.8676 + 0.0018
0.8695 + 0.0016
NONFAT
0.7444 0.0042
0.7460 = 0.0032
0.7645 = 0.0061
0.7553 0.0029
THIS-ODE
0.7511 + 0.0052
0.7539 = 0.0011
0.7614 + 0.0024
0.7620 + 0.0032
DEMOTE
0.6327 + 0.0119
0.6109 = 0.0056
0.6172 = 0.0075
O. 0.6344 0.0085
CA Traffic
CP-DTLD
0.6498 + 0.0257
0.6424 0.0266
0.6436 = 0.0268
0.6405 - 0.0262
GP-DTLD
0.6309 0.0167
0.6290 = 0.0185
0.6383 = 0.0204
0.6496 0.0193
NN-DTLD
0.6528 0.0230
0.6545 + 0.0244
0.6401 = 0.0282
0,6136 +: 0.0338
CP-DTND
0.6497 + 0.0245
0.6456 = 0.0265
0.6431 = 0.0263
0.6419 + 0.0259
GP-DTND
0.6544 + 0.0213
0.6559 = 0.0224
0.6604 = 0.0243
0.6674 + 0.0214
NN-DTND
0.6578 0.0248
0.6528 0.0256
0.6519 = 0.0249
0.6482 0.0261
CP-CT
0.9858 0.0120
0.9972 + 0.0056
0.9816 1 0.0136
0.9991 0.0120
GP-CT
0.6610 0.0207
0.6668 1 0.0191
0.6756 1 0.0190
0.6768 1: 0.0196
NN-CT
0.9804 +: 0.0017
0.9815 + 0.0015
0.9791 + 0.0012
0.9802 1 0.0017
NONFAT
0.4461 + 0.0247
0.4610 0.0231
0.5031 = 0.0155
0.6307 = 0.0847
THIS ODE
0.6603 0.0230
0.6536 = 0.0212
0.6838 = 0.0193
0.6378 0.0142
DEMOTE
0.3601 0.0334
0.2972 = 0.0099
0.3174 = 0.0118
0.3269 - 0.0162
Server Room
CP-DTLD
0.4211 + 0.0029
0.4209 = 0.0031
0.4208 = 0.0028
0.4208 = 0.0028
GP-DTLD
0.0914 + 0.0020
0.0791 + 0.0010
0.0739 + 0.0014
0.0753 0.0013
NN-DTLD
0.4213 0.0032
0.4213 + 0.0032
0.4212 + 0.0034
0.4205 + 0.0030
CP-DTND
0.2835 0.0160
0.1751 + 0.0020
0.1174 + 0.0011
0.0829 0.0044
GP-DTND
0.0925 0.0013
0.0784 + 0.0011
0.0739 1 0.0009
0.0774 1: 0.0009
NN-DTND
0.4213 1 0.0032
0.4212 + 0.0030
0.4211 + 0.0032
0.4205 1: 0.0030
CP-CT
0.9919 0.0096
0.9951 + 0.0050
0.9862 1 0.0109
1.0121 1: 0.0070
GP-CT
0.1385 0.0020
0.1223 = 0.0016
0.1275 0.0014
0.1365 0.0014
NN-CT
0.1193 0.0030
0.1140 0.0015
0.1113 0.0027
0.1149 - 0.0028
NONFAT
0.1468 + 0.0026
0.1407 = 0.0023
0.1396 = 0.0022
0.1409 + 0.0030
THIS-ODE
0.1412 + 0.0021
0.1312 = 0.0013
0.1304
1 0.0016
0.1350 + 0.0019
DEMOTE
0.0536 0.0031
0.0403 0.0014
0.0393 + 0.0018
0.0403 + 0.0027
Table 1: Normalized Root Mean-Square Error (nRMSE). The results were averaged from five runs.
As shown in Table DEMOTE consistently achieves the best prediction accuracy and in many cases outperforms the competing methods by a large margin. Although learning an embedding trajectory much more challenging than learning fixed-value embedding, the experimental results have demonstrated the ad advantage of our method in predictive performance. To investigate the effect of the two processes our model, with the reaction process only on all the datasets. Their performance, as compared with DEMOTE, is shown in Fig. 3 We can see that each individual component can lead to good prediction accuracy. However, each component is worse than DEMOTE that synergizes the two components together. Therefore, the results show that each process effective, and more important, the two processes can bolster each other to further improve the performance when they are combined.
also examined our method with the diffusion process only and
6.3 Learning Result Investigation
Next, we looked into the learned embedding trajectories and checked if they exhibit patterns. To do so, we set and ran DEMOTE on Server Room dataset. In Fig. 4 we show the learned embedding trajectories for the first location (a-c), the first air condition mode (d-f) and the first power usage level (g-i). we can see, even for the same object, e.g., a particular location, the corresponding embedding trajectories vary quite differently, implying the evolution of different underlying properties, such as the workload, memory usage, and network latency. We further found there are underlying structures within the embeddings during their evolution. We listed the results in Appendix (Sec. A).

---

0.5
1.5
Time
Time
Time
(a) (t)
(b) (())
(c) (())())
-2
10
Time
Time
Time
(d) uf.1 (1)
(e) (1)
(1) (1)
1.5
0.5
Time
Time
Time
(g) (t)
(h) (t)
(i) (t)
Figure 4: The learned embedding trajectorics for location (a-c), air conditional mode (d-f), and power usage
level (g-i) in Server Room datasct.
We showcase the temporal predictions for tensor entries, As we can see from Fig. 5 given only a few training points (blue), our method can predict the test points (green) much more ac- curately, as compared with GPCT, and the predictive uncertainty (reflected by the noise variance o²) is much smaller. This might be due to that via the diffusion- reactionn process, and our method can more effectively extract the temporal knowledge from sparse data. For example, DEMOTE successfully captured the periodic nature in the first entry (Fig. (D)) while GPCT treated the fluc- tuation as noises and ended up with much inaccurate predictions and larger predictive variances. Computational Efficiency. We compared the per- epoch/iteration running time DEMOTE with the other methods. We tested all the methods in workstation with one NVIDIA GeForce RTX 3090 Graphics Card, 10th Generation Intel Core 19- 10850K Processor, 32 GB RAM. and TB SSD. The results shown in Table in Appendix. We can see that the running speed of DEMOTE is com- parable to NONFAT and other NN decomposition methods. We also compared with running DE- MOTE with naive sampling (DEMOTE-NS). The stratified sampling led to 4x to 22x speed-up.
DEMOTE GF Trainino Data Testing
(a) (2, 1, 1)
(b) (3, 1, I)
Figure Entry value prediction on Server Room.
Conclusion
We have presented DEMOTE, neural diffusion-reaction process model to learn dynamic embeddings for dynamic tensor decomposition. The predictive performance encouraging, and the learned embedding trajectories exhibit interesting patterns. Currently, our method is limited to small number of entities since it has to integrate the entire multi-partite graph construct the diffusion process. In the future work, we plan to develop graph cut algorithms to partition the graph into set of small sub-graphs so that we can construct multiple diffusion processes parallel so as to scale up our model to big graphs and to large tensors.
Acknowledgments
This work has been supported by NSF CAREER Award IIS-20466 10-2006..55.
10

---

References
Ahn, D., Jang, J.-G., and Kang, U. (2021). Time-aware tensor decomposition for sparse tensors. Machine Learning, pages 1-22.
Al-Kababji, A., Bensaali, F., and Dakua, S. P. (2022). Scheduling techniques for liver segmentation: Reducelronplateau vs onecyclelr. arXiv preprint arXiv:2202.0637.
Atwood, J. and Towsley D. (2016). Diffusiion-comvlutiional neural networks. Advances in neural information processing systems, 29.
Chamberlain, B., Rowbottom, J., Gorinova, M. L., Bronstein, M.. Webb, S., and Rossi, E. (2021). Grand: Graph neural diffusion. In International Conference Machine Learning, pages 1407- 1418. PMLR.
Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K. (2018). Neural ordinary differential equations. Advances neural information processing systems, 31.
Dormand, J. R. and Prince, P. J. (1980). A family of embedded runge-kutta formulae. Journal of computational and applied mathematics, 6(1):19-26.
Du, Y., Zheng, Y., Lee, K.-c,, and Zhe, S. (2018). Probabilistic streaming tensor decomposition. In 2018 IEEE International Conference on Data Mining (ICDM), pages 99-108. IEEE.
Fang. S., Kirby, R. M., and Zhe, S. (2021a). Bayesian streaming sparse Tucker decomposition. In Uncertainty in Artificial Intelligence, pages 558-567. PMLR.
Fang, S., Narayan, A., Kirby, R., and Zhe, S. (2022). Bayesian continuous- time Tucker decomposition.
In International Conference Machine Learning, pages 6235-6245. PMLR.
Fang. S., Wang, Z., Pan, Z., Liu, J., and Zhe, S. (2021b). Streaming Bayesian deep tensor factorization. In International Conference Machine Learning, pages 3133-3142. PMLR.
Harshman, A. (1970). Foundations of the PARAFAC procedure: Model and conditions for an" exxxplaaator"" "astoooy"mmuiiiiii mmulti-mod factor analysis. UCLA Working Papers in Phonetics, 16:1-84.
Hensman, J., Fusi, N., and Lawrence, N. D. (2013). Gaussian processes for big data. In Proceedings of the Twenty-Ninth Conference on Uncertainty in Artificial Intelligence, pages 282-290. AUAI Press.
Huang. Z., Sun, Y., and Wang, W. (2021). Coupled graph ode for learning interacting system dynamics. In 27th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, KDD 2021, pages 705-715.
Kang, U., Papalexakis, E., Harpale, A., and Faloutsos, C. (2012). Gigatensor: scaling tensor analysis up by 100 times-algorithms and discoveries. In Proceedings of the 18th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 316-324. ACM.
Ketchen, D. and Shook, C. L. (1996). The application of cluster analysis in strategic management research: an analysis and critique. Strategic management journal, 17(6):441-44
Kingma, D. P. and Ba, J. (2014). Adam: A method for stochastic optimization. cite Published as a conference paper at the 3rd International Conference for Learning Representations, San Diego, 2015.
Kolda, T. G. (2006). Multilinear operators for higher-order decompositions, volume 2. United States. Department of Energy.
Li, S., Kirby, R., and Zhe, S. (2022). Decomposing temporal high-order interactions via latent odes. In International Conference Machine Learning, pages 12797-12812. PMLR.
Liu, H., Li, Y., Tsang, M., and Liu, Y. (2019). CoSTCo: Neural Tensor Completion Model for Sparse Tensors, page 324-334. Association for Computing Machinery, New York, NY, USA.

---

Moosavi, S., Samavatian, M. H., Nandi. A., Parthasarathy, S., and Ramnath, R. (2019). Short and long-term pattern discovery over large-scale peo-spatiotemporal data. In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, pages 1905-29133
Pan, Z.,, Wang, Z.,, Phillips, J. M., and Zhe, S. (2021). Self-adaptable point processes with nonpara- metric time decays. Advances in Neural Information Processing Systems, 34:4594-4606.
Pan, Z., Wang, and Zhe, (2020a). Scalable nonparametric factorization for high-order interaction events. In International Conference Artificial Intelligence and Statistics, pages 0225-333. PMLR.
Pan, Z., Wang, Z., and Zhe, S. (2020b). Streaming nonlinear Bayesian tensor decomposition. In Conference on Uncertainty Artificial Intelligence, pages 190-499. PMLR.
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein,
N., Antiga, L., et al. (2019). Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32:8026-8037.
Pontryagin, L. S. (1987), Mathematical theory of optimal processes. CRC press.
Rai, P., Wang, Y., Guo, S., Chen, G., Dunson, D., and Carin, L. (2014). Scalable Bayesian low- rank decomposition of incomplete multiway tensors. In Proceedings of the 31th International Conference on Machine Learning (ICML).
Rasmussen, C. E. and Williams, C. K. (2006). Gaussian Processes for Machine Learning. MIT Press.
Rogers, M., Li, L., and Russell, S. J. (2013). Multilinear dynamical systems for tensor time series. Advances in Neural Information Processing Systems, 26:2634-244 2642..
Rubanova, Y., Chen, R. T. Q., and Duvenaud, D. K. (2019). Latent ordinary differential equations for
iregularly-sampler time series. In Wallach, H., Larochelle, H., Beygelzimer, A., d' l'Alché-Bo Buc, F, Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.
Schein, A., Paisley, J.. Blei, D. M., and Wallach, H. (2015). Bayesian poisson tensor factorization for inferring multilateral relations from sparse dyadic event counts. In Proceedings of the 2 ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 1045-1054, ACM.
Schein, A.. Zhou, M., Blei, D. M., and Wallach, H. (2016). Bayesian poisson tucker decomposition for learning the structure of international relations. In Proceedings of the 33rd International Conference on International Conference on Machine Learning Volume 48, ICML'16, pages 2810-2819. JMLR.org.
Song, Q., Huang, X., Ge, H., Caverlee, J., and Hu, X. (2017). Multi-aspect streaming tensor completion. In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 435-443.
Tillinghast, C., Fang, S., Zhang, K., and Zhe, S. (2020). Probabilistic neural-kernel tensor de- composition. 2020 IEEE International Conference on Data Mining (ICDM), pages 531-540. IEEE.
Tillinghast, C., Wang. Z., and Zhe, S. (2022). Nonparametric sparse tensor factorization with hierarchical Gamma processes. In International Conference on Machine Learning, pages 21432- 21448. PMLR.
Tillinghast, C. and Zhe, S. (2021). Nonparametric decomposition of sparse tensors. In International Conference Machine Learning, pages 10301-10311. PMLR.
Tucker, L. (1966). Some mathematical notes on three-mode factor analysis. Psychometrika, 31:279- 311.
12

---

Wang, Z., Chu. X., and Zhe, S. (2020). Self-modulating nonparametric event-tensor factorization. In International Conference on Machine Learning, pages 9857-9867. PMLR.
Wang, Z., Xu, Y., Tillinghast, C., Li, S., Narayan, A., and Zhe, S. (2022). Nonparametric embeddings of sparse high-order interaction events. In International Conference on Machine Learning, pages 23237-23253 PMLR.
Wu, X., Shi, B., Dong, Y., Huang. C., and Chawla, N. V. (2019). Neural tensor factorization for temporal interaction learning. In Proceedings of the Twelfth ACM International Conference Web Search and Data Mining, pages 537-545.
Xiong, L.. Chen, X., Huang, T.-K., Schneider, J., and Carbonell, J. G. (2010). Temporal collabora- tive filtering with bayesian probabilistic tensor factorization. In Proceedings of the 2010 SIAM International Conference on Data Mining, pages 211-222. SIAM.
Xu, Z., Yan, F,, and Qi, Y. A. (2012). Infinite tucker decomposition: Nonparametric bayesian models for multiway data analysis. In ICML.
Yang. Y. and Dunson, D. (2013). Bayesian conditional tensor factorizations for high-dimensional classification. Journal of the Royal Statistical Society B. revision submitted.
Zhang, Y., Bi, X., Tang, N., and Qu, A. (2021). Dynamic tensor recommender systems. Journal of Machine Learning Research, 22(65):
Zhe, S. and Du Y. (2018). Stochastic nonparametric event-tensor decomposition. In Advances in Neural Information Processing Systems, pages 6856-6866
Zhe, S., Qi, Y., Park, Y., Xu, Z.. Molloy, I., and Chari, S. (2016a). Dintucker: Scaling up Gaussian process models on large multidimensional arrays. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 30.
Zhe, S., Xu, Z., Chu, X., Qi, Y., and Park, Y. (2015). Scalable nonparametric multiway data analysis. In Artificial Intelligence and Statistics, pages 125-1134. PMLR.
Zhe, S., Zhang, K., Wang, P., Lee, K.-c,, Xu, Z., Qi, Y., and Ghahramani, Z. (2016b). Distributed flexible nonlinear tensor factorization. Advances in neural information processing systems, 29.
13

---

Appendix
A Investigation of Embedding Dynamics
We investigated if there are underlying structures within the embeddings during their evolution. To this end, looked into the embeddings of the 34 locations Server Room dataset five time points (l = 1, 20, 50, 80, 100). At each time point, we ran the k-means algorithm over the embeddings extract the clustering structures. We used the elbow method (Ketchen and Shook T996) to select the number of clusters. We can see that at carlier time (t < 50), the clusters are more compact, while the later stages, the clusters become more scattered. This reflects how the structure of those entities (i.e., locations) evolves along with time. It is interesting to see that some locations are in the same cluster all the time, like location (5.7) and location (16, 32). implies that their underlying properties might have quite similar (or correlated) evolution patterns. Some locations are grouped in the cluster at the beginning, e.g., location (32, 34] (at 1), but later moves to different clusters (t >> 1). It implies their evolution patterns can vary significantly. leading to the change of the cluster memberships.
CA Weather
CA Traffic
Server Room
CP-DTLD
0.037
0.086
0.023
GP-DTLD
0,246
0,247
0.248
NN-DTLD
2,400
4,730
1.080
CP-DTND
0.038
0.087
0.025
GP-DTND
0.119
0.242
0.080
NN-DIND
2,360
4,701
1.060
CP-CT
0.025
0.052
0.018
GP-CT
0.068
0,216
0.105
NN-CT
2.310
3.885
1.030
NONFAT
0.952
1,925
0.571
THIS-OD
58.710
136.100
7.190
DEMOTE
1.390
1.895
0.309
DEMOTE-NS 6.12
10.42
7.06
Table 2: Per-epoch/ierration running time (in seconds). DEMOTE-NS means running DEMOTE with naive sampling of min-batches rather than the stratified sampling.
14

---

219
#28
63
X29
#30
31 1
29
xT028
*17
#32
15
*
Cluster
X32
*16
25
O Cluster 2
X36
Cluster
15
C27
25
C24
13
*2
22
95
27
x12
X24
13
C23
23
X26
7
(a) ===1
(b) f 20
9
28
C7
34
17
22
23 24
C25
C609
O1
10
3
C5
29
#3066
4
C27
X32
155
C1333
24
C83
$4428
* 28
44
#23
CR
1
C5
29
25
219
C6
(c) f 50
(d) f 80
#3
*31
*17
*29
*15
BB
L3
19
11
*4
X28
*44
6
25
C24
9
C8
5
26
(c) t 100
Figure 6: Evolution of the clustering structure within the 34 locations on Server Room dataset.
15


Reviews:
[{"review_id": "o9All78gCE", "reviewer": "Reviewer_9boz", "summary": "This paper proposes a dynamic tensor decomposition (TD) model, in which the latent factors are functions of time stamps described by neural ODEs. Moreover, they encode an adjacent matrix into the diffusion process to better model the correlations among different tensor modes. Finally, they use an additional neural network to map latent factors to tensor entries. Using differentiable ODE solvers, all the parameters can be learned by optimizing the joint likelihood with SGD. For empirical evaluations, the authors first conduct a simulation study to show that the proposed model can learn latent factor trajectories. Then, they demonstrate the completion performance on several dynamic tensor datasets.", "questions": "1. In the equation between Line 178 and 179, does $\\mathcal{L}$ equal to $\\log (Prior)$ **plus** the likelihood?", "limitations": "None", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The authors propose to use a diffusion process to learn tensor factor trajectories, which is new in the field.\n\n2. The experiments show that the proposed model can learn factor and entry trajectories. Besides, it achieves a significant improvement on the completion performance.", "weaknesses": "1. In the paper, the authors construct a complex TD with ODEs and NN to parameterize the latent relation. Besides, the model parameters are learnt from MAP estimates (Eq. 3). I am wondering if such a model suffers from overfitting for sparse tensor data?\n\n2. In Equation 1, the authors define the ODE of the factor trajectories. Since $\\mathcal{F}(\\mathcal{U}, t)$ is already parameterized by flexible neural networks, I am wondering if it is necessary to add the multi-partite graph term, i.e., $(\\mathcal{W} - \\mathcal{A})\\mathcal{U(t)}$? Can we construct a specific network structure that is equivalent to Eq. 1?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their valuable comments. Here are our responses. C: comments; R: response.\n\n>C1: \"the authors construct a complex TD with ODEs and NN to parameterize the latent relation. Besides, the model parameters are learnt from MAP estimates (Eq. 3). I am wondering if such a model suffers from overfitting for sparse tensor data?\"\n\nR1: Great question. We do agree solely using the NN components can bring the risk of overfitting. However, our work also constructs a multi-partitie graph, and integrates the graph into the embedding trajectory learning, namely, the diffusion process. As one major contribution of our work, the diffusion process not only incorporates the structural knowledge of the data, it also uses **the structural knowledge to regularize the learning of the embedding trajectories**. Specifically, it encourages homophily of the trajectories for neighboring nodes (as they are correlated) --- the change rate of $\\mathbf{u}^k_j(t)$ is modeled to correspond to the difference with the neighborhood's trajectories; see Line 140-144 and the equation in between. Therefore, we believe the diffusion process can effectively inhibit overfitting for sparse data. Our experimental results have shown that our method performs much better than just learning each trajectory independently (with the same NN components, i.e., no diffusion, just reaction process); see Fig. 3.\n\n>C2: \"In Equation 1, the authors define the ODE of the factor trajectories. Since $\\mathcal{F}(\\mathcal{U},t)$ is already parameterized by flexible neural networks, I am wondering if it is necessary to add the multi-partite graph term, i.e., $(\\mathcal{W}-\\mathcal{A})\\mathcal{U}(t)$ Can we construct a specific network structure that is equivalent to Eq. 1?\"\n\nR2: Thanks for the great question. Please see R1 for our explanation about why the multi-partite graph term is important --- we believe it not only leverages the structural knowledge to facilitate the trajectory learning, but also induces a proper regularization to avoid overfitting.\n\nRegarding ``a specific'' equivalent ''network structure'', yes, we can view the RHS of Eq1 as modeled by a Residual Network. That is, the output of this network is a linear transformation of $\\mathcal{U}$ (i.e., the multi-partite graph term) **plus** an NN transform, namely $\\mathcal{F}(\\mathcal{U},t)$.  \n\n\n>C3: \"In the equation between Line 178 and 179, does $\\mathcal{L}$ equal to $log(Prior)$ plus the likelihood?\"\n\nR3: Yes, exactly. We appreciate the reviewer catching this typo. We will fix it in the paper."}, {"review_id": "ebRLf7ixAN", "reviewer": "Reviewer_gF6D", "summary": "In this paper, to capture higher-order interaction among different modes of tensor data and their evolution process over time, the authors proposed an ODE-based method called DEMOTE to capture the dynamics of latent embeddings.   In detail, the authors propose to use a diffusion-reaction process to model the evolution of these embeddings.  What's more, structure knowledge within tensor data is considered in this diffusion-reaction process, which enhances the model's representation ability.  Experiments show the effectiveness of the proposed method.", "questions": "1. I think the contribution should be further clarified.", "limitations": "Yes, the authors have adequately addressed the limitations and potential negative societal impact of their work.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "1. The experimental results show an advantage compared with related work.\n\n2. The way of using the structure knowledge is novel and effective in modeling tensor data. It is particularly interesting to embed this structured knowledge with the framework of ODE.", "weaknesses": "1. The main concern lies in its novelty since there are some efforts studying the interaction results of higher order continuous time data, as clarified in the paper (that is [1]). Actually, in [1], they use ODE to model the entry values' trajectories and depend on $m_{i}$. In this work, the authors use ODE to model $m_{i}$ directly. Although these two papers correspond to different ways of handling problems and the authors have stated the advantage of the proposed method, they do share many similar intuitions towards the same task, which hinders their novelty.\n\nTypo: Page 3, last paragraph: $(v_{j}^{k}, v_{j'}^{k'})\\in \\mathbb{E}$ $\\rightarrow$  $(v_{j}^{k}, v_{j'}^{k'})\\in E$.\n\n[1] Li S, Kirby R, Zhe S. Decomposing Temporal High-Order Interactions via Latent ODEs[C]//International Conference on Machine Learning. PMLR, 2022: 12797-12812.", "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: responses\n\n> C1: \"The main concern lies in its novelty since there are some efforts studying the interaction results of higher order continuous time data, as clarified in the paper (that is [1]) ... Although these two papers correspond to different ways of handling problems and the authors have stated the advantage of the proposed method, they do share many similar intuitions towards the same task, which hinders their novelty\" \"I think the contribution should be further clarified.\"\n\nR1: We appreciate the reviewer's concern. We do agree that both our work and [1] use ODE as  the modeling framework to handle temporal interaction data. While we have clarified the difference in our paper (Line 226-228), here we would like to give more detailed and comprehensive clarification. In fact, the task, model design, algorithm and evaluation of [1] are all very different from our work.\n\nFirst, the task/goal of [1] is **NOT** to estimate dynamic embeddings; it still estimates **static** embeddings (time-invariant), like most tensor decomposition methods do. The latent ODE model of [1] is only designed to predict the temporal interaction result. By contrast, the task of our work is to learn the embedding trajectories to capture the temporal evolution of the representations; we consider the embedding as a function of time. We believe this is scarce in existing works and our work fills an important gap in tensor decomposition.\n\nCorrespondingly, a large portion of our experiment is to examine if the embedding trajectory recovered by our method is accurate (see Section 6.1 and Fig. 2 for Simulation Study) and meaningful/interpretable (see Section 6.3 Fig. 4, and Fig. 1 of Appendix). Such evaluations have never shown up in [1]. \n\nThird, [1] has never leveraged the structural knowledge underlying the observed interactions; [1] constructs a neural ODE model in a straightforward manner. The key hypothesis of our work is that such knowledge is important for trajectory learning. We proposed a novel k-partite graph to represent the structure, based on which we proposed a neural diffusion-reaction process to integrate the structural knowledge into the embedding trajectory learning and interaction result prediction.\n\nFourth, [1] uses forward sensitivity for model training. That is, [1] augments the ODE model with the sensitivity (i.e., gradient), and performs a forward solve to jointly compute the ODE states and state gradient. This is feasible, because [1] models each entry value as an ODE, and the state is just a scalar (one dimensional). However, the learning algorithm of [1] does not apply to our work, because our ODE state includes all the embedding trajectories and is high-dimensional. For efficient learning, we devised a stratified mini-batch sampling method, and combined with the adjoint state method (backward) or computational graph to compute the gradient for training.\n\nThanks for bringing up the question. We will add more discussions about [1] and clarification about the difference of our work."}, {"review_id": "DBYNI9fEPa", "reviewer": "Reviewer_3waT", "summary": "In this work, the authors introduce DEMOTE which can both consider both structure information and temporal information via dynamic tensor decomposition. Structure information is captured by multi-partite graph diffusion process which in substance is a special ODE, whereas temporal information is analyzed using another ODE. Experimental results demonstrate DEMOTE in general outperforms one baseline for synthetic data and many baselines for real-world data.", "questions": "1. For the statement in Page 3 Section 3 Lines 113-115, why can learning dynamics produce more reliable and long-term trajectory predictions? The time derivative can provide the general dynamics of the whole time series and mainly focus on the global sequential information. Please check this paper \"Neural Rough Differential Equations for Long Time Series\" [1].\n\n2. Could the authors provide any explanation and examples on why K should be involved and represent structure information using multi-partite graphs is necessary? It seems that including the multi-partite graphs is very complicated, and in what circumstances the multi-partite graphs should be applied to describe what types of intercorrelation among multi-dimensional data, without incurring overfitting?\n\n3. In Figure 2(c), why do the red lines decline? In the long-term, the performance of DEMOTE is worse than the baseline model NONFAT. Could the authors articulate it?\n\n4. For the results of DEMOTE in Figure 5(a) and (b), why are the entry value predictions with high and sharp spikes across time? It seems that DEMOTE can capture the periodic fluctuations, but DEMOTE by nature is based on two ODEs. In specific, ODE and neural ODE assume that the modeled trajectory should be smooth, or at least Lipschitz continuous [2,3]. Could the authors also provide their code?\n\nReferences\n\n[1] Morrill, J., Salvi, C., Kidger, P., & Foster, J. (2021, July). Neural rough differential equations for long time series. In International Conference on Machine Learning (pp. 7829-7838). PMLR.\n\n[2] Chen, R. T., Rubanova, Y., Bettencourt, J., & Duvenaud, D. K. (2018). Neural ordinary differential equations. Advances in neural information processing systems, 31.\n\n[3] Süli, E., & Mayers, D. F. (2003). An introduction to numerical analysis. Cambridge university press.", "limitations": "The limitations are sufficiently addressed by the authors.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. It is interesting to model the spatial interactions for values at the same time point using multi-partite graphs through graph diffusion, for the small-scale datasets.\n\n2. Addressing the sparsity in tensorial data is another strength of this work. It has been an important problem in tensorial data and appeal to the tensor learning community.", "weaknesses": "1. There have been some existing work on dynamic tensor decomposition using neural ODEs and relational information in the recent years, whereas not discussed and compared in this work, for example, TANGO [1].\n\n2. The presentation should be improved. The notations should be defined in a clear manner, for example, the adjacency tensor W needs more clarification on its order and how it is computed, from its indices. \n\nReference\n\n[1] Han, Z., Ding, Z., Ma, Y., Gu, Y., & Tresp, V. (2021, November). Learning neural ordinary equations for forecasting future links on temporal knowledge graphs. In *Proceedings of the 2021 conference on empirical methods in natural language processing* (pp. 8352-8364).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your insightful comments and the great reference TANGO. We will add discussion and comparison with this nice, related work. \n\n>C1: why can learning dynamics produce more reliable and long-term trajectory predictions? \n\nR1: Great question. To better illustrate it, we can consider an alternative popular model, Gaussian process (GP), to directly learn the trajectories. GP is known to work well when making prediction near the training inputs. However, at places far from the training inputs, the prediction accuracy degrades severely, because the interpolation on the training points becomes ineffective (the similarity to the training inputs drops quickly). Hence, GP cannot provide reliable/long-term trajectory prediction. However, if we use ODE modeling, we learn the dynamics rather than try to interpolate the prediction. As the reviewer pointed out,  \"the time derivative can provide the general dynamics of the whole time series\". Hence, as long as we can well capture the time derivative --- this is what ODE modeling is good at --- we can use it to infer long-term trajectories reliably, no matter if the places are close to or far away from the training inputs.\n\nIn fact, our intuition has been verified by the previous work (Rubanova et al., 2019)  that builds latent ODE models for time series (cited at Line 42). In Fig. 4a and 5b of their paper, they show their ODE models can use a small amount of training observations to well extrapolate to a long time horizon; also see Table 3 their much better quantitative results in extrapolation (Enc-Dec).\n\nThird, thanks for the great reference (Morrill et al., 2021). But this work is about **Neural Control Differential Equation models, i.e., Neural CDE rather than Neural ODE**. Due to a different model design (Neural CDE is a continuous version of RNN while Neural ODE is a continuous version of ResNet), Neural CDEs have an additional, significant learning challenge: one has to densely approximate the time derivative of the observed sequence, which is also integrated into the dynamics. According to their paper, \"Neural CDEs, as with RNNs, begin to break down for long time series. Loss/accuracy worsens, and training time becomes prohibitive due to the sheer number of forward operations within each training epoch\". However, Neural ODEs do not have to face such challenge, and the bottleneck might not apply to them.\n\nWe will add more detailed discussions about the model choice in the paper.\n\n>C2: \"any explanation and examples on why K should be involved and represent structure information using multi-partite graphs is necessary? ...\"\n\nR2: Great questions. First, tensor data are about multiway interactions, e.g., *(customer, seller, product)*, and the entry value is the interaction result, e.g., payment. The presence of these interactions indicate intercorrelations between the involved entities (in each mode). We are motivated to leverage these correlations to enhance learning embedding trajectories in our ODE framework. One elegant way is to design a diffusion process, which *needs to be on a graph*. We therefore convert the observed multiway interactions into a K-partite graph. We believe this graph is reasonable to encode various correlations in the data. For example, if any two entities appear in the same observed interaction (tensor entry), we add an edge between them to indicate their inter-correlation. Note that we do not add edges between the entities in the same mode, because a tensor entry always represents the interaction between different modes. That's why our graph is a multi-partite graph. See Line 122-135 for more detailed explanation/rationale for our graph construction. With the graph, we are then able to design our diffusion-reaction process (Eq1).\n\nSecond, we can choose to remove the graph and hence remove the diffusion process in learning. However, this will treat the embedding trajectory of each entity independently and bring the risk of overfitting. The K-partite graph and the associated diffusion process not only bring additional structural knowledge to facilitate learning, but also induces a proper regularization (homophily) to alleviate overfitting. In fact, our experiments (see Fig. 3) have shown that only learning with the reaction process (no diffusion) performs much worse. Please see R1 for Reviewer 9boz for more detailed discussion.\n\nThird, we agree that the K-partite graph might not be the only choice. There might be other possible graph representations of the observed data. This is an interesting open equation. We will investigate it in the future study.\n\n>C3: \"In Fig. 2(c), why do the red lines decline? In the long-term, the performance of DEMOTE is worse than the baseline model NONFAT. Could the authors articulate it?\"\n\nR3: Actually, Fig. 2c is **NOT** about the prediction accuracy. The decline of the red line does **NOT** mean worse long-term performance. The whole Fig. 2 shows the embedding trajectories learned by different methods (see the caption). Fig. 2a & 2d show the ground-truth, Fig. 2b & 2e show the trajectories learned by NONFAT while Fig. 2c & 2f the learned trajectories by DEMOTE. Comparing Fig. 2c and 2a, we can see that the learned trajectories by DEMOTE is much closer to the ground-truth than those learned by NONFAT (see Fig. 2b). Hence, the performance of DEMOTE actually is better than NONFAT.\n\n>C4: \"For the results in Figure 5(a) and (b), why are the entry value predictions with high and sharp spikes across time?...\"\n\nR4: Great question. Our method does NOT model the entry value as an ODE. We model the embedding trajectories in each mode as an ODE; see Eq1. Given the trajectory values, the entry value at a specific time point is then predicted by a neural network (rather than an ODE); see Eq 2. Therefore, it might not be surprising to see fluctuations at the final outcomes, i.e., entry value prediction, since our model only learns smooth trajectories in the latent embedding space (see Fig. 4)."}, {"review_id": "eBcWcLhQ2n", "reviewer": "Reviewer_wffj", "summary": "The paper presents the framework to generate dynamic embeddings for multi-entity (e.g., traffic type, severity, latitude, and longitude) interactions with timestamps. The temporal high-order interactions are modeled via neural diffusion-reaction processes, and to tack the temporal drift of underlying node, authors propose the dynamic embeddings that follows a diffusion-based ODE model. The diffusion process uses a multi-partite graph to coevolve the trajectories, and the reaction process models the individual evolution of each object. The combination of these two processes can capture both the commonalities and personalities in learning the embedding trajectories for different objects. Authors demonstrate that the proposed model can capture the temporal dynamic well through both synthetic and real-world datasets with comparison to a large number of existing methods.", "questions": "1. [Line 6] It seems that \"I\" in \"EMbedIngs\" should be lowercase.\n\n2. [Lines 92-100] It seems that there are not enough conclusions to this paragraph. What are the advantages and disadvantages of the indicated \"continuous CP decomposition\" approach?\n\n3. [Line 106] \"typical\" -> \"typically\".\n\n4. [Line 270] \"We compared <our approach>...\"", "limitations": "The proposed method is limited to a small number of entities, however, further development of the approach using the graph cut algorithms could potentially improve its scalability.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 3, "strengths": "1. The proposed ODE-based model looks like a reasonable way of representing the dynamics of embeddings.\n\n2. The dynamic embeddings for entities are generated in continuous, rather than discrete time, unlike most existing approaches.\n\n3. Numerical experiments demonstrate that the proposed method outperforms baselines on 3 real-world datasets.", "weaknesses": "1. It seems that the text can be further improved to present better. Some phrases are repeated several times unnecessarily, and some points, on the contrary, are covered too briefly. I would recommend the authors to carry out additional proofreading of the text.\n\n2. The source code is not provided, and this complicates the reproducibility of the results.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your valuable comments and suggestions. Here are our responses. C: comments; R: responses\n\n> C1: \"It seems that the text can be further improved to present better. Some phrases are repeated several times unnecessarily, and some points, on the contrary, are covered too briefly. I would recommend the authors to carry out additional proofreading of the text.\"\n\nR1: Great suggestion and we do agree. We will polish our paper to further improve the readability and clarify.\n\n> C2: \"The source code is not provided\"\n\nR2: We appreciate the reviewer's concern. We have provided the code via an anonymous link (see the global response R1). \n\n> C3: \"[Lines 92-100] It seems that there are not enough conclusions to this paragraph. What are the advantages and disadvantages of the indicated \"continuous CP decomposition\" approach?\"\n\nR3: Great question. The continuous-time CP decomposition integrates the continuous time information into the CP framework, and hence will not lose information like other methods that simply bin the timestamps into coarse steps (see the references in Line 94-95). However, the disadvantage of this approach is that it still estimates a static embedding and cannot capture the evolution of the representations (see the motivation of our work presented in Line 102-109). In addition, CP is a simple multilinear decomposition model, and is not flexible enough to capture nonlinear interactions of the entities. Hence, our work uses a neural network decomposition model instead (see Eq2 and Line 160-166). We will supplement this discussion into the paper."}]

Author rebuttal:
We thank all the reviewers for their careful, insightful and constructive comments! Here we post a few global responses. C: comments; R: responses

>C1: provide the code

R1: Yes, we are happy to provide the code. We have created an anonymous repository that includes the implementation of our model and algorithm and scripts to run the experiments. According to the rule of author response, `"If you were asked by the reviewers to provide code, please send an anonymized link to the AC in a separate comment"`, we have sent the anonymized link to the AC. 

>C2: improve the presentation 

R2: We appreciate the reviewer's suggestion. We will improve the notations and presentation to further enhance the readability and clarity of our paper.

Meta-review: This paper proposes a dynamic tensor decomposition (TD) model, in which the latent factors are functions of time stamps described by neural ODEs. It demonstrates that the proposed model can capture the temporal dynamic well through both synthetic and real-world datasets with comparison to a large number of existing methods. This work proposes a new trajectory learning model for tensors, which could be interesting to the field. The authors' rebuttal have successfully addressed some questions and minor issues.
Decision: Accept (spotlight)