E$^2$M: Double Bounded $\alpha$-Divergence Optimization for Tensor-based Discrete Density Estimation
Authors: Kazu Ghalamkari
Jesper Løve Hinrich
Morten Mørup
Venue: Accepted by TMLR ()

Abstract: Tensor-based discrete density estimation requires flexible modeling and proper divergence criteria to enable effective learning; however, traditional approaches using α-divergence face analytical challenges due to the α-power terms in the objective function, which hinder the derivation of closed-form update rules. We present a generalization of the expectation-maximization~(EM) algorithm, called the E2M algorithm. It circumvents this issue by first relaxing the optimization into the minimization of a surrogate objective based on the Kullback–Leibler (KL) divergence, which is tractable via the standard EM algorithm, and subsequently applying a tensor many-body approximation in the M-step to enable simultaneous closed-form updates of all parameters. Our approach offers flexible modeling for not only a variety of low-rank structures, including the CP, Tucker, and Tensor Train formats, but also their mixtures, thus allowing us to leverage the strengths of different low-rank structures. We evaluate the effectiveness of our approach on synthetic and real datasets, highlighting its comparable convergence to gradient-based procedures, robustness to outliers, and favorable density estimation performance compared to prominent existing tensor-based methods.

Published in Transactions on Machine Learning Research (4/2026)
M2M: Double Bounded a- a-Divergence Optimization for Tensor-based Discrete Density Estimation
Kazu Ghalamkari
kezfu dk
Deportment of Applied Mathematics and Computer Science Technicol University of Desmark
Jesper Love Hinrich
jehi@dtu. dk
Deportment of Applied Mathematics and Computer Science Technical University of Denmork
Morten Mørup
mmorūdtu. dk
Department of Applied Muthematics and Computer Science Technical University of Denmark
Reviewed on OpenReview: https: openreview. net/ formm? id- 964C)hXXXI
Abstract
Tensor-based discrete density estimation requires flexible modeling and proper divergence criteria to enable effective learning; however, traditional approaches using a-divergence face analytical challenges due to the a-power terms in the objective function, which hinder the derivation of closed-form update rules. We present generalization of the expectation- maximization (EM) algorithm, called the E algorithm. It circumvents this issue by first relaxing the optimization into the minimization of a surrogate objective based on the Kullback Leibler (KL) divergence, which is tractable via the standard EM algorithm, and subsequently applying tensor many-bodt approximation in the M-step enable simultaneous closed-form updates of all parameters. Our approach offers flexible modeling for not only varicty of low-rank structures, including the CP, Tucker, and Tensor Train formats, but also their mixtures, thus allowing us to leverage the strengths of different low-rank structures. We evaluate the effectiveness of our approach synthetic and real datasets, highlighting its comparable convergence to gradient-basee proceduros, robustnoss to outliers, and favorable density estimation performance compared to prominent existing tensor-based methods.
1 Introduction
Tensors are versatile data structures used in broad range of fields, including signal processing Sidiropoulos et al. 2017), computer vision (Panagakis et al., 2021), and data mining (Papalexakis et al.) 2016). It is an established fact that features can be extracted from consor-formatted data by low-rank decompositions that approximate the tensor by a linear combination of a few bases (Cichooki al. 2016 Liu et al., 2022 There are numerous variations of low-rank structures used for the decomposed represcntation, such as the CP [Hitchoock] 1927; Carroll & Chang 1970), Tucker (Tucker) 1966), and Tensor Train (TT) (Oseledets) 2011 and the user typically selects evaluates each choice of structure exhaustively in order to identify suitable decomposition.
A series of recent sludies (Glasser et al. 2019 Novikov et al. 2021 show that low-rank tensor decompositions are also useful for discrete density estimation by taking advantage of the discreteness the tensor indices. Specifically, given observed discrete samples (1)),.. (N), the empirical distribution p(x:) can be regarded normalized non-negative tensor and its low-rank reconstruction P estimates the distribution underlying the samples, rated in Figure [a).

---

Published in Transactions on Machine Learning Research (4/2026)
a.
T
(x40)
22
xx(()
p
I (F, (F ②)
Discrete tabular data
12
U uper bound of upper bound
x1
(1)
x
(1)
x5²
E2 step (update Φ)
1/N
a-div.
L, P) Upper bound
x1
x2
1/M
step
1
I/N
I/N
E1 step (update f)
/N
x1
x2
x
(N)
(
La E Empirical distribution
bjective function
M Mixture low-rank density
(n) Figure 1: (a) A discrete density estimation by samples x(1), ,(.),..., (N) for x(1) = (α) 22 23
(r)) and
(n) [Za]. The empirical distribution p(z) is identical to normalized non-negative tensor T, and the true distribution is estimated by its low-rank approximation P. (b) The E algorithm includes two F-steps. The El-step makes the upper bound tight w.r.t. the objective function, and the E2-step makes the upper bound of the upper bound tight w.r.t. the upper bound.
To learn an appropriate P from samples, it is required to optimize natural measures between distributions, such the a-divergence (Rényi 1961 Amari 2012 van Erven & Harremos 2014), which is rich family
of measures that includes the KL divergence, reverse KL. divergence, and the Hellinger distance as special
cases. The a-divergence thereby enables to adjust the modeling to be mass-covering or mode-seeking by
controlling the scalar a (Li & Turrner 2016). Although the usefulness of the a-divergence is well recognized
Lobatt in the machine learning community (Hernandez- Lobato et al. 2016 Li & Gal 2017 Cai et al., 2020 Gong
ctt al., 2021 Wang et al. 2021), its optimization is more challenging compared to KL divergence optimization due to the a-power lerms in the objective function.
The a-EM algorithm proposed by Matsuyama 2000 2003) provides generalized EM algorithm based on the a-log likelihood to optimize the a-divergence. However, this approach rarely yields closed-form update rules because of the non-linearily of the a-power dependent terms. Similarly, the multiplicative update (MU) for the a-divergence optimization (Kim et al., 2008) has no guarantee to provide closed-form update rules when imposing normalization constraints on the reconstructions in order to produce valid densities. As a result, optimization of the a-divergence often relies on gradient-base methods (Daudel et al., 2021 Daudel
& Doue 2021 Sharma & Pradhan 2023
Daudel et al. 2023b), requiring careful learning rate tuning.
To establish tensor-based discrete density estimation with closed-form updates for a-divergence optimization, we propose an alternative generalization of the EM algorithm, termed the E algorithm, consisting of two E-steps denoted E1 and E2 followed by an M-step. Specifically, our approach constructs double-bound on the a-divergence by use of Jenson's incquality twice. The first bound relaxes the a-divergence- basedd based low-rank approximation into the corresponding KI-divvrgence--baee second bound further relaxes the decomposition into a many-body approximation (Ghalamkari et al.] (002)), which represents tensors with reduced interactions among modes to be solved in the M-step. Importantly, derive exact closed-form solutions for the many-body approximations appearing in the M-step for the Tucker which
low-rank approximation while the
and TT decompasitions, which eliminates the need for learning-rate tuning in various low-rank approximations by combining these two solutions. By subsequently maximizing the expectations, the second bound (E2-step) and first bound (E1-step) become tight, as illustrated in Figure [b).
In addition, our double-bounding scheme naturally integrates the E1- and E2-steps into a unified E2-step, enabling memory-efficiet updates in the sparse density estimation setting. We also show that the tightness of the bounds ensures convergence of the E algorithm regardless of the choice of low-rank structure assumed in the model and establish the connection between approach and the recent auxiliary function minimization framework for a-divergence optimization in the continuous and variational inference setting (Daudel et al. 2023b). We further explore how our double bounds enable the proposed method to handle mixtures of low-rank structures. Our flexible mixture modeling automatically determines appropriate weights not only for mixed low-rank structures, climinating the need for the user to predefine a single low-rank structure advance, but also for adaptive background uniform density levels to make the learning stable. In summary, make the following contributions:
2

---

Published in Transactions on Machine Learning Research (4/2026)
We derive double bound of the a-divergence, which relaxes the low-rank approximation into a many-body approximation.
We establish the E framework to optimize the a-divergence for discrete density estimation based on non-negative low-rank mixtures with convergence guarantee and closed-form updates for various low-rank structures.
We extensively evaluate the developed E procedure for tensor-based discrete density estimation, considering its i) convergence properties compared to conventional gradient- based procedures, ii) stabilized learning induced by background density term, iii) ability to enhance robustness outliers by varying a, and iv) performance compared existing non-negative tensor- based methods.
1.1 Related work
This work focuses on the a-divergence, rich family of divergences, including the KL divergence. It has been reported that adjusting the scale parameter a controls the sensitivity to outliers (Miri Rekavandi et al.) 2024
The utility of the a-divergence has been validated in various contexts, including variational inference (Li &
Turner 2016), autoencoders (Daudel et al., 2023ab, generative adversarial networks (GANs) (Kurri et al.
2022), dictionary learning (Iqbal & Seghouanc 2019). and Gaussian processes (Villacampa-Catee Calvo et al. 2022
These works often rely on gradient-based methods, requiring careful learning rale Luning.
Several studies have shown that decomposing a normalized non-negative tensor constructed from observed samples be used to estimate the underlying distribution behind the data (Kargas et al. 2018 Glasser
et al. 2019 Vora et al.) 2021). In contrast to the well-established SVD-based methods for real-valued
tensor networks (Cheng et al. 2019 Tblisdir et al. 2007 Romén 2014), unified framework for normalized
non-negative tensor decompositions that supports various low-rank structures and optimizes natural measures
between distributions is lacking. Consequently, some existing studies of tensor- density estimation
have been performed via least. squares minimization (Kargas et al., 2018 Dolgov et al.) 2020 Novikov ct al. 2021 Vora et al. 2021) RRccently, density ostimation methods using second and third-order marginals have been developed (Kargas & Sidiropoulos, 2017 Ibrahim & Ful (202), which are so far limited to considering the CP decomposition- based model and typically require hyperparameters for the associated gradient-based optimization. In contrast, our approach directly applics to various low-rank structures with hyper-parameter--free optimization.
Although some studies employ the MU method to op- timize the divergence (Kim & Choil [0077] Kim et al.) 2008 Cichocki et al. 2009) (aMU), cxisting works do not take normalization into account, making them un- suitable for density estimation. We note that imposing normalization on the MU method is nontrivial, since the a-power in the auxiliary function makes it diffi- cult, to obtain the Lagrange multipliers in closed-form. Even if this limitation were to be addressed, the MU method still lacks flexibility for mixtures of low-rank models because the auxiliary functions include nonlinear terms such as a-power or a-logarithms, which prevent the decoupling of the components and thus hinder the derivation of closod-form update rules.
Table Low-rank- based discrete density learning for a-divergence.
E
aEM
aMU
GD
Normalization Step-size- free Closod-form update
X
/
/
/
V
X
X
X
Low-rank mixtures
X
X
/
While the EM algorithm natural approach to optimize the KL divergence (Dempster et al. 1977). EM-based normalized negativee tensor factorization has not been established except for the CP decomposition (Huang
& Sidiropoulos 2017
Yeredor & Haardt, 2019 Choge et al. 2022
The a-EM algorithm is known as a natural extension of the EM algorithm for a-divergence oplinization. However, the a-power term in the a-divergence prevents the decoupling of factors, and thus closed-form update rules are rarely available (Matsuyama) 2000 2003). These propertics are summarized in Table 1 Although the a-divergence is a variant of the f-divergence (Csiszár 1967] generalized EM algorithms for f-divergences have not been thoroughly investigated. This is likely because their naive formulations share
3

---

Published in Transactions on Machine Learning Research (4/2026)
the same limitation as a-EM: the lack of closed-form updates in the M-step, which requires pradient-based optimization at each iteration.
DNN-based density estimators While the effectiveness of deep neural network (DNN)-based density estimators in continuous settings has been demonstrated over the past decade (Rezende & Mohamed 2015 Papamakarios et al. 2021). an extension of these methods to discrete settings typically requires encoding categorical variables into numeric representations (Kobyzev et al.) 2021 For example, the naive extensions using simple one-hot encodings substantially increase the number of parameters; alternatively, embedding layers introduce additional hyperparameters (e.g., the embedding dimensionalityy.. On the other hand, tensor-based approaches naturally capture high-order interactions among discrete variables by leveraging the discreteness of tensor modes, and have repeatedly shown competitive performance against DNN-based methods (Novikov et al. discrete normalizing flows
2021 Amiridi et al. 2022 Wu et al.) 2025 Loconte et. al. 2025). More recently,
compositions of invertible discrete mappings that model probability mass directly without dequantization have attracted attention as they avoid the continuous relaxxttton step (Tran et al.
2019: Hoogeboom et al.) 2019 However, these mappings involve dineeifeernti operations, such as argmax, which complicates optimization and imposes constraints on network depth. As a result, discrete normalizing flows often require tuning an additional temperature hyperparameter, in addition to learning rates and architectural choices. Moreover, autoregressive discrete normalizing flows are order-dependent: the estimated discrete density is not invariant under permutations of the feature order. In contrast, one of the key advantages of tensor-based methods is that order dependence can be flexibly controlled by choosing a low-rank structure; for example, CP decompositions are order-indeeendent. whereas TT decompositions inherently depend on the ordering of the modes (Kolda & Bader 2009] Tichavský & Straka 2025). Additionally, depending on the optimization method, DNN-based methods require further assumptions
the gradient, such Lipschitz contimuity (Goodfellow et al., 2016) to guarantee convergence. In contrast, our proposed algorithm offers lyperparameter--ieee optimization and guarantees monotonic decrease and convergence without requiring any assumptions on the gradient.
Probabilistic tensor models for density estimation In probabilistic tensor decomposition, tensor elements are assumed to be sampled from a distribution pe. and the model parameters are typically estimated using the FM algorithm (Yilmaz Cemgil 2010: Rai et al.] [015 Hinrich Mørup 2019) under a specified prior, minimizing the KL-divergence. These approaches generally rely on conjugate model specifications, in which the posterior distribution shares the same distributional form as the prior. In contrast, our framework treats the entire non non-negative normalized tensor a discrete probability distribution, rather than modeling each tensor element as an independent sample from an underlying distribution.
1.2 Problem setup
We describe the problem setup here. Although all notations and abbreviations used in this paper are defined in the main text, tables summarizing the notations and abbreviations are provided in Tables and respectively, in the Appendix for convenience.
We consider N discrete samples (1)) (N), each consisting of D categorical features. Each feature value is encoded as natural number so that (n) E [71] .... x [Tn], where Ia denotes the number of categories of the d-th feature. The empirical distribution of the samples is represented by normalized D-th order tensor T E x..xx...xx Specifically, the tensor element T is regarded as the value of the distribution 20 p(#1 = .,,..., Φn = in). We denote the domain of é = (ii,..., in) by Ωr i.e., i 22 = [T]] x x.. x [Tn]. Each element of the tensor T is defined as T = L/N E 8(i, (")), where ò(ē, (()) is the Kronecker product, i.., 5(i,z(")) = 1 if i = ("), O otherwise. To estimate the discrete probability distribution underlying the samples, we consider the problem of finding a convex linear combination of K tensors plu),..., plk] that minimizes the a-divergence from the tensor T,
(N),
D. (T, a(1 - a)
1
Σ EEa
s.t. P: Σ 1lp!p and P!N
Σ rEA
olk
(1)
4

---

Published in Transactions on Machine Learning Research (4/2026)
where each plk] can be represented structures such as the CP, Tucker, or TT. Specifically, the normalized tensor Qlki and the index set Sgz will be introduced in this section. The mixture ratio = (1,")...,,11 satisfies Σx 1 = and 171h 2 0 for all k. The sctup for D = 3 is shown in Figure Ia).
normalized non-negative low-rank tensor using decomposition
Tt follows from l'Hôpital's rule (Brdlee et al. 2015) that Equation [ reduces to the KT. divergence in the limit as a + 1 and to the reverse KL divergence in the limit as a + 0. Accordingly, we define D (T, P) = (T, and P) DR. (P, T) where the KL divergence is defined DKi. (T. paper focuses on the range KL divergence to the mode-sccking KL divergence, including the Hellinger distance at a = 0.5. In this problem setup, the value of a is regarded as a hyperparameter controlling the model behavior, and we provide guidance on its selection in Section 2.3
Ti log (T//P.). This 1), which represents important spectrum from the crass-covering reverse
Low-rank structures on tensors Any non-negative low-rank tensor P.M. E 2220 R4Ix-II can be representod /xxxxxpxxxx x-xx Ryp ie., by the summation over indices r = (ri,...,Ty') of a higher-order tensor @l 2R
it holds that p!! = Erer, O!6 for 0g= = [Ri] x .... x [Ry*]. The integer V and the structure of Q depends on the low-rank structure in plk where k indicates a specific low-rank structure. The subscript ir denotes the concatenation of i E St and r E Agg, Le., ir = (fi....,in) in, (il..,.rry))) In the following, we consider the three prominent examples given by k E {CP, Tucker, TT}:
D = II d=1
D II
D II giu) ra 1ara d=1
(d) igr
of Mucker] = Gr
( lara d=1
TT ID
(2)
where the integer V is 1 for CP, D for Tucker, and D - for TT. We suppose that To rD - for (())
QITTI, On the hggt- side of Q.CP), the matrices A ..... AD) correspond to the factor matrices of R In the Tucker decomposition, the tensor is reconstructed ER Re In the TT decomposition, the tensor is reconstructed using sequence of D tensors ((1), ,.... "..,,....... G(D), where each gid) E R2 1xlaX Ri, We collectively refer to the matrices A(d) and the tensors g and g(d) as factors. Note that while pla is assumed to be normalized, the factors themselves are not required to be normalized. The tuple of integers Ryyx), representing the degrees of freedom for the index r E 2gz= = [R x... ..... [Ryx], corresponds to the CP rank, Tucker rank, or Tensor Train rank. In the following, square brackets are used to indicate low-rank structures, such as [CP), [Tucker], and [TT], whereas parentheses (.) are used mainly in Appendix to specify tensor factors.
the CP decomposition, where each A(d) E 20
using a core tensor 9 E R 20
x... RD and factor matrices A ...,, (D), where each A(d)
1.3 Remarks for tensor- based discrete density estimation
We here provide a remark regarding the formulation of tensor-based discrete density estimation and its scalability induced by the property of the a-divergence.
Tensor ranks as hidden variables Assuming the normalization of the tensor QikI, i.c., Eir 2 = 1, the low-rank tensor pikl can be regarded as the marginalization of a joint distribution over indices r, such that pe(ii,..., = E,
where pe(fi,...,$p) and ps(fi...., ii,,,,,.,. ry.) correspond to p!h and ///1, respectively. While the variable i directly corresponds to the features of the samples and appears explicitly as a tensor index of the model pla,, the discrete variable r is invisible in the model because it is marginalized out. In this sense, r can be regarded as a hidden or latent variable.
Scalability by exploiting non-zero entries in the a-divergence We highlight a key property of the a-divergence: it inherently ignores contributions from tensor entries that are equal to zero. Thus, we replace the summation over all indices Eien, in Equation ones Eieny where 227 = {i T + 0}. As we further explore in Section 3.4 this property justifics updating only the tensor elements (Pi)sen; and enables scalable inference when the tensor 'T is sparse, a common assumption in density estimation as T has at most N non-zero entries. This property does not hold in general for other divergences, including the 3-divergence (Basu et al. 1998 Eguehi & Komori 2022).
with summation over only the observed
5

---

Published in Transactions on Machine Learning Research (4/2026)
2 Tensor-based density estimation via a-divergence optimization
We develop the E M algorithm for discrete density estimation via a-divergence optimization. Our framework offers the flexibility to incorporate various low-rank structures, such as CP, Tucker, TT, their mixtures, and adaptive background terms. The key idea behind the E²M algorithm is as follows. (1) We first relax the a-divergence minimization to the KL divergence minimization using Jensen's inequality, (2) Then, we further relax the KL divergence optimization into many-body approximation, which can be solved exactly either by the convex optimization or by the dlseed-form solutions derived in Section 3.2
2.1 Two upper bounds of the a-divergence forming the E algorithm
For given D-th order normalized non-negative tensor T, we find mixture of tensors p(u,... ...., p[k] minimize the following objective function:
[[k]
K 1 log Σ P.P: where P: = Σ Hp!h and P!R
(P)
(3)
iESt
TE 22
where the function (P) ) is often referred to as the Rényi a-divergence (Rényi 1961 Póczos & Schneider 201 By the monotonicity of the logarithmic function, the above optimization is equivalent to the minimization of the a-divergence Da (T, P) given in Equation (I) for a E 10, Setting K = and a - 1 provides conventional KL-divergence basod low-rank tensor dccomposition. Each tensor pl is normalized, i.c., Eieo, p(k] = For simplicity, we introduce the weighted Qlki as 2 = 1 o and refer to the tuple of them (Q"), ..... (IK)) Q In the following, we denote the cross-entropy H between tensors A and B by H(A, B), that is, H)A, Eien, A log Bg We also refer to the non-negative vector (n!"),...,,))" defining the mixture weights 77. The symbol represents the element-wise product between tensors, and T denotes the element-wise exponentiation of T to the power of a. The difficulty in the optimization in Equation 3 arises from the two summations, E; and Et,r, inside the logarithmic function. We address this issuc using a double bounding strategy: the first bound moves E. outside the logarithm, and the socond bound moves Σer outside the logarithm as seen in the following proposition.
-.Ixx--IIpp let the set of tensors C and D be Proposition For any a 10 and non-negative tensor T E 2
c F Σ
ign, FTFt= =
plix---XIn 20
(4)
D = {(00))
l)k1 1E rE
! 1, (1
R200
1 xxxx
(5)
ir
respectiucly. Then, for any tensors F EC and +- = (T11) {bkk) E DD it holds that
L. < La P) < (F,
where the upper bounds are given as
o F,P) ha(F),
(6)
= ()
+ ha (F) + Σxe H(M HMM ((H)),
(7)
and the function J, the function ha, and each element of the tensor Ml E R as
1ix-xXlp x-. are given
J(n) = Σ
K
M! I1, ha(F) = H(T oFFF)
M ir = F.
(8)
1 1
iED
1= PER
respectively-

---

Published in Transactions on Machine Learning Research (4/2026)
Algorithm 1: ECM algorithm for non-negative tensor mixture learning
input : Tensor T, the number of mixtures K, ranks (R' ..... RK), and (0, 1]
Initialize Q! and 17() for all k;
repeat
Pi +- Σx "MMlp! where P!M - ErER 0
3
4
M
for all k;
1/ E2-step
@lk Many-body approximation of M Eqs. [44.
15 and 16
// M-step
6 Update mixture ratio 7 using Equation for all k;
// M-step
7 until Convergence;
8
return P
Proof. See Appendix A
In the standard EM algorithm, a bound is typically derived using Jensen's incquality (Jensen 1906 However, since the lensor T is not guaranteed be normalized, this approach is not directly applicable. The essence of the first bound is introducing a tensor to define an alternative normalized distribution W T and enabling the application of Jensen's inequality.
Notably, monotonic a-divergence optimization method has recently been developed based on bound derived using auxiliary function (Daudel et al., 2023b) exploring the property log(u") S u² - whereas we explore Jensen's 's inequality. In Appendix [.77. we derive the correspondence of our first bound and their bound by extending our above framework to the Gaussian Mixture Model considered in their work. We presently use the above double bound to relax the a-divergence based low-rank approximations into corresponding many-body approximations, which thereby enables closcd-form updates, as described below.
2.2 EaM-algorithm for a-divergence optimization of low-rank tensor structures
The proposed method monotonically decreases the objective function La in Equation 33) by iteratively performing the El-, E2-, and M-steps described below. Notably, the El- and E2-steps, can be jointly executed single E2-step, which enables memory-efficient updates by leveraging the sparsity of the data.
E1-step minimizes the upper bound (J,P) for the tensor F that is, F = arg minzec (FF, (FFP), where C is the set of tensors satisfying the condition Eieo, F as defined in Equation 44) This subproblem is equivalent to maximizing the expectation Ew (PI-a/ for W = = oF, which is detailed in Appendix B.1. As shown in Proposition 2 in Appendix A the optimal solution is
Pl-a
F
(9)
Σ Eien, T!P!!-
By simply putting the above optimal tensor F into Equation 6 it can be verified that the tensor F
makes
the bound tight, i.e., (P) = La (J=, P), which is also ratenstraaed in the proof of Proposition 2 Thus, we can optimize the objective function La by iteratively updating F and P to minimize the bound L, (F, (F,P). P). However, the naive update of P to minimize Lo requires non-convex optimization due to the summation Σx,r in the bound I.p Thus, we focus on the second bound T., in the following E2-step and M-step. We highlight that the bound (F, successfully relaxes the a-divergence optimization into the KL divergence optimization since it holds that arg minp I., (F. argmin, (TL F.P)). (ol..... (IKI)., i.e., $* arg m [. &b, ②), where the set D is defined in Equation 5 For the tensor M defined in Equa- tion 88, this subproblem is equivalent to maximizing the expectation EM (②/b)), which is detailed in
E2-step minimizes the bound (F, Q) for tensors
7

---

Published in Transactions on Machine Learning Research (4/2026)
Appendix B.1 As shown in Proposition 3 in Appendix A
the optimal solution ((A) is
6!H Ek=1 Ererm
(10)
By putting the above optimal d into Equation also demonstrated in the proof of Proposition 3 Thus, we can optimize the upper bound Taa by iteratively
(7 it can be shown that (=, = P), which is
updating & and Ô to minimize the bound La.
M-step minimizes the upper bound (F, Φ Q) for tensors Q = (Q @(k]) and non-negative weights 7 (r)").... ,,....,,,,()),,, We recall that the weighted 2 defined as 2 n D Since the bound can be dccoupled for each k as shown in Equation 7 the required optimizations in the M-step are follows:
arg min H (M )
arg min J(n).
(11)
0gn(t)g1.> 1/1=1
where the model space b is defined as the set of tensors Qlk) that can be represented as a form in Equation and whose sum over the index r yields the low-rank structure k; i.e.,
B - {OHH
off P!
Oikl E R2
/xxxxxxxxR..... - RIxxxRx Ryp
In the above optimization for the tensor @ in Equation 11; the logarithmic term no longer contains E; or Ex. making the optimization much easier compared the original problem for Lo. In fact, this optimization is known the tensor many-body approximation (as detailed in Section 3), which is convex optimization problem with guaranteed global optimum. Furthermore, since the lensor O(Ai can be given in product form seen in Equation (2), the properties of the logarithm allow the objective to decouple each factor, resulting in all factors being simultaneously optimized in closed-form providing a globally optimal solution for the M-step in many cases. The closed-form update rules for & E {CP, Tucker,TTI are introduced in Section We also provide the optimal update rule for the mixture ratio ny closed-form as follows:
EvER, Ererpt M! ir
(12)
Ek=E Eieos Erea
M 8
ir
which is derived in Proposition 4 in Appendix A
To avoid memory inefficiency by handling the dense Lensors and the above El and E2-steps, we substitute the updated expressions from Equations 9 and 10) into the definition of the tensor M Equation i.e.,,, MM! = F. T,
Σ, T. P.: step. Thus, instead of F and Φ, we only store the tensor M, which is proportional to T, resulting in significant. memory savings when T is sparse. The complete algorithm is presented in Algorithm [ The normalization is guaranteed at each iteration since the updated tensor Qihl] E BBk]] satisfies Eien, EreR e!4 1, which
We refer to this integrated E-step as an E
immediately implies Eien, P: = Eien, E. (pllt = =x Etea, ErER 1 =
When a + 1, it holds that F3 = 1 in Equation (9 and consequently, the entire algorithm reduces to the standard EM algorithm since minimizing the upper bound L in the M-step is equivalent to maximizing the lower bound of the likeliood.o...
Convergence guarantee of the E algorithm Notably, both the E1 and E2-steps tighten the upper bounds as shown in Figure [b), which makes the objective function not increase in each iteration and converge regardless of the low-rank structure, as formally proven in Theorem [[] The statement in the following theorem consistent with existing results on convergence of the objective function in the general theory of the EM-algorithm and Najoriattior-Miinnizzaooo
Minimization methods (Jeff Wu 1983 Sun et al. 2016).
Theorem 1. The objective function of the mixture of tensor factorizations using the E² M M-procedure always connerges regariless of the choice of low-rank structure and mixtures.
8

---

Published in Transactions on Machine Learning Research (4/2026)
Proof. See Appendix
A
We note that Theorem establishes only the convergence of the objective function sequence; it does not assert convergence of the iterates (the tensor P) nor convergence to a stationary point. These stronger forms of convergence are left for future work and are expected to require additional assumptions on the objective function.
Adaptive background term for stable learning Our approach deals not only with typical low-rank structures and their mixtures but also with learnable uniform background term. We here define model P the mixture of low-rank tensor ppllow-rank] and a normalized uniform tensor plbg, whose elements all are 1/|0/] as P = (1 (pbgllppllowwee + pllbs, for [Sty] = 1122....... Ip. The learnable parameter rylbal indicates the magnit ude of global background density in the data. It is a well-known heuristic that learning can be stabilized by adding a small constant to the tensor (Cichocki & Phan, 2009) Gillis & Glineur 2012). Our mixture modcling with such background density enables this constant to be learned from data instead of being manually tuned, promoting stable learning. We note that we do not need any modification of the algorithm for the background term since all elements in p are fixed, and we just update its ratio 17 1718
other mixture ratios by Equation 12
In the following, we denote the proposed methods as B + structures, where 'B' indicates an adaptive background term (e.g., E*MCPTTB is a mixture of CP, TT, and background). We also denote the proposed methods with fixed a 1.0 as EM structures, e-g-, EMTTB is mixture of TT and background optimizing the KL divergencc.
1D Gaussian Mixture (2 components)
0.5
Estimated dist. (a=0.9)
Estimated dist. (a-D.9)
Estimated dis., (a-0.7)
Estimated dist. (a=0.8)
Density
Empiniaa dist.
0.3
2.3 How to choose the a value?
0.2
0.1
As we demonstrate in Section 4 reconstruction becomes less sen- sitive to outliers, noise, and mislabeling for smaller a, and more sensitive for large a. This observation is also supported by the results of fitting Gaussian distributions using the a-divergence, as shown in Figure 2 for smaller a, the fitting emphasizes the dominant modes while largely ignoring minor modes. Thus, we suggest using a smaller when data contamination is expected, and a larger otherwise. If no prior knowledge about contam- ination is available, a can be tuned to maximize the validation score. The continuous extension and the setup for the E algorithu for the Gaussian mixture model fitting
-4
-2
2
4
6
Figure 2: Fitting Gaussian mixture model by E*M algorithm with various values. Lower a values make the model fit the dominant modes well, providing robustness against data contamination.
in Figure 2 are provided in Appendix B.7]
3 Many-body approximation meets the E'M algorithm
The naive EM-based formulation typically relies on an iterative gradient method during the M-step (Lange 1995 Wang et al. 2015). Importantly, we point out that the M-step in Equation (II) can be regarded many-body approximation. We derive closed-form exact solutions of the many-body approximation for various low-rank structures such CP, Tucker, and TT. Consequently, the gradient-based M-step for these low-rank structures and their mixtures is no longer required. In the following, we review the many-body approximation in Section 3.1 We then derive closed-fore form solutions for the M-step in Section 3.2 and sketch
their derivations in Section 3.3 Finally, with the closed-fon form updates in hand, we discuss the computational
complexity of the proposed E M algorithm in Section 3.4
3.1 Review of the body-body approximation
Low-rank tensor dccompositions, such as Tijes ~ Qijut = E. based on the CP structure, have had great success in various applications. However, the rank-based modeling causes the summation Σx to appear inside the logarithm in the KI-divergence objective DRI(T, 2), resulting in a non-convex
9

---

Published in Transactions on Machine Learning Research (4/2026)
a. b.
C.
d.
TD
k
470
!!D
A
A
...
A
9
g
1D
Lp
2
Figure 3: Interaction diagram for (a) Qiskt = Aijk CapCa (b) [
A
fir Apr (c) o...ckerr]
(1) 4...A (orD' and (d) o = 9 ... G .. The summation over the indices recovers the CP, G,A in
Tucker, and TT structures, respectively, as p [k| Er in the factorized representation in Equation
!h! and each interaction,
corresponds to a factor
optimization problem. The many-body approximation provides an alternative tensor modeling that avoids sunmation the logarithun function to keep the optimization convex. In particular, it represents a tensor a purely element-wist multiplicative form, for example, 2iyu a,bycadi, Qyyut CaDDD and Qijki = Zoyllis The factor tensors, i.c., the lower-order tensors the right-hand side of the above examples, are interpreted as interactions among tensor modes. For example, AG represents the two-body interaction between the indices i and j, and Kiyk represents the three-body interaction among the indices i, j, and k. Rather than tensor rank, this modeling requires specifying an interaction by an interaction diagram, which visualizes the presence or absence of interactions among tensor modes by a factor graph where cach node corresponds to a tensor mode and interacting modes are connected via an edge labeled
We provide here an example of the diagram for a fourth-order tensor in Figure 3a). For given tensor M, the approximation corresponding to the diagram can be written Mijui Likki A. where matrices B, C, and D define two-body interactions, whereas the tensor A defines a three-body interaction. The many-body approximation always provides globally optimal non-negative normalized tensor Q that minimizes the cross-enlropy from the lensor M as described in Section 3.2
3.2 Many-body approximation with exact solution
To discuss the closed-form updates of the M-step in the E algorithm, we here consider the following three kinds of many-body approximations for normalized non-negative tensors, namely a (D + 1)-th order tensor and a (2D - 1)-th order tensor Their interactions are illustrated in Figure b), (c), and (d) respectively, and these tensors can therefore be factorized as shown in Equation ( In the following, the symbol Q with upper indices refers to the index set for all indices other than the upper indices, e.g.:
QiCP| E R x
2D-th order tensor [[Tucker] E 1 xx,.xxxpxxxxx..xxx xxxx-xxxxXx x
@ E 12
[LD] R.raira- [R,] .... x [Ra2]] x (ra-1) {ra} [Re+a] x x [Ry].
2id [41] x [[--1] x {id) [La+1]
2
Ad,d-1
From the general theory presented in {CP, Tucker, TT), the many-body approximation finds the globally optimal tensor Q! that minimizes the cross-entropy,
Sugiyama al.
2019 for given tensor M and k E
H(M, = Σ Σ Mer log 0/4 iESI rEfR
(13)
assuming the normalization Eier, Erern 0 1. This optimization is guarantced to be convex problem regardless of the choice of interaction. It has been shown in (Huang & Sidiropoulos, 2017 Yeredor & Haardt 2019) that the following matrices globally maximize H(M, Q(CPI))
Eienya Mir
(14)
1/D (Ere,, Mar) 1-1/D
10

---

Published in Transactions on Machine Learning Research (4/2026)
a.
b.
(1)
(D)
8(1)
(D-1) g (D)
A00
A (2)
A
9
9
Figure 4: The sketch of absorbing the scaling redundancy of the Tucker (a) and the TT structure (b). Vectors a!d) and g( correspond to the scaler to make the factors normalized and the root factor (bold boxed) absorb them. The definition of the scaler g(d) can be found in Equation (29) in Appendix
for j Σ tEn, Erean Mir: which is consistent with mcan-field approximation Ghalaamkari & Sugiyama 2021: 2023): As a generalization of the above result, presently provide the optimal solution of the
many- many-body approximation of Q TTuckerr and Q in closed-form. The following tensors globally maximizes H(M, Q [ucker))):
Eienpy Erente Mir (d) Mara = Eier, Ereah Mir
Eieo, Mir
Gr Σieo, Ererg Mir
(15)
and the following tensors globally maximizes H(M, QITTI):
Eienyy Ereage Riraves
Mir
(16)
Eien, Erenga Mir Rira
We note that both the denominators and numerators are real values in Equations 14] 15 and [6]. We formally derive the above closod-form solutions Propositions 5 and Appendix Although we do not guarantee that the many-body approximation admits a closed-form solution for arbitrary interactions as the decoupling of the normalizing condition is nontrivial, more complicated many-body approximations by combining these solutions, discussed in Appendix [.5]
obtain closed-form solutions even for
Above Equations 14 16 obviously provide the closed-form update for the M-step in Equation 111 where
the model space B K is defined as the set of tensors with interaction described in Figure 3 ()) (d),
accordingly to the low-rank structure k.
We also note that, although the original many-body approximation assumes that both the input tensor M and the reconstructed tensor Ql6i are normalized, the above discussion, the normalization constraint is imposed only on (!A), This is because the tensor M that appeared the M-step is not guaranteed satisfy the normalization condition, i.e., Eir Mit 1, when the number of mixture components K > 1. Therefore, the above discussion with K leads to the connection of the M-step to an extended many-body approximation that does not rely on the normalization the input tensor.
3.3 Sketch of the proof of closed-form updates
A we outline the proof of the closed-form updates in Although the formal derivation is provided in Appendix
Equations [5] and [16] here. As mentioned in Section 2 the bound L,, has summation inside the logarithm function, which enables us to decouple the objective function into independent factors leveraging the basic property of the logarithm function, i.c., logzy = log x + log However, the normalizing condition is exposed on each low-rank reconstruction Pik) but not on each factor. This renders the decoupling of the Lagrange function nontrivial. To address this difficulty, as sketched in Figure 4 we eliminate the scale redundancy of low-rank tensors p(k), i.c., the invariance under scaling one factor by v and another by 1/v, by transforming the factors and absorbing the normalization constraint into single factor, referred to as the rool. factor.
11

---

Published in Transactions on Machine Learning Research (4/2026)
For the M-step of the Tucker decomposition particular example, the Lagrange function with the normalization constraint is given as
-- E- Mär log Gr ... A - ( BEDI PERR
GnA ... A -1)
ΣΣ iE rella
(17)
However, the closed-form solution satisfying 0C/04 al/09, 0 is non-trivial since the normalization
constraint is not decoupled yet. Then, we define the rescaled factors as A / /ard where (d) Eia Aigras and we rewrite the Lagrange function in Equation 17
D (a)
== E E RED, rEia
Mir log 4
A - a > --- -1-
Σ id
1 kark -1)
(18)
TM
We note that the core tensor absorbs the normalizers a1)).... ........e (D)) Consoquently, we also define the trans-
formed core lensor 9r- = Grar: .rr where the normalization Σien, pf [TTucker) = requires Ererg G. = 1.
Since the normalization condition is also decoupled into that of each factor in Equation (18 the critical
condition /C/0A TTT ac/0û, 0 leads to the forem update rules. Appendices A and B further illustrate these decoupling techniques for TT and Tensor Tree structures.
3.4 Scalability of the E'M-algorithm for tensor learning
Given closed-form update rules, we here discuss the computational complexity of the E algorithm. Since the tensor to be optimized in the M-step is proportional to the input tensor T, i.c., Mi!i ir ir: we can leverage the sparsity of to make the proposed framework scalable with respect the dimension D of the data, except in the case of MTTcker,,, which involves the dense D-th order core tensor G. Specifically, in Equations
Low-rank model Complexity
EºMCP
O(DNR)
ErMTucker Е?
O(DNRD)
O(DNR²)
Table 2: Computational complexity of each iteration for the tensor dimension D, number of samples N, and rank R.
[14] T6) we can replace Σien, with Σien where 92 is the set of indices of nonzero values the tensor T. The complexity of the E2-step is O(V[9R= since this step involves marginalizing the resulting tensor @ik) in the M-step, i... P. Er O and updating M!k ir for all indices ir E 22 x 2px. The resulting computational time complexity is ODyDN for E2MCP and O(>DNRD) for MTMTucker and E?MTT, respectively, where y is the number of iteration, assuming ranks are (R,..., R) for all low-rank models. Furthermore, we reduce the complexity of E2MTT to (() by merging cores from both edge sides to obtain TT format, detailed in Appendix B.2 These complexities summarized in Table 2 The complexity for K > is simply the sum of the complexity of each low-rank structure in the mixture. We note that when applied to data with many samples, i.c., a very large N, the proposed method can be scaled using a batched FM approach (Cappé & Moulines 2009 Chen et al. 2018). Notably, the overall computational complexity of the E procedure equivalent to standard gradient- based mothods that rely on the chain rule, as we demonstrate in Appendix B.3
4 Results and Discussion
To verify the practical benefits of our framework, we conduct four experiments that evaluate: i) optimization performance, ii) lcarning stability offered by the background term, iii) robustness to outliers controlled by the a parameter, and iv) generalization in density estimation lasks. While the proposed framework can explore a variety of low-rank structures and their mixtures, we here mainly focus on the performance of E*MCPTTB due scalable propertics. The experimental results for other low-rank structures as ablation studics are available in Appendix C We also provide the experimental setup, dataset description, bascline selection, implementation details, and hyperparameter tuning in Appendix [)
12

---

Published in Transactions on Machine Learning Research (4/2026)
a.
o:: 0.7
b.
16-5
QC 0.7
d: 0.7
d.
10
Neggetiv Iogs- Baegg-BBii Booled
DMCRS PMTTR
Adagrad
Adagrad
EMMPPTB
RMSppop
5
RMSSprop
10
EMCTTTB
S60
S6D
Adam
MCPTTBB
Change in Alpha-diveess Allha-diverrgecc
Adaam
(3MCPTTT
10
0:3
Alpha-divergence
10
3
Change in Alpha-divergent
8.0
10
500
1000
1500
2000
1-
10
Nogadee thodd
-00
0-
10
10
20
155
-2
-0
10
10
10a
10
10
10a
10a
10
500
1000 3500 2000
Iteration
Iteration
Iteration
# parameters
Figure 5: The objective function in each iteration (a) and its iteration- change (b) (c) to reconstruct color image (i.e., SIPI house 4). The rank set to R[CP] 35 and R(TT] (8,8) so that the mumber of model parameters is approximately half the size of the input data. The weight decay of the baseline methods is set to the default value. In the above, (a), (b), and (c) are plotted using all results obtained from five random initializations. (d) Density estimation for the CarEvaluation data, varying the number of model parameters with (lop) and without (bottom) adaptive background modeling, respectively. The results show the mean and standard deviation over five random initializations.
i) Performance on optimization To evaluate the optimization performance of the E algorithm, we approximate the normalized SIP1 House image as a tensor T E 12 using mixture of CP, TT, and the adaptive background term. Although the SIPI image is not categorical tabular dataset, we intentionally used it to evaluate the potential usefulness of the proposed method beyond density estimation. With fixed at 0.7 and five random initializations, we observe the values of the objective function at each iteration in Figure 5 a) and its iteration-wise change in Figure 5(b) b) and (c), where (c) shows a log-log plot. These results are reported with the gradient-based baselines SGD Robbins & Monro, 1951] D, RMSprop (Tieleman & Hinton 2012 Adagrad (Duchi et al., 2011), and Adam (Kingma & Ba 2015 which also allow for a-divergence minimization of mixtures of low-rank structures as discussed in Table 1 To ensure that these bascline methods satisfy the non-negativity and normalization conditions, we used softmax transformation as detailed in Appendix [D..1]] Learning rates for the baseline methods are tuned to achieve the lowest a-divergence. The proposed method consistently decreases the objective function monotonically, which is consistent with our theoretical result in Theorem and achieves costs comparable to those of the baselines.
We include only the proposed method Figure 5 c), since the baseline methods do not guarantee a monotonic decrease and may produce a negative change in the a-divergence objective, which cannot be plotted on a log scale. As observed in the full results in Figure I0 in Appendix C which vary the learning rates of the baselines and the values of a; gradient-based methods often oscillate or fail to converge when the learning rate is poorly tuned, resulting in larger variance across random initializations when evaluating the objective after sulficient iterations. The E algorithm requires no learning-rale tuning, guarantees monotonic decrease, and reliably converges to solutions of comparable quality to those obtained by gradient-baser methods.
ii) Effect of the adaptive background term We verify that the adaptive background term stabilizes learning considering density estimation for the CarEvaluation dataset 2 Specifically, we construct an empirical tensor T from the training samples in the dataset, factorize it by the proposed methods with and without the adaptive background term, and evaluate the reconstruction using negative log-likelihood for validation samples, varying tensor ranks. As demonstrated in Figure Bid), the validation negative log- ligeellldoood often becomes large when the model has no adaptive background term, caused by numerical instability due to
usc. http:////ssp..ucceee//tttttataet https://archive. ics. wee...edd/daaaae adu/dataset/19/
13

---

Published in Transactions on Machine Learning Research (4/2026)
KL divergence
Hellinger distance
Reverse KL
divergence
Jeffreys divergence
7:00
Tist emer Thaining
0.60-
7.00
6.00
800
6.00
0:50-
6.00
Traning Tist erro
Training eror
5.00
Number noise
Treining fest
4.00-
5.00
Trainin: Tust
0:45
0.50
5.00
0.20
0.00
600
0.40
Traingg
4.00
0.15-
3.00
3.00
2.00
0.35
2.00-
0.30
2.00
200
0.30
1.00-
0.25
1.00
0.03
1.00
0.00
0.00
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
100
0.0
0.5
1.0
0.0
0.5
1.0
Figure 7: Reconstruction error of noisy half-moon data using E2MCPTTB, varying the noise level and a. The optimal value of a depends on the noise level; for example, in the case of the Jeffreys divergence, the test error is minimized around a = 0.3 under high noise, whereas a = 0.7 becomes optimal in the low-noise regime.
extremely small values in the logarithm function. The background term eliminates this problem and prevents the rapid increase in validation error when the number of model parameters is large.
iii) Mass covering property induced by the a parameter half-moon dataset, discretize the v) coordinates, and generate a tensor T E 1 The first and second modes correspond to the discretized coordinates (I, whereas the third mode corresponds to the class label. The 25 randomly selected class labels are flipped, and a further 25 random outliers are added to T and obtain Tuoiny, Further details regarding the data synthetization can be found in Section D.1.2 in the Appendix. The reconstruction P obtained by E²MCPTTB is shown in Figure 6 for a = 0.1 and a 0.9, where red indicates higher probability of class and blue indicates higher probability of class 2 at each coordinate. See Appendix D.1.2 for details of the visualization. Full results, varying and the number of noisy samples, can be found in Figure [3 in Appendix C The CP-rank and TT-rank are set as RCPP] = 20 and R[TT] = (20, 2), respectively. outliers and mislabeled samples when a is larger, while it more accurately estimates the true distribution when a is small. Specifically, in Figure G the mislabeled samples in the yellow box and outliers in the dashed cyan boxes adversely affect the reconstruction when a 0.9, resulting in reduced confidence in the yellow box and noisy prediction in the dashed cyan boxes. In contrast, the reconstruction with a = 0.1 is robust to outliers mislabeled samples. As a result, the proposed E method leverages the a-divergence advantage of controlling the sensitivity to outliers through suitable tuning of a. The above experimental results are consistent with existing works that demonstrate the robustness of the a-divergence
We sample 5000 points from the
We observe that the reconstruction is affected by
against noise and outliers (Cichocki & Amari
2010
Cichocki et al. 2011 Iqbal & Seghouane 2019
Noisy empirical dist.
0:0.1
a:0.9
Furthermore, Figure 7 shows how the reconstruction quality changes quantitatively as a varies. Specif- ically, evaluated how well the mixed low-rank approximation P of the input distribut ion approxi- mates the noisy imput Tuoisy (training error) and the noiseless data T (test error) using five measures: KL divergence (a = 1.0), Hellinger distance (a = 0.5), re-
Figure 6: Reconstructions of the noisy empirical distri-
verse KT. divergence (a- = 0.0), Consen-Shannor (JS)
divergence (Lin 1991p. and Jeffreys divergence (Jef- bution by E2MCPTTB with 0.1 and 0.9.
freys. 1946). The definitions of JS divergence and Jeffreys divergence can be found in Appendix D.1.2 We observe the following: 1) The optimal value of a that minimizes the test error depends on the noise level, indicating that noisc robustness can be achieved by appropriately selecting a. 2) Optimizing the KL divergence does not necessarily lead to optimal performance under other divergence measures. 3) The test error sometimes decreases as the noise level increases, pointing to the noise acting as a regularization that makes the models less prone to overfitting to the training data. As such, this phenomenon was not observed when considering a model with low-rank and therefore less prone to overfitting (see Figure 14a in Appendix Cp.
The similar results with noisy real dataset be found in Figure 14b in Appendix C
14

---

Published in Transactions on Machine Learning Research (4/2026)
Table 3: Accuracy of the classification task (top lines, higher is better) and negative log- oiglliielihood (bottom lines, lower is better) for the first datasets (in alphabetical order) out of a total of 34. The full results of all 34 datasets are available in Table c in Appendix. The bottom row reports the mean accuracy across all 34 datascts. Error bars are given standard deviation of the mean across 5 randomly initialized runs and for the bottom row the 34 datasets.
Dataset
BM
KLTT
CNMF
EMCP
ЕĄМCРТТB
Optimal a
0.571 (0.00)
0.589 (0.02)
0.571 (0.00)
0.571 (0.00)
0.500 (0.00)
0.90
AsiaLung
2.240 (0.00)
2.295 (0.06)
3.130 (0.09)
2.262 (0.00)
2.480 (0.16)
0.95
0.777 (0.00)
0.830 (0.00)
0.649 (0.07)
0.868 (0.01)
0.872 (0.00)
0.90
B.Scale
7.317 (0.03)
7.480 (0.15)
7.418 (0.02)
8.979 (0.49)
7.102 (0.00)
1.00
0.901 (0.01)
0.923 (0.01)
0.956 (0.00)
0.912 (0.00)
0.934 (0.00)
1.00
BCW
13.120 (0.12)
12.63 (0.004)
12.776 (0.15)
12.815 (0.00)
12.623 (0.00)
1.00
0.847 (0.04)
0.892 (0.00)
0.736 (0.01)
0.942 (0.02)
0.950 (0.02)
1.00
CarEval.
7.964 (0.15)
7.752 (0.02)
8.186 (0.06)
7.777 (0.04)
7.847 (0.07)
1.00
0.500 (0.00)
0.500 (0.00)
0.500 (0.00)
0.500 (0.00)
0.500
(0.00)
0.75
Chess
12.45 (0.30)
12.45 (0.17)
14.739 (0.04)
10.306 (0.05)
10.886 (0.01)
0.95
0.157 (0.00)
0.277 (0.01)
0.243 (0.01)
0.439 (0.01)
0.573
(0.00)
0.95
Chess2
13.134 (0.00)
12.561 (0.12)
12.772 (0.03)
11.948 (0.08)
11.637 (0.00)
0.95
0.537 (0.02)
0.585 (0.07)
0.573 (0.04)
0.561 (0.00)
0.585 (0.00)
0.95
Cleveland
6.379 (0.04)
6.353 (0.27)
6.144 (0.03)
6.139 (0.00)
6.235 (0.00)
0.95
0.864 (0.02)
0.909 (0.00)
0.909 (0.00)
0.909
(0.00)
0.909 (0.00)
1.00
ConfAd
6.919 (0.04)
6.979 (0.01)
6.918 (0.03)
6.846 (0.00)
6.859 (0.00)
1.00
Mean Acc.
0.708 (0.04) 0.698 (0.04) 0.649 (0.04)
0.743 (0.04) 0.776 (0.04)
iv) Density estimation on real datasets We use 34 real-world discrete datasets obtained from the repositories listed in Table D in Appendix D to perform density estimation. We emphasize that our experiments include large-scale tensors such as the Chess and Mushroom datasets, whose numbers of modes are 36 and 22, respectively, and total numbers of elements are on the order of one trillion. We construct an empirical tensor T from the training samples, factorize it, and evaluate the reconstruction using negative log-likelihood and classification accuracy. The hyperparameters are tuned to maximize the validation score. Further details are provided in Appendix D We compared our approach to the four lensor-based baselines: pairwise marginalized
method (CNMF) [lbrahim & Fu, 2021), KL-based TT (KLTT) (Glasser et al., 2019), EMCP (Huang &
Sidiropoulos 2017 Yeredor & Haardt 2019), and Born machine (BM) (Fan et al. 2008). See Appendix D.1.4
for details of the baseline selection. Due to the page limitation, only results from the first eight datasets (sorted alphabetically) are displayed in terms of the accuracy of classification and negative log-likelihodd of the estimated density in Table 3 The full results on all 34 datasets are available in Table G in Appendix C whereas the mcan accuracy across all 34 datasets is displayed at the bottom of Table 3 Overall, E²MCPTTB outperforms single low-rank models. This improvement cannot be attributed solely to the expressiveness of TT, as ECMCTT also outperforms KLTT and E²MTTB as seen in Table 7 in Appendix C Finally, we conduct Wilcoxon signed-rank tests to evaluate whether our proposed method, E2MCPTTB, statistically outperforms the baseline methods (BM, KLTT, CNMF, and EMCP) in terms of classification accuracy. The resulting p-vauues after the Holm correction for multiple comparisons are (0.0017, 0.0003, 0.0003, 0.0299). We also provide experimental results comparing non-tensor-based methods in Table in Appendix C
5 Conclusion
We establish the E algorithm, a general franework for tensor- based discrete density estimation via a- divergence optimization. Importantly, our double bounding of the divergence relaxes the optimization problem
15

---

Published in Transactions on Machine Learning Research (4/2026)
into tensor many-body approximation, thereby eliminating the traditional difficulty raised by a-power terms in the objective function. Our framework not only uniformly handles various low-rank structures but also supports mixtures of low-rank structures and an adaptive background term without losing convergence guarantecs. The essence of the E algorithm lies removing all summations (c.g., Ex Ex. E,) from inside the logarithmic function thereby leveraging the basic properly log(ry) = log(a)) + log(y), which decouples the optimization problem into solvable independent parts. This approach guarantecs closed-form updates for all steps whenever the standard EM algorithm admits a closed-form M-step. Although this work focuses on discrete density estimation
simple example of non-negative tensor factorization, tensor factorization also be interpreted as a probabilistic circuit. (Loconte et al. 2025 graphical model (Robeva & Seigal 2018), hidden Markov model (Glasser et al. 2019). Therefore, we can naturally expect fut-ure applications of the E algorithm in these related fields, in addition to existing applications of non-negative lensor faclorizations
beyond density estimation (Nickel et al. 2011 Panagakis et al. 2021 Cichocki et al. 2016
Limitations
While all steps in E algorithm are guaranteed to be globally optimal, the entire procedure remains non-convex, in the standard EM procedure. For large N, the E algorithm requires mini-batched learning techniques. The thcoretical analysis supporting the generalization performance of the mixture models remains for future work. The number of ranks that need to be tuned is larger in mixture models than in the non-mixture low-rank model.
Ethical Statement
The aim of the developed E2M algorithm is to advance the use of low-rank decompositions, exploring properties of the a-divergence. As for all machine learning technologies, there are many societal concerns about the applications of the developed algorithm. This includes the misuse of the developed methodology for unethical purposes, sensitive data, and applications that raise privacy concerns. Furthermore, low-rank decompositions may exacerbate biases and represent aspects in terms of stereotypical structures, missing important nuances. As the associated optimization problems are non-convex, the results and conclusions drawn can potentially change based on the obtained solutions. Consequently, when making interpretations of the model representations, the reliability has to be systematically assessed, and care has to be taken not. to over-interpret model results.
Acknowledgments
This work was supported by Danish Data Science Academy, which is funded by the Novo Nordisk Foundation, Grant Number NNF21SA0069429 (KG), NII Open Collaborative Research 2025 Grant Number 24FP07 (KG), JST, CREST Grant Number JPMJCR22D3, Japan (KG). JLH and MM were supported by the Independent. Research Fund Denmark, Grant Number 10. 00.46540/2035500220 Nordisk Foundation, Grant Number NNF230C0083224.
MM was further supported by the Novo
References
S. Amari. Differential-geometriia methods in statistics, volume 28. Springer Science & Business Media, 2012.
M. Amiridi, N. Kargas, and N. D. Sidiropoulos. Low-rank characteristic tensor density estimation part Foundations. IEEE Transactions on Signal Processing, 70:2654 226654-668,, 2022.
H. Bao, S. Sakane, and Y. Takezawa. sty-stepsize gradient descent for separable data under Tenchel-Yourd losses. arXiv preprint or Xrv::2200 04889, 2025.
A. Basu, I. R. Harris, N. L. Hjort, and M. C. Jones. Robust and efficient estimation by minimising a density power divergence. Biometrika, 85(3):549 559, 09 1998.
L. E. Baum and T. Petrie. Statistical Inference for Probabilistic Functions of Finite State Markov Chains. The Annals of Mathematical Statistics, 37(6):1554 1563, 1966.
16

---

Published in Transactions on Machine Learning Research (4/2026)
C. M. Bishop. Pattern Recogrition and Machine Learning (Information Science and Statistics). Springer- Verlag, Berlin, Heidelberg, 2006. ISBN 0387310738.
R. E. Bradley, S. J. Petrilli, and C. E. Sandifer. L L'Hópital's Analyse des infiniments petits: An Annotated Translation with Source Material by Johann Bernoulli, volume 50 of Science Networks. Historical Studies. Springer International Publishing, 2015. ISBN 978-3-319-17115-9.
S. Bubeck. Convex Optimization: Algorithms and Complexity. Found. Trends Mach. Learn., 8(3 (33-4)231-335) 357, November 2015. ISSN 1935-8237.
L. Cai, Y. Chen, N. Cai, W. Cheng, and H. Wang. Utilizing amari-alpha divergence stabilize the training of generative adversarial networks. Entropy, 22(4):410, 2020.
O. Cappé and E. Moulines. On-line expectation maximization algorithm for latent data models. Journal of the Royel Statistical Society Series B: Statistical Methodology, 71(3):593-611 2009.
J. D. Carroll and J-J Chang. Analysis of individual differences in multidimensional scaling via an n-way genoralization of "Eckart- Young" docomposition. Psychometrika, (5(3):283-311 1970.
J. Chege, M. J. Grasis, A. Manina, Yeredor, and M. Haardt. Efficient probability mass function estimation from partially observed data. In 2022 56th Asilomar Conference on Signals, Systems, and Computers, pp. 256-262. TEEE, 2022.
J. Chen, J. Zhu, Y. W. Tch, and T. Zhang. Stochastic expectation maximization with variance reduction. Advances in Newral Information Processing Systems, 31, 2018.
S. Cheng, L. Wang, T. Xiang, and P. Zhang. Tree tensor networks for generative modeling. Physical Review B, 99:155131, Apr 2019. doi: 10.1103/PMM (PhyseeeB.99.1500
C. Chow and C. Liu. Approximating discrete probability distributions with dependence trees. TEEE Transactions on Information Theory, 14(3):462 467, 1968.
A. Cichocki and S. Amari. Familics of alpha- beta- and gamma- divergences: Flexible and robust measures of similarities. Entropy, 12(6):1532-1168) 2010. ISSN 1099-4300.
A. Cichocki and A. Phan. Fast local algorithms for large scale nonnegative matrix and tensor factorizations. IEICE Transactions Pundamentals of Electronics, communications and computer sciences, 92(3):708 721, 2009.
A. Cichocki, S. Cruces, and S. Amari. Generalized alpha-beta divergences and their application to robust nonnegative matrix factorization. Entropy, 13(1):134-170, 2011. ISSN 1099-4300.
Cichocki, N. Lee, I. Oseledets, A. Phan, Q Zhao, and D. P. Mandic, Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions. Foundations and Trends@ in Machine Learning, (4-5)::249 429, 2016.
R. Cichocki, A. amd Zdunek, A. H. Phan, and S. Amari. Nonnegative Matrix and Tensor Factorizations: Applications to Exploratory Multi- woy Data Analysis and Blind Source Separation, volume 194. John Wiley & Sons, 2009.
I. Csiszár. On information-type measure of difference of probability distributions and indirect observations. Studia Sci. Math. Hungar., 2:299-31 318, 1967.
Daudel and R. Douc. Mixture weights optimisation for alpha-divergence variational inference. Advances Neural Information Processing Systems, 34:4397-4408, 2021.
Daudel, R. Douc, and Portier. Aufnite-dimensionn gradient- based descent for alpha-divergence mininisation. The Annals of Statistics, 49(4):2250-222) -2700, 2021.
17

---

Published in Transactions on Machine Learning Research (4/2026)
K. Daudel, J. Benton, Y. Shi, and A. Doucet. Alegaeeivergence variational inference meets importance weighted auto-encoders: Methodology and asymptotios. Journal of Machine Learning Research, 24(243): 1 138 2023a.
K. Daudel, R. Douc, and F. Roueff. Monotonic alpha-divergence minimisation for variational inference. Journal of Machine Learning Research, 24(62):1-76, 2023b.
A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, 39(1):1-22 1977.
J. Deng, W. Dong, R. Socher, Li-Jia Li, L. Kai, and F. Li. Imagenet: A large-scale hierarchical image database. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pp. 248 255, 2009.
L. Dinh, J. Sohl-Dickstein, and S. Bengio. Density estimation using real NVP. In International Conference on Learning Representations, 2017.
S. Dolgov, K. Anaya-Izquirdo,, C. Fox, and R. Scheichl. Approximation and sampling of multivariate probability distributions in the tensor train decomposition. Statistics and Computing, 30:603-624 625, 2020.
J. Duchi, E. Hazan, and Y. Singer. Adaptive subgradient methods for online learning and stochastic oplimization. Journal of Machine Learning Research, 12(7), 2011.
C. Durkan, A. Bekasov, I. Murray, and G. Papamakarios. Neural spline flows. In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
S. Eguchi and O. Komori. Minimum Divergence Methods in Statisticol Machine Leurnóngg: From on Information Geometric Viewpoint. Springer, 2022.
R. E. Fan, Kai-Wei Chang, Cho-Jui Hsich, Xiang- Rui Wang, and Chih-Jen Lin. LIBLINEAR: library for large linear classification. the Journal of Machine Learning Research, 9:1871-1874 2008
M. Germain, K. Gregor, I. Murray, and H. Larochelle. Made: Masked autoencoder for distribution estimation. Tn Francis Bach and David Blei (eds.), Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 881 889, Lille, France, 07 09 Jul 2015. PMLR.
K. Ghalamkari and M. Sugiyama. Fast Tucker rank reduction for non-negative tensors using mean-field approximation. In Advances in Neural Information Processing Systems, volume 34, pp. 443 443-454, Virtual Event, December 2021.
K. Ghalamkari and M. Sugiyama. Non-negative low-rank approximations for mulii-dimennsioall arrays on statistical manifold. Information Geometry, 6:257-292, 2023.
K. Ghalamkari, M. Sugiyama, and Y. Kawahara. Many-body approximation for non-negative tensors. In Advances in Neural Information Processing Systems, volume 36, pp. 257 292, New Orleans, US, December 2023.
N. Gillis and F. Glineur. Accelerated multiplicative updates and hicrarchical als algorithms for nonnegative matrix factorization. Neural Compulation, 24(4):1085 1105, 2012.
I. Glasser, R. Sweke, N. Pancotti, J. Eiscrt, and I. Cirac. Expressive power of tensor-network factorizations for probabilistic modeling. Advances in Newral Information Processing Systems, 32, 2019.
C. Gong, D. Wang, and Q Liu. Alphamatch: Improving consistency for semi-supervised learning with divea-eivergence. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 13683 -3692, 2021.
I. Goodfellow, Y. Bengio, A. Courville, and Y. Bengio. Deep learning, volume 1. MIT Press Cambridge, 2016.
W. Grathwohl, R. T. Q Chen, J. Bettencourt, and D. Duvenaud. Scalable reversible generative models with free-form continuous dynamics. In International Conference Learning Represendations, 2019.
18

---

Published in Transactions on Machine Learning Research (4/2026)
J. Hernandez-Lobato, Y. Li, M. Rowland, T. Bui, D. Hernández-Lobato and R. Turner. Black-box alpha divergence minimization. In International Conference on Machine Learning, pp. 1511-1520. PMLR, 2016.
J. L. Hinrich and M. Mørup. Probabilistic tensor train decomposition. In 2019 27th European Signal Processing Conference (EUSIPCO), pp. 1-5. IEEE, 2019.
F. L. Hitchcock. The expression of a tensor or a polyadic as sum of products. Journal of Mathemalics and Physics, 6(1-4):164-189, 1927.
E. Hoogeboom, J. Peters, R. van den Berg, and M. Welling. Integer discrete flows and lossless compression. In Advances in Neural Information Processing Systems, volume 32, Curran Associates, Inc., 2019.
K. Huang and N. D. Sidiropoulos. Kullback-Leible principal component for tensors not NP-hard. In 2017 51st Asilomar Conference on Signals, Systems, and Computers, pp. 503-697. IEFE, 2017.
S. Iblisdir, R. Orus, and I. Latorre. Matrix product states algorithms and continuous systems. Physical Review B, 75(10)1104305, 2007.
S. Ibrahim and X. Fu. Recovering joint probability of discrete random variables from pairwise marginals. IEEE Transactions on Signal Processing, 69:4116-4131, 2021.
A. Iqbal and A. K. Seghouane. An a-divergence-basse based approach for robust dictionary learning. 1EEE Transactions on Image Processing, (5722 55399, 2019.
C. F. Jeff Wu. On the convergence properties of the EM algorithm. The Anmals of Statistics, pp. 95-103, 1983.
H. Jeffreys. An invariant form for the prior probability in estimation problems. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 186(1007):433440. 1946.
J. Jensen. Sur les fonctions convexes les inégalités entre les valeurs moyennes. Acta mathematicn, 30(1):
175 193, 1906.
K. Kargas, N. D. Sidiropoulos, and X. Fu. Tensors, learning, and "Kolmogorov extension" for ainite-alphabet random vectors. IEEE Transoctions on Signal Processing, 66(18):4864-4888, 2018.
N. Kargas and N. D. Sidiropoulos. Completing a joint PMF from projections: low-rank coupled tensor factorization approach. In 2017 Information Theory and Applications Workshop (ITA), pp. 6. IEEE, 2017.
Y. D. Kim and Choi. Nonnegative Tucker decomposition. In 2007 IEEE conference on computer vision and pallern recognilion, pp. 1-8. IEFE, 2007.
Y. D. Kim, A. Cichocki, and S. Choi. Nonnegative Tucker decomposition with alpha-divergence. In 2008 IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 1829- 1832,, 2008.
D. P. Kingma and J. Ba. Adam: method for stochastic optimization. In International Conference on Learning Representations (TCLR), 2015.
N. K. Kitson, A. C. Constantinou, Z. Guo, Y. Liu, and K. Chobtham. A survey of bayesian network structure learning. Artificial Intelligence Review, 56(8))8721-855) 8814, 2023.
I. Kobyzev, S. J.D. Prince, and M. A. Brubaker. Normalizing flows: An introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, 43(11):3964-35)) 3979, 2021.
T. G. Kolda and B. W. Bader. Tensor decompositions and applications. STAM review, (3):45455 (000)000,,000 2009.
A. Krizhevsky. Learning multiple layers of features from tiny images. Technical report, University of Toronto, 2009.
G. R. Kurri, M. Welfert, T. Sypherd, and L. Sankar. o-gan: Convergence and estimation guarantees. In
2022 IEFE International Symposium on Information Theory (ISTT), pp. IEFF, 2022.
19

---

Published in Transactions on Machine Learning Research (4/2026)
K. Lange. A gradient algorithm locally equivalent to the EM algorithm. Journal of the Royal Statistical Sociely: Series B (Melhodological), 57(2):425 437, 1995.
Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner. Gradient-based learning applied document recognition. Proceedings of the IEFE, 86(11):2278 2324, 1998.
Y. Li and Y. Gal. Dropout inference in bayesian neural networks with alpha-divergences. In International Conference on Machine Learning, pp. 2052-2061. PMLR, 2017.
Y. Li and R. E. Turner. Rényi divergence variational inference. In Advances in Neural Information Processing Systems, volume 29, 2016.
J. Lin. Divergence measures based the shannon entropy. IEEE Transactions Information Theory, 37 (1):145-151, 1991.
Y. Liu, F. Shang, L. Jiao, J. Cheng, and H. Cheng. Trace norm regularized CANDECOMP/PARAL decomposition with missing data. IEEE Transactions Cybernetics, 45(11):2437-2448 2015.
Y. Liu, Z. Long, and C. Zhu. Image completion using low tensor tree rank and total variation minimization. IEEE Transactions on Multimedia, 21(2):338-35) 35000 2018.
Y. Liu, J. Liu, Z. Long, and C. Zhu. Tensor computation for data analysis. Springer, 2022.
L. Loconte, A. Mari, G. Gala, R. Peharz, C. de Campos, E. Quacghebeur, G. Vessio, and A. Vergari. What is the relationship between tensor factorizations and circuits (and how can we exploit it)? Transactions Machine Learning Research, 2025. ISSN 2835-8856.
Y. Matsuyama. The a-EM algorithm and its applications. In 2000 IEEE International Conference on Acoustics, Speech, and Signal Processing. Proceedings (Cat. No. 00CH37100), volume pp. 292-595, IEEE, 2000.
Matsuyama. The a-EM algorithm: surrogate likelihood maximization using a-logarithmic information measurcs. IEEE Transactions on Information Theory, 4933):692-706, 2003.
A. Miri Rekavandi, A. K. Seghouane, and R. J. Evans. Learning robust and sparse principal components with the divieergence. IEEE Transactions on Image Processing, 33:3441-3455, 2024.
M. Nickel, V. Tresp, and H.-P. Kriegel. A three-way model for collective learning on multi-relational data. In Proceedings of the 28th International Conference on International Conference on Machine Learning, pp. 809 816, Madison, WI, USA, 2011. Omnipress. ISBN 9781450306195.
G. S. Novikov, M. E. Panov, and I. V. Oseledets. Tensor-train density estimation. In Uncertainty in Artificial Intelligence, pp. 1321 1331. PMLR, 2021.
R. S. Olson, W. La Cava, P. Orzechowski, Ryan J. Urbanowicz, and Jason H. Moore. PMLB: large benchmark suite for machine learning evaluation and comparison. BioData Mining, 10(1):36, Dec 2017. ISSN 1756-0381.
I. V. Oseledets. Tensor-train decomposition. SIAM Journal on Scientific Computing, (3(5):2295 2317, 2011.
Y. Panagakis, J. Kossaifi, G. G. Chrysos, J. Oldfield, M. A. Nicolaou, A. Anandkumar, and S. Zafeirion. Tensor methods in computer vision and deep learning. Proceedings of the IEEE, 109(5):863 890, 2021.
E. E. Papalexakis, C. Faloutsos, and N. D. Sidiropoulos. Tensors for data mining and data fusion: Models, applications, and scalable algorithms. ACM Trans. Inlell. Syst. Technol, 8(2), oct 2016. ISSN 2157-6904.
G. Papamakarios, E. Nalismick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1 64, 2021.
A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and
A. Lerer. Automatic differentiation in pytorch. 2017.
20

---

Published in Transactions on Machine Learning Research (4/2026)
J. Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988. ISBN 1558604790.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825 2830, 2011.
B. Póczos and J. Schneider. On the estimation of alpha-divergenoes. In Proceedings of the Fourteenth International Conference Artificial Intelligence and Statistics, pp. 609-617. JMLR Workshop and Conference Proceedings, 2011.
P. Rai, C. Hu, M. Harding, and L. Carin. Scalable probabilistic tensor factorization for binary and count data. In International Joint Conferences Artificial Intelligence, 2015.
Rényi. On measures of entropy and information. In Proceedings of the fourth Berkeley symposium on mathematical statistics and probability, volume contributions to the theory of statistics, volume 4, pp. 547-562. University of California Press, 1961.
D. Rezende and S. Mohamed. Variational inference with normalizing flows. In Francis Bach and David Blei (eds.), Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 1530 1538, Lille, France, 07 09 Jul 2015. PMLR.
H. Robbins and S. Monro. stochastic approximation method. The Annals of Mathematical Statistics, pp- 400-407, 1951.
E. Robeva and A. Seigal. Duality of graphical models and tensor networks. Information and Inference: A Journal of the IMA, 8(2):273 288, 06 2018. ISSN 2049-8772.
O. Román. A practical introduction to tensor networks: Matrix product states and projected entangled pair stales. Annals of Physics, 149:117-158, 2014. ISSN 0003-4916.
J. Schreiber. Pomegranate: fast and flexible probabilistic modeling in python. Journal of Machine Learning Research, 18(164):1-6, 2018.
P. Sharma and P. M. Pradhan. Development of amari alpha divergence- based gradient-descent least mean square algorithm. IEEE Transactions Cireuits and Systems II: Express Briefs, 70(8):319443111 2023.
N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos. Tensor decomposition for signal processing and machine learning. IEFF Transactions Signal Processing, 65 (13):3551 3582, 2017.
Q Song, H. Ge, J. Caverlee, and X Hu. Tensor completion algorithms in big data analytics. ACM Trans. Knowl. Discov. Data, 13(1), 2019. ISSN 1556-4681.
M. Sugiyama, H. Nakahara, and K. Tsuda. Legendre decomposition for tensors. Journal of Statistical Mechanics: Theory and Experiment, 2019(12):124017, 2019.
Y. Sun, P. Babu, and D. P. Palomar. Majjoriation-mmninizaaoo algorithms in signal processing, communica- tions, and machine learning. JEEE Transactions on Signal Processing, 65(3):794 816, 2016.
P. Tichavský and O. Straka. Optimizing the order of modes in Tensor Train decomposition. IEEE Signal Processing Letlers, 32:1361 1365, 2025.
T. Tieleman and G. Hinton. Lecture 6.5 6.5-RMSprop, coursera: Neural networks for machine learning, University of Toronto, Technical Report, 6, 2012.
G. Tomasi and R. Bro. PARAFAC and missing values. Chemometries and Intelligent Laboratory Systems, 75 (2):163-16) 180, 2005.
D. Tran, K. Vafa, K. Agrawal, L. Dinh, and B. Poole. Discrete flows: Invertible generative models of discrete data. In Advances in Neural Informalion Processing Systems, volume 32. Curran Associates, Tne., 2019.
21

---

Published in Transactions on Machine Learning Research (4/2026)
L. R. Tucker. Some mathematical notes threc-mode factor analysis. Psychometrika, 31(3):279 --1, 1966.
van Erven and P. Harremos. Rényi divergence and Kullback-Leibler divergence. 1EEE Transactions on Information Theory, 60(7):3797-3820, 2014.
C. Villacampa-Calvo, G. Hernandez-Munoz, and D. Hernandez-Lobato. Alpha-divergence minimization for deep gaussian processes. International Journal of Approzimate Reasoning, 150:139-171 2022.
J. Vora, K. S. Gurumoorthy, and A. Rajwade. Recovery of joint probability distribution from one-way marginals: Low rank tensors and random projections. In 2021 TEEE Statistionl Signal Processing Workshop (SSP), pp. 481 485. IEEE, 2021.
D. Wang, C. Gong, M. Li, Q. Liu, and V. Chandra. Alphanet: Improved training of supernets with alpha-divergence. In International Conference Machine Learning, pp. 10760-10771. PMLR, 2021.
Z. Wang, Q. Gu, Y. Ning, and H. Liu. High dimensional EM algorithm: Statistical optimization and asymptotic normality. Advances Advonces in Neural Information Processing Systems, volume 28, 2015.
R. Wu, J. Liu, B. Li, Anh-Huy Phan, Ivan V. Oseledets, C. Zhu, and Y. Liu. TERM model: Tensor ring mixture model for density estimation. IEEE Transactions on Big Data, pp. 1-16, 2025.
A. Yeredor and M. Haardt. Maximum likelihood estimation of low-rank probability mass tensor from partial obscrvations. IEEE Signal Processing Letters, 26(10):1551-1555, 2019.
Y. K. Yılmaz and A. T. Cemgil. Probabilistic latent tensor factorization. In International Conference on Latent Variable Analysis and Signal Separation, pp. 346-353 Springer, 2010.
22

---

Published in Transactions on Machine Learning Research (4/2026)
Supplementary Material
Table 4: List of the notations used throughout the paper.
Symbol
Description
The label of structures, e-g., {CP, Tucker, TT, bg} The dimension (the number of modes) of the tensor T, P, plki, and F The number of hidden indices for the structure k, e.g,, V(CP) = =1 and V = -
D
Vk
D vk
The dimension of the tensor M, @ik), @1h,, and Φ
[T]
The set of natural number smaller than I, i.c., [T] = (1,2..... 1
T E 12-00
Given empirical non-negative normalized tensor
DE RR Iip
Estimated (mixture) low-rank lensor, P Σx plk]
Pl
phxxun paxxxx x..xxxx
k-th low-rank tensor component in the mixture model P Mixture ratio for component (Ex r = 1 and r/ 2 0) Higher-order tensor before summation over indices r, i.e., E Eren 2 P!!!!
20
= (1
[kl)
Q
QIM
The weighted @ik), i.e. @lki 2
F EC
The tensor used in the first upper bound (E1-step)
The tensor normalized over r and h; i.e., Σk Eren $
14 ir
=
The tuple of 114, i.e., + = - T(K])
c
The of tensor JF satisfying E Eien, F Fi =1
M
Intermediate tensor M P [g
[k]
a E (0,1)
a-divergence paramcter controlling cass-covering vs mode-secking Objective function based on the a-divergence
La
La
The first bound of the I., The second bound of the La
La
H(A,B)
Cross-entropy between tensors A and B, H(A,B) -E log Bi Standardized cross-entropy for the tensor F, ha(F) = H(T o F, FF(1((11 a a-divergence from T to P
hF(F)
Da(T,P)
(RI(T,P)
The KL divergence from T to P Factor matrices matrices/tensors for CP, Tucker, and TT structures
A(d), G, g(d)
The index sets of and P, Le., Q2 [h]
[Tp]
Index set except mode d, i.e., 2.d [7] x x [Ta-1] (ia) [T4t1] ... [TD]
27 C St
The set of indices of non-zero value of T, i.e., 229 = {i E Ω2 T # 0}
Oge
The partial index sets of Qill and Φlk), i.e., Oxx = [R]] ... [Rya] The tensor rank associated with & structure, c.g., Ca-rank RCP and TT-rank RTT Uniform background tensor (P) bg]
R
plbk)
IR. xpp
1/192/]
Weight of the background term (0 < 1/ [bb] 1)
5(,,j)
Kronecker delta, i.e., 5(i.j) 1 if Landau symbol for time computational complexity
5(i.j) 0 otherwise
O
Ep [X]]
Expectation of X' with respect the distribution P, i.c., Ep [2] = E; PiKk Real numbers
IR
12
Non-negative real numbers
The index of the tensor T, F, and P, i.c., i (f.....,ip)
2ge
A part of the index of the tensor of Mil QIN), and कlA,, i.c., r (FI....
23

---

Published in Transactions on Machine Learning Research (4/2026)
Table 5: Abbreviations used throughout the paper.
Abbreviation
Meaning
E algorithm
Proposed double-bounded EM algorithm optimizing the a-divergence Expectation-Maximiaalli algorithm optimizing the KI-divergence First expectation step in E2M tightening the upper bound La Second expectation slep in E tightening the second upper bound Lo Integrated step combining the E1- and E2-steps Minimizing the second upper bound La via many-body approximation Kullback-Leible divergence
EM algorithm
El-step
E2-step
E²-step
M-step
KT. divergence
JS divergence
Jensen-Shanmnn divergence
CP
SANDECOMP/PARAPCA decomposition
TT
Tensor Train dccomposition
CNMF
Coupled Non-negative Matrix Factorization
BM
Born Machine
GD
Gradient Descent
MU
Multiplicative Update
DNN
Deep Neural Network
A Proofs
A.1 Proofs for E*M-algorithm
We prove the propositions used to derive the E algorithm for non-negative tensor mixture learning
Section 2
Proposition 1. For any (0,1) and negateee tensor T, let the set of tensors C and D be
C
F
T 1, Fe
R
XID
D
(T11) ,..., xlk))
Σ
Σ
T
$
R xxxxxxxxx X Rayk
k
respectively. Then, for any tensors F E C and $ = (011) (Dl'!...,,$11]
E D. it holds that
(P) < La P) < (F,
where the upper bounds are given as
(T, H(T o F, P) + ha(F).
T.a O + ha (F) + Exe H(M ((k)),
and the functions J and ha and each. element of the tensor M
<R RR.
Ixx Ruy
given
M [8 hT(F) H(T F,F))
J(n) - ΣΣ Σ
M!K TT F. [t]
F.
a
iED; 1=1 r
respectively-
24

---

Published in Transactions on Machine Learning Research (4/2026)
Proof. The first bound (P) s La(F (F, P)
can be shown follows:
1
(P( =
log tEll
Tipp!
log
Fi
tE
p pl-a
5
FiT.
log
(19)
Fi
iE
Σ
FIT: T log Pi --1
Σ
FiT:" log
iESt
ieft
H(T o, (F) (F,P)
where the following relation, as defined by Jensen's 's incquality (Jensen 1906).
used:
M Am
M
(xmms
(20)
m=
valid for any concave function f:R. + and real numbers At,..., Am that satisfies Em= Am = 1. This
inequality can be applied in Equation 19 for (0, observing that the logarithm function concave function and it holds that E.
by the construction F FeC C.
The second bound (F,P) - (F ②)
follows by defining the any tensor + = (Tl1)... satisfies Σ EreR,e o [M] = 1. Using this tensor, we can transform the upper bound (F, P) as follows:
k(K)) E D that
La(FP))= -- Σ F.T log Ps + ha(FF
ig
K FiT. log
[M] ir
=
Σ
ir
+ (F)
(21)
iE0i
k=1 EE
ir
[
-
M F.T. log
$ R
+ (F)
tEDi 1= rER
Σ E Σ
FiTs 4 14
log 6.h
HMM!), (nh)) (F(FF)
iED 1=1 rER
(() + hF) H(M o 11 (M l(kl) LL 4 2)
where we again use Jensen's inequality 20 in Equation 21 As a result, we obtain the double bound (P) s Lo (F,P) L. ②)
Remark 1. In Equation (3, the summation Er inside the (1 - a) power function makes the optimization challenging. As seen in the above proof, our upper bound La eliminates the problematic power term by leveraging the properties of logarithms, i.e., log P! a1- This trick motivates optimizing Equation 3 instead of Equation ① directly.
Proposition 2. For the E1-slep, the optåmal update for F E EC === {F EiEn, T. = =1 that minimizes the upper bound La (F, P) is given as
p1-
F
Eien,
T.P!
25

---

Published in Transactions on Machine Learning Research (4/2026)
Proof.
We put Equation 9 into the upper bound La in Equation 6
LF P) H o F¹ P) + (F*)
P!-a
F th log
Fi
iE
P Pl-a
T:
log iE
Σ
T.
Q
Σ Liefi
T:pp--a
iE
1
log Σ !:
a
iE
(P)).
Proposition ] shows that
La(P) = LF S (P(
(22)
for any tensors F E C. Thus, the tensor F is optimal.
Proposition 3. For the E2-step, the optimal update for + = (Φ11), "oll),..,,,.A ,bk1)) E D that minimizes the upper
bound [F, T ②) is given
2 1 Ck=1 Erer
bb
[k] ir
gk
Proof.
We put Equation
10 into the upper bound La in Equation
Σx K
(F 2-), J(n + ha (F) +
H(M
@ - HM (+))
k=1
[ [k]
Fr [K]
log
+ (F)
(23)
iEfi k= g
K
Σ ΣEΣ
ir
log Σ 1=
2
hF)FF)
ief; R= 1 rEst gk
2k=1
K
E Ereal np
rEs
Σ T F log Σ
SlF + (F)FF
iED;
rEQ Rh
Σ iER;
T:e F; log P; + (F)
HTT FF hha =
Proposition D shows that
[F(F = (F, TTa +, (1
(24)
for any tensors & D where D (TDl1) ,.., plK1) Ex Ereg T1 =1. Thus, the tensors Φ* (p-lIl..... (x1 are optimal.
Proposition 4. For the M-step, the optimal update for the miature ratio n = (p......p)k)) that optimizes ((n in Equation given as
Eien, Erer gt ir EiEr, Σk Erer M.
M!
n
26

---

Published in Transactions on Machine Learning Research (4/2026)
Proof. We optimize the decoupled objective function
J(n) Σ Σ
M
[k] log 1g
iE k=
re
with conditions E 1/ = and n 2 0. Thus we consider the following Lagrange function
K
K
M
log log n -
Σ
1 [kl
iEft
R= rER
The condition 0L/0m/ 0 leads to the optimal ratio
Σ
M
iest rE
where the multiplier À is given
A = Σ Σ M! iEtl k=
Theorem The objective function of the misture of tensor factorizations using the E² M M-procedure always converges regardless of the choice of low-rank structure and mindures.
Proof. We prove the convergence of the E algorithm using the following two observations: (i) the objective
function La is bounded below, and (ii) the objective function L. is non-increasing along the iterations. The proof of (i) can be found in Theorem The El-step in iteration updates J by optimal FF11 to minimize the upper bound such that
(van Erven & Harremos 2014). We show observation (ii) below.
La (P())
Eq. 22
La ( pt) & La(F P*).
The E2-step in ileration updates df by optimal t g++ to minimize the upper bound for t such that
(F pt)
Eq. 24
(a ++ ②') s L. I (F++1 4 ②³).
Then, the M-step updates Q¹ by Qu+l to minimize the second upper bound for å such that
[Ft+1 +1 Ø++1) < La(tt 4
2¹).
In the next update in E2-step, x++1 is replaced by the optimal 4'+2
(F p(+1)= ( . Q++1) Ia ( $141, Q
Again, the El-step updates F++ by optimal F++ to minimize the upper bound such that Combining the above three relations, we obtain
L. (P#+1) Eq. L (T++2 p(t1) & (Ft+l pttl).
La (pt+l) = La p(+1) < LF (Itt1,ptt)) = Lo (Ftt1, git2 ,422,0041) (t+1)
s (Fttl 16+1 ( -60) x
(Ft+l ,p') La (Pt).
Thus, it holds that (pit1) S L. (Pt). This concludes the proof, since bounded monotone sequence converges.
For simplicity of both proof and notation, consider the E2-step rather than the E1-step following the M-stcp. This does affect the algorithm, El-step and F22step independent of each other. fact, Algorithm performs these steps jointly F stepp
27

---

Published in Transactions on Machine Learning Research (4/2026)
A.2 Proofs for exact solutions of many-body approximation
First, show the known solution formulas for the best CP rank-1 approximation that globally minimizes the KL divergence from the given tensor.
Theorem 2 (Optimal M-step in CP dccompasition (Huang & Sidiropoulos 2017
For given non-negative
tensor M E 12 (24xx-xXXXX ils many-body approcimation with interactions as described in Figure 3/b) b) is given by
Eienty, Mir
100 bur
== E iEDi refla
Σ Mir
WIDD (Eien, 1-1/D #=
Proof. Please refer to the original paper by Huang & Sidiropoulos 2017
In the following, we provide proofs of the closed-form solutions of many-body approximation in Figure 3 and (d). When a factor in low-body tensor multiplied by v, the value of the objective function of many-body approximation remains the same if another factor is multiplied by 1/V, which we call the scaling redundancy. The key idea in the following proofs is reducing the scaling redundancy and absorbing the normalizing conditions the entire tensor into single factor. This enables the docoupling of the normalizod condition of the entire lensor into independent. conditions for each of the factors. This trick is also used in Section B.5 to decouple more complicated low-rank structures into combination of CP, Tucker, and TT dccompositions.
Proposition 5 (The optimal M-step in Tucker decomposition). For given tensor M R R4x.../D x 20
its mody-body approzimation with the interaction described in Figure 3/e) is given by
Eienyy Erenhy Mir Eieni Ereche Mir
EiEr, Mir Gr Eien, Erena Mir'
A(d) Lara
Proof.
The objective function of the many-body approximation is
HiM, Q (Tucke))) = Σ
Σ
Mir log 2
(Tcker
(25)
tEN, refta
where
(1) ... ( W..diom...ro Gr....r
Tucker]
Since the many-body approximation parameterizes tensors as discrete probability distributions, we optimize the above objective function with the normalizing condition Σien, Erean Tuckke Then, we consider the following Lagrange function:
L= L= E E-E
Σ Mir log (1 ... ((0) - (
Σ E iE rella
A A..A ADrD
-1) -1)
RED2 PELLR
To reduce the scaling redundancy and decouple the normalizing condition, introduce scaled factor matrices ( given by
A (d) Sara (d)
((d)
where
(d) Σ
Alara
(26)
ar
for -= 1,2,.... D and the scaled core tensor,
Gr G,a ... ..)
28

---

Published in Transactions on Machine Learning Research (4/2026)
The normalizing condition Eien, Ereag ir
Tucker]
guarantees the normalization of the tensor 9
Σ
(27)
rett
The tensor Q can be represented with the scaled matrices A(d) and core tensor 9 introduced above
Tuckeer Gr A = G Ji A
We optimize 9 and A iara
instead of & and Alard Thus the Lagrange function can be written as
Σ Gr -1
-1)
C= Σ Σ
A - x Mer log 9
A ( Σ (d)
-1
(28)
iE rEit
d=1
nd
id
The condition
OC
00,
θ A
iara
leads the optimal core tensor and factor matrices given by
Gr
Mar.
1
Σ Ara
E
Σ
Mir.
EE
iEn),
ren
The values of the Lagrange multipliers are obtained by the normalizing conditions 26 and
27 as
== Σ E iEN, rein
E
Mers
1 = Σ Σ
Mar.
iEf
ren
Proposition (The optimal M-step in Train decomposition). For given tensor M
R 1 xxxx ils many-body approzimation with interactions described in Figure 3(d) is given by 20
Eieni Drenkera
Mir
(d) gra
Eien, Erenge Mir Rira
for d = I,..., D, assuming Tn = =D 1.
Proof.
The objective function of the many-body approximation is
((TT) = Σ Σ
Mir log ,TT
iel refla
where
(( g(2) GDDD
Since the many-body approximation parameterizes tensors as discrete probability distributions, we optimize the above objective function with the normalizing condition EieR, Erean Q.TT] 1. Then, we consider the following Lagrange function:
-= Σ Σ
Mar log 9,10 Griarz
c(2)
((D) P
>
Σ Σ
(0 g.2) ... D -1)
iED: refl
29

---

Published in Transactions on Machine Learning Research (4/2026)
To decouple the normalizing condition and make the problem simpler, we introduce scaled core tensors g())
G(D-1 I that are normalized over rd-1 and ia
(d-1)
g
gra-1
(d)
Pa-riars
(d) ra-riaara
gra
where define
(1 = ΣΣ F--1 bd
- (d-1)
Pa-riara 9ra-1
(29)
with gro 1. We assume
-TD = Using the scaled core tensors, the tensor Q can be written as
Train io....o
- g(2) Prears ... 0 ITp-isp PiD
((1))
0(2)
Prears ... ((D)
with
1 (D)
(D-1) ID-
(30)
The matrix g(P) is normalized, satislying Σrp-1 Eip -rD- written
(D)
Thus, the Lagrange function can be
D-1
C E E
Mar log ((1) ((2) inf F122r2
GLD)
Σ
(a)
ΣΣ
g(l) Pa-ifare -1)
TD-
iEs, PEOR
d=1
Pd-1 3o
(DD
Σ Σ
(DD)
(31)
TD-1
FD
The critical condition
Oc 0
(d)
1 1tara
leads to the optimal core tensors
g(d)
Σ
Σ
M
Xra
VER).,
rEa
1d.-d1 R.ra-ra-1
where the values multipliers are identified by the normalizing conditions Equation 30
as
(d) = Σ Σ Ard
Mir-
iE
rEn
R.ra
B Additional Remarks
B.1 Expectation maximization and E-steps
Assuming a < 1, the El-step is equivalent to maximizing the expectation of log log(pl-a// apl-a/F] with respect to the distribution W = T o F since the objective can be reformulated as follows:
arg min (F, FEC
p1-a FiT.
arg min
EEAT
log
Fec
Fi
iEn,
pl-a
arg max Ew
log
Fec
F
30

---

(IT)) such that Published in Transactions on Machine Learning Research (4/2026)
where used Equation 19 and the definition of the expectation over the discrete domain Ep|X] = E, Pad for any non-negative normalized tensor P. Here, PI-0/F denotes the element-wise (1 a-) power of P, followed by element-wise division by F. The logarithm function is element-wisely operated. We note that the normalization of the tensor W is immediately guaranteed by the definition of tensor F in Equation 4 Similarly, the E2-step is equivalent lo maximizing the expectation of log (Ω/D) with respect to the distribution M since the objective can be reformulated as follows:
arg min L (F, 2, 2) EED
K
0 log
= arg min RED
Σ
T:
to
[k]
① ir
iEli 1=1 rERRt
arg max EM BED
log
where we used Equation 23
We also note that the normalization of the tensor M, i.c E, =1,
is guaranteed by its definition in Equation
B.2 Scalable E E-M-algorithm for tensor train
Assuming TT-rank R[TT]- (R....: : R), we reduce the computational cost of the E2MTT to Q(yDN by computing the sum over the indices r E RRTT as follows. Firstly, we introduce the following tensors
Σ
( -
gldtt P
Σ
gldt1) Sraiepiaas glutt=)) Trab (i42amppp.
(32)
ra+1
with ((+11) = g(), -(0-1-)='' and g g(D+) 1. Each complexity is O(Ra) in order to get ((+d) and g(d+) when compute them in the order of g(+2), g((+3) ..... ,...,g(-D)) and G(D-2-), G(D --- ....., g(++), respectively. The tensor P[TT] can be written p(TT) -((=D) g(n-), Consequently, the update rule in Equation 16 can be written
Eierple Eieng P. (
TH P.
gid) Pia-zian Pist1....p.r iD.r
Gldt)
G Iiar piar
(33)
Te (-od) gl(t)
for 22 = -29 - We used the relation Mir TilirFa and Dir = Pir/P; to get the above update rule. tensors; however, as shown in Equation 33 we only require the elements on indices i in the subsets 529. Therefore, under the sparsity assumption on the given tensor T, the merged cores remain sparse and their cost of storing, proportional to the rank R, is negligible compared to that of the updated core tensor g(d), which scales by the square of R. Since the number of elements in U2 N, the resulting complexity O(>DNR²). We provide the above procedure of E2MTT in Algorithm 2 As seen above, the developed algorithm climinates the complexity of 1o achieving scalability that is proportional to the number of non-zero values N, thus making it suitable for nigh-dimensional data.
As D increases, the merged cores in Equation 332 become higher-dimensions
B.3 Complexity of gradient-based method
We consider a low-rank structure P: P==Di Pr=Ere Qir and corresponding T interactions g!!)....,!!!
Q+r 0,10 p p @ The indices of the tensor Q can be given as
(in...,, in) E [4] x ... [/ = Ut, r == (ri,...,.y) E [Ri] x .... x [Rv] QR.
assuming P E RR XID and Q R x Rv The indices of each interaction ((t) can be written
it = (ik, " i (kg, [Tx]] x [Txy] x ... [n,y,, r' r-= Tea FEm.. E [Re,) x [Ri,] x ... x (Reach }
31

---

Published in Transactions on Machine Learning Research (4/2026)
Algorithm 2: Effficient E algorithm for 'T decomposition
input : Non-ncgative tensor and train rank R!TT] (Ri.... RD-1).
Initialize core lensors g())...,
repeat
Update F; for all E S2% using Equation 3
9
// El-step
4
for d 1 to D do Obtain G(-+d) and G(D-d-) using Equation ,
5
32
Update ga), ,.... G(D) using Equation 33 for d x-1 1 to D do
// Combined E2 and M-step
7
L Obtain ((-+) using Equation 32
8
P!TT! + g ( for i SED2 27;
9
10 until Convergence;
11 return gi.,,..,I
assuming ky < kut1 and ty < fut1 without loss of generality, where M, and S are the number of visible and hidden variables in the interaction ((0), respectively. Technically, indices k2 and by should be labeled with the index of interaction I, such as kn and e, respectively; however, we omitted the label to simplify the notation. The derivative of the objective a-divergence (a P) with respect each element of the interaction @ can be given
@D.
Σ P PER;
OP, 00
Σ d'ER, r'E
Σ
0Qr To P
(34)
00 (1) p'
00
(1)
j'pd
j'pe
The summation over i' E Qr can be replaced by the summation over i' E 92 as discussed in Section 3.4 Thus, the naive complexity is O(N(RR)). In addition, thanks to the derivative
dQr
2 (e)
077 007
0 - Oftip piiree O i (0) r"),
00 00(0) (8) pt
ae
2
where 8(a,b) b) is if b, 0 otherwise, the summation survive only when 3' it and r² r" Thus, we need consider summation over only
Ω N2..tgmmn.n i' E 22 n 2 1 A and r' E Rire Feg
ka....ky
The resulting complexity is consistent with Table 2 When we use the softmax transformation to enforce non-negativity and normalization constraints, applying the chain rule in Equation 34 yields the same result.
B.4 Low-rank approximation meets many-body approximation in the EM algorithm
In Section 2 the tensor T constructed by observed samples is approximated by a low-rank tensor, i... = P; = Krr 0.4l We naturally interpret i as visible variable because ia represents the d-th feature of the observed data. In contrast, the variable r in the tensor Q(k) does not have a direct correspondence to any obscrved feature, and the model P implicitly includes r by summing over r, i.c., Σ, Thus, we can regard as hidden variable and its summation Ex as marginalization, following standard conventions in the machine learning community (Huang & Sidiropoulos 2017 Ibrahim & Fu 2021). The degrees of freedom of the hidden variables (R1.... Ry*) correspond to the tensor rank of pMa]
The non-convex nature of low-rank tensor decomposition arises from the hidden variables r in the model. On the other hand, in the many-body approximation optimizing Equation (13). the model Q contains no indicos, explicitly or implicitly, that are not present in the given tensor M. In this sense, the many-body approximation corresponds optimization without hidden variables, which leads convex optimization. The essence of this work lies in relaxing a low-rank approximation with hidden variables into many-body
32

---

Published in Transactions on Machine Learning Research (4/2026)
approximation without hidden variables, thereby allowing find approximate solutions via iteratively applying the global optimal solutions of the many-body approximation obtained by convex optimization. In this interpretation, the tensor @ik] corresponds to the joint distribution pi iD "',...,,'') the tensor (l corresponds to the conditional distribution p(k, F1...., TV (i,...;ip), and the model P corresponds to the marginalized distribution p(; Ex, (ii dD ,
B.5 Closed-form updates for more general low-rank structures
We now discuss how to find the solution for the many-body ap- proximation required in the M-step when a more complex low-rank structure is assumed in the model. Our strategy is to decouple the cross-entropy H to be optimized in the M-step and the normaliza- tion condition into components that admit closed-form solutions for structures such as CP, Tucker, and TT.
Tucker
Train
G
H
12
A
D
E
In the following, we discuss the decomposition as described Figure 8 which known as a typical tensor tree structure (Liu et al.) 2018), while the generalization to arbitrary
13
1.
low-rank structures is straightforward. For given tensor M[ the objective function of many-body approxination in the M-step can be decoupled as follows:
Figure 8: An example a tensor tree structure represented by the tensor net- work.
H HMM [Teel, Q [reel)
E Σ
M log Grirers Bixrs Croro Hygrers Eraia br
iESt rEiR (Tucker)
O ([Tukker)) + H(MITT) Q(TT)),
(35)
where we define the tensor Treel as
Q
[Treel
Grivars Aisri Bigra Crove Disra Hagerro Brais
(36)
ir
and
[Tucker]
Σ [Tree]
M TT igAerSrATE
Σ
M
Tree
Vaierara
$1f2r1r275
[Tucker]] Cistrrr2rste
Grirers A Bigry Cryrer
TT Mgtarstare
Disry Hiyners Erais
Since the tensor TTTe))) needs to satisfy the normalized condition, we consider the following Lagrange function
(>1
L == (MM((M( @[Tucker))- H (MIIT) QITTI) - E
2rr
Tree]
(37)
iEft
refa
Although the objective function has been decoupled as seen Equation the normalizing condition
(35
we still need a treatmentt for
Σ
Treel
(38)
dr
ig
reta
in order to apply the closed-form solutions in Equations Equation (35). respectively.
15 and
16 to the first and second terms in
To decouple the normalizing condition, we reduce scaling redundancy by scaling each factor and decouple the
Lagrange function into independent parts as explained at the beginning of Section A.2 More specifically, we define a single root tensor and introduce normalized factors that sum over the edges that lic below the root. Although the choice of the root lensor is not unique, we let tensor G be the root lensor and introduce
1
1
lree
Any On
Air Bigra
b Bigra: Craro
Crarer
(39)
Cra
1
1
drg Hrararor hro
D Diyrs drs
Digras
Eagrs
Paraa Hrggare
(40)
eri
33

---

Published in Transactions on Machine Learning Research (4/2026)
where each normalizer is defined as
d7 Airrs
brz
Σ
B Bigrzs
Crs = >x=Econ
Crsrohre
1/
drg
D Digron
Cra Σ
Birrs
hroe Σ
Crgers Hestavo
1374
As a result, it holds that
Σ Ann
Σ Bigra =
E Digrs
Eears Fears=
Σ Hrgraro 1.
We define the tensor & as Grivars ar, brg Grivers and putting Equations we obtain the normalizing condition for the root tensor G as
39 and 40 into Equations 36
and
38
Σ PIP2P8
Inrzro = 1.
Then, the tensor can be written as
Q [Teeel = Grears Aigrz Bigrz Orgro Disry Hergors Brgis.
(41)
The above approach to reduce scaling redundancy is illustrated in Figure 9 Finally, the original optimization problem with the Lagrange function in Equation 37 is equivalent to the problem with the Lagrange function
C L Tucker] + C(T
where
C [Tucker]
Σ
M
Tucker] iaierrranooe
log Grizro Aiir Bigry Orgro + 2
Σ Gnprars -1)
Σ
Aum
Σ A
Birra -
Σ x
Σ Crara
-1
+
This is equivalent to the Lagrange function for the Tucker decomposition given in Equation 28 and
CTTT]_
Σ
M [T
Vaiersaree log Hergrare Digry Bigra
iglarara7a
Σ H Hergaro -1).
Σ D iars -
-1+
+ + X Σ Eiars
Σ
H A
which is also equivalent to the Lagrange function for the TT decomposition given in Equation 31 assuming G(D) is a normalized uniform tensor. Then, we solve these independent many-body approximations by the closed-form solution given by Equations (15) and (16) for tensors M[Tucker] and MITTI, respectively, and multiply solutions to get optimal tensor Q as Equation 41h, which satisfied the normalizing condition given Equation
B.6 Selection of mixture components
Our framework supports mixtures of low-rank tensors. This flexibility raises the question of how to determine the components of the mixture. One approach to manually select the low-rank structures to leverage the strengths of different structures to improve generalization in density estimation. For instance, the CP decomposition offers intuitive representation with fewer parameters; the TT decomposition captures latent.
34

---

Published in Transactions on Machine Learning Research (4/2026)
Tucker
Train
h
C
g
H
a
b
d
e
A
B
D
E
Figure 9: We normalize all tensors except for the root tensor, which is enclosed in bold line. We then push the normalizer of each lensor, a, b, c, e,d,e d, e, and on the root lensor. The root lensor absorbs scaling redundancy. This procedure decouples the Lagrangian into two independent problems, Tuccker and LTT
pairwise interactions using more parameters; while the Tucker decomposition accounts for all combinations of latent. interactions across the attributes, resulting in an exponential increase in the number of parameters. A Alternatively, we can employ cross-validatio to automatically select suitable structures. We also emphasize that it not always possible uniquely define the optimal mixture components for given dataset, low-rank structures can overlap for example, rank-R CP decomposition be expressed using rank-Rt TT decomposition.
B.7 EMM-algorithm for continuous distributions
A generalization of the E algorithm to contimuous distributions is straightforward, as shown below. We consider optimizing the a-divergence from given distribution T(x) to a model P(z),
(T(T,P) all all- - 1 Jax T(x) dx
{
(42)
where distributions T(x) and P(x) are defined on the continuous sample space Qx of the random variable x E Qx. For a E (0, 1), this optimization is equivalent to minimizing the following Rényi divergence:
La(P) = a 1 log
P(x)
Jüx
We assume the model P is defined mixture of K distributions .i....,Qu Qx with hidden variables =, i... the model P can be given as follows::
P(x) = Ef -ro)=E/ee
②(k, x z)dz where =
& SEB
for the discrete index & E { 1.2...... K} and the weight ryk) satisfying Σg (kk) and 17(k) [0, 1]. The double bounds in Equations 6 and 7 and the update rules in Equations (10) and 12 hold by replacing the sum Σx E, with Ex dzz, the tensors & ($121, ...., ..., ΦIK])E D with the conditional distribution 4(k: x), the tensors = M[K]] with the joint distribution M(k, a, z), and the lensor Fe C with any continuous function satisfying the condition f F(x))"((x)) = 1. When the given T(x) is general continuous probability [= dx, is nontrivial. However, in typical density estimation scenario where T represents observed data Ω2 (ar ..... (e("),..,,!!N))) (N)) C Ωxx we have T(x) = I/N for x E 9% and T(x) = 0 otherwise. Consequently, the denominator becomes MN 1-a
distribution, the integral in the denominator Equation
35

---

Published in Transactions on Machine Learning Research (4/2026)
B.7.1 Example: Gaussian Mixture Model
As an example of the continuous extension of the E algorithm, we consider a one-dimensional Gaussian mixture model
K Po(x) = Σ (k)) k=
(a- - mx))
where (o(k,z)= = V2T0
cxp
(43)
20
for æ E f2x = R. We denote the collection of all parameters by 8, i.c., 0 = (0,,..., 0x) and 0g = (px, ok)
We fit given one-dimensional N samples (2), #2,..., IN) to the model in Equation 43p. The model has
no hidden variable in each mixed component (k. Given Equations (9 (10p
12 and the known
closed-form update rule in the standard EM algorithm for Gaussian mixture (Bishop following update rules:
2006
we obtain the
E²-step:
Po(xn)) Σn Po(
n0k) for all and n.
Mkn t Na
(44)
Po(xn)
M-step:
()k) + Na
MEn for all k
(45)
n=1
ΣN
Mania for all ky
1k +
(46)
m== Mkn
ΣN
Man(In - m&)²
o
for all k
(47)
EN Mkn
where M is K x N matrix. When updating the variance of in Equation 47 we use the newly updated
mean A4 in Equation 446). After each M-step, we update the distributions Po(x) and (k(k using the
updated parameters Pk. T% and mh by Equation (43
To verify the above update rules, we generated 10,000 samples from one-dimensional Gaussian mixture with two components, where 0, (0.0, 0.8), 0₂ (4.0, 1.0), and n(1) = 0.6, 772) 0.4, and fitted Gaussian mixture model with K = 2 by optimizing the a-divergence for various values of The true distribution and the resulting estimated distributions are visualized in Figure 2, where we can confirm consistency with the known fact that a smaller a focuses on the dominant mode [L & Turner 2016
We remark on the connection between our framework and the recent work by Daudel et al. (002bb) Their approach derives bound for the a-divergence by exploiting the incquality log ud S u and introducing an auxiliary function g(0, 8') such that g(0,0) S g(0,0' S 9(0', 0)) in contrast our Jensen's 's s-inequality-bue based double bound. Although these two bounds different, algorithmically, the resulting update rule for the one-dimensional Gaussian mixture is identical.
C Additional experimental results
In this section, we discuss additional experimental results that could not be included in the main text due Lo page limitations. We note that the existing FMCP (Huang & Sidiropoulos, 2017 Yeredor & Haardt, 2019 is equivalent to the proposed method for k = CP with a = 1.0 and K 1.
Performance on optimization Figure 10 revisits Figure aa by changing the value of the learning rate of the bascline methods and the a. Despite requiring no hyperparameter tuning for optimization, our method achieves better convergence performance than modern maddient-based approaches with carefully tuned learning rates. These results demonstrate that the proposed EPM algorithm guarantees monotonic
36

---

Published in Transactions on Machine Learning Research (4/2026)
0.3
0.5
0.7
0.9
1
10 1
10"2
10ah
10-2
10
a-div.
a-div.
a-div.
a-div
10-3
10-3
10-3
10-3
102
10
102
10
102
10*
102
10
Iteration
Iteration
Iteration
Iteration
Adagrad (Ir:0.0001)
RMSprop (Ir=0.0001)
SGD (Ir=0.1)
SGD 1000)
Adam (Ir=0.01)
Adagrad (Ir=0.001)
RMSprop (Ir=0.001)
SGD (Ir=1)
Adam (Ir=0,0001)
Adam (Ir=0.1)
Adagrad (lr=0,01)
RMSprap (Ir=0.01)
SGD (Ir=10)
Adam (Ir=0.001)
ECMCPTTB
Adagrad (Ir=0.1)
RMSprop (Ir=0.1)
SGD (Ir=100)
Figure 10: Learning curves with varying a values and different optimizers for reconstructing the SIPI House color image (log-log scale).
C1: 0.3
0: 0.5
Q:: 0.7
0.9
0.10
0.10
0.10
0.10
0.06
0.06
a-div.
a-div.
a.d.v.
a.div.
0.04
0.04
0.04
0.04
0.02
0.02
0.02
0.02
0.00
0.00
0.00
0.00
0
200
400
600
200
400
600
200
400
600
0
200
400
600
Iteration
Iteration
Iteration
Iteration
Adagrad (Ir=0.1)
RMSprop (Ir=0.1)
(Im=1000)
(Ir=0.1)
Adagrad (Ir=0.01)
RMSprop (Ir=0.01)
SGD (lr=100)
Adam (Ir=1)
(Ir=0.001) EMCPTTB
Adagrad (Ir=0.001)
RMSprop (Ir=0.001)
Adam (Ir=0.1)
Adagrad (Ir=0.0001)
RMSprop (lr=0.0001)
(Ir=1)
Comparison of the number of iterations required to reconstruct the SIP1 House color image by the mixtures of CP, TT, and background term, with varying a values.
Figure 11:
decrease in the objective function, whereas existing gradient- based methods may suffer from oscillation divergence. In addition, we investigate the behavior of the objective function during the carlier iterations of oplimization, as shown in Figure []] Interestingly, we observe that convergence is slower for smaller values of a. To demonstrate the reproducibility of this phenomenon, Figure [2] presents an additional comparison of optimization performance, including the Tucker structure in the model, showing that convergence is again relatively slower for smaller a. We hypothesize that this slowdown is related the increasing gap between Ia and T., for larger values of It has been reported that in gradient-based a-divergence optimization, the
37

---

Published in Transactions on Machine Learning Research (4/2026)
α: 0.3
0:
0.5
0.7
0.9
0.10
0.10
0.10
0.10
0.08
0.08
0.08
0.06
0.06
a-div.
a-div.
a-div.
a.div.
0.04
0.04
0.04
0.04
0.02
0.02
0.02
0.02
0.00
0.00
0.00
0.00
0
200
400
600
200
400
600
200 400
600
0
200
400
600
Iteration
Iteration
Iteration
Iteration
Adagrad (Ir-0.1)
RMSprop (Ir-0.1)
SGD (I--1000)
SGD (Ir-0.1)
Adam (r-00001
Adagrad (Ir=0,01)
RMSprop (I==0.01)
SGD (Ir=100)
Adam (Ir=1)
Adam (Ir=0.001)
Adagrad (Ir:0,001)
RMSprop (I=00.001)
SGD (Ir=10)
Adam (I==0.1)
EPMCPTucke/TTB
Adagrad (Ir-0.0001)
RMSprop (lr00.0001) SGD (Ir-1)
Comparison of the number of iterations required to reconstruct the SIPI House color image by the mixture of CP, Tucker, TT, and background term, with varying a values.
Figure 12:
convergence rale depends on the value of a (Bao el al. 2025). A theoretical investigation of the convergence rate within the proposed E'M algorithm, particularly in relation to the value of a, would be an interesting direction for future work.
Mass covering property induced by the a paramcter To investigate the dependence of noise sensitivity on the value of a, we inject mislabeled samples and p outliers into the empirical distribution, and observe the resulting changes in Figure 6. Figure 13 shows the results for varying values of and a. In the visualization of reconstructions and empirical distributions, each pixel (i,j) is plotted with a color indicating its label assignment: red represents label one and blue represents label two, which is more detailed in Section D.1.2 with the description of how we supply the outliers and the mislabeled samples. For p 0 in the top row, the reconstruction is accurate regardless of the value of However, as p incrcases, the reconstruction becomes more sensitive to outliers and mislabeled samples for large a, while the reconstruction is still accurate for small a, which is also confirmed in Figure 6
We also conducted the quantitative study shown in Figure uusin a smaller low-rank model with R[CP] and RTTT] = (2, 2), and the results are presented in Figure 14a When the model size is reduced, overfitting does not occur, and the model cannot learn noise, resulting in weaker dependence of the test error on noise compared to Figure 7 In addition, conducted the above experiments the noisy CarEvaluatoo dataset with increasing noise levels. We augmented the dataset with random noise samples drawn from uniform distribution and randomly flipped the class labels of a subset of samples. The number of added noise samples matches the number of label-flipped samples. The rank was tuned so that the number of parameters of the low-rank model is set to about 40% of the tensor size. In Figure 14b] we consistently observe that minimizing the KT, divergence does not necessarily yield optimal performance with respect to other divergence measures, whereas our E algorithm, with an appropriatcly chosen a, improves performance under alternative divergences.
Density estimation with real datasets and ablation study The full version of Table 3 is available in Table 6 We also provide tuned
a of E2MCPTTB that maximizes the score for validation datasets on the rightmost column in Table ( In the density estimation evaluated by negative log-likllhhod, it is natural that 1.0 yields the best validation score a setting without distribution shift because optimizing the KT. divergence is equivalent to optimizing the negative log-likelihood. On the other hand,
38

---

Published in Transactions on Machine Learning Research (4/2026)
Empirical dist.
a:0.1
a: 0.3
a: 0.5
a: 0.7
a: 0.9
Figure 13: Reconstructions of the empirical distribution with outliers by MCCTTТ using different values. The further down the rows, the greater the number outlicrs and noise contained in the empirical distribution. Specifically, each row shows the rection the empirical distribution with 2p contaminations (consisting of p outliers and p random noise points) for p 50, 75, 100, and 200 from top to bottom.
39

---

Published in Transactions on Machine Learning Research (4/2026)
KL divergence
Hellinger distance
Reverse KL
J5 divergence
Jeffreys divergence
0.80
12.00
0.50
12.00
15.00
0.45
10.00
Test 6100
0.70
Test smar
210
0.40
800
10.00
0.35
0.50
6.00
6.00
0.25
600
4.00
Number of notse
0.20
0.00
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.8-
0.5
15.0
0.7
8.0
Trahing error
Training orror
8.0
0.4
0.6
6.0
Trining
0.3
Training
6.0
5.0
0.6
4.0-
0.2
0.5
0.0
0.5
1.0
0.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
(a) Noisy half-moon data with smaller ranks [CCP] and R!T (2.2).
KL divergence
Hellinger distance
Reverse KL
Js divergence
Jeffreys divergence
2,00
3.00
0.35
2.00
0.08
2.00
1.50
1.00
1.00
0.25
1.00
0.04
Number noise
0.00
0.20
0.50
4.00
2.50
0.35-
2.00
0.10
0.08
Train
2.00
0.30
1.00
1.00
0.25
Fa:::
0.20
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
1.0
(b) Noisy CarEvaluation data with ranks [CCP]] 120 and TT (2, ...,....,))
Figure 14: Reconstruction error of noisy data using E²MCPTTB, varying the noise level and the value of a.
the fact that classification performance validation datasets is maximized at + 1.0 suggests that the KL divergence is not necessarily the most appropriate measure for classification tasks, depending on the dataset. In the last row of Table G we report. the average classification accuracy to evaluate the overall performance across multiple datasets. In contrast, for density estimation tasks, such an aggregate evaluation is not feasible due differences in the scale of the sample space across datasets.
For the ablation study, also compare the generalization performance of FºMCPB, F²MTTB, and E?MCPTTB for density estimation in Table 7 Interestingly, the TT model cxhibits lower generalization performance compared to the CP models. However, our mixture low-rank modeling enables elfective lybridization of the representation of CP and the TT, thereby achieving better generalization performance.
Comparison to non-tensor-bassed methods for density estimation We also provide experimental results comparing the proposed E2M algorithm to non-tensor-basee methods, including a tree-structured Bayesian network (BaysNet) (Pearl 1988 Kitson al. 2023), a hidden Markov model (HMM) (Baum & Petrie 1966), and discrete normalizing flow (DiscreteF low) (Tan al.) 2019h, in Table & Restricting its network structure to tree, BaysNet. enables efficient. and fully automated structure learning from
40

---

Published in Transactions on Machine Learning Research (4/2026)
data (Schreiber 2018 However, because it cannot model latent variables, its performance is expected to degrade on datasets where latent representations are essential. DiscreteFlow requires tuning of the number of units, the number of layers, the temperature parameter, and the learning rate. Both the HMM and the DiscreteFlow are not invariant to permutations of feature axes; in other words, they exhibit order dependency. In contrast, the lensor-based approaches do not suffer from order dependency when using CP decompositions, whereas the TT decomposition does (Tichavský & Straka, 2025). A benefit of tensor-basod methods is that the user can flexibly choose whether to incorporate order dependencies by selecting an appropriate low-rank structure. Although the most suitable method depends on the characteristics of the dataset, such whether linear or nonlinear interactions dominate, whether order dependency is present, and whether latent variables are needed, our tensor-bass based framework consistently achieves performance comparable to existing approaches.
41

---

Published in Transactions on Machine Learning Research (4/2026)
Table 6: Accuracy of the classification task (top lines, higher is better) and negative log-likelioood (bottom lines, lower is better) for all datasets. The bottom row reports the mean accuracy across all 34 datasets. Error bars are given as the standard deviation of the mean across five randomly initialized runs, and for the bottom row, the 34 datasets. The rightmost column shows the value of for E²MCPTTB that maximizes the score for validation data.
Dataset
BM
KLTT
CNMF
EMCP
ECMCTTТB
Optimal
0.571 (0.000)
0.589 (0.018)
0.571 (0.00)
0.571 (0.000)
0.500 (0.000)
0.90
AsiaLung
2.240 (0.003)
2.295 (0.055)
3.130 (0.091)
2.262 (0.000)
2.480 (0.163)
0.95
0.777 (0.000)
0.830 (0.000)
0.649 (0.074)
0.868 (0.014)
0.872 (0.000)
0.90
B.Scale
7.317 (0.033)
7.480 (0.153)
7.418 (0.022)
8.979 (0.491)
7.102 (0.000)
1.00
0.901 (0.011)
0.923 (0.011)
0.956 (0.000)
0.912 (0.000)
0.934 (0.000)
1.00
BCW
13.120 (0.115)
12.63 (0.004)
12.776 (0.146)
12.815 (0.000)
12.623 (0.000)
1.00
0.847 (0.044)
0.892 (0.003)
0.736 (0.006)
0.942 (0.018)
0.950 (0.020)
1.00
CarFxal.
7.964 (0.145)
7.752 (0.017)
8. 186 (0.062)
7.777 (0.037)
7.847 (0.067)
1.00
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.75
Chess
12.45 (0.300)
12.45 (0.170)
14.739 (0.010)
10.306 (0.050)
10.886 (0.011)
0.95
0.157 (0.000)
0.277 (0.01)
0.243 (0.012)
0.439
(0.007)
0.573
(0.003)
0.95
Chess2
13.134 (0.001)
12.561 (0.123)
12.772 (0.029)
11.948 (0.081)
11.637 (0.004)
0.95
0.537 (0.024)
0.585
(0.073)
0.573
(0.037)
0.561
(0.000)
0.585
(0.000)
0.95
Cleveland
6.379 (0.048)
6.353 (0.272)
6.144 (0.033)
6.139 (0.000)
6.235 (0.000)
0.95
0.864 (0.015)
0.909 (0.000)
0.909 (0.000)
0.909
(0.000)
0.909
(0.000)
1.00
ConfAd
6.919 (0.035)
6.979 (0.006)
6.918 (0.032)
6.846 (0.000)
6.859 (0.000)
1.00
0.633 (0.011)
0.644 (0.000)
0.644 (0.000)
0.644
(0.000)
0.644
(0.000)
0.50
Coronary
3.510 (0.000)
3.525 (0.003)
3.568 (0.008)
3.511 (0.004)
3.527 (0.000)
1.00
0.773 (0.058)
0.547 (0.000)
0.677 (0.131)
0.547
(0.000)
0.872 (0.020)
1.00
Credit
7.117 (0.013)
7.334 (0.062)
8.518 (0.167)
8.371 (0.000)
7.034 (0.020)
1.00
0.167 (0.014)
0,158 (0.032)
0.203 (0.032)
0.234
(0.000)
0.153 (0.000)
1.00
DMFT
7.313 (0.190)
7.212 (0.010)
7.565 (0.091)
7.165 (0.000)
7.213 (0.017)
0.95
0.906 (0.010)
0.896 (0.000)
0.854 (0.021)
0.625
(0.000)
0.844 (0.010)
1.00
DTCR
10.013 (0.275)
10.111 (0.100)
10.921 (0.177)
12.065 (0.000)
10.239 (0.158)
1.00
0.767 (0.057)
0.814 (0.043)
0.792 (0.042)
0.831
(0.000)
0.814 (0.000)
0.90
GermanGSS
6.561 (0.015)
6.514 (0.130)
6.507 (0.042)
6.312 (0.000)
6.456 (0.040)
0.95
0.639 (0.028)
0.639 (0.028)
0.306 (0.028)
0.689
(0.027)
0.792
(0.173)
0.95
Hayesroth
6.640 (0.093)
6.762 (0.136)
6.457 (0.187)
5.887 (0.109)
5.879 (0.100)
1.00
0.517 (0.000)
0.515 (0.000)
0.514 (0.003)
0.509
(0.000)
0.545
(0.003)
0.95
Income
10.510 (0.306)
10.686 (0.018)
11.646 (0.056)
11.408 (0.000)
9.435 (0.004)
0.95
0.273 (0.043)
0.191 (0.032)
0.245 (0.014)
0.446
(0.000)
0.378
(0.004)
1.00
Led7
5.584 (0.238)
5.936 (0.237)
5.741 (0.014)
4.646 (0.026)
4.870 (0.006)
1.00
0.750 (0.000)
0.750 (0.000)
0.750
(0.000)
0.750
(0.000)
0.750 (0.000)
0.85
Lenses
2.945 (0.000)
2.945 (0.000)
3.028 (0.026)
2.945 (0.000)
3.407 (0.035)
1.00
0.705 (0.023)
0.705 (0.023)
0.523 (0.023)
0.500 (0.000)
0.659 (0.023)
0.95
Tymphography
12.855 (0.093)
13.117 (0.125)
14.028 (0.092)
14.080 (0.000)
12.894 (0.074)
1.00
1.000 (0.000)
0.974 (0.015)
0.803 (0.003)
0.957 (0.019)
1.000 (0.000)
0.90
Mofn
7.065 (0.010)
7.053 (0.048)
7.343 (0.018)
7.212 (0.005)
7.192 (0.011)
1.00
0.985 (0.027)
0.769 (0.237)
0.604 (0.073)
1.000 (0.000)
0.996 (0.007)
0.75
Monk
6.404 (0.048)
6.654 (0.277)
6.828 (0.041)
6.272 (0.059)
6.527 (0.065)
0.90
42

---

Published in Transactions on Machine Learning Research (4/2026)
0.945 (0.007)
0.832 (0.034)
0.817 (0.075)
0.975
(0.010)
0.978
(0.014)
1.00
Mushroom
13.196 (0.515)
16.613 (0.109)
17.439 (0.566)
10.118 (0.284)
10.756 (0.000)
1.00
0.730 (0.006)
0.700 (0.018)
0.572 (0.075)
0.991 (0.003)
0.991
(0.001)
1.00
Nursery
10.004 (0.041)
10.022 (0.014)
10.495 (0.029)
9.642 (0.018)
9.689 (0.014)
1.00
0.575 (0.246)
0.866 (0.232)
0.477 (0.020)
1.000 (0.000)
1.000 (0.000)
0.90
Parity5p5
7.539 (0.316)
7.201 (0.326)
7.658 (0.007)
7.254 (0.026)
7.265 (0.068)
1.00
0.654 (0.192)
0.692 (0.000)
0.615 (0.000)
0.615
(0.000)
0.692 (0.077)
0.50
PPD
7.072 (0.009)
7.072 (0.009)
7.073 (0.039)
7.038 (0.000)
7.196 (0.000)
0.95
0.795 (0.000)
0.727 (0.000)
0.659 (0.000)
0.739
(0.011)
0.705 (0.000)
0.75
PTumor
7.766 (0.204)
7.734 (0.016)
7.731 (0.036)
7.391 (0.000)
7.737 (0.035)
1.00
0.552 (0.006)
0.657 (0.017)
0.669 (0.006)
0.721
(0.000)
0.686 (0.023)
0.95
Sensory
11.794 (0.207)
11.554 (0.003)
12.613 (0.124)
11.721 (0.077)
11.349 (0.180)
1.00
0.988 (0.004)
0.992 (0.000)
0.992
(0.000)
0.992
(0.000)
0.992
(0.000)
0.15
SolarFlare
6.778 (0.146)
6.921 (0.125)
7.725 (0.176)
6.700 (0.155)
6.459 (0.063)
1.00
0.846 (0.137)
0.949 (0.000)
0.949
(0.000)
0.936
(0.022)
0.910 (0.043)
0.90
SPECT
13.085 (0.270)
13.130 (0.373)
13.028 (0.408)
12.828 (0.335)
12.658 (0.262)
1.00
0.987 (0.013)
0.987 (0.013)
0.760 (0.032)
0.987 (0.008)
1.000 (0.000)
0.90
ThreeOfNine
6.639 (0.034)
6.504 (0.008)
6.822 (0.012)
6.652 (0.031)
6.717 (0.011)
0.95
0.326 (0.000)
0.239 (0.000)
0.239 (0.000)
0.304 (0.000)
0.402 (0.033)
0.90
Tumor
7.942 (0.054)
8.578 (0.418)
7.808 (0.036)
7.739 (0.000)
7.610 (0.003)
0.95
0.714 (0.143)
0.286 (0.143)
0.929 (0.071)
1.000 (0.000)
1.000 (0.000)
1.00
Vehicle
4.993 (0.000)
4.993 (0.000)
4.615 (0.102)
4.523 (0.000)
5.305 (0.146)
0.95
0.758 (0.016)
0.806 (0.000)
0.823 (0.048)
0.839 (0.000)
0.968 (0.000)
1.00
Votes
8,416 (0.159)
7.964 (0.078)
7.680 (0.031)
8.026 (0.000)
8.058 (0.000)
1.00
1.000 (0.000)
1.000
(0.000)
0.741 (0.015)
1.000 (0.000)
1.000 (0.000)
0.95
XD6
6.357 (0.022)
6.335 (0.042)
6.818 (0.005)
6.484 (0.014)
6.497 (0.003)
1.00
Mean Acc.
0.708 (0.040)
0.698 (0.042)
0.649 (0.038)
0.743 (0.039)
0.776 (0.039)
Table 7: Accuracy of the classification task (top lines, higher is better) and negative log-likelihood (bottom linos, lower is better) for all datasets. The bottom row reports the mean accuracy across all 34 datascts. Error bars are given as the standard deviation of the mean across five randomly initialized runs, and for the bottom row, the 34 datasets.
Datasel
EMCP
EMCPB
EMTTB
EMCPTТB
E*MCPB
E.MTBB
E2MCPTTB
0.571 (0.000)
0.464 (0.000)
0.464 (0.000)
0.464 (0.000)
0.464 (0.000)
0.464 (0.000)
0.500 (0.000)
AsiaLung
2.262 (0.000)
2.620 (0.001)
2.619 (0.000)
2.580 (0.000)
2.591 (0.000)
2.594 (0.000)
2.480 (0.
0.856 (0.005)
0.835 (0.005)
0.843 (0.018)
0.872 (0.000)
0.849 (0.016)
0.843 (0.018)
0.872 (0.000)
B.Scale
8.979 (0.491)
7.158 (0.002)
7.106 (0.020)
7.102 (0.000)
7.158 (0.002)
7.106 (0.020)
7.102 (0.000)
0.912 (0.000)
0.879 (0.000)
0.549 (0.000)
0.934 (0.000)
0.879 (0.000)
0.549 (0.000)
0.934 (0.000)
BCW
12.815 (0.000)
12.936 (0.000)
13.908 (0.000)
12.623 (0.000)
12.936 (0.000)
13.908 (0.000)
12.623 (0.000)
0.942 (0.018)
0.905 (0.010)
0.928 (0.010)
0.950 (0.020)
0.905 (0.010)
0.928 (0.010)
0.950 (0.020)
CarEval.
7.777 (0.037)
8.022 (0.001)
7.947 (0.063)
7.847 (0.067)
8.022 (0.001)
7.947 (0.063)
7.847 (0.067)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
Chess
10.306 (0.046)
10.229 (0.023)
13.902 (0.016)
10.908 (0.050)
10.229 (0.023)
13.902 (0.016)
10.886 (0.009)
0.439 (0.010)
0.646 (0.013)
0.280 (0.012)
0.562 (0.008)
0.658 (0.005)
0.280 (0.012)
0.573 (0.003)
Chess2
11.948 (0.081)
11.549 (0.004)
12.647 (0.016)
11.681 (0.010)
11.549 (0.004)
12.647 (0.016)
11.637 (0.004)
(continued on next page)
43

---

Published in Transactions on Machine Learning Research (4/2026)
Dataset
EMCP
EMCPB
EMTTB
EMCPTTB
E*MCPB
E.MTTB
E*MOPTTB
0.585 (0.000)
0.610 (0.000)
0.585 (0.000)
0.537 (0.024)
0.610 (0.000)
0.517 (0.082)
0.585 (0.000)
Cleveland
6.139 (0.000)
6.270 (0.000)
6.290 (0.000)
6.233 (0.000)
6.268 (0.000)
6,290 (0.000)
6.235 (0.000)
0.909 (0.000)
0.879 (0.000)
0.800 (0.036)
0.909
(0.000)
0.879 (0.000)
0.806 (0.024)
0.909 (0.000)
ConfAd
6.846 (0.000)
6.957 (0.000)
7.017 (0.079)
6.859 (0.000)
6.957 (0.000)
7.017 (0.079)
6.859 (0.000)
0.644 (0.000)
0.644 (0.000)
0.644 (0.000)
0.644
(0.000)
0.644 (0.000)
0.644 (0.000)
0.644 (0.000)
Coronary
3.511 (0.004)
3.563 (0.002)
3.581 (0.000)
3.527 (0.000)
3.563 (0.002)
3.574 (0.007)
3.527 (0.000)
0.547 (0.000)
0.860 (0.000)
0.802 (0.007)
0.872
(0.020)
0.814 (0.013)
0.814 (0.010)
0.872 (0.020)
Credit
8.371 (0.000)
7.080 (0.035)
7.102 (0.039)
7.034 (0.020)
7.043 (0.021)
7.102 (0.039)
7.034 (0.020)
0.234 (0.000)
0.176 (0.014)
0.216 (0.008)
0.153
(0.000)
0.176 (0.014)
0.202 (0.019)
0.153 (0.000)
DMFT
7.165 (0.000)
7.211 (0.000)
7.203 (0.018)
7.186 (0.011)
7.210 (0.000)
7.203 (0.018)
7.213 (0.017)
0.625 (0.000)
0.792 (0.000)
0.838 (0.008)
0.844
(0.010)
0.792 (0.000)
0.838 (0.008)
0.844 (0.010)
DTCR
12.065 (0.000)
10.235 (0.100)
10.530 (0.072)
10.239 (0.158)
10.235 (0.100)
10.530 (0.072)
10.239 (0.158)
0.835 (0.007)
0.729 (0.051)
0.797 (0.037)
0.814
(0.000)
0.800 (0.022)
0.793 (0.027)
0.814 (0.000)
GermanGSS
6.312 (0.000)
6.493 (0.000)
6.493 (0.000)
6.461 (0.053)
6.493 (0.000)
6.493 (0.000)
6.456 (0.040)
0.694 (0.028)
0.306 (0.028)
0.500 (0.035)
0.875 (0.024)
0.300 (0.044)
0.344 (0.022)
0.792 (0.173)
Hayesroth
5.802 (0.046)
5.623 (0.073)
5.648 (0.075)
5.879 (0.100)
5.778 (0.068)
5.594 (0.025)
5.879 (0.100)
0.509 (0.000)
0.547
(0.000)
0.529 (0.003)
0.544 (0.002)
0.547
(0.000)
0.530 (0.001)
0.545 (0.003)
Income
11.408 (0.000)
9.365 (0.001)
9.569 (0.013)
9.446 (0.021)
9.365
(0.001)
9.569 (0.013)
9.435 (0.00-1)
0.446 (0.000)
0.432
(0.000)
0.249 (0.030)
0.378 (0.004)
0.422
(0.017)
0.276 (0.015)
0.378 (0.004)
Led7
1.646 (0.026)
4.823 (0.100)
5.453 (0.053)
4.870 (0.006)
4.756 (0.044)
5.417 (0.049)
4.870 (0.006)
0.750 (0.000)
0.750
(0.000)
0.750 (0.000)
0.625 (0.125)
0.750
(0.000)
0.750 (0.000)
0.750 (0.000)
Lenses
2.945 (0.000)
3.098 (0.000)
3.198 (0.201)
3.407 (0.055)
3.277 (0.073)
3.198 (0.201)
3.407 (0.055)
0.500 (0.000)
0.682
(0.000)
0.555 (0.018)
0.591 (0.000)
0.664 (0.022)
0.564 (0.022)
0.659 (0.023)
Lymphography
14.080 (0.000)
13.641 (0.000)
13.330 (0.571)
12.894 (0.074)
13.641 (0.000)
13.330 (0.571)
12.894 (0.074)
0.967 (0.003)
0.923 (0.026)
0.986 (0.018)
0.995 (0.005)
0.923 (0.026)
0.999 (0.002)
1.000 (0.000)
Mofn
7.211 (0.005)
7.420 (0.000)
7.393 (0.021)
7.192 (0.011)
7.420 (0.000)
7.392 (0.023)
7.192 (0.011)
1.000 (0.000)
1.000 (0.000)
0.982 (0.023)
1.000 (0.000)
1.000 (0.000)
0.982 (0.023)
0.996 (0.007)
Monk
6.309 (0.088)
6.436 (0.008)
6.585 (0.085)
6.465 (0.054)
6,436 (0.008)
6.585 (0.085)
6.527 (0.065)
0.975 (0.010)
0.994 (0.006)
0.950 (0.008)
0.978 (0.014)
0.994 (0.006)
0.950 (0.007)
0.978 (0.014)
Mushroom
10.118 (0.284)
9.566 (0.125)
12.317 (0.331)
10.756 (0.000)
9.566 (0.125)
12.317 (0.331)
10.756 (0.000)
0.991 (0.003)
0.994 (0.002)
0.980 (0.013)
0.991
(0.001)
0.994
(0.002)
0.980 (0.013)
0.991 (0.001)
Nursery
9.642 (0.018)
9.800 (0.001)
9.732 (0.046)
9.689 (0.014)
9.800 (0.001)
9.732 (0.046)
9.689 (0.014)
0.954 (0.037)
0.973 (0.027)
1.000
(0.000)
0.995
(0.008)
1.000 (0.000)
1.000
(0.000)
1.000 (0.000)
Parity5p5
7.431 (0.026)
7.395 (0.003)
7.164 (0.012)
7.265 (0.068)
7.395 (0.003)
7.164 (0.012)
7.265 (0.068)
0.615 (0.000)
0.615 (0.000)
0.615 (0.000)
0.692 (0.077)
0.615 (0.000)
0.615 (0.000)
0.615 (0.000)
PPD
7.038 (0.000)
7.314 (0.000)
7.314 (0.000)
7.496 (0.000)
7.719 (0.064)
7.359 (0.000)
7.524 (0.000)
0.739 (0.011)
0.659 (0.000)
0.732 (0.009)
0.682 (0.000)
0.686 (0.056)
0.745 (0.027)
0.705 (0.000)
Pluumor
7.391 (0.000)
7.664 (0.222)
7.745 (0.015)
7.737 (0.035)
7.664 (0.222)
7.745 (0.045)
7.737 (0.035)
0.721 (0.000)
0.703 (0.006)
0.635 (0.037)
0.669 (0.006)
0.709 (0.019)
0.644 (0.038)
0.686 (0.023)
Sensory
11.721 (0.077)
11.928 (0.006)
11.003 (0.395)
11.349 (0.180)
11.928 (0.006)
11.003 (0.395)
11.349 (0.180)
0.992 (0.000)
0.992 (0.000)
0.992 (0.000)
0.992 (0.000)
0.992
(0.000)
0.992 (0.000)
0.992 (0.000)
SolarFlare
7.303 (0.071)
6.510 (0.000)
6.486 (0.026)
6.459 (0.063)
6.458 (0.046)
6.486 (0.026)
6.459 (0.063)
0.936 (0.022)
0.949 (0.000)
0.949 (0.000)
0.904 (0.033)
0.949
(0.000)
0.949 (0.000)
0.910 (0.043)
SPECT
12.630 (0. (13)
12.576 (0.090)
13.450 (0.161)
12.658 (0.262)
12.576 (0.090)
13.450 (0.
12.658 (0.262)
(continued on next page)
44

---

Published in Transactions on Machine Learning Research (4/2026)
Dataset
EMCP
EMCPR
EMTTB
EMCPTTB
E*MCPB
Е*MTBB
E*MOPTTB
0.961 (0.013)
0.968 (0.019)
0.919 (0.082)
0.981
(0.006)
0.968 (0.019)
0.919 (0.082)
1.000 (0.000)
OfreeOfNine
6.657 (0.002)
6.811 (0.013)
6.802 (0.101)
6.685
(0.006)
6.811 (0.013)
6.802 (0.101)
6.717 (0.011)
0.304 (0.000)
0.239 (0.000)
0.274 (0.029)
0.391
(0.000)
0.261 (0.000)
0.274 (0.029)
0.402 (0.033)
Tumor
7.739 (0.000)
7.669 (0.058)
7.970 (0.069)
7.576
(0.044)
7.669 (0.058)
7.970 (0.069)
7.610 (0.003)
1.000 (0.000)
0.857 (0.143)
0.629 (0.114)
1.000
(0.000)
1.000 (0.000)
0.886 (0.140)
1.000 (0.000)
Vehicle
4.523 (0.000)
4.952 (0.000)
4.952 (0.000)
5.239 (0.000)
4.952 (0.000)
4.639 (0.067)
5.305 (0.146)
0.839 (0.000)
0.968 (0.000)
0.774 (0.000)
0.968
(0.000)
0.968 (0.000)
0.871 (0.020)
0.968 (0.000)
Votes
7.867 (0.000)
8.092 (0.000)
8.175 (0.000)
8.058 (0.000)
8.092 (0.000)
8.175 (0.000)
8.058 (0.000)
1.000 (0.000)
0.978 (0.022)
0.788 (0.057)
1.000 (0.000)
0.978 (0.022)
0.775 (0.113)
1.000 (0.000)
XD6
6.487 (0.008)
6.731 (0.010)
6.644 (0.040)
6.497 (0.003)
6.731 (0.010)
6.644 (0.040)
6.497 (0.003)
Mean Acc.
0.743 (0.039)
0.741 (0.042)
0.698 (0.038)
0.767 (0.041)
0.748 (0.041)
0.703 (0.043)
0.776 (0.039)
45

---

Published in Transactions on Machine Learning Research (4/2026)
Table 8: Accuracy of the classification task (top lines, higher is better) and negative log- log-likelihood (bottom lines, lower is better) for all datasets, comparing non-tensor-bart based methods (BayNet, HMM, and DiscreteFlow to the proposed sensor-based methods (EMCPB and ECPTTT). The bottom reports the mean accuracy across all 34 datasets. Error bars are given as the standard deviation of the mean across five randomly initialized runs, and for the bottom row, the 34 datasets.
Dataset
BaysNet
HMM
DiscreteFlow
EMCPB
E?MCPTTE
0.817
0.768 (0.032)
0.580 (0.018)
0.464 (0.000)
0,500 (0.000)
AsiaLung
5.455
2.622 (0.081)
2.431 (0.046)
2.620 (0.001)
2,480 (0.163)
0.941
0.784 (0.066)
0.934
(0.016)
0.879 (0.000)
0.934 (0.000)
BCW
27.489
14.577 (0.195)
15.143 (0.478)
12.936 (0.000)
12.623 (0.000)
0.883
0.753 (0.028)
0.896 (0.013)
0.835 (0.005)
0.872 (0.000)
BalanceScale
7.126
7.374 (0.008)
7.021 (0.323)
7.158 (0.002)
7.102 (0.000)
0.846
0.710 (0.000)
0.874 (0.109)
0.905 (0.010)
0.950 (0.020)
Car CarFvaluation
7.851
8.197 (0.024)
8.020 (0.218)
8.022 (0.001)
7.847 (0.067)
0.522
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
0.500 (0.000)
Chess
12.567
14.8745 (0.002)
15.625 (0.209)
10.229 (0.023)
10.886 (0.009)
0.345
0.175 (0.004)
0.313 (0.089)
0.646 (0.013)
0.573 (0.003)
Chess2
12.386
13.031 (0.011)
13.799 (0.304)
11.549 (0.004)
11.637 (0.004)
0.600
0.678
(0.029)
0.549 (0.047)
0.610 (0.000)
0.585 (0.000)
Cleveland
7.054
6.211
(0.027)
6.868 (0.223)
6.270 (0.000)
6.235 (0.000)
0.892
0.892
(0.000)
0.909 (0.000)
0.879 (0.000)
0.909 (0.000)
ConfAd
10.728
6.825 (0.046)
7.290 (0.077)
6.957 (0.000)
6.859 (0.000)
0.888
0.888 (0.000)
0.644 (0.000)
0.644 (0.000)
0.644 (0.000)
Coronary
4.425
3.667 (0.024)
3.580 (0.004)
3.563 (0.002)
3.527 (0.000)
0.911
0.705 (0.086)
0.802 (0.171)
0.860 (0.000)
0.872 (0.020)
Credit
7.605
7.878 (0.018)
9.665 (0.320)
7,080 (0.035)
7.034 (0.020)
0.175
0.190 (0.049)
0.171 (0.019)
0.176
(0.014)
0.153 (0.000)
DMFT
7.527
7.316 (0.042)
7.642 (0.108)
7.211
(0.000)
7.213 (0.017)
0.947
0.908 (0.022)
0.906 (0.040)
0.792
(0.000)
0.844 (0.010)
DTCR
17.194
10.547 (0.091)
12.660 (0.165)
10.235 (0.100)
10.239 (0.158)
0.881
0.763 (0.000)
0.881 (0.081)
0.729 (0.051)
0.814 (0.000)
GermanGSS
6.388
6.558 (0.000)
6.759 (0.101)
6.493 (0.000)
6.456 (0.040)
0.600
0.537 (0.125)
0.639 (0.195)
0.306 (0.028)
0.792 (0.173)
Hayesroth
8.044
6.174 (0.037)
5.871 (0.063)
5.623 (0.073)
5.879 (0.100)
0.747
0.561 (0.119)
0.315 (0.086)
0.432 (0.000)
0.378 (0.004)
Led7
6.300
6.004 (0.025)
6.997 (0.283)
4.823 (0.100)
4.870 (0.006)
0.750
0.750 (0.000)
0.812 (0.125)
0.750 (0.000)
0.750 (0.000)
Lenses
9.612
2.945 (0.000)
2.793 (0.070)
3.098 (0.000)
3,407 (0.055)
0.818
0.602 (0.078)
0.670 (0.120)
0.682 (0.000)
0.659 (0.023)
Lymphography
25.344
12.984 (0.093)
16.037 (0.598)
13.641 (0.000)
12.894 (0.074)
0.849
0.774 (0.000)
1.000 (0.000)
0.923 (0.026)
1.000 (0.000)
Mofn
7.200
7.450 (0.015)
7.008 (0.041)
7.420 (0.000)
7.192 (0.011)
0.723
0.685 (0.029)
0.977 (0.046)
1.000 (0.000)
0.996 (0.007)
Monk
6.684
6.843 (0.031)
6.553 (0.232)
6.436 (0.008)
6.527 (0.065)
0.983
0.885 (0.045)
0.999 (0.000)
0.994 (0.006)
0.978 (0.014)
Mushroom
14.329
18.179 (0.143)
20.914 (1.093)
9.566 (0.125)
10.756 (0.000)
(continued next page)
46

---

Published in Transactions on Machine Learning Research (4/2026)
Dataset
BaysNet
HMM
DiscreteFlow
EMCPR
E*MCPTTR
0.901
0.728 (0.009)
0.813 (0.324)
0.994 (0.002)
0.991 (0.001)
Nursery
9.768
10.468 (0.033)
10.243 (0.588)
9.800 (0.001)
9.689 (0.014)
0.615
0.615 (0.000)
0.615 (0.000)
0.615 (0.000)
0.615 (0.000)
PPD
11.729
7.038 (0.000)
7.406 (0.040)
7.314 (0.000)
7.524 (0.000)
0.740
0.700 (0.000)
0.688 (0.034)
0.659 (0.000)
0.705 (0.000)
PTumor
7.812
7.558 (0.016)
8.084 (0.038)
7.664 (0.222)
7.737 (0.035)
0.450
0.482 (0.022)
0.753 (0.213)
0.973 (0.027)
1.000 (0.000)
Parityōp5
7.694
7.641 (0.000)
7.533 (0.201)
7.395 (0.003)
7.265 (0.068)
0.600
0.950 (0.000)
0.949 (0.000)
0.949 (0.000)
0.910 (0.043)
SPECT
11.282
13.189 (0.119)
12.347 (0.062)
12.576 (0.090)
12.658 (0.262)
0.709
0.674 (0.000)
0.695 (0.024)
0.703 (0.006)
0.686 (0.023)
Sensory
11.586
12.731 (0.035)
12.700 (0.486)
11.928 (0.006)
11.349 (0.180)
0.995
0.995 (0.000)
0.992 (0.000)
0.992
(0.000)
0.992 (0.000)
SolarFlare
7.936
7.549 (0.120)
9.013 (0.226)
6.510 (0.000)
6.459 (0.063)
0.831
0.653 (0.055)
1.000 (0.000)
0.968 (0.019)
1.000 (0.000)
ThreeOfNine
6.681
6.949 (0.024)
6.509 (0.064)
6.811 (0.013)
6.717 (0.011)
0.471
0.216 (0.000)
0.245 (0.011)
0.239
(0.000)
0.402 (0.033)
Tumor
11.746
8.062 (0.057)
10.982 (0.299)
7.669 (0.058)
7.610 (0.003)
1.000
0.321 (0.071)
0.893 (0.071)
0.857
(0.143)
1.000 (0.000)
Vehicle
4.593
4.997 (0.000)
5.124 (0.342)
4.952 (0.000)
5.305 (0.146)
0.943
0.814 (0.049)
0.968 (0.000)
0.968 (0.000)
0.968 (0.000)
Votes
8.270
9.497 (0.150)
8.741 (0.238)
8.092 (0.000)
8.058 (0.000)
0.801
0.717 (0.007)
1.000 (0.000)
0.978
(0.022)
1.000 (0.000)
XD6
6.687
6.821 (0.015)
6.413 (0.032)
6.731 (0.010)
6.497 (0.003)
Mean Acc.
0.756 (0.040)
0.658 (0.043)
0.725 (0.007)
0.748 (0.041)
0.776 (0.039)
47

---

Published in Transactions on Machine Learning Research (4/2026)
Table 9: Datasets used experiments.
# Feature
# Non-zero values
Tensor size
Sparsity N/121
# Classes
D
N
|221]
ID
AsiaLung
8
40
2.56e+-02
1.56e-01
2
BCW
10
321
1.80c++09
1.78c-07
2
BalanceScale
5
437
1.88c+03
2.33c-01
3
CarEvaluation
7
1209
6.91e+003
1.75e-01
4
Chess
36
6266
1.03e+11
6.08e-08
Chess2
7
19639
1.18c+06
1.66o-02
18
Cleveland
8
139
0.76
2.41e-02
ConfAd
7
129
0.88++03
2.19c-02
2
Coronary
6
61
6.40e+01
9.53e-01
2
Credit
10
309
1.09c+05
2.84c-03
2
DMFT
5
425
2.27e+03
1.87e-01
DTCR
16
179
9.68c+07
1.85c-66
GermanGSS
6
280
0.
3.50e-01
Hayesroth
5
61
02
1.06e-01
Income
12123
1.85e+088
6.56e-05
4
Led7
8
281
003
2.20e-01
10
Lenses
4
8
+-
3.33e-1
3
Lymphography
18
103
0.13
9.10e-07
4
Mofn
11
777
+.
3.79c-01
Monk
7
302
8.64++02
3.50e-01
Mushroom
22
5686
1.
1.17e-10
Nursery
9
9072
6.48e+04
1.40e-01
PPD
9
53
0.
6.82c-03
2
PTumor
16
167
9.83e+04
1.70e03
2
Parity5p5
11
735
00
3.59e-01
2
SPECT
23
163
0.
1.94e-05
2
Sensory
12
403
0.
9.11e-04
SolarFlare
13
416
0.97c+06
6.97c-05
ThreeOfNine
10
358
1.02e+03
3.50e-01
2
Tumor
13
195
1.29e+005
1.51e-03
21
Vehicle
5
33
9.60e+011
3.44e-01
2
Votes
17
123
1.31c+05
9.38c-04
2
XD6
10
455
1.02e+03
4.44e-01
2
Experimental details
D.1 Experimental setup
D.1.1 Experiments for optimization
We downloaded the SIPI House color image of size 512 x 512 x 3 and resized it to 170 x 170 x 3. Each pixel value in the resized tensor is then rescaled to the range [0, 1], and the resulting tensor is normalized to obtain the tensor T. We reconstruct the tensor using a mixture CP, TT, and background term with the CP-rank 35 and TT rank (8,8). The ranks are chosen such that the number of parameters in the CP and TT structures is approximately equal, Since the normalization and non-negative conditions are not satisfied for naive baseline methods, employ parameterivation via softmax functions to satisfy these conditions. As
48

---

Published in Transactions on Machine Learning Research (4/2026)
Table 10: Datasets source and license.
Dataset
License
URL
AsiaLung
BCW
CC BY 4.0
BY http://tttiiittaaa//aa//aaaaaaaa
BalanceScale
CarEvaluation
CC
Chess
Chess2
CC BY 4.0
Cleveland
Public CCO Public
httpe: http:///aoaannaaaaaa////caaaaaaa
ConfAd
Coronary
Credit
CC BY 4.0 http://aaaaalllaaaad///////////r Public hrw.a///aaaaaaaaaaaaa/aaaaaan..n
DMFT
attps: http://////////ooo////////aaaaaa
DTCR
GermanGSS
at http://///aaa//////////o//ooaaaa
Hayesroth
CC http:
Income
CC BY 4.0 ht
Lod7
Public
attps://wwe http://///../////////////////aaa
Lenses
CC BY
bttps /ttp:////aaaaaeeeaa.aaaa///////0
Lymphography
CC BY http:////eoaoiiiiaaaaa/a/a/aaaaa
Mofn
Public
https: http:///////oo...aaa.a///o//asaa
Monk
at tataa///aaaaa../o//////////////r
Mushroom
https: o o
Nursery
CC BY http:////s//oooo//o/////aaaaaaaa
PPD
Public Public MIT
PTumor
https: http:////ccannnnnoooaoaaacaaaaaa
Parity5p5
SPECT
CC BY 4.0
Sensory
Public
https: http://///////////////////aaaaaa
SolarFlare
CC /ttp://///o/aaaa.....o.aaaaa//a0 hutp
ThreeO
Public
Tumor
CC BY http:///arriiiio.c..//a//a////ar
Vehicle
Public
attps: http:////so//ooooooaaaaaaaaaaaaa
Votes
CC BY 4.0 Public
XD6
at http://///a//a///aaooaaaaaeaaaaa
example, we transform the CP model
R Σ
p(CP) tjk
A By,
=
into its softmax representation
PICPI tyk
ur
Xir
Nr
der
En ewi Ep exte
Ey En clur Ele
and optimize each element of w, X, Y and Z without any constraints. We note that the additional parameters w are required for fair comparison with the proposed E algorithm, as the above softmax transformation reduces the number of free parameters by column-wise normalization of the factor matrices, which is an additional constraint on the proposed models. We also conduct a softmax transformation for the TT-structure, We plot the value of the objective function each itcration while varying the value of a. Each method is repeated five times, and all five results are plotted. The learning rate varied from 0.001 10,000 multiplicative steps of 10, and all converged results are plotted
49

---

Published in Transactions on Machine Learning Research (4/2026)
D.1.2 Experiments with synthetic Moon datasets
We use the make moons function in the scikitlearn library (Pedregosa et al. 2011] in Python to obtain 5000 samples of (x.y) coordinates. The noise value set to 0.07. Each sample has a label 2 (1, By discretizing cach sample (z. y), we create the tensor E x0x00 We add mislabeled samples and outliers as follows. First, we randomly select p points from the discretized points (Iz, y1) 2 and flipped the labels. Second, define the region Y A U B LUCU D where
A [0, 20) x 20); B
[70, 90] x [0, 20). C C= [0, 20) [70, 90),
D = [70, 90) (70,90) 90,,
and we randomly sample additional p points in Y, and the label of each is randomly determined. These p samples are regarded as outliers, which are added to the tensor T to obtain Tnoisy We reconstruct the tensor Tholay and obtain mixture low-rank tensor P by E*MCPTТB with CP-rank of 20 and TT-rank of (20, 2), and observe the changes in reconstruction while varying the value of a. In the visualization of reconstructions and empirical distributions, each pixel (i.j) is plotted with color indicating its label assignment: red represents label one and blue represents label two. More specifically, when Paja 0 and Pisa = 1/N, the pixel is colored red, whereas it is colored blue when Piy1 = 1/N and Payz = 0 for the number of samples N. If both Piyi and Pijz are zero, the pixel appears black. When Piyı = Pijz > 0, the pixel is colored purple.
We also evaluated the test error D(T, P) and training error (T P)) using KL divergence (a = 1.0), Hellinger distance (a = 0.5), and reverse KL divergence (a = 0.0), as well as Jensen-Shamie (JS) divergence and Jelfreys divergence, defined as follows:
1 DJs(T.P) 2 DKL (T, RR)
1 (P,R), where R= 2
T+P 2
and
Djer = DKL(T,P) + DKL(P,
respectively.
D.1.3 Density estimation with real datasets
We download categorical tabular datasets from the UCI database OpenML database 5 and Penn Machine Learning Benchmarks (Olson et al. 2017). Since there are missing values in Credit, Lymphography, and Votes, used only the features with missing values in the experiment. Fach dataset contains N categorical .N. samples ((() = ip) for n I...,,N. Both the number of samples, N, and the sample dimension, D, vary across datasets as seen in Table 9 In the original datasets, each i represents a categorical quantity such as color, location, gender, etc. By mapping these to natural numbers, each feature is is converted to a natural number from 1 to Ia, where Ia is the dogree of froedom in the d-th feature. We randomly select 70% of the N samples to create the training index set Dtrain 15% of samples to create the validation index set Qualid and the final 15% of the samples form the test index set Ωtest, Some datasets may contain exactly the same samples. To deal with such datasets, we suppose that these index sets Dtrain Qualid and Qrest may contain multiple identical elements.
We create empirical tensors Ttrain Tvalid, and Tlest, where each value Téd is defined as the number of i in the set o2 for E { train, valid, test ). They are typically very sparse tensors. We normalize each tensor by dividing all the elements by the sum of the tensor to map them to a discrete probability distribution. The above procedure to create empirical tensors consistent with the discussion in Section 1.2
We optimize the a-divergence D. (Ttrann P) for (mixture of) low-rank tensor P by the proposed methods. In the same way, we optimize the KL divergence (T P)) for P by baseline methods. We adjust hyperparameters such as tensor ranks, the value of a, and learning rates to minimize the KT. divergence vaT P). Finally, we evaluate the generalization error by negative likelihood (T P), where the
https://archive.. ics. uci. edi.
https://www.oeenl.cc
50

---

Published in Transactions on Machine Learning Research (4/2026)
tensor P is the reconstruction approximating Terain with tuned rank. We note that optimizing the KL divergence from the test datasets PKI.(Tus) P) is equivalent to optimization for the above negative likelihood.
We further evaluated the quality of the obtained density P by classification task. Specifically, the final feature iD in the density corresponds to the class label of the discrete item (f.....;i--1). Thus, for a given test sample .i....,jj-1), we see the estimated density -ii,;....--1,, for each class label c. The sample is then classified into the class с that maximizes the density. We adjust the hyperparameters to maximize the accuracy of the classification for validation datasets Tvalid.
The proposed method and the baseline method have initial value dependence. Hence, all calculations are repeated five times to evaluate the mean and standard deviation of the negative log-likelihoodd and classification accuracy.
D.1.4 Baseline selection for discrete density estimation with real datasets
Since this study proposes a novel framework for estimating a non-negative low-rank mixture tensor by optimizing the a-divergence, our experiments compare the proposed method with existing non-negative tensor-based approaches for discrete density estimation [Glasser et al. 2019 Tbrahim & Ful 2021 Huang & Sidiropoulos 2017 Yeredor& & Haardy, 2019). The tensor train-based density estimation (TTDE) (Novikvv et al. ) does not optimize the KL, divergence and assumes that the squared value of each element of a real- valued tensor corresponds to a probability value, rather than imposing non-negative constraints on the factors of the docomposition. Therefore, instead of TTDE, we use KL divergenco- basedd basod non-negative decomposition (KLTT) (Glasser et. al. 2019 as a baseline. Although the lensor ring-based density estimation optimizes the KL divergence (Wu et all 2025). its formulation also relies on squared real values, and the source code is not publicly available.
Although we evaluate the obtained density by classification accuracy, our scope is not merely classification but density estimation, which is essentially an unsupervised learning problem. Indeed, the classification accuracy is evaluated based on the ability of conditional distributions to impute class labels, as discussed in Appendix D.1.3 Therefore, comparisons with supervised classification methods are outside the scope of this work. We also note that prior studies (Novikov et al., 2021 Amiridi et al.] 2022 Wu et al., 2025 Loconte
ct al. 2025 have systematically shown that tensor-based continuous density estimation methods perform well
compared to DNN-based continuous density estimatiio (Germain et al. 2015 Dinh et al. 2017 Grathwohl
et al.
2019 Durkan et al. 2019
Notably, existing methods such as a-MU (Kim et al. 2008) and a-EM (Matsuyama 2000), as discussed in Section 1.] do not accommodate mixtures of low-rank decompositions or account for normalization to produce valid densities, making them unsuitable for density estimation. According to the existing literature, EMCP (Huang & Sidiropoulos) 2017 Yeredor & Haardt 2019) is the only EM-based non-negative normalized tensor decomposition reported to date and this method is special case of the proposed framework (with K 1 and a fixed to 1.0). EMCP was employed as baseline. The FM algorithm is also commonly used for tensor completion (Tomasi & Bro 2005 Liu et al. 2015 Song et al. 2019). which assumes the observed tensor T has missing values, such that P is the reconstructed low-rank tensor, optimized considering only indices corresponding to the observed entries of T. However, density estimation does not impose this constraint, and thus these methods are not suitable as basclines. To the best of our knowledge, our method is the first discrete density estimation approach that feasibly optimizes a-divergence and naturally captures high-order interactions among discrete variables through tensor decompositions.
Finally, datasets such as CIFAR (Krizhevsky) 2009), MNIST (Lecun et al. 1998). or ImageNet (Dong et al. 2009) are not appropriate for discrete density estimation because they assume continuous-valued features. Although images are discrete objects, their pixel intensities are integers in the range [0, 255]. not categorical variables. Technically, one could treat cach integer as a distinct category; however, for an I color image, this would imply a sample space of size 25531J, making a continuous formulation far more natural given the data characteristics.
51

---

Published in Transactions on Machine Learning Research (4/2026)
D.2 Implementation details and ppeer-paramtter tuning
We describe the implementation details of the proposed and baseline methods in the following. All experiments are conducted using Python 3.12.3. We provide our source code for all experiments in the supplementary matcrial.
Rank tuning In the experiments for density estimation with real datasets, rank tuning for all methods is carried out according to the following policy. For pure models, we consider eight equally spaced ranks so that the number of parameters in the model does not exceed N/2, where N is the number of samples. To make the tuning easier, we let R1 R2 validation score these candidates. For mixture models, grid search is performed by selecting five of the smallest numbers from the eight candidates for each mixed structure.
Ry- We evaluate the test scores with the rank that maximizes the
Proposed method The pscudocodes for the proposed methods are described in Algorithms 1 and 2 All lensors and mixture ratios were initialized with a uniform distribution bet ween 0 and 1 and normalized as necessary. The algorithm was terminated when the number of iterations of the EM step exceeded 1200 or when the difference of the log-likelihood from the previous itcration was below 10c-8. Hyperparameters a are tuned in the range a E f0.15, 0.5, 0.75, 0.85, 0.9, 0.95, 1.003.
SGD, RMSprop, Adam, Adagrad We used the torch.optim module from PyTorch 2.7 (Paskee et al.) 2017 The learning rate is varied by powers of ten from 0.0001 to 10,000. All other hyperparameters, such
weight: gggh_dccay, are set to their default values. In Figures TO 11 and 12 diverged loss curves were omitted to improve readability.
We downloaded the source code for Positive Matrix Product State (KLTT) and Real Born Machine (BM) from the official repository. ( The license the code is described in the repository as MIT License. They optimize real-valued lensors and square each element to obtain non-negative lensors. Each element of these real-valued tensors is initialized with the standard normal distribution. We varied the learning rates from 1.0c-4, 1.0c-3, ... to 1.0 for the training data. We then used the learning rate that provided the smallest validation score. According to the description in the README file, the batch size was fixed at 20, and the number of iterations was set to 10,000. We performed each experiment. five times with each rank and evaluated the mean and standard deviation. We varied the rank of the model as 1, 2, .... 8, and evaluated the test data with the rank that best fits the validation data.
BM and KLTT
CNMFOPT We implemented CNMFOPT according to the original paper [lbrahim & Fu 2021 , We iteratively optimize the factor matrices A AD)) and the weights À one after the other. We use the exponentiated gradient method (Bubeck 2015) to optimize each factor matrix and weight. Although this optimization be performed by closed-form update rules, all parameters cannot be updated simultancously. Thus, loop is required for each update. This loop is called the inner iteration, which is repeated 100 times. (("),....A"") A A) described above are called outer iteration. In the update rule for the exponential gradient method, the product of the derivative and the learning rate a is included in the exponential function. The learning rate a is selected from 0.005, 0.001, 0.0005, 0.0001, and 0.00005 to minimize the validation error. Learning with a larger learning rate was not feasible because caused an overflow of the exponential function. The initial values of the parameters follow a uniform random distribution. The algorithm terminates when one of the following conditions is met: (1) We compute the KL divergence from the input data to the reconstruction after updating all factor matrices and weights. The change is less than 1.0e- compared to the previous outer loop. (2) The number of outer iterations exceeds 600. When the condition (2) is met, the algorithm performs up to 60000(D) + 1) inner iterations in total.
The iterations of updating (A(1), .....
Bayesian network We used the pomegranate library (version 0.1.4) to learn a Bayesian network using the Chow Liu tree method (Schreiber) 2018). The network was learned automatically from the training data only, assuming no hidden variables, so that the structure maximizes the sum of mutual information
52

---

Published in Transactions on Machine Learning Research (4/2026)
between connected variables, which is equivalent minimizing the KL divergence from the empirical distribution (Chow & Liu, 1968). Thus, we did not use validation data to train the Bayesian network. The algorithm does not depend on initial values and therefore only run once.
Hidden Markov model We evaluated layered single Hidden Markov Model using the CategoricalHMN in hmmlearn library We used the default values of maximum iterations and convergence tolerance. Each feature treated as one time step, with K states per layer. We tume the hyper-parameter K by maximizing K 21, 3, 12, 16, 32, 50, 64, 80) with five random initializations for each K.
validation score
Discrete normalizing flow We downloaded the source code from the official repository of discrete flow Tran et al. 2019) and trained auto-regressive flow models following the demo code provided as Jupyter notebook Training was run for 2,000 epochs on CPU. We performed a grid search over learning rates {0.1, 0.01, 0.001 temperature parameters {0.05, 0.1, 0.15, 0.2), number of flow layers {1, 2, 4}, and hidden units {10, 20, 40, 80), to select the hyperparameter combination that maximized validation performance for each dataset. To account for initialization dependence, each configuration was executed five times, and we chose the setting with the highest mean validation score. Following the default configuration, we used the Adam optimizer in PyTorch, set the batch size equal to the number of training samples, fixed the mumber of layers to three, and initialized the base distribution in the demo code.
D.3 Additional information for reproducibility
The license and source of each dataset used in the experiments are described in Table 10 The imported data were directly converted to tensors by the procedure described in Section D.I Experiments were conducted Ubuntu 20.04.1 with a single core of 2.1GHz Intel Xeon CPU Gold 5218 and 128GB of memory. This work does not require GPU computing. The total computation time for all experiments, including tuning the lcarning rate of the basolincs, was less than 320 hours, using 88 throads of parallel computing. The source code used in experiments provided in the supplementary materials.
53


Reviews:
[{"review_id": "RMUMhjqtXE", "reviewer": "Reviewer_CdUi", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes a new EM-type algorithm for nonnegative tensor factorization under alpha-divergence objectives, with application in tensor-based discrete density estimation. The authors derive two upper bounds on the alpha-divergence objective using Jensen’s inequality, resulting in an algorithm that involves two expectation (E) steps, hence the name E2M. The authors also establish a closed-form update for the M step by formulating the problem as a many-body approximation, which is more efficient and requires less hyperparameter tuning compared to the iterative gradient-based methods used by prior approaches. The authors demonstrate the effectiveness of the proposed method on synthetic and real data and across multiple low-rank structures, including CP, Tucker, Tensor Train, and their mixtures.", "claims_and_evidence_yn": "No", "claims_explanation": "- My major concern is about Theorem 1. The theorem states that the objective always converges. However, the proof only shows that the objective value is non-increasing in each iteration. This is not sufficient, since the objective value may stay constant across iterations, which is not convergence. To state that the objective converges, a sufficient decrease condition is needed. Can the authors clarify this?\n\n\nFurthermore, the clarity of this paper can be improved. The current writing and use of notations can be relatively hard-to-follow for readers not familiar with tensor factorization or many-body approximation. A notable amount of methodological and experimental details seem missing. Please see my specific questions and comments below.\n\n- It would be helpful to clarify why the subproblems in the E1 and E2 steps are equivalent to their respective expectation maximization problems, i.e., to provide more details on the statements \"This subproblem is equivalent to maximizing the expectation...\" above both (9) and (10).\n\n- Can the authors explain why, below (9), the optimal tensor F* satisfies $L_\\alpha (P) = Lbar_\\alpha (F*,P)$? Prop. 1 only shows that $Lbar_\\alpha$ is an upper bound of $L_\\alpha$, but not that it is tight at F*. Similar question regarding $\\Phi^*$ below (10).\n\n- It would be helpful to define the concept of an interaction diagram in Figure 2. Since it plays a key role in many-body approximation, a brief explanation would improve clarity for readers unfamiliar with this matter. It would also be helpful to provide more background on many-body approximation and why (13) is related to this problem.\n\n- The statement \"Above Equations (14)–(16) obviously provide the closed form update for the M-step in Equation (11)\" is unclear. Eq (11) is constrained optimization with feasible set $B^{[k]}$, whereas (14)-(16) are solutions to unconstrained optimization problems. Can the authors clarify how they relate?\n\n- Figure 3a: Can the authors provide more details on the baseline methods? Do they apply gradient-based optimization on the original objective (1)? How are the nonnegativity constraints enforced?\n\n- In Figure 5, it seems that the optimal alpha is mostly dependent on the evaluation metric. This is reasonable as the optimal alpha corresponds to the alpha used in the evaluation metric. However, the text states that the optimal alpha depends on noise level, which is unclear from the figure. Can the authors clarify this?", "audience_yn": "Yes", "audience_explanation": "The paper presents a new approach for tensor factorization that could have broad applications.", "requested_changes": "- Please address the major concern on the proof for Theorem 1.\n- Please consider improving the clarity of the paper, regarding the questions above.\n\n\nOther minor comments:\n- Section 2.2, Paragraph \"Adaptive background term for stable learning\": Can the authors clarify how Algorithm 1 should be modified to incorporate the background term?\n\n- (14) (15) (16): Can the authors clarify if the division is element-wise?\n\n- Section 3.2: Since the text only discussed the computation of M step, does the iteration complexity shown here (and in Table 2) include both the E^2 and M steps, or just the M step?\n\n- Figure 3c only provides the result of the proposed method. Can the authors clarify why other methods were excluded?\n\n- Figure 3a: Multiple curves are shown for each method. Can the authors clarify what each curve represents?\n\n- Section 4, i) and Figure 3 a-c: These figures compare the result of different methods at the same number of iterations. Can the authors also comment on how the **per-iteration complexity** compares among different methods? Does the proposed method have the lowest per-iteration complexity as well?\n\n- Section 4, ii): It would be helpful to refer to the relevant figure in the main body.\n\n- \"The 25 randomly selected class labels are flipped, and a further 25 random outliers are added to T and obtain T noisy.\"\nCan the authors clarify how the \"outliers\" are constructed?", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "E4Y1oOC6Ws", "reviewer": "Reviewer_dknd", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper proposes the E2M algorithm as a generalization of the standard EM algorithm. It is proposed to optimize the alpha divergence in tensor-based discrete density estimation.  The method is evaluated on the synthetic and real datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The authors derive a double upper bound on the alpha divergence objective using Jensen's inequality. The first bound relaxes the alpha divergence minimization into a surrogate KL-divergence minimization. The second bound further relaxes it into a many-body approximation. \n\nThe approach supports different tensor decompositions. The paper provides theoretical derivations to support the claims. \n\nThe empirical claims like robustness and perforamnce on real data are not fully verified from the paper.", "audience_yn": "Yes", "audience_explanation": "The machine learning researchers focused on probabilistic modeling, optimization algorithms, tensor decompositions, or density estimation would likely be interested in the findings of the paper.", "requested_changes": "Provide a table that explains the notations and terminology. Make the paper more readable.\n\nAdd complexity analysis and optimization for large-scale tensors.\n\nInclude ablation studies on the impact of double bounds. \n\nCompare to non-tensor methods and discuss when E2M is preferable.", "broader_impact_concerns": "Tensor-based density estimation could be applied to sensitive categorical data.\n\nLow-rank approximations might amplify biases in datasets.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "aMUSk82ddg", "reviewer": "Reviewer_Gpwu", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors present E$^2$M, a generalized Expectation-Maximization algorithm designed for discrete density estimation using non-negative tensor factorization. The paper addresses the difficulty of optimizing the $\\alpha$-divergence (which generalizes KL, Reverse-KL, and Hellinger distances) in tensor models, where non-linearities usually prevent closed-form updates. The key contributions include a theoretical framework that uses a double-bound strategy, derivation of closed-form M-step updates for CP, Tucker, and Tensor Train (TT) decompositions, and flexible mixture modeling.\n\nKey strengths:\n1. The double-bound derivation cleverly decouples the mode interactions, allowing for the application of convex optimization results to complex tensor networks.\n2. By avoiding gradient-based optimization, the method guarantees monotonic convergence and eliminates the hyperparameter sensitivity.\n3. The ability to seamlessly mix CP, Tucker, and TT structures within a single probabilistic framework is a significant advancement over prior works.\n\nWeaknesses:\n1. The mathematical formulation, particularly regarding the indices for the many-body approximations, is extremely dense and required significant effort to parse.\n2. The closed-form update for the Tucker decomposition involves a dense core tensor, scaling exponentially with the dimension, which may limit its applicability to low-order tensors compared to CP or TT.\n3. The evaluation is restricted to other tensor-based methods. While valid for the paper's specific scope, a comparison against standard discrete density estimators like autoregressive neural networks would provide better context for the method's general utility.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The claims are well-supported by both rigorous theoretical derivations that prove the closed-form updates and monotonic convergence, and extensive empirical evidence. The experiments convincingly demonstrate the method's superior stability over gradient-based baselines by eliminating learning rate tuning and its ability to handle outliers by tuning the $\\alpha$-parameter.", "audience_yn": "Yes", "audience_explanation": "The paper addresses fundamental challenges in probabilistic modeling and non-convex optimization, specifically solving the intractability of \n$\\alpha$-divergence minimization for tensors without requiring gradient descent tuning. This is of significant interest to researchers working on tensor methods, unsupervised learning, robust estimation, and optimization stability.", "requested_changes": "1. Improve Derivation Accessibility (Strengthen): The mathematical derivations for the closed-form updates of Tucker and TT models involving scaling redundancy removal are dense and difficult to parse. Please include a simplified schematic, flow-chart, or a high-level \"walkthrough\" of the variable dependencies and update steps. This would significantly aid readers in reproducing the implementation.\n\n2. Contextualize Baselines (Strengthen): While the comparison to gradient-based tensor methods is exhaustive, the paper would benefit from a broader discussion of the density estimation landscape. Please add a brief section (in the main text or Appendix) contrasting tensor-based approaches with:\n   - Autoregressive Neural Networks: Explain the trade-offs, e.g., interpretability, exact marginals vs. expressivity, to justify the focus on tensors.\n   - Probabilistic Tensor Factorization: Briefly mention how E$^2$M relates to or differs from broader Bayesian tensor decomposition methods beyond just EMCP.\n\n3. Hyperparameter Guidance (Strengthen): The paper demonstrates that $\\alpha$ controls robustness, but practical guidance is limited. Please provide a brief discussion or heuristic on how a practitioner might select a starting $\\alpha$ based on dataset characteristics, e.g., sparsity, expected noise level, if extensive validation sets are not available.\n\n4. Initialization Sensitivity (Strengthen): Since tensor factorization optimizes a non-convex objective, please briefly comment on the sensitivity of E$^2$M to initialization compared to the baselines. Does the method exhibit lower variance in final performance across the 5 random seeds compared to gradient-based methods like Adam?\n\n5. Scalability Clarification (Strengthen): Regarding Algorithm 2 (Scalable TT), please clarify if the \"merging cores\" strategy introduces significant intermediate memory overhead for large ranks, effectively trading compute for memory, or if the impact is negligible.", "broader_impact_concerns": "No broader impact concerns. The paper presents a fundamental optimization algorithm for density estimation. This specific algorithmic contribution does not introduce new or unique ethical risks that would necessitate a dedicated impact statement. In fact, by eliminating the need for extensive hyperparameter tuning, the method likely reduces the computational energy footprint compared to baseline gradient-based approaches.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)