Undirected Probabilistic Model for Tensor Decomposition
Authors: Zerui Tao
Toshihisa Tanaka
Qibin Zhao
Venue: NeurIPS 2023 poster (2023)
Keywords: Tensor decomposition
tensor completion
probabilistic methods
TL;DR: We introduce an innovative approach to undirected probabilistic tensor decomposition, characterized by its exceptional flexibility in accommodating various structures and distributions.

Abstract: Tensor decompositions (TDs) serve as a powerful tool for analyzing multiway data. Traditional TDs incorporate prior knowledge about the data into the model, such as a directed generative process from latent factors to observations. In practice, selecting proper structural or distributional assumptions beforehand is crucial for obtaining a promising TD representation. However, since such prior knowledge is typically unavailable in real-world applications, choosing an appropriate TD model can be challenging. This paper aims to address this issue by introducing a flexible TD framework that discards the structural and distributional assumptions, in order to learn as much information from the data. Specifically, we construct a TD model that captures the joint probability of the data and latent tensor factors through a deep energy-based model (EBM). Neural networks are then employed to parameterize the joint energy function of tensor factors and tensor entries. The flexibility of EBM and neural networks enables the learning of underlying structures and distributions. In addition, by designing the energy function, our model unifies the learning process of different types of tensors, such as static tensors and dynamic tensors with time stamps. The resulting model presents a doubly intractable nature due to the presence of latent tensor factors and the unnormalized probability function. To efficiently train the model, we derive a variational upper bound of the conditional noise-contrastive estimation objective that learns the unnormalized joint probability by distinguishing data from conditional noises. We show advantages of our model on both synthetic and several real-world datasets.

Undirected Probabilistic Model for Tensor Decomposition
Zerui T
Toshihisa Tanaka 1.2
Qibin Zhao
zerui. tetuiittt@riiin tao@riken. tanakat@cc. tuat. jc. jp
qibin. zhao@riken -jp
'Tokyo University of Agriculture and Technology
PRIKEN AIP
Abstract
Tensor decompositions (TDs) as powerful tool for analyzing multiway data. Traditional TDs incorporate prior knowledge about the data into the model, such directed generative process from latent factors observations. In practice, selecting proper structural or distributional assumptions beforehand is crucial for obtaining promising TD representation. However, since such prior knowledge is typically unavailable in real-world applications, choosing appropriate TD model can be challenging. This paper aims address this issue by introducing a flexible TD framework that discards the structural and distributional assumptions, in order to learn as much information from the data. Specifically, we construct TD model that captures the joint probability of the data and latent tensor factors through deep energy-based model (EBM). Neural networks are then employed to parameterize the joint energy function of tensor factors and tensor entries. The flexibility of EBM and neural networks enables the learning of underlying structures and distributions. In addition, by designing the energy function, our model unifies the learning process of different types of tensors, such as static tensors and dynamic tensors with time stamps. The resulting model presents doubly intractable nature due the presence of latent tensor factors and the unnormalized probability function. To efficiently train the model, derive variational upper bound of the conditional coise-contrastive estimation objective that learns the unnormalized joint probability by distinguishing data from conditional noises. We show advantages of our model on both synthetic and several real-world datasets.
1 Introduction
Tensor decompositions (TDs) serve as powerful tools for analyzing high-order and high-dimensional data, aiming to capture the inter-dependencies among different modes by utilizing multiple latent factors. TDs have demonstrated remarkable success various machine learning tasks, including data imputation [53, 9], factor analysis [4], time-series forecasting [27], model compression [28, 41], generative models [10, 20] among others.
Existing TDs typically incorporate predefined directed graphical models into the generative process. These models specify the priors latent factors and the conditional probabilities of observations, following specific contraction rules associated with the latent factors. Traditional contraction rules predominantly employ multi-linear products, like CP [15]. Tucker [42], tensor train [29] and other variants [19, 6). However, selecting an appropriate contraction rule for specific datasets is often chal- lenging in real-world applications. Recent research, known as tensor network structure search (TNSS, 22, 23], has demonstrated that selecting appropriate TN contraction rule significantly enhances the factorization performance, Another promising approach involves learning non-linear mappings from the data, utilizing techniques like nonparametric models [5, 45, 53] and deep neural networks
"Corresooningg author
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

[25.9] 9]. Empirical results demonstrate that non non-linear TDs exhibit superior performance compared to traditional multi-linear TDs in various upplications, attributed to their enhanced expressive power.
Despite the success of non-linear TDs in reducing structural assumptions, they often rely on simplistic distributional assumptions. Typically, specific directed graphical model is adopted to model the generative process from latent factors to tensor entries, represented as p(x) = J p(x =)p(z) dz, where z denotes tensor factors and represents observations. Additionally, the distributions are usually selected from exponential families for tractability, such as Gaussian and Bernoulli distributions. For instance, a Gaussian prior can be assigned to latent factors, and observed entries can be modeled using Gaussian distribution [32, 49] or Gaussian process [45, 53]. However, these prior assumptions regarding the probabilistic model can introduce model bias and reduce the effectiveness of TD models. In real-world applications, the latent factors might originate from unknown distributions, and the observations can exhibit complex multi-modal generative processes. Without knowing the underlying generative process, these simplistic assumptions can lead to inaccurate estimations.
To address these issues, this paper proposes to construct an undirected graphical model of TD. More specifically, a TD model that captures the joint probability of the data and latent tensor factors is con- structed through deep energy-based model (EBM), represented as p(x, exp(-f 2)). Neural networks (NNs) are then employed to parameterize the joint energy function f(x z). The flexibility of EBM and NNs facilitates the learning of underlying structures and distributions. Furthermore, our model unifies the learning process in the presence of side information, such as dynamic tensors with time stamps, by designing the energy function. The resulting model presents a doubly intractable nature due to the presence of latent tensor factors and the unnormalized probability density function (pdf). For efficient model training, we derive a variant of conditional noise-contrastive estimation [3] algorithm that learns the unnormalized joint probability by distinguishing data from conditional noises. The proposed model offers several advantages: (1) features a flexible structure that can adapt to different distributions; (2) the undirected nature allows us to learn more general correlations than traditional directed TDs; (3) it can handle diverse tasks and encode auxiliary information by adjusting the energy function.
Experiments are conducted on synthetic and real-world datasets showcase the advantages of our model. Through simulation studies, we demonstrate the capability of our model to handle data generated from diverse distributions, in contrast to traditional Gaussian-basee models that yield unfaithful and biased estimates. Subsequently, experiments are performed on multiple real-world datasets to evaluate sparse and continuous- time tensor completion. Our model outperforms various baselines across multiple metrics and settings, highlighting the generality of the proposed model.
2 Backgrounds
We adopt similar notations with [19]. Throughout the paper, we use lowercase letters, bold lowercase letters, bold capital letters and calligraphic bold capital letters to represent scalars, vectors, matrices and tensors, e.g., x, x, X and X. Tensors refer to multi-way wayli-wayyaee arrays which generalize matrices. For a D-order tensor X E TR xxxxxxxppp we denote its (i.,...,úi)) ip), entry as Ti.
Notations
2.1 Tensor decomposition
Given a D-order tensor X E B tensor decomposition (TD) aims to factorize X into D smaller latent factors Z d E Rax by using some predefined tensor contraction rules. The classical Tucker decomposition [42] assumes X = Wx1 Z xx. DD where W E R. x...xxx. RD is the coefficient and Xa denotes the matrix-tensor contraction [19]. Equivalently, each entry can be written as x; = Σ ... ER Wr....r = ... >.prrr where the tuple Rp) is the Tucker rank of tensor X. The latent factors zd can capture information of each tensor mode and W represents the weight of each factors. CP decomposition [15] is a restricted form of Tucker by assuming W is super-diagonal, i.e., 3; = D Wr W =P where we simplify w, wy In this paper, we focus on probabilistic version of TDs, which serves as generalizations of traditional ones. The standard approach is to formulate TDs as directed graphical model, x(Xx) = [(x ))p( 2) dZ, where Z denotes [ ...., z for simplicity. For continuous data, the x Z) is usually assumed to be Gaussian and TDs are used to parameterize the mean of corresponding Gaussian distribution [32, 49, 50].

---

Despite the elegant form of these multi-lineas contraction rules, they have limited flexibility that can be mitigated by extending TDs their non non-linear counterparts. We can think of TD as a function that maps the multiway latent factors to tensor entries. One extension is to add Gaussian process (GP) priors on the function to obtain a nonparametric model, which resembles a GP latent variable model. In particular, [53] proposed to stack the latent factors as m, [z], ........ ...
=| E RDD and then assign GP prior the functional mapping. In specific, for continuous data, it assumes I; N M(4, o²), where the mean function is a GP fا; = fm)) ~ (O k(m,))) associated with kernel k(:;)) :). Since GP has large computational complexity and designing the kemel requires ad ad-hoc expert domain knowledge, [25] proposed to parameterize the function using neural networks (NNs), Ti ~ N NN44, o²) where 14 = JNN (m:). However, NNs easily overfit due the high-dimensional and sparse nature of tensor data. To address this issue, [9] proposed to use Bayesian NN with spike-and-slab prior for sparse weights. All these models are based on directed graphical model, that assumes there exists direct mapping from the latent factors to tensor entries and use simplistic distributions.
2.2 Energy-based model
Energy-basee based model [EBM, 21] is class of undirected probabilistic models, which uses an energy function to characterize the data distribution. Given observed data z, the basic idea of EBM is to approximate the data distribution by Boltzmann distribution, (x) exp(-f(z;0)) where
2((0)
Z(0) exp(-f(z; 0)) dr is the normalization constant (a.k..,, the partition function) and f(; 0) is the energy function. One classical example of EBM is the restricted Boltzmann machine (RBM) where the energy function has bi-linear form for tractability. In deep EBMs, the energy function is typically parameterized by deep neural networks. The main difficulty for training EBMs is to deal with the intractable normalization constant [36]. There several ways to train EBMs, including contrastive divergence [14]. score matching [16]. noise-contrastive estimation INCE, 11] and so on. In this paper, we focus on NCE, due to its efficiency and ability of tackling different data types.
We denote the unnormalized pdf as 60x; 0) = exp(-J(x; 0)). NCE consider the normalization constant 7()) as trainable parameter. However, maximum likelihood estimation (MLF) does not work for this case since (0) can be arbitrarily small and the log- likelihood goes to infinity. Instead, the NCE can be obtained by maximizing the following objective,
(0) = F.g log h(r;; 0) + viy log(1 hyy; 0)),
where x denotes observed data and y denotes noises generated from some known distribution Pn, v is the ratio between noise and sample sizes, i.e., #y/ #y/#r. #x. And h(-) can be regarded as classifier )(u:0)
It has been shown that that distinguish data from noises, defined as follows, h(u;; 0) (vp:)))pp.."
NCE is consistent with MLE [11].
Although the noise distribution is essential for training efficiency, the selection of noises is currently limited to heuristics. A common intuitive is that the noise should be similar with the data. Indeed, [3] proposed to conditional noises y Py ar) and minimizing the following loss function,
LCNCE(f) ()) = 2Ezy log[1 + exp(-G(x, y))l.
where G(u1, W2; 0) = log ((u1:))pe( (e wa) with y drawn from Po(y x))
3 Proposed model
3.1 Energy-based tensor decomposition
Even though many non-linear tensor decompositions (TD) have been proposed to enhance flexibility, existing methods typically adopt simplistic distributions such Gaussian. This can be problematic for complex real-world data. To address the issue, we propose energy-based tensor decomposition SinergyTD), by integrating EBMs in TD framework. Given order-D - tensor X' of shape I x
TD we aim to factorize it into D smaller latent factors zd E TR. l Vd 1,..., D. We denote the latent factor associated with the (f1,..., ip ) entry as m; = where 211 Rd is the ih row of z Unlike traditional directed TDs trying to parameterize the conditional expectation E[[z] m:] f(m)), we model the joint distribution using EBM,
=P E DR
p(a;, m;; 0) exp( (f m4; 0))
(2)
2()
3

---

where - ((:-;)) ) is the energy function and Z f eexp(_f(2) m.: f)) dz; dm; is the partition function to make it a valid pdf. We further assume the joint probability of all entries are independent. Le., p(X, m) lien P(I;, 2/, ...,.,,, = ). where f denotes the set of observed entries. This is standard setting TDs and the dependence of tensor entries can he captured by sharing latent factors.
The expressive nature of the energy function enables us to easily handle diverse data types. For example, we can deal with discrete data by plugging one-hot codings into Eq. (2) to represent categorical probabilities. Additionally, the flexibility of NNs allows us to model tensors with side information, where each tensor entry incorporates additional features [34]. Specifically, in this paper, we focus on a particular case of dynamic tensors with continuous time stamps [48]. In this case, we consider an observed tensor as time series Xt where the time stamp continuous, and each entry has its own specific time stamp 4. To model the tensor time series, we assume that each entry follows the same distribution and construct the time-dependent energy function, p(Zi, m4; 0, (4) exp(- (2;, mi, 4;; f)), where the time stamp 4 is considered as an auxiliary feature. The flexibility of NNs allows this function to learn general patterns across continuous time stamps. Experimental results demonstrate that this simple treatment can achieve good performances.
Network architecture The network architecture plays a crucial role in learning accurate probabilis- tic manifolds. Specifically, we define the energy function as [[1 m,) = 91 (92 (2)()1), 91 (m;))) where J3 and 91 are MLP layers that encode information from zj and m., respectively. 92 is a summation concatenation layer that induce coupling between tensor values and latent factors, and 91 is the output layer. Although we currently utilize only MLPs, it is worth noting that con- volutional architectures, demonstrated by [39, 25], can also be employed, which is topic for future research. To handle dynamic tensors, we incorporate an extra sinusoidal positional encod- ing layer [37] denoted as (t) to capture temporal information. This embedding utilizes random Fourier features as proposed by [31]. Consequently, the energy function can be expressed as [ai mi 44 = 9i (#2(g3 (Ii), gs (m.) (44)). This architecture is commonly employed to capture temporal information and has been demonstrated to effectively learn high-frequency information when combined with MLPs [37].
Posterior sampling A significant application of TDs is to estimate the posterior of missing entries. Unlike traditional TDs, direct predictions cannot be obtained even after learning the latent factors due to the utilization of an undirected probabilistic model. Instead, we need to seek for sampling methods of p(z; m4). One choice is score-based samplers by utilizing the score function vz, log p(æ) p(z).m))
mi) = m log p(mm) - v., ( mi). such as Langevin dynamics [44]. Score-based samplers are not suitable for handling discrete data. However, in our case, we model the one-dimensional pdf for each entry, enabling us to directly sample the discrete data. Consequently, for continuous data,
the use of grid search is viable approach to obtain maximum a posteriori (MAP) estimations.
3.2 Learning objective
Despite the flexibility of the proposed model in Eq. (2), obtaining maximum likelihood estimation (MLE) becomes doubly intractable, as both the partition function (0) and the marginal distribution p(z;) are intractable. Therefore, the CNCE loss Eq. (1) cannot be directly applied. In this section, we extend the variational approach [33] to construct upper bound that addresses the challenge posed by intractable marginal distributions.
Denote the unnormalized pdf as d(1,, m;; 0) = exp(=f(xi, = m;; 0)) and the unnormalized marginal pdf as o( 0) I mx;; 0) dvn;, For clarity, we omit the index in the subsequent context of this subsection. We follow the idea of CNCE [3] to distinguish data x from conditional noises y ~ (y 2)) Firstly, Eq. (1) can be rewritten as
LCNCE 2.Ezy log(1 + 1/r(2, 0)}
(3)
where
r(x, 0) = أ<x; 0) f) x)
(4)
:(y: () y)
However, for our problem, the unnormalized marginal probability o(x;;) 0) is unknown. An additional variational distribution q(m;;) up) is used to approximate the true posterior p(m X;; 2). Note that
in TDs where the data size is static, there is no need for amortized inference, which is different

---

from previous ones like [2]. Equipped with the variational distribution, the unnormalized marginal distribution can be computed using importance sampling.
6(x, m. ()q(m: dm m((p)))))
[( m, 01
6(x; 0)
(5)
p(m;
g(m sp)
Plugging Eq. (5) into Eq. (4), we have
r(z,y: y: E By(ms (x, m; ()/((m////////mmmmmmmm)))))))// (>)p(()
(6)
oly: px [x] y)
Since Eq. (3) is convex function w.r.t. r(z, yi 0). plugging Eq. (6) into Eq. (3) and applying the Jensen's inequality, we have the upper bound,
(CNCE(0) = my log[1 + 1/r(z, (y, 0)]
I (((p)) log 1+ O)pe<> () LVCNCE(O.))
(7)
Φ(z, m; p) x)
Following [33], we have the theorem about the tightness of the bound.
Theorem The difference between the VCNCE loss Eq. (7) and CNCE loss Eq. (1) is the expectation of the f-divergence,
LVONCE :p) LCNCE(O) = Ezy Da,, (p(m x;; [;)]]]((m;)))
where (uy(u) = llog(
Maotu Nzv
with Kizy
p(x:0)p. ((.:))).(()))
The proof can be found in appendix. Based on the theorem, we have the following corollaries to justify the optimization process.
Corollary When q(m; :p) equals to the true posterior, the CVNCE bound is tight, i.e.,
LVCNCE = LeNcE > q(m;; up) p(m x;; 0).
Corollary The following two optimization problems are equivalent,
min (C) = min min EvCNCe(0,9)) ((mpp)))
In practice, we need seek for the sampled version of Eq. (7). Supposing we have N observed samples [ noises {yn,) for each sample 2; and using importance sampling for duy 0), the sampled objective function is,
E
Oon m2);88 g(mp;p)
Pe(ai V1.a)a((m vp)
(VCNCE(0,9))
E Em(mp))) log 1
(1 my; 0 Op (4))
Specifically, we formulate g(m; vp) as diagonal Gaussian and use reparameterization trick [18] to compute the expectation. When dealing with continuous data, we typically select conditional Gaussian noises, represented (1 z) N2 x, o²). This choice entails only one hyper- parameter a that needs to be tuned. Another benefit is the symmetry of the conditional distribution for Gaussian noise, expressed (y æ) (x(x y). Hence, the objective function can be further reduced. For binary or categorical data, such symmetric noises can also be derived [3].
The time complexity of the proposed objective is OlyB(DRH + LH?)), where B is the batch size, is the number of conditional noises, H/ the number of hidden units per layer, layers and D is the tensor order. The time complexity of our model is times greater than traditional TDs, since we need to compute forward passes for u particles. However, as we only use small networks, the computational speed is still very fast (See Appendix C.3 for an illustration).
number of
5

---

4 Related work
Traditional tensor decompositions (TDs) are based on multi-linear contraction rules, such as CP [15], Tucker [42], tensor networks [29, 51] and their variations [19, 6]. In this paper, we mainly focus on probabilistic TDs, which extend traditional methods by providing uncertainty estimates about both observations and latent factors [32, 49, 50, 26, 38]. These models build directed mapping from latent factors to tensor entries using multi-linear contraction rules, resulting in limited flexibility when dealing with complex datasets. An alternative approach involves replacing the multi-linen relations with non-lineat ones. [5, 45, 52] introduced the use of tensor- Gaussian processes (GPs) for achieving nonparametric factorization. [53] further expanded this concept by incorporating GP prior on the function that maps latent factors to tensor entries, resulting in a nonparametric TD for sparse tensors. GP-based TDs are further extended using hierarchical priors [40], stochastic processes [43, 8J. Despite the success of GP-based TDs, nonparametric approaches can encounter computational challenges and may still lack sufficient flexibility. Recently, neural networks (NNs) are also applied to TDs. [25] suggested the utilization of convolutional NNs to map latent factors to tensor entries. Besides, [7] built hierarchical version of the Tucker model and introduced non-linear mappings within each hierarchy of latent factors. To mitigate overfitting, [39] suggested the adoption of deep kernels in GP-based TD rather than using NNs directly. On the other hand, [9] proposed to use Bayesian NN with spike-and-slad prior to prevent from overfitting and obtain probabilistic estimates. More recently, [24] adopt neural ODEs to capture dynamic tensor trajectories. Other works regarding more flexible exponential families [13] or mixture of Gaussians [12] employ linear structures. While all these methods using directed mapping from latent factors to tensor entries, our model is fundamentally different from them, in that we construct much more flexible undirected probabilistic model of TD that can deal with diverse distributions and structures.
Another related direction is the energy-based model (EBM). To address the intractable pdf of EBMs, various training methods have been proposed, including contrastive divergence [CD, 14]. score matching [SM, 16]. noise-contrastive estimation INCE. 11]. CD requires large steps of Monte Carlo samples, which can be computationally expensive for high-dimensional tensors. SM cannot handle discrete data, and learning latent variables with SM requires complex bi-level optimization [2]. Therefore, we focus on NCE in this paper. Learning energy- -baed TD even more challenging because it involves multiple coupled latent factors that cannot be analytically marginalized. [33] proposed VNCE to handle unnormalized models with latent factors. We enhance their algorithm by using conditional noises [3] to improve the learning efficiency. One fundamental distinction between our model and traditional EBMs is that, through TD construction, we only need to learn one-dimensional distributions for each scalar entry, instead of the original high-dimensional tensor. Hence, our model avoids performance degradation when learning high-ding -igh-dimensional data using NCE.
5 Experiments
We demonstrate the proposed energy-based TD (Energy TD) on synthetic data and several real- world applications. All the experiments are conducted on a Linux workstation with Intel Xeon Silver 4316 CPU, 256GB RAM and NVIDIA RTX A5000 GPUs (24GB memory each). The code is implemented based on PyTorch 1.12.1 [30]. More experimental details can be found in the appendix. The code is available at https: // //g/tttttubu/ com/ comn/toozerui/energy
5.1 Simulation study
Tensors with non-Gaussian distributions Traditional TDs commonly assume that tensor entries follow a Gaussian distribution. However, in real-world applications, we often encounter highly complex distributions. In this experiment, we evaluate the capability of our model to learn distributions that deviate from the Gaussian assumption.
We consider two-mode tensor of shape 7 x 7, where we set / = 8. Firstly, two latent factors of shape / x R are generated, where the rank R is set to 5. Then, conditioned on the latent factors, we generated tensor observations from particular distributions. For each entry, we generate N 200 samples. Three types of distributions are considered: (1) Beta distribution; (2) Mixture of Gaussians (MoG) and (3) Exponential distribution. For Beta distribution, we generate latent factors from
uniform distribution i=12 Uni(0.0, 1.1). Then, sample the observed tensor from Beta

---

2.0
Tuth
2.0
1.5
Density
1.5
Density,
Density
Enegy TD
1.0
1.0
0.5
0.5
0.0
0.0
0.0
0.0
0.5
1.0
1.5
-1.0
0.0
1.0
2.0
-1.0
0.0
1.0
2.0
3.0
(a) Beta distribution
(b) Mixture of Gaussians
(c) Exponential distribution
Figure 1: Simulation results for different distributions. The blue line is the ground truth pdf. The yellow line is the kernel density estimation (KDE) plot of observed samples. The red line is the GPTF estimation, which is Gaussian pdf. The green line is method, computed by evaluating the unnormalized pdf on grids and calculating the partition function using Gaussian quadrature.
2.0
2.0
Test
1.5
1.5
1.5
Truhh
NUNFAT
EnergyTD
0.5
0.5
0.0
0.0
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.5
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Time
Time
Time
(a) MR 0.3
(b) MR 0.1
(c) Fully observed
Figure 2: Simulation results for continuous-time tensor decomposition. The blue regions are observed and the red regions missing. The trajectories of ground truth, BCTT, NONFAT and our model drawn in black, blue, red and green lines, respectively.
distribution Tij Beta I 22 1.2). For MoG distribution, draw latent factors from uniform distribution 2 iid Uni(0.00 1.0). The tensor entries are then drawn from MoG T.y iid 0.6 N/ (cooss( z (zy))) + 0.4 (N(sin((2 Z² 0.25?). For Exponential distribution, generate latent factors from uniform distribution zi=1, lid Uni(O.O, 1.0). The tensor entries iid
then sampled from Exponential distribution Vij ~ factorization [GPFT, 53], which assumes the entries follow Gaussian distribution.
Exp(( 2 (TT)))) We compare with GP tensor
The results of probability density function (pdf) estimation are presented in Fig. 1. We display the learned pdf of one single tensor entry. This reveals that GPTF is limited capturing only the 1st moments (mean values) and overlooks higher-order information of more complex distributions. Our model exhibits greater flexibility handling non- oonGaassiia distributions.
Continuous-ting tensors We then consider dynamic tensor, where each entry is time series. We follow similar setting with [9] and use the same data size with the previous simulation, i.e., two-mode tensor of shape 8 x 8 with each entry being time series of length 200. We firstly generate latent factors of shape 2, with each rows drawn from 2 ~ (O0 2- 2-T) and 2? ~ ([1, 1], 2- T). Then generate N 200 observed entries from time E [0. 1]. The tensor entries computed by (1 - Σr,=1 =,,,, 22 Prir (l), where wn (1) sin(2 (12 cos(2m/(), w21 (4) 2 (2rt) and wez (t) = cos(5xt)) sin (5xt). Finally, the data are normalized to range [0, [0,2]. The synthetic data consist of low-free low-frequency trends and high-frrquency fluctuations. Apart from fully observed case, we also test with missing rates (MR) 10% and 30%. In specific, for each entry, we randomly select starting time and set the following consecutive 10% 30% time stamps missing.
compare with two methods that are designed for dynamic tensors, the Bayesian continuous-time Tucker decomposition [BCTT, 8] and the nonparametrie factor trajectory learning [NONFAT, 43].
7

---

Table 1: Sparse tensor completion
RMSF.
MAE
Alog
Rank 3
Rank 5
Rank
Rank 10
Rank 3
Rank
Rank
Rank 10
CP.WOPT 1.486 1: 0.282 1,386 0,043 1,228 1: 0.063 1.355 1: 0.079 006994 0.098 0.664 0.018 0.610 1: 0.027 0.658 1: 0,026
GPTF
0,911 + 0.008 0.867 + 0.008 0.878 + 0.009 0.884
+ 0.009
0.511 = 0.005 0.494 = 0.004 0.530 + 0.004 0.554
+ 0.006
HGP-GPTF
0.896
0.009 0,850 + 0.011 0.844
0.006
0.479
0.473
0.474 + 0.004 0.480 0.004
POND
0.885 + 0.010 0,871 1: 0.013 0.838 1: 0.009 0,857 1 0.011
0.463
+ 0.004 0.454 : 0.005 0.444 1 0.005 0.443 + 0.006
CoSTCo
0.999 0.007
0.017 0,930 1: 0.024 0.909
(0.523
0.007 0.514 1 0.031 0.481 0.008
EnergyDD 0.864 + 0.011 0,835 + 0.011 0.840 + 0.013 0,833
+ 0.016
0.450
0.006 0,433
= 0.006 0.424 + 0.005 0.409
+ 0.004
ACC
CP-WOPT 0.533 + 0.039 0.592 0.037 0.603 +: 0.028 0,589 +: 0.022
00.138 0.004 0.147 + 0.005
0.148 +: 0.003 0.147 + 0.004
GPTF
0.367 1: 0.001 0.357 1: 0.001 0.359 1: 0.001 0.368 1 0.001
001152
0.002 0.150 1: 0.001 0.167 1 0.002 0.182 1 0.001
HGP-GPTE 0.355 + 0.001 0.344 +: 0.001 0.341 + 0.001 0.338 + 0.001
0.125 + 0.003
0.129 1 0.001
0.139
+ 0.000 0. 145 + 0.002
CoSTCo 0.385 0.003 0.376 0.018 0,363
+ 0.004
0.348
+ 0.002
0.111 0.004 0.137 0.020
0.107
- 0.004
0.101 0.004
EnergyID 0,348 + 0.005 0.336 + 0.004 0.328 + 0.003
0.328
+ 0.003
0.110
= 0.008
0.101
= 0.006
0,094
+ 0.006
0.101 = 0.009
BCTT treats Tucker core tensors as functions of time and NONFAT treats all GPTF factors as time series. Unlike BCTT with GP prior on the time domain, NONFAT uses GP prior on the frequency domain through inverse Fourier transform of original time series.
Fig. displays the completion results. The learned trajectory of single tensor entry is plotted. Higher missing rates result in the inability of BCTT to capture accurate trajectories, particularly in missing regions. NONEAT achieves more stable predictions, yet it tends to favor over-smoothed trajectories while disregarding high-frequency fluctuations. This behavior may be attributed to its unique construction, which introduces GP prior in the frequency domain. Our utilization of flexible neural networks allows us to adapt to complex situations encompassing both low-frequency and figh-free information.
5.2 Tensor completion
We evaluate our model on two sparse tensor and two dynamic tensor completion applications. For real datasets, the energy function be difficult learn, when the pdf has very sharp curve. Motivated by the idea of noise-perturbed score estimation [35], we add small i.i.d. Gaussian noises on the data during the training of Energy TD form of smoothing technique. The results reported on clean test data. For EnergyTD, we use MAP estimates as described in Section 3.1.
5.2.1 Sparse tensor completion
We test our model on two sparsely observed tensors: (1) Alog, file access log dataset [52] of shape 200 users x 100 actions x 200 resources with about 0.33% nonzero entries; (2) ACC, a three-way tensor generated from code repository management system [52] of shape 3k users 150 actions x 30k resources with about 0.009% nonzero entries. We use the same dataset split as in [52] and report the 5-fold cross validation results.
Competing methods We compare with five baselines: (1) CP-WOPT [1], CP decomposition with stochastic optimization; (2) GPTF [53], a GP-based tensor factorization using stochastic variational inference; (3) HGP-GPTF [40], a GPTF equipped with hierarchical Gamma process prior; (4) POND [39], a probabilistic non-linear TD using deep kernels with convolutional NNs (CNNs); (5) CoSTCo [25]. a non-linear TD that uses CNNs to map latent factors to tensor entries. CP-WOPT is provided
in Matlab Tensor Toolbox [1]. We implement GPTF based on PyTorch by ourselves and use official implementations for HGP-GPTF2, POND and CoSTCo".
Experimental settings and results We set batch size 1000 and run 1000 epochs for Alog, 100 epochs for ACC. For our model, we use Adam [17] optimizer. Learning rates of all models are chosen from (1e-2, 2, 1e-3, 1e 1e-4]. For all methods, we evaluate with rank R E {3, 5, 8, 10). All methods are evaluated by 5 runs with different random seeds.
http://///gghabcccccacllll////// gegpthab..cm/cttillin///nn "nttps://githu" hatpa//o///---M-aaaa///1111----0
8

---

Table 2: Continuous-ting time tensor completion
RMSF
MAE
Air
Rank 3
Rank 5
Rank 8
Runk 10
Runk 3
Rank 5
Runk
Rank 10
CTCP
1.020 1: 0,002 1,022 1: 0.002 1,022 1: 0.002 1,022 1: 0.002
00.784 4: 0.002 0.785 0.002 0.787 0.002 0.787 1: 0.002
CTGP
0.475 + 0.000 0.463 + 0.000 0.459 + 0.000 0.458 + 0.000
0.318 = 0.000 0.304 0.000 0.301 = 0.000 0.299
= 0.000
CTNN
1.013
0.001 1,005
0.005
1,013 + 0.002
0.780
0.777
0.776 0.003 0.780 0.001
NNDTN
0.377
0.00 0.364 + 0.002 0.334 1: 0.004 0.328 1: 0.004
0.244 = 0.003 0.239 1 0.002 0.217 = 0.003 0.212 = 0.004
NONEAT
0.339 0.002 0.335
1: 0.002
(0.224
0.219
0.221
0.001
THIS-ODE
0.569 + 0.001 0.566
+: 0.004 0.542
+ 0.005
0.541
+ 0.002
0,415 0.002 0.409 0.004
0.395 0.004
0.391
= 0.001
EnergyTD 0.302 + 0.008
0.291 0.006
0.300
0.012
0.283
+ 0.004
0.184 0.006
0.177
0.003
0.172 0.006 0.184 0.003
Click
CTCP
2.063 0,009
2,020 1: 0.025
2,068
0.012
2.009 1 0.023
1.000
0.009
0.977 0.021
1.005 0.010
0:969 1 0.012
CTGP
0.002
1.423
0.004
1.404
0.004
1.392
1: 0.002
0.880 = 0.003 0.877 0.003 0.856
0.002 0.849 + 0.001
CTNN
1,820
0.005
1,820
+
1.820
0.005
1,820
+ 0.005
1.077
1.053 0.027
1.083
0.016
1,071 0.024
NNDTN
1.488
+ 0.005
1.409
+ 0.004
1.407
+ 0.002
1.410
+ 0.004
0.858
= 0.002 0.856 = 0.002 0.859 = 0.003 0.863 = 0.006
NONEAT
1,400 0.008
1.351
1: 0.002
00855
0.812 0.002
THIS-ODE
1.421 1: 0,004
1.413
4: 0.002
1.408
1: 0.002
1,395
1: 0.003
0.836
0.004 0.836 0.003
0.832 0.002 0.829 1 0.003
EnergyTD
1.39
+ 0.003
1,385
0.003
1.356
0.001
1.357
+ 0.001
0.777 0.003 0.775 0.003
0.772 0.002 0.773 - 0.001
The completion results are shown in Table 1. The mean and standard deviation of the root mean square eπror (RMSE) and mean absolute eπror (MAE) are reported. Results from POND are not included for ACC due to the slow code execution and the substantial memory requirements for this dataset. Our model outperforms the baselines in nearly all cases on both RMSE and MAE. Moreover, the improvement is significant (p < 0.05) for most cases. All non-linear methods exhibit substantial superiority over CP-WOPT, thereby highlighting the advantage of employing more flexible structures. Furthermore, the enhanced performance of our model may be attributed to the utilization more flexible undirected probabilistic distributions. It is worth noting that our model, despite employing MLP layers, outperforms POND and CoSTCo, which utilize convolutional layers. We believe that model can be further improved by designing more suitable network architectures.
5.2.2 Continuous-time tensor completion
In this subsection, we evaluate our model on two continuous-time tensor datasets: (1) Air, the Beijing air quality dataset [47] of shape 12 sites x pollutants with about 1 x 10 observations different time stamps; (2) Click, an ads click dataset [43] of shape 7 banner positions x 2842 site domains x 4127 mobile APPs with about 5 x 10 entries different time stamps. We use the same dataset split as in [43] and report the 5-fold cross validation results.
Competing methods We compare with (1) Nonparametric factor trajectory learning INONFAT, 43], which is the SOTA method for continuous-time tensor completion, and other baselines including: (2) Continuous-time CP [CTCP, 48], which uses polynomial splines to model dynamics of CP coefficients; (3) Continuous-time GP (CTGP), which extends GPTF [53] by injecting time stamps into GP kernels; (5) Continuous- time NN decomposition (CTNN), which directly uses time stamps as inputs in CoSTCo [25] to learn continuous dynamics; (5) Discrete-tima NN decomposition with non-linear dynamics INNDTN, 43], which employs RNN dynamics for time steps; (6) Tensor high- order interaction learning via ODEs [THIS-ODE ODE, 24], where continuous-time trajectories of tensor entries are captured by neural ODEs. For baseline models (1-5), we use the official implementation" provided by [43]. For (6) THIS-ODE, its official implementation' is used.
Experimental settings and results We set batch size 128 and run 400 epochs for Air, 200 epochs for Click. Notably, in the case of THIS-ODE, we only run 200 epochs for Air and 75 epochs for Click due to the slow speed of ODE solvers (See Appendix C.3 for an illustration). For our model, we use Adam [17] optimizer with learning rate chosen from f1c-2, 3). The network architecture is described in Section 1. For baseline models, provided settings in their codebases are used. For all methods, we evaluate with rank R E {3, 5, 8, 10). All methods are evaluated by 5 runs with different random seeds.
hhetpe://gichabb.com//ee

---

Table 2 presents the completion results. It should be noted that the results presented here differ from those reported in [43] as we adhere the standard definition of RMSE and MAE (see appendix for detail). Our model surpasses the baseline methods in almost all cases for both RMSE and MAE, with statistically significant improvements (p < 0.05) observed in most cases. Particularly, we observe that the improvements in MAE are notably more significant. One possible reason is that NONFAT is trained implicitly by minimizing the square loss (with regularization) as it adopts Gaussian assumption about the data. However, our model does not make such assumptions about the data distribution and the loss function. Hence, can adapt to the data distribution more flexibly, Additionally, we observe that directly injecting time stamps into neural networks, done in CTNN, is ineffective, thus highlighting the advantage of our model in learning more informative structures.
Conclusion
We introduce an innovative approach undirected probabilistic tensor decomposition (TD), character- ized by its exceptional flexibility in accommodating various structures and distributions. Specifically, our model integrates deep EBMs in TD relax both structural and distributional assumptions, enabling it to handle complex real-world applications. To efficiently learn the doubly intractable pdf, derive VCNCE objective that is the upper bound of the CNCE loss. Experimental results demonstrate that our model can handle diverse distributions and outperforms baseline methods in multiple real-world applications. One limitation is that our final loss function not fully variational upper bound of CNCE, since we have to use importance samples to approximate the pdf of noise samples in Eq. (7). In the future, we aim to derive a fully variational bound as in [46]. Finally, we did not delve into the interpretability of learned factors in this work. However, exploring the interpretability of these factors represents a promising avenue for future research in the realm of tensor decompositions.
Acknowledgments
Zerui Tao was supported by the RIKEN Junior Research Associate Program. This work was supported by the JSPS KAKENHI Grant Numbers JP20H04249, JP23H03419.
References
[1] Brett W Bader and Tamara G Kolda. Efficient matlab computations with sparse and factored tensors. SIAM Journal on Scientific Computing, 30(1):205-231, 2008
[2] Fan Bao, Chongxuan Li, Kun Xu, Hang Su, Jun Zhu, and Bo Zhang. Bi-level score matching for learning energy-based latent variable models. Advances in Neural Information Processing Systems, 33:181 18121122, 2020.
[3] Ciwan Ceylan and Michael U Gutmann. Conditional noise-contrastive estimation of unnor- malised models. In International Conference on Machine Learning, pages 726-734. PMLR, 2018.
[4] Rong Chen, Dan Yang, and Cun-Hui Zhang. Factor models for high-dimensional tensor time series. Journal of the American Statistical Association, 117553))94-116, 2022.
[5] Wei Chu and Zoubin Ghahramani. Probabilistic models for incomplete multi-dimensional arrays. Artificial Intelligence and Statistics, pages 99. PMLR, 2009.
[6] Andrzej Cichocki, Namgil Lee, Ivan Oseledets, Anh-Huy Phan, Qibin Zhao, Danilo Mandic, et al. Tensor networks for dimensionality reduction and large-scale optimization: Part low-rank tensor decompositions. Foundations and Trends& in Machine Learning, 94455:2449-429, 2016.
[7] Jicong Fan. Multi-mode deep matrix and tensor factorization. In international conference on learning representations, 2021.
[8] Shikai Fang, Akil Narayan, Robert Kirby, and Shandian Zhe. Bayesian continuous-time tucker decomposition. In International Conference on Machine Learning, pages 6235-6245. PMLR, 2022.
10

---

[9] Shikai Fang, Zheng Wang, Zhimeng Pan. Ji Liu, and Shandian Zhe. Streaming bayesian deep tensor factorization. In International Conference on Machine Learning, pages 3133-3142. PMLR, 2021.
[10] Ivan Glasser, Ryan Sweke, Nicola Pancotti, Jens Eisert, and Ignacio Cirac. Expressive power of tensor-network factorizations for probabilistic modeling. Advances in neural information processing systems, 32, 2019.
[11] Michael Gutmann and Aapo Hyvärinen. Noise-contrastive estimation: A new estimation principle for unnormalized statistical models. In Proceedings of the thirleenth international con- ference on artificial intelligence and statistics, pages 297-304. JMLR Workshop and Conference Proceedings, 2010.
[12] Zhi Han, Yao Wang, Qian Zhao, Deyu Meng, Lin Lin, Yandong Tang, et al. A generalized model for robust tensor factorization with noise modeling by mixture of gaussians. JEEE transactions on neural nelworks and learning systems, (9(11):3380-333, 2018.
[13] Kohei Hayashi, Takashi Takenouchi, Tomohiro Shibata, Yuki Kamiya, Daishi Kato, Kazuo Kumieda, Keiji Yamada, and Kazushi Ikeda. Exponential family tensor factorization for missing- values prediction and anomaly detection. In 2010 IEEE International Conference Data Mining, pages 216-225. IEEE, 2010.
[14] Geoffrey E Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 14(8):1771-180, 2002.
[15] Frank L Hitchcock. The expression of tensor or a polyadic as a sum of products. Journal of Mathematics and Physics, 6(1-4):164-189 1927.
[16] Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005.
[17] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint 1rXiv:1412.6980, 2014.
[18] Diederik P Kingma and Max Welling, Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114 2013.
[19] Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3):555-500 5000 2009.
[20] Maxim Kuznetsov, Daniil Polykovskiy, Dmitry Vetrov, and Alex Zhebrak. A prior of a googol gaussians: tensor ring induced prior for generative models. Advances in Neural Information Processing Systems, 32, 2019.
[21] Yann LeCun, Sumit Chopra, Raia Hadsell, M Ranzato, and Fujie Huang. A tutorial on energy- based learning. Predicting structured data, 10)), 2006.
[22] Chao Li and Zhun Sun. Evolutionary topology search for tensor network decomposition. In International Conference on Machine Learning, pages 5947-5957. PMLR, 2020.
[23] Chao Li, Junhua Zeng, Zerui Tao, and Qibin Zhao. Permutation search of tensor network structures via local sampling. In International Conference on Machine Learning, pages 13106- 13124. PMLR, 2022.
[24] Shibo Li, Robert Kirby, and Shandian Zhe. Decomposing temporal high-order interactions via latent ODEs. In Proceedings of the 39th International Conference on Machine Learning, volume 162 of Proceedings of Machine Learning Research, pages 12797-12812. PMLR, 2022.
[25] Hanpeng Liu, Yaguang Li, Michael Tsang, and Yan Liu. Costco: A neural tensor completion model for sparse tensors. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 224-334, 2019.
[26] Zhen Long, Ce 7hu, Jiani Liu, and Yipeng Liu. Bayesian low rank tensor ring for image recovery. IEEE Transactions on Image Processing. 50:3568-3588 2021.

---

[27] Jacob Miller, Guillaume Rabusseau, and John Terilla. Tensor networks for probabilisticc sequence modeling. In International Conference on Artificial Intelligence and Statistics, pages 3079-3087.
PMLR, 2021.
[28] Alexander Novikov, Dmitrii Podoprikhin, Anton Osokin, and Dmitry P Vetrov. Tensorizing neural networks. Advances in neural information processing systems, 28, 2015.
[29] Ivan V Oseledets. Tensor-train decomposition. SIAM Journal on Scientific Computing,
33(5):2295-2311 2011.
[30] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang. Zachary Vitio,, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: An imperative style, high-performance deep learning library. In Advances in Neural Information Processing Systems 32, pages 8024-8035. Curran Associates, Inc., 2019.
[31] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007.
[32] Piyush Rai, Yingjian Wang. Shengbo Guo, Gary Chen, David Dunson, and Lawrence Carin. Scalable bayesian low-rank decomposition of incomplete multiway tensors. In International Conference on Machine Learning, pages 1800-1808. PMLR, 2014.
[33] Benjamin Rhodes and Michael Gutmann. Variational noise-contrastive estimation. In The 22nd International Conference Artificial Intelligence and Statistics, pages 2741-2750. PMLR, 2019.
[34] Qingquan Song, Hancheng Ge, James Caverlee, and Xia Hu. Tensor completion algorithms in big data analytics. ACM Transactions Knowledge Discovery from Data (TKDD), 13(1):1-48, 2019.
[35] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances neural information processing systems, 32, 2019.
[36] Yang Song and Diederik Kingma. How to train your energy-based models. arXiv preprint arXiv:2101.03288, 2021.
[37] Matthew Tancik, Pratul Srinivasan, Ben Mildenhall, Sara Kriil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan Barron, and Ren Ng. Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Processing Systems, 33:7537-7547, 2020.
[38] Zerui Tao, Xuyang Zhao, Toshihisa Tanaka, and Qibin Zhao. Bayesian latent factor model for higher-order data. In Asian Conference on Machine Learning, pages 1285-1300. PMLR, 2021.
[39] Conor Tillinghast, Shikai Fang, Kai Zhang, and Shandian Zhe. Probabilistic neural-kernel tensor decomposition. In 2020 IEEE International Conference on Data Mining (ICDM), pages 531-540. IEEE, 2020.
[40] Conor Tillinghast, Zheng Wang, and Shandian Zhe. Nonparametric sparse tensor factorization with hierarchical gamma processes. In International Conference on Machine Learning, pages 21432-21448. PMLR, 2022.
[41] Andros Tjandra, Sakriani Sakti, and Satoshi Nakamura. Compressing recurrent neural network with tensor train. In 2017 International Joint Conference on Neural Networks (IJCNN), pages 4451-4458. IEEE, 2017.
[42] Ledyard R Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika,
31 31(3):279-311, 1966.
[43] Zheng Wang and Shandian Zhe. Nonparametric factor trajectory learning for dynamic tensor decomposition. In International Conference on Machine Learning, pages 23459-23469. PMIR, 2022.
12

---

[44] Max Welling and Yee W Teh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of the 28th international conference on machine learning (ICML-11), pages 681-688, 2011.
[45] Zenglin Xu, Feng Yan, and Yuan Qi. Infinite tucker decomposition: nonparametric bayesian models for multiway data analysis. In Proceedings of the 29th International Coference on International Conference Machine Learning. pages 1675-1682, 2012.
[46] Christopher Zach. Fully variational noise-contrastive estimation. In Image Analysis: 23rd Scandinavian Conference, SCIA 2023, Sirkka, Finland, April 18-21, 2023, Proceedings, Part II, pages 175-190. Springer, 2023.
[47] Shuyi Zhang, Bin Guo, Anlan Dong, Jing He, Ziping Xu, and Song Xi Chen. Cautionary tales on air-quality improvement in beijing. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2205):2017045,, 2017.
[48] Yanqing Zhang, Xuan Bi, Niansheng Tang, and Annie Qu. Dynamic tensor recommender systems. The Journal of Machine Learning Research, 22(1):3032-3066, 2021.
[49] Qibin Zhao, Liqing Zhang. and Andrzej Cichocki. Bayesian ep factorization of incomplete tensors with automatic rank determination. IEEE transactions on pattern analysis and machine intelligence, 37(9):1751-177 2015.
[50] Qibin Zhao, Liqing Zhang, and Andrzej Cichocki. Bayesian sparse tucker models for dimension reduction and tensor completion. arXiv preprini arXiv: 1505.02353. 2015.
[51] Qibin Zhao, Guoxu Zhou, Shengli Xie, Liqing Zhang, and Andrzej Cichocki. Tensor ring decomposition. arXiv preprint arXiv:1606.000 05535, 2016.
[52] Shandian Zhe, Zenglin Xu, Xinqi Chu, Yuan Qi, and Youngja Park. Scalable nonparametric multiway data analysis. In Artificial Intelligence and Statistics, pages 1125-1134, PMLR, 2015.
[53] Shandian Zhe, Kai Zhang, Pengyuan Wang, Kuang-chih Lee, Zenglin Xu, Yuan Qi, and Zoubin Ghahramani. Distributed flexible nonlinear tensor factorization, Advances in neural information processing systems, 29, 2016.
13

---

A Proof of Theorem 1
Firstly, we give the definition of f-divergence.
Definition (f (f-divergence) The f f-divergence between two probability density functions (pdf) p and is defined as,
Dy(pllg) =E, ( (
where f: : |0 oc) I. is a convex function and f(1)
As shown in [33], since partition functions for (b(r, m: 0) and o(2; 0) are the same, have the following factorization,
Φ(z, m: 0) o(x; pm x; The difference between the two objective becomes, LVENCE(O, s)) - LONCE(0)
1 o(y; O)p v)g(m;
- log + ol: 0)p (x y)
log
Φ(2,, mi (1) x)
Ø(x; Py x)
( rr m. 0)0(;; 0 P x)) + d(y; 0)p ())(;; (){(m;
-2F Ey(mm)) log
Φ xr m; m; ())(); p
x) + o(v; 0)pe(z 1))(( m;
p(m m x:; x;)))(
x)) o(y; 0 () v)g(m: 4)
(15m:Byll(m))))) log
pmm xz; 0);((;;;;) (y
x) oly; ()pe (x y)p(m x;
B2 [D, Jay (p(m)
x; (m)))
where
Kay + w-
fay(u) log
Kizy +
with Kzy
((vl)))e(vl It straightforward to verify that f(1) 0. The derivatives of f is (( w)
2um +
('()) u + ("(u)
(u2x + u)²
Since x and are positive, f is a convex function. Therefore, f satisfy the requirements of f-
divergence.
B Proof of Corollaries and
Corollary is straightforward consequence of Theorem Since the f-divergence becomes zero and only if the two distributions are identical, we have,
LVENCE :p) = LoNCE(0) + ((m;))))= = p(m x;; 0).
Moreover, since the f-divergence is positive and Theorem we have
p(m x; 0) = arg min LVeNCE(O, ((m;; ((va:p))
Then, plugging the optimal distribution gives the tight bound, have,
min LONCE(0) = min min LVCNCE(O, p). ((m)p))
C Experimental details
C.1 Simulation study
Tensors with non-Gaussian distributions For both GPTF and our model, we set batch size to 1000 and run 500 epochs with Adam optimizer. The initial learning rate is 1e-3 and subsequently reduced
14

---

by 0.3 at 60%, 75% and 90% of the maximum epochs. Moreover, the rank is set to 3 for both models. For GPTF, radial basis function (RBF) kernel with band width 1.0 used, where 100 inducing points is adopted for approximation. For the conditional distribution p(æ; m;) Vr; [mm)). o?) in GPTF, a fixed and chosen as the sample standard variance. For our model, we use hidden layers of width 64 for both 91. 93 and Is defined in Section 3. 92 is summation layer. We use ELU activation for non-linearity. For the VCNCE loss, the conditional noise distribution is set as Po(y z) = Nyy 2, 0.3²) and 10 noise samples are used for each data point.
Continuous-ting tensors The data sizes and optimization parameters are the same with the previous simulation. The rank of all models are set to 3. For NONFAT, 100 inducing points are used approximate the kernel function. We run the NONFAT model for 5000 epochs because we find that the algorithm converges very slowly. Other hyper-parameters settings. For BCTT, we do not modify their code and settings. For our model, we use 3 hidden layers of length 64 with ELU activation. The conditional noise distribution in the VCNCE loss is set to Poly x) = Ny y x, 1) and 20 noise samples are used for each datum.
chosen by their default
C.2 Tensor completion
For all datasets, when training our model, we scale the data to [0. 1] based on the training data. For testing, we multiply the scale statistic computed by the training data and evaluate the performance on the original domain. We do not employ such data normalization for baselines models, because that will influence their default settings.
C.2.1 Sparse tensor completion
For both Alog and ACC, the batch size is set to 1000. We run 1000 epochs for Alog and 100 epochs for ACC due to their different sample numbers. For Alog dataset, we add i.i.d. Gaussian noises from (O0, 0.05²) during training, while for ACC, the standard variance is set to 0.02. The Adam optimizer is used with learning rate chosen from (1e-2, 1e-3, le-4]. We also use gradient clip with maximum infinity norm of 2.0 for training stability. Moreover, we use learning rate scheduler by reducing the initial learning rate by 0.3 at 40%. 60%, and 80% of the total iterations. For both datasets, use 2 hidden layers of length 50 with ELU activation for 91, 93 and 94 for our model. For the VCNCE loss, we set // = 20 noise samples with noise variance tuned from {0.3², 0.5², 0.82, 1.0"). In practice, we find that the noise variance is influential to the final performance, even using conditional noises. However, with VCNCE, there is only one hyper- parameter for the noise distribution. While for CNCE, one may need to tune both mean and variance of the noise.
C.2.2 Continuous-time tensor completion
For Air and Click datasets, we set batch size to 128. We run 400 epochs for Air and 200 epochs for Click due their different data sizes. For Alog dataset, we add i.i.d. Gaussian noises from (0,0 0.05?) during training, while for ACC, the variance is set as 0.15?. To encode the temporal information into the energy function, we use the sinusoidal positional encoding, described in Section 3. Other settings are the same with Appendix C.2.1.
It should be noted that we use the standard definition of root mean square error (RMSE) and mean absolute error (MAE), namely,
En (z+ 2,22 N
N I 2z/
RMSE
MAE
where I; is the ground truth and i;j is the estimate. Therefore, the results are different from those presented in [43], where the authors used relative versions of RMSE and MAE,
N
[1;
22)2
N
T
RMSE
MAE
1=1
We modify the evaluation part of their code? and report the results.
hhttps://iithub http://iitth...ooo//utuutttt con/vzhut/NONFA.
15

---

C.3 Computational time
The time complexity of our model is v times greater than traditional TDs, since we need compute forward passes for particles. However, as we only use small networks, the computational speed still very fast. To illustrate, we compare the runtime performance of several baselines and our model on single RTX A5000 GPU. We conduct tests on the Air dataset, with batch size of 128 and tensor rank 5. The reported running time the average of the first 10 epochs. For our model, we set = 20. The default settings are used for other baselines. Table 3 lists the computational time of CTGP, NNDTN, NONFAT and THIS-ODE, all of which perform better than other baselines. The results show that our model achieves better performances within reasonable time, especially compared to THIS-ODE.
Table 3: Computing time
CTGP
NNDTN
EnergyTD 2.51+0.13 464.+131. 5.30++0.37
NONFAT
THIS-ODE
Time//ipoch (in seconds)
1.17++0.30
2.184
C.4 Ablation study on the objective function
We conduct additional ablation study to show the advantage of VCNCE over the variational noise- contrastive estimation [VNCE, 33] objective. The main difference between the VNCE and VCNCE that VNCE uses noises from fixed Gaussian distribution, e.g., y ~ Pv (y) = Ny μ, o²)²), while VCNCE uses conditional noises, e.g., Pe(y x) N(w xr (.). Hence, these two strategies yield different objective functions. The objective function of VNCE defined as
Ó(x, me; 0)
CYNCE Egggg(m log
Φ(z m; + vm(m x; (x)
vpn(y)
+ MEy log
vpn(yy + (mm/))
<(y.m:0) g(mly)
where Pn (-) is the fixed noise distribution. For VNCE, choosing inappropriate noise distributions may result in bad performances.
We test the proposed model on the Air dataset, training on the VCNCE loss and VNCE loss, respectively. We set the batch size 128 and run 400 epochs. Adam optimizer with initial learning rate 1e-2 is adopted. The initial learning rate is subsequently reduce by 0.3 at 20%, 50% and 80% of the total epochs. For VNCE, set μ 0, which is common practice in relevant literature. To show how the noise variance affects the learning process, we test different noise variances, egg., {0.3, 0.5, 0.7] for both VNCE and VCNCE. Other settings are the same with Appendix C.2.2.
Fig. 3 depicts the RMSE and MAE on the test data when optimizing VNCE and VCNCE objective functions. We test five runs, plot mean values in lines and standard deviations in shadowed areas. It shown that VCNCE gets better and more stable results on both RMSE and MAE.
16

---

0.7
0.7
0.7
mmm
mma
RMSE
RMSE
RMSE
0.5
0.5
0.5
0.3
0.3
0.3
200
400
200
Epoch
Epoch
Epoch
VNCE VCNCE
0.5
M
0.4
0.4
MAE
MAE
0.3
0.3
0.3
0.2
0.2
0.2
100
200 Epoch
300
400
100
200
300
400
100
200
300
400
Epoch
Epoch
(a) Noise 0.3
(b) Noise 0.5
(c) Noise 0.7
Figure 3: Learning process of optimizing the VNCE and VCNCE loss. The first row is RMSE and the second row is MAE.
17


Reviews:
[{"review_id": "yM8lbgiT7l", "reviewer": "Reviewer_buha", "summary": "The paper proposes an innovative approach for non-linear tensor decomposition. It utilizes the deep energy-based model (EBM) to model the joint energy function of tensor observations and latent factors. This design enables a more flexible decomposition without the need for explicitly defining the interaction between the latent factor and the commonly used Gaussian prior. Furthermore, the expressive nature of the energy function allows for the modeling of side information, such as timestamps.", "questions": "See weakness", "limitations": "See weakness", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper presents a clear motivation and introduces innovative and sophisticated solutions. To my knowledge, this is the first work to adopt Evidence-Based Medicine (EBM) for tensor decomposition without the need for an explicit design of the decomposition form.\n\nFurthermore, the adaptability of the proposed method for dynamic tensors is impressive. The experimental setup is robust and comprehensive, demonstrating the thoroughness of the research.", "weaknesses": "1. Instead of using NON-FAT or BCTT for experiments on dynamic tensors' decomposition and visualization, it would be more appropriate to choose THIS-ODE[1] as the baseline. While NON-FAT and BCTT are designed for modeling temporal dynamics of factors in a tensor, THIS-ODE aligns with the proposed model by having an entry-wise dynamic design.\n \n2. It appears that the objective function involves entry-wise importance sampling, and completing the posterior will require another round of entry-wise sampling. I am unsure if this will lead to training and inference costs becoming problematic, especially in scenarios with a large number of entries. It would be beneficial to have a more in-depth discussion about the scalability aspect.\n\n\n[1] Decomposing Temporal High-Order Interactions via Latent ODEs, Shibo Li, Robert M. Kirby, and Shandian Zhe, The 39th International Conference on Machine Learning, 2022", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## W1: Comparison with THIS-ODE.\n\nThanks to the reviewer for pointing out this reference, which is very relevant to our model. We would like to add some discussions about it in the manuscript and compare it with our model in experiments. However, its official implementation runs very slow and we are not able to complete very sufficient experiments due to the time limit. Here is a brief comparison about the computational time of several models, running on a single RTX A5000 GPU. Apart from THIS-ODE, we also list the computational time of CTGP, NNDTN and NONFAT, which perform better than other baselines. We test on the *Air* dataset, setting batch size to 128 and tensor rank to 5, and report the running time of one epoch by averaging on 10 runs. For our model, we set $\\nu = 20$. For NONFAT and THIS-ODE, the default settings are used.\n\n|                         | CTGP            | NNDTN           | NONFAT          | THIS-ODE        | EnergyTD        |\n| ----------------------- | --------------- | --------------- | --------------- | --------------- | --------------- |\n| Time/Epoch (in seconds) | 1.17 $\\pm$ 0.30 | 2.18 $\\pm$ 0.04 | 2.51 $\\pm$ 0.13 | 464. $\\pm$ 131. | 5.30 $\\pm$ 0.37 |\n\nAs we can see, our model is slower than NONFAT, but much faster than THIS-ODE. While NONFAT is faster for each epoch, it converges much slower than our model. As suggested by the paper of NONFAT, we run 10,000 epochs in order to get a good performance. For our model, only 400 epochs is sufficient to converge.\n\nFor THIS-ODE, it takes more than 30 hours to run 100 epochs on the *Air* dataset, and more than 40 hours to run 50 epochs on the *Click* dataset. Due to the time limit, we do not run more epochs. Moreover, the learning process shows that the model almost converges at this point, so the results are also convincing. The preliminary results are listed below and we will definitely conduct more sufficient experiments later. Our model outperforms THIS-ODE with much lower computational cost.\n\n\n|          | RMSE              |                   |                    |                   | MAE               |                   |                   |                   |\n| -------- | ----------------- | ----------------- | ------------------ | ----------------- | ----------------- | ----------------- | ----------------- | ----------------- |\n| *Air*    | Rank 3            | Rank 5            | Rank 8             | Rank 10           | Rank 3            | Rank 5            | Rank 8            | Rank 10           |\n| THIS-ODE | 0.588 $\\pm$ 0.003 | 0.612 $\\pm$ 0.000 | 0.586 $\\pm$ 0.000  | 0.578 $\\pm$ 0.000 | 0.434 $\\pm$ 0.002 | 0.451 $\\pm$ 0.000 | 0.426 $\\pm$ 0.000 | 0.424 $\\pm$ 0.000 |\n| EnergyTD | 0.302 $\\pm$ 0.008 | 0.291 $\\pm$ 0.006 | 0.300 $\\pm$ 0.012  | 0.283 $\\pm$ 0.004 | 0.184 $\\pm$ 0.006 | 0.177 $\\pm$ 0.003 | 0.172 $\\pm$ 0.006 | 0.184 $\\pm$ 0.003 |\n| *Click*  |                   |                   |                    |                   |                   |                   |                   |                   |\n| THIS-ODE | 1.408 $\\pm$ 0.009 | 1.409 $\\pm$ 0.008 | 1.405 $\\pm$ 0.005. | 1.405 $\\pm$ 0.008 | 0.846 $\\pm$ 0.009 | 0.835 $\\pm$ 0.007 | 0.846 $\\pm$ 0.002 | 0.843 $\\pm$ 0.007 |\n| EnergyTD | 1.396 $\\pm$ 0.003 | 1.385 $\\pm$ 0.003 | 1.356 $\\pm$ 0.001  | 1.357 $\\pm$ 0.001 | 0.777 $\\pm$ 0.003 | 0.775 $\\pm$ 0.003 | 0.772 $\\pm$ 0.002 | 0.773 $\\pm$ 0.001 |\n\nCurrently, we have not got very faithful trajectory estimates using the simulation data with THIS-ODE\n\n## W2: Training and inference cost.\n\nFirstly, for the training, the computational cost mainly comes from the evaluation of the energy function. In particular, we need to evaluate the energy function for $\\nu$ times, where $\\nu$ is the number of conditional noises and is typically 20 in our experiments. The time complexity is $\\mathcal{O}(\\nu B ( D R H + L H^2))$, where $B$ is the batch size, $\\nu$ is the number of conditional noises, $H$ is the number of hidden units per layer, $L$ is the number of layers and $D$ is the tensor order. Since we are using very small NNs, e.g., several MLP layers with hundreds of hidden units each, the computational cost is not large. From the table above, we can see that our model is slower than NONFAT, but much faster than THIS-ODE. More importantly, since the forward passes for each particle can be computed in parallel, we believe the computational speed can be easily improved using parallel programming libraries such as `JAX` and `functorch`.\n\nThen, for the inference, apart from traditional MCMC samplers, we can also perform a very effective grid search to find MAP estimates in parallel. This can be much faster than the THIS-ODE model, which uses time-consuming diffusion processes for training and inference."}, {"review_id": "bTI5arDBNM", "reviewer": "Reviewer_32MU", "summary": "This paper uses the Energy Based Model (EBM) framework to capture the joint probability of the data and latent tensor factors to learn as much information from data as possible, which discards the structural and distributional assumptions and thus avoid picking an inappropriate TD model.  To further flexibly learn the unnormalized probability function, the authors derive a variational upper bound of the conditional noise-contrastive estimation objective.  The main contribution lies in inserting the EBM framework into modeling probabilistic Tensor Decompositions.", "questions": "Questions:\n\n1. As mentioned in the last part, I am confused about why not directly model $\\phi(x;\\theta)$ and $p_{\\theta}(x)$.", "limitations": "Yes, the authors have adequately addressed the limitations and potential negative societal impact of their work.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Pros:\n\n1. This paper first introduces undirect graphic models into modeling probabilistic tensor decomposition. In detail, the authors consider latent factor matrics as latent factors in the Variational Noise Contrastive Estimation framework, leading to flexible modeling of joint distribution over observations. To the best of my knowledge, this is interesting and novel to the Tensor Decomposition community.\n\n2. This paper is well-written and the technological parts are sound and valid.", "weaknesses": "Cons:\n\n1. If what I understand is correct,  one important goal of this paper is to find the $\\theta$ and $p_{\\theta}(x)$. The authors use the framework of variational noise contrastive estimation to obtain the unnormalized function $\\phi(x;\\theta)$. However, it seems to be easier to directly model $\\phi(x;\\theta)$ and thus minimize $\\mathcal{L}_{CNCE}(\\theta)$. Furthermore, it seems that these $\\textbf{m}$s are unnecessary and can have arbitrary dimensions.\n\n2. Most technical contributions are in the form of a combination of existing theory and framework. This actually decreases this paper's contributions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## W1: Directly modeling $\\phi(x; \\theta)$ and the choice of $m$.\n\nThanks for the question. Apart from the data distribution $\\phi(x; \\theta)$, in tensor decompositions, we aim to infer the latent factors, i.e., $z$ in the manuscript. If we directly model $\\phi(x; \\theta)$, we cannot obtain the information of these latent factors. Therefore, we need to model the joint energy function $\\phi(x, m; \\theta)$, where $m$ consists of $z$. Equipped with this joint energy function, we can sample posterior of the latent factors from $p(z \\mid x)$. However, the data energy function $\\phi(x; \\theta) = \\int \\phi(x, m; \\theta) dm$ becomes intractable, since we adopt non-linear mappings. This is one main difficulty in the learning process. Variational inference (VI) has become one of the most popular tools to deal with such intractable integrations of latent variables. By using VI, we finally obtain the lower bound of the CNCE loss in Eq. (5), which is tractable to compute.\n\n$m$ actually consists of latent tensor factors (as described in Line 132) and the size of $m$ is $D R$, where $D$ is the tensor order and $R$ is the tensor rank. The rank serves as a trade-off between model complexity and expressive power, which is usually treated as a hyper-parameter. In tensor decomposition, these tensor factors are crucial to learn latent correlations of the high-order data. Taking tensor completion applications as an example, once we learn these tensor factors, we are able to predict missing entries given very sparsely observed data.\n\n## W2: Contribution.\n\nOur main contribution is to establish a new tensor decomposition model. As **Reviewer buha** suggested, there is no previous literature trying to use EBMs to model tensor data. Our model shows the potential to deal with diverse types of applications in tensor decompositions. Although both tensor decompositions and EBMs are existing models, it is generally hard to learn EBMs for such high dimensional data in the presence of latent factors. To enable efficient learning processes, we establish the variant of CNCE loss, as well as ad hoc network architectures, that finally result in a powerful TD model.\n\n## Q1: Why not directly model $\\phi(x; \\theta)$ and $p_\\theta(x)$.\n\nPlease see our response to *W1*."}, {"review_id": "mV3hozfFS7", "reviewer": "Reviewer_Duah", "summary": "This paper proposes a probabilistic tensor decomposition model called EnergyTD that integrates a deep energy-based model (EBM) in tensor decomposition. The EnergyTD does not model the values in a tensor as the conditional probability conditioned on the latent factors based on the predefined models but models them as the joint probability of the data and the latent factors. It provides a learning algorithm that extends the CNCE loss using the variational approach. The experimental results show that EnergyTD outperforms conventional methods in the tensor decomposition tasks assuming diverse distributions including continuous-time tensors, tensor completion, and continuous-time tensor completion.", "questions": "What is the time/space complexity of the proposed method? How is the proposed method computationally effective compared to the conventional methods?", "limitations": "As the authors describe in the Conclusion, there is room to improve theoretical analysis. Moreover, the architecture of the network structure should be discussed more in the future.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "This paper proposes a probabilistic tensor decomposition model that integrates EBM in tensor decomposition along with a theoretical guarantee. The experimental results show that the proposed method outperforms conventional methods in the tensor decomposition tasks assuming diverse distributions. This paper is well-organized and easy to follow.", "weaknesses": "The time/space complexity of the proposed method is only poorly discussed, and there is no evaluation regarding the computational cost. This is one of the major concerns in the tensor decomposition research community.\n\nA minor problem: The \"difficultly\" in line 109 should be \"difficulty\".", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## W1: Time/space complexity.\n\nThe time complexity of training should be $\\mathcal{O}(\\nu B ( D R H + L H^2))$, where $B$ is the batch size, $\\nu$ is the number of conditional noises, $H$ is the number of hidden units per layer, $L$ is the number of layers and $D$ is the tensor order. The space complexity should be $\\mathcal{O}(D R + LH^2)$ to store tensor factors and NN parameters. The space complexity is linear to the tensor order and rank, which is similar to many classical TDs such as CP and tensor train decomposition. We additionally use an NN. However, the network is small, e.g., several MLP layers with hundreds of hidden units each. The time complexity is $\\nu$ times more than traditional ones, since we need to compute forward passes for $\\nu$ particles. However, since these forward passes can be computed in parallel, we believe the computational speed can be easily improved using parallel programming libraries such as `JAX` and `functorch`. Moreover, unlike many traditional TDs performing ALS updates on the whole datasets, we use stochastic optimization, which is highly scalable to large datasets.\n\nTo illustrate the computing time, we post an example of training the *Air* dataset, with rank 5. We compare with CTGP, NNDTN, NONFAT, which are better than other baselines, and THIS-ODE suggested by **Reviewer buha**. All experiments are conducted on a single RTX A5000 GPU. For more details, please check our response to **Reviewer buha**.\n\n|                         | CTGP            | NNDTN           | NONFAT          | THIS-ODE        | EnergyTD        |\n| ----------------------- | --------------- | --------------- | --------------- | --------------- | --------------- |\n| Time/Epoch (in seconds) | 1.17 $\\pm$ 0.30 | 2.18 $\\pm$ 0.04 | 2.51 $\\pm$ 0.13 | 464. $\\pm$ 131. | 5.30 $\\pm$ 0.37 |\n\nOur model is slower than NONFAT, but much faster than THIS-ODE. While NONFAT is faster for each epoch, it converges much slower than our model. As suggested by the paper of NONFAT, we run 10,000 epochs to get a good performance. For our model, only 200 epochs is sufficient to converge.\n\n\n## W2: Typo\n\nThanks for pointing out the typo. We will fix it.\n\n\n## Q1: What is the time/space complexity?\n\nPlease see our response for *W1*."}, {"review_id": "vTAVApqYs8", "reviewer": "Reviewer_7dHC", "summary": "The authors propose a new Probabilistic Tensor Decomposition (TD) method. By modeling the joint probability of the data and latent tensor factors using an Energy Based approach, they make no (possibly restrictive) structural and distributional assumptions on the generative process linking latent and observations. The model is trained using an upper bound of the Conditional Noise-Contrastive Estimation Loss. It is then evaluated on simulations using non Gaussian generative models and on tensor completion task on real world datasets.\n", "questions": "1) The paper only discusses 2 and 3-way tensors. How does the method scale with higher order tensors ?  \n\n2) Beside Tensor completion, an interesting aspect of TD methods is often to explain a dataset $\\mathcal{X}$ using simpler yet unobserved factors $\\mathcal{Z}$. How robust are the inferred factors depending on the initialization of the model ? Can you measure it using e.g. a similarity metric and compare it to existing methods ?\n\n3) The authors argue that the use of overly simple generative model $p(x|z)$ can bias the TD decomposition. They provide convincing illustration on their toydataset but it would be interesting to better understand why it is the case on the tested the real world dataset.  \n\n4) The approach developed by the authors seems extremely flexible. Yet, by allowing almost arbitrary complex link between the latent and the observation, one might fear that the inferred latent distribution is \"non-unique\", which might hinder the interpretation of the discovered factors. Can you say more about identifiability ?", "limitations": "There is no discussion on the limitation of the method when the posterior $p(z|x)$ itself is non Gaussian. The approximation $q(m; \\phi)$ is set to be multivariate normal with diagonal covariance, and, so, it cannot model, for example, multi-modal distribution on $z$. \n\n", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The possibility to bypass the specification of both structural and distributional dependencies is very interesting. It has many applications beyond the tensor completion tasks, for example when analyzing datasets with complicated dependencies on significantly simpler latent variables. \n\nThe simulated datasets used in the paper are simple but illustrate this point well. \n\nResults on real world tensor completion tasks are convincing.", "weaknesses": "Major Points:\n\n1) The code is not provided. \n\n2) In the introduction, it would be worth motivating the use of probabilistic methods over traditional ones.\n\n3) I did not find the Energy Based Model section (2.2) very well written. In particular, it is not clear how to choose the conditional $p_c(y|x)$  to improve efficiency. \n\n4) The method is only evaluated on tensor completion tasks. Do the inferred factors have any explanatory value ? How robust is the decomposition ?\n\n5) There is no discussion on the compute time of the proposed method.\n\nMinor Points:\n\n1) l.43 It is worth mentioning that Generalized CP decomposition exist but that most of them do not allow a probabilistic treatment of the data. Moreover, beyond Gaussian and Binomial, tractable decomposition have been developed for Negative Binomial Distributions.\n\n2) l.93 I think the authors implicitly assume $R_d = R$ for all $d$ without explicitly saying it.\n\n3) Experimental details in the main paper and supplementary are not consistent (experiment 1: main paper $R=5$, supplementary $R=3$).\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "## W1: Code availability.\nWe send an anonymous link to the AC via the official comment.\n\n## W2: Motivation of using probabilistic methods over traditional ones.\nThere are several advantages of using probabilistic models, such as,\n\n1. Using probabilistic models, we can deal with different data types.\nFor examples, by adopting different distributions, it provides a principled way\nof designing proper loss functions.\n2. Probabilistic models can give uncertainty estimates about the observations and latent factors.\n3. Bayesian models have other potential benefits, such as continual learning, less likely to over-fit and so on.\n\nWe will add more discussions about this part in the introduction.\n\n## W3: How to choose conditional noises to improve efficiency.\n\nThe basic idea of NCE is to identify the noises from the data points, which shares some similarities with density ratio estimation such as in GAN. Unfortunately, there is no theoretical guidance about how to choose noise distributions yet. A common belief is that the noises should be similar to data. This is the motivation of using data-dependent noises, rather than the same noise for each data point. We will rearrange the presentation to make it more clear.\n\n## W4: Explanatory value of factors? How robust is the decomposition?\n\nCurrently, we do not focus on learning explanatory factors. To learn such factors, one may need to add additional constraints on the tensor decomposition, e.g., identifiability of the factors. This is generally hard for non-linear models using neural networks. However, we conduct some simple investigation on the *Alog* dataset following the setting of [1]. In specific, the tensor latent factors are firstly mapped to the two-dimensional plane through PCA and then the $k$-means clustering results are plotted for each data point. The figure can be found in the attached PDF. We find hidden patterns similar to Figure 4 in [1]. Moreover, we investigate the results for different ranks and initializations, showing that this pattern is not due to random effects. Although we initially did not consider this problem, our model seems to have better clustering results compared with other baselines. This indicates the potential of our model to discover explanatory latent factors.\n\nCurrently, we do not conduct experiments about robustness. It would be a very interesting direction since our model is capable of handling different data distributions. Moreover, the adversarial robustness would also be an interesting direction to investigate.\n\n[1]. Tillinghast, C., Wang, Z., & Zhe, S. (2022). Nonparametric sparse tensor factorization with hierarchical Gamma processes. In ICML.\n\n## W5: Compute time.\n\nDue to the space limit, please check our response to *W1* raised by **Reviewer Duah*.\n\n## Minor1: Tractable tensor decomposition for other distributions\n\nWe agree that there are many tractable TDs for other distributions. We will add several references about this point, such as [1-3].\n\n[1]. Schein, A., Zhou, M., Blei, D. & Wallach, H.. (2016). Bayesian Poisson Tucker Decomposition for Learning the Structure of International Relations. In ICML.\n\n[2]. Hong, D., Kolda, T. G., & Duersch, J. A. (2020). Generalized canonical polyadic tensor decomposition. SIAM Review.\n\n[3]. Soulat, H., Keshavarzi, S., Margrie, T., & Sahani, M. (2021). Probabilistic tensor decomposition of neural population spiking activity. In NeurIPS.\n\n## Minor2: Assume $R_d = R$.\n\nWe will highlight this assumption.\n\n## Minor3: Inconsistent setting of rank.\n\nWe checked experimental records and the rank should be 3. We will fix it.\n\n## Q1: Scalability with higher order tensors\n\nIn principle, both the time and space complexity of our model is linear with the tensor order $D$. Given the same number of observations $N$, our model is scalable with the tensor order $D$.\n\n## Q2: Robustness of the learning process of latent factors.\n\nGenerally, we think it is hard to learn  identifiable latent factors for non-linear models, especially using neural networks. However, we do find robustness of these latent factors in some sense. For example, in our response to *W4*,\nwe plot the latent factors for different ranks and different initializations.\nThey show similar patterns in the clustering.\nThis would be an interesting direction to study further, including new metrics to evaluate these similarities.\n\n## Q3: Different distributions in real applications.\n\nDemonstrating this phenomenon is hard in real applications, since we do not know the underlying real distributions. However, we can find some evidence from current results. In Table 1 of the manuscript and the additional results in our response to **Reviewer buha**, we observe that the improvements in MAE are notably more significant. One possible reason is that other baselines are trained implicitly by minimizing the square loss as they adopt Gaussian assumptions about the data. However, our model does not make such assumptions about the data and the loss function. Hence, it can adapt to the data distribution more flexibly. Additionally, in Table 2 of manuscript, while simply using NN to model the trajectory is hard, as shown by NNDTL and CTNN, our model outperforms them and NONFAT, which is the SOTA in the field. This may indicate that our model efficiently learns how the data distribution shifts with time.\n\n## Q4: Identifiability about the tensor factors.\n\nUnfortunately, we think the model cannot learn identifiable latent factors currently, since we are using non-linear neural networks. However, as our response to *W4*, our model possibly learns some sharing latent patterns under different settings and different initializations.\n\n## Limitation: Non-Gaussian posteriors.\n\nWe agree with the reviewer about this limitation. However, due to the computational efficiency, mean-field Gaussian assumption is widely used in variational inference. It could be a potential direction to improve our model, by using more expressive distributions."}]

Author rebuttal:
We thank all the reviewers for their constructive comments and suggestions. Below we respond to them respectively.

Meta-review: Looks like a good paper. Accept.
Decision: Accept (poster)