figuret SMReferences for Supplementary Materials

##

Kazu Ghalamkari
RIKEN AIP
<EMAIL_ADDRESS>
Jesper Løve Hinrich
Technical University of Denmark
<EMAIL_ADDRESS>
&Morten Mørup
Technical University of Denmark
<EMAIL_ADDRESS>

## Non-negative Tensor Mixture Learning
for Discrete Density Estimation

Kazu Ghalamkari
RIKEN AIP
<EMAIL_ADDRESS>
Jesper Løve Hinrich
Technical University of Denmark
<EMAIL_ADDRESS>
&Morten Mørup
Technical University of Denmark
<EMAIL_ADDRESS>

###### Abstract

top.-1top.-1Tensor Mixture LearningTensor Mixture Learningtop.-1 We present an
expectation-maximization (EM) based unified framework for non-negative tensor
decomposition that optimizes the Kullback-Leibler divergence. To avoid
iterations in each M-step and learning rate tuning, we establish a general
relationship between low-rank decomposition and many-body approximation. Using
this connection, we exploit that the closed-form solution of the many-body
approximation can be used to update all parameters simultaneously in the
M-step. Our framework not only offers a unified methodology for a variety of
low-rank structures, including CP, Tucker, and Train decompositions, but also
their combinations forming mixtures of tensors as well as robust adaptive
noise modeling. Empirically, we demonstrate that our framework provides
superior generalization for discrete density estimation compared to
conventional tensor-based approaches.

###  1 Introduction

In machine learning, a key aim of unsupervised learning is to perform accurate
density estimation of high-dimensional data. The estimated density can then be
used for multiple purposes such as predicting new data points, inferring
missing values, performing outlier detection, or generating new data points
for data augmentation. Density estimation has typically been addressed by the
use of parametric and non-parametric methodologies, including mixture modeling
approaches and kernel density estimation procedures [22]. Recently, deep
learning-based models have been found to accurately characterize the
underlying data distribution by use of variational auto-encoders, auto-
regressive, as well as flow and diffusion-based generative models; for an
overview, see [26].

In this paper, we focus on discrete density estimation which seeks to
characterize the joint distribution $P(x_{1},...,x_{D})$ where each discrete
random variable $x_{d}\in\\{1,2,\dots,I_{d}\\}$ corresponds to a categorical
feature. Such discrete distributions are naturally identical to normalized
non-negative _tensors_ , or multidimensional arrays. A series of recent
studies estimate the discrete distribution behind categorical data [10, 19] by
approximating a tensor constructed from the observed data, i.e., the
_empirical tensor_ , with a low-rank tensor described by a linear combination
of a small number of bases; however, two challenges remain in these current
works.

The first challenge is to develop a unified formulation of nonnegative tensor
decomposition that works with various kinds of low-rank structures, optimizing
the Kullback–Leibler (KL) divergence, which is a natural measure of similarity
between probability distributions. In contrast to the well-established SVD-
based methods for real- and complex-valued tensor networks [5, 13, 20], a
general framework for nonnegative tensor decompositions optimizing the KL
divergence is not well developed. Consequently, some existing studies of
tensor-based discrete density estimation have been performed via optimization
of the Frobenius norm, which imposes unnatural modeling assumptions when
applied to probability distributions [19]. In addition, due to the lack of a
unified KL-divergence-based framework, users can try only popular low-rank
structures that have been developed piecemeal or differentiate the cost
function for each low-rank structure by themselves. A framework that allows
users to try various low-rank structures more intuitively and freely is,
therefore, desirable.

The second challenge is scalability, as the size of the tensor — which
corresponds to the sample space size of the discrete distribution — increases
exponentially with the number of features $D$, namely a tensor has $I^{D}$
elements for the degree of freedom of each feature $I$. However, the number of
samples available for training $N$ — typically indicates the number of nonzero
values of the empirical tensor — is often limited, also deemed the curse of
dimensionality. Therefore, it is desirable to develop scalable tensor
decomposition methods that can operate well for high-dimensional data despite
a limited number of samples available.

To address these two challenges, this paper proposes a unified non-negative
tensor factorization method based on the expectation-maximization (EM)
algorithm [7]. The EM algorithm is an iterative framework for maximum
likelihood estimation that repeatedly maximizes the lower bound of the log-
likelihood function in two steps, the E-step and the M-step. The naive
formulation of the EM-based general non-negative tensor decomposition requires
either an alternating optimization bounding the log-likelihood (E-step) and
maximizing (M-step) each factor at a time or gradient methods in each M-step.
It requires either tailoring methods to each type of decomposition or time-
consuming double iterative methods. To overcome this issue, we exploit that
the optimization required in each M-step coincides with a tensor many-body
approximation [9] that decomposes tensors by a reduced representation of the
interactions among tensor modes. We derive the exact closed-form solution for
many-body approximation that appears in EM algorithms for Tucker and Train
decomposition, and thereby successfully remove the gradient method in M-step
for various kinds of low-rank approximation by combining these two closed-form
formulas.

Our framework inherits the properties of the EM algorithm [27], and therefore,
it always converges regardless of the choice of the low-rank structure assumed
in the model. Furthermore, the tensor to be approximated in the M-step is
sparse as it is defined in terms of elementwise multiplications by the
empirical tensor. As a result, the computational complexity of the proposed
method is proportional to the number of observations $N$. Notably, as the EM
algorithm is frequently used for maximum likelihood estimation of mixture
models, the proposed method allows for density estimation using mixtures of
low-rank tensor decompositions providing a flexible modeling framework that
has not been studied previously. We empirically show that learning such tensor
mixtures provides better generalization than pure low-rank tensor models for
discrete density estimation. We summarize our contribution as follows:

  * •

We reveal a relationship between tensor many-body approximation and
traditional low-rank approximation.

  * •

Using this relationship, we provide a unified EM-based framework for non-
negative low-rank approximation optimizing the KL-divergence, notably
providing simultaneous closed-form updates for all parameters in the M-step
while exploiting the sparsity structure of the empirically estimated joint
distribution.

  * •

Based on the proposed framework, we develop a mixture of low-rank tensor
modeling approach that empirically demonstrates inferential robustness and
improved generalization.

###  2 Low-rank approximation and many-body approximation

Figure 1: (a,b,c) Interaction diagrams for tensors
$\mathcal{Q}^{\mathrm{CP}}$,$\mathcal{Q}^{\mathrm{Tucker}}$, and
$\mathcal{Q}^{\mathrm{Train}}$, respectively, in Equation (1). Each node
represents a tensor mode, and the black square, $\blacksquare$, represents the
interaction between modes. (d,e,f) Tensor networks for
$\mathcal{P}^{\mathrm{CP}}$, $\mathcal{P}^{\mathrm{Tucker}}$, and
$\mathcal{P}^{\mathrm{Train}}$, respectively, in Equation (6). Nodes represent
factor tensors, and edges represent mode products. The gray-filled $\Lambda$
is a tensor whose hyper-diagonal elements are $1$ otherwise $0$. These low-
rank tensors correspond to marginalizations of $\boldsymbol{r}$ modes in the
many-body approximation.

The naive direct EM-based formulation for non-negative low-rank approximation
bounding the likelihood in the E-step and optimizing all parameters in the
M-step relies on an iterative gradient method in the M-step. We presently
explore the relationship between low-rank decomposition and many-body
approximation and leverage that the many-body approximation provides exact
closed-form solutions to eliminate the gradient method in the M-step. We
describe in Section 2.2 that the marginalization along the tensor modes
connects many-body approximation to low-rank decomposition.

####  2.1 Many-body approximation with exact closed-form solution

Many-body approximation factorizes tensors by given interactions among tensor
modes. These interactions are described through an interaction diagram, as
shown in Figure 1. Each node (circle) corresponds to a tensor index and edges
through a black square, $\blacksquare$, mean the existence of interaction. In
the following, we denote by $[d]=\\{1,2,\dots,d\\}$ for a positive integer
$d$.We provide three examples, a ($D+1$)-th order tensor
$\mathcal{Q}^{\mathrm{CP}}\in\mathbb{R}_{\geq 0}^{I_{1}\dots\times\dots
I_{D}\times R}$, a $2D$-th order tensor
$\mathcal{Q}^{\mathrm{Tucker}}\in\mathbb{R}^{I_{1}\times\dots\times
I_{D}\times R_{1}\times\dots\times R_{D}}_{\geq 0}$, and a ($2D-1$)-th order
tensor $\mathcal{Q}^{\mathrm{Train}}\in\mathbb{R}^{I_{1}\times\dots\times
I_{D}\times R_{1}\times\dots\times R_{D-1}}_{\geq 0}$. When their interactions
are described as Figure 1(a), (b), and (c) respectively, they can be
factorized as

$\displaystyle\mathcal{Q}^{\mathrm{CP}}_{i_{1}\dots
i_{D}r}=\prod_{d=1}^{D}A^{(d)}_{i_{d}r},\quad\mathcal{Q}^{\mathrm{Tucker}}_{i_{1}\dots
i_{D}r_{1}\dots r_{D}}=\mathcal{G}_{r_{1}\dots
r_{D}}\prod_{d=1}^{D}A^{(d)}_{i_{d}r_{d}},\quad\mathcal{Q}^{\mathrm{Train}}_{i_{1}\dots
i_{D}r_{1}\dots
r_{D-1}}=\prod_{d=1}^{D}\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}},$ (1)

where $i_{d}\in[I_{d}]$ and $r_{d}\in[R_{d}]$ for $d=1,\dots,D$. For
simplicity, we suppose $r_{0}=r_{D}=1$ for $\mathcal{Q}^{\mathrm{Train}}$. Let
tensor indices $i_{1},\dots,i_{D}$ and $r_{1},\dots,r_{V}$ be $\boldsymbol{i}$
and $\boldsymbol{r}$, respectively. Here, the integer $V$ is $1$ for
$\mathcal{Q}^{\mathrm{CP}}$, $D$ for $\mathcal{Q}^{\mathrm{Tucker}}$, and
$D-1$ for $\mathcal{Q}^{\mathrm{Train}}$. We denote the domains of
$\boldsymbol{i}$ and $\boldsymbol{r}$ by $\Omega_{I}$ and $\Omega_{R}$,
respectively, i.e.,
$\boldsymbol{i}\in\Omega_{I}=[I_{1}]\times\dots\times[I_{D}]$ and
$\boldsymbol{r}\in\Omega_{R}=[R_{1}]\times\dots\times[R_{V}]$. The symbol
$\Omega$ with upper indices refers to the index set for all indices other than
the upper indices, e.g.,

$\displaystyle\Omega_{I}^{\backslash d}$
$\displaystyle=[I_{1}]\times\dots\times[I_{d-1}]\times[I_{d+1}]\times\dots\times[I_{D}],$
$\displaystyle\Omega_{R}^{\backslash d,d-1}$
$\displaystyle=[R_{1}]\times\dots\times[R_{d-2}]\times[R_{d+1}]\times\dots\times[R_{V}].$

Many-body approximation parameterizes tensors as discrete distributions, where
the random variables correspond to the tensor modes, and the sample space
corresponds to the index set of a tensor. Maximum likelihood estimation then
finds the globally optimal tensor that minimizes the cross-entropy from the
data tensor $\mathcal{M}$ to the tensor in the model space
$\boldsymbol{\mathcal{B}}$, which is the set of tensors with specific
interactions. Thus, for a given tensor $\mathcal{M}$, the many-body
approximation based on the above three interactions maximizes

$\displaystyle
L_{\mathrm{MBA}}(\mathcal{M};\mathcal{Q}^{k})=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{k}},\quad
k\in\\{\mathrm{CP},\ \ \mathrm{Tucker},\ \ \mathrm{Train}\\}.$ (2)

This optimization is equivalent to optimizing the KL divergence from tensor
$\mathcal{M}$ to $\mathcal{Q}^{k}$ and is guaranteed to be a convex problem
regardless of the choice of interaction. While the conventional method finds a
numerical solution of general many-body approximation by the Natural gradient
method, it has been shown in [12] that the following factors globally maximize
$L_{\mathrm{MBA}}(\mathcal{M};\mathcal{Q}^{\mathrm{CP}})$:

$\displaystyle
A^{(d)}_{i_{d}r}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}r}}{\mu^{1/D}\left(\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}r}\right)^{1-{1/D}}},\quad\mu=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{r\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}r},$
(3)

which is consistent with mean-field approximation [8]. As a generalization of
the above result, we presently provide the optimal solution of the many-body
approximation of $\mathcal{Q}^{\mathrm{Tucker}}$ and
$\mathcal{Q}^{\mathrm{Train}}$ in closed-form. The following tensors globally
maximizes $L_{\mathrm{MBA}}(\mathcal{M};\mathcal{Q}^{\mathrm{Tucker}})$:

$\displaystyle\mathcal{G}_{\boldsymbol{r}}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}},\quad
A^{(d)}_{i_{d}r_{d}}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}},$ (4)

and following tensors globally maximizes
$L_{\mathrm{MBA}}(\mathcal{M};\mathcal{Q}^{\mathrm{Train}})$:

$\displaystyle\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d,d-1}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}},$ (5)

where we assume $r_{0}=r_{D}=1$. We formally derive the above closed-form
solutions in Theorems 2 and 3 in the supplementary material. Notably, we can
obtain exact solutions also for more complicated many-body approximations by
combining these solutions, which we discuss in Section 3.3.

####  2.2 Many-body approximation and low-rank approximation

When summing the tensors in Equation (1) over indices $\boldsymbol{r}$, these
models are identical to the traditional low-rank models based on CP, Tucker,
and Tensor Train formats, respectively.

$\displaystyle\mathcal{P}^{\mathrm{CP}}_{\boldsymbol{i}}=\sum_{r}\mathcal{Q}^{\mathrm{CP}}_{\boldsymbol{i}r},\quad\mathcal{P}^{\mathrm{Tucker}}_{\boldsymbol{i}}=\sum_{r_{1}\dots
r_{D}}\mathcal{Q}^{\mathrm{Tucker}}_{\boldsymbol{i}r_{1}\dots
r_{D}},\quad\mathcal{P}^{\mathrm{Train}}_{\boldsymbol{i}}=\sum_{r_{1}\dots
r_{D-1}}\mathcal{Q}^{\mathrm{Train}}_{\boldsymbol{i}r_{1}\dots r_{D-1}}.$ (6)

Since many-body approximation treats the tensor indices as discrete random
variables, we consider these low-rank tensors as models in which the random
variables $\boldsymbol{i}$ and $\boldsymbol{r}$ represent visible and hidden
variables, respectively. The degree of freedom of hidden variables
$(R_{1},\dots,R_{V})$ corresponds to _CP rank_ , _Tucker rank_ , and _train
rank_ with $V=1$ for $\mathcal{Q}^{\mathrm{CP}}$, $V=D$ for
$\mathcal{Q}^{\mathrm{Tucker}}$, and $V=D-1$ for
$\mathcal{Q}^{\mathrm{Train}}$, respectively. Since any low-rank factorization
decomposes a tensor by summing over its ranks, _any low-rank approximation can
be regarded as a many-body approximation with hidden variables_. When a low-
rank tensor $\mathcal{P}$ is obtained by marginalization of a tensor
$\mathcal{Q}$ with appropriately selected modes, we refer to tensor
$\mathcal{Q}$ as the _low-body tensor_ corresponding to $\mathcal{P}$. For
example, $\mathcal{Q}^{\mathrm{CP}}$, $\mathcal{Q}^{\mathrm{Tucker}}$ and
$\mathcal{Q}^{\mathrm{Train}}$ are low-body tensors corresponding to a low CP-
rank tensor, low Tucker-rank tensor, and low Train-rank tensor, respectively.
In this paper, hidden variables — corresponding to modes summed according to
the low-rank structure — are denoted by $\boldsymbol{r}$ and visible variables
— variables other than hidden variables — by $\boldsymbol{i}$.

The many-body approximation involves only visible variables in the model,
making maximum likelihood estimation a convex optimization problem, while a
low-rank approximation becomes a nonconvex optimization problem due to the
marginalized hidden variables in the model. Therefore, although finding an
exact solution for the low-rank approximation is challenging, the optimization
remains tractable through the EM algorithm [7], a well-known general framework
for maximum likelihood estimation that accommodates hidden variables.

Whereas it is already known that many-body approximation coincides with the
non-negative tensor network using Dirac tensors on the nodes [9], it is our
contribution to clarify the relationship between general low-rank
approximation and many-body approximation.

###  3 Discrete density estimation via non-negative tensor learning

Based on the exact solutions of many-body approximation, we develop a novel
framework _EM non-negative tensor learning_ for discrete density estimation.
We construct a normalized histogram, or empirical distribution, from tabular
data with $D$ categorical features. This histogram is identical to a
normalized $D$-th order tensor $\mathcal{T}\in\mathbb{R}_{\geq
0}^{I_{1}\times\dots\times I_{D}}$ where $I_{d}$ is the degree of freedom
(i.e., categories) of the $d$-th categorical feature. The number of non-zero
values in the tensor is equivalent to the number of observed samples $N$ if
there is no identical sample in the data. To estimate the discrete probability
distribution underlying the data, we approximate such a tensor $\mathcal{T}$
with a low-rank mixture $\mathcal{P}$, which is a convex linear combination of
$K$ low-rank tensors.

Our framework has two advantages: (1) it achieves linear computational
complexity relative to the number of nonzero elements in the input tensor, as
described in Section 3.2; and (2) it has the flexibility to generalize
complicated low-rank structures and adaptive regularization terms without
losing the convexity of the E and M-step, as detailed in Sections 3.3 and 3.4.
Neither of these advantages has been explored in previous studies [10, 12].

####  3.1 EM-method for mixture of low-rank approximations

Here, we provide a unified EM-based method for a mixture of low-rank
approximations. We maximize the negative cross-entropy from a given $D$-th
order tensor $\mathcal{T}$ to a mixture of $K$ low-rank tensors
$\mathcal{P}^{1},\dots,\mathcal{P}^{K}$,

$\displaystyle
L(\hat{\mathcal{Q}})=\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{T}_{\boldsymbol{i}}\log{\sum_{k=1}^{K}\eta^{k}\mathcal{P}^{k}_{\boldsymbol{i}}},\quad\quad\mathcal{P}^{k}_{\boldsymbol{i}}=\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{Q}^{k}_{\boldsymbol{i}\boldsymbol{r}},$
(7)

where the mixture ratio satisfies $\sum_{k}\eta^{k}=1$ and $\eta^{k}\geq 0$
for all $k\in[K]$ and setting $K=1$ yields a conventional low-rank tensor
decomposition. Each low-rank tensor $\mathcal{P}^{k}$ is normalized, i.e.,
$\sum_{\boldsymbol{i}}\mathcal{P}_{\boldsymbol{i}}^{k}=1$. We let the model
$\mathcal{P}$ be the convex linear combination of low-rank tensors, that is,
$\mathcal{P}=\sum_{k}\eta^{k}\mathcal{P}^{k}$. We denote the number of hidden
variables in the tensor $\mathcal{Q}^{k}$ by $V^{k}$, and
$\boldsymbol{r}\in\Omega_{R^{k}}=[R^{k}_{1}]\times\dots\times[R^{k}_{V^{k}}]$.
For simplicity, we introduce
$\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}=\eta^{k}\mathcal{Q}^{k}_{\boldsymbol{i}\boldsymbol{r}}$,
and refer to the tensors $\hat{\mathcal{Q}}^{1},\dots,\hat{\mathcal{Q}}^{K}$
as $\hat{\mathcal{Q}}$. We apply Jensen’s inequality [14] to the objective
function $L(\hat{\mathcal{Q}})$ in order to move the summation over hidden
variables $\boldsymbol{r}$ outside the logarithm function thereby obtaining
the lower bound,

$\displaystyle L(\hat{\mathcal{Q}})$
$\displaystyle\geq\overline{L}(\hat{\mathcal{Q}},\Phi)=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}\log{\frac{\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}}}$
(8)

for any $K$ tensors $\Phi=(\Phi^{1},\dots,\Phi^{K})$ where $\Phi^{k}$ is a
$(D+V^{k})$-th order tensor satisfying
$\sum_{k}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}=1$.
The derivation of the inequality is described in Proposition 1 in the
supplementary material. The above lower bound can be decoupled into
independent multiple many-body approximations for tensors
$\mathcal{M}^{1},\dots,\mathcal{M}^{K}$ where each tensor $\mathcal{M}^{k}$ is
defined as
$\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}^{k}=\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}$,
and an optimization problem for the mixture ratio
$\eta=(\eta^{1},\dots,\eta^{K})$. More specifically, we decouple the lower
bound as follows:

$\displaystyle\overline{L}(\hat{\mathcal{Q}},\Phi)=\sum_{k=1}^{K}L_{\mathrm{MBA}}(\mathcal{M}^{k};\mathcal{Q}^{k})+J(\eta),\quad
J(\eta)=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}\log{\eta^{k}},$
(9)

where the objective function of many-body approximation $L_{\mathrm{MBA}}$ is
introduced in Equation (2). The EM algorithm iteratively optimizes the lower
bound for tensors $\Phi$ during the E-step and tensors $\hat{\mathcal{Q}}$
during the M-step until convergence. Each step is a convex optimization while
Equation (7) is a non-convex function. This procedure is guaranteed to
converge, and each iteration increases the objective function (7)
monotonically, which we prove in Theorem 4 in the supplementary material.

E-step  We maximize the lower bound $\overline{L}(\hat{\mathcal{Q}},\Phi)$ for
$\Phi$, that is,

$\displaystyle\Phi=\operatorname*{arg\,max}_{\Phi\in\boldsymbol{\mathcal{D}}}\overline{L}(\hat{\mathcal{Q}},\Phi),$

where the solution space $\boldsymbol{\mathcal{D}}$ is a set of $K$ tensors
such that each tensor is normalized by hidden variables, i.e.,
$\boldsymbol{\mathcal{D}}=\left\\{\left(\Phi^{1},\dots,\Phi^{K}\right)\mid\sum_{k}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}=1\right\\}$.
The optimal solution $\Phi^{k}$ is given as

$\displaystyle\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}=\frac{\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}},$
(10)

as shown in Proposition 2 in the supplementary material. The denominator in
(10) is equivalent to the definition of the model $\mathcal{P}$.

M-step  We maximize the lower bound for tensors
$\mathcal{Q}=(\mathcal{Q}^{1},\dots,\mathcal{Q}^{K})$ and non-negative values
$\eta$. Since the lower bound can be decoupled as shown in Equation (9), the
required optimizations in the M-step are as follows:

$\displaystyle\mathcal{Q}^{k}=\operatorname*{arg\,max}_{\mathcal{Q}^{k}\in\boldsymbol{\mathcal{B}}^{k}}L_{\mathrm{MBA}}(\mathcal{M}^{k};\mathcal{Q}^{k}),\quad\eta=\operatorname*{arg\,max}_{0\leq\eta^{k}\leq
1,\sum_{k}\eta^{k}=1}J(\eta),\quad$ (11)

where the model space $\boldsymbol{\mathcal{B}}^{k}$ is the set of low-body
tensors corresponding to the $k$-th low-rank tensor $\mathcal{P}^{k}$. The
naive implementation of the M-step requires a gradient method to solve the
many-body approximation for $\mathcal{M}^{k}$. The existing algorithm for
many-body approximation is based on the Natural gradient method [3], which
requires cubic computational complexity for the number of parameters in the
low-body tensor $\mathcal{Q}^{k}$. Thus, repeating the gradient method in each
M-step is computationally expensive. However, we introduced the closed-form
solution of the many-body approximation in Section 2.1, which eliminates the
iterative gradient method in the M-step for Tucker and Train decomposition. We
discuss more complicated low-rank structures in Section 3.3.

We also provide optimal updates for the mixture ratio $\eta$. By the condition
$\partial J(\eta)/\partial\eta^{k}=0$ and the normalizing condition
$\sum_{k}\eta^{k}=1$, the optimal $\eta^{k}$ is simply given as follows:

$\displaystyle\eta^{k}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}^{k}}{\sum_{k=1}^{K}\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}^{k}}$
(12)

which are shown in Proposition 3 in the supplementary material. It is
straightforward to check that the M-step is a convex optimization problem
whatever the low-rank structure assumed on $\mathcal{P}^{k}$, given that we
can decouple the lower bound into multiple independent convex many-body
approximations.

While the EM-based method for CP decomposition optimizing KL divergence has
already been developed in [12], our proposed framework addresses more general
low-rank decompositions and their mixtures. We provide the entire algorithm in
Algorithm 1. Notably, the framework does not require a learning rate that
needs to be carefully tuned when relying on gradient-based methods.

##### EM tensor factorization for general non-negative tensors

Non-negative tensor factorization optimizing KL divergence or cross-entropy is
frequently used beyond density estimation and in various fields such as sound
source separation [17], computer vision [15, 21], and data mining [6, 24].
Although the given tensor $\mathcal{T}$ is not necessarily normalized in such
applications, EM tensor factorization can be used for them as follows. First,
we obtain the total sum of the input tensor
$\mu=\sum_{\boldsymbol{i}}\mathcal{T}_{\boldsymbol{i}}$ and then perform EM
tensor factorization on the normalized tensor $\mathcal{T}$ by dividing all
elements of $\mathcal{T}$ by $\mu$. Finally, all elements of the resulting
mixture low-rank tensor are multiplied by $\mu$. This procedure is justified
by the property of the KL divergence,
$D_{KL}(\mu\mathcal{P},\mu\mathcal{T})=\mu D_{KL}(\mathcal{P},\mathcal{T})$,
where $\mu$ is any non-negative value.

input : Non-negative tensor $\mathcal{T}$, the number of mixtures $K$, and
ranks $(R^{1},\dots,R^{K})$

Initialize $\mathcal{Q}^{k}$, $\eta^{k}$ for all $k\in[K]$;

repeat

$\mathcal{P}_{\boldsymbol{i}}\leftarrow\sum_{k}\eta^{k}\mathcal{P}_{\boldsymbol{i}}^{k}$
where $\mathcal{P}^{k}_{\boldsymbol{i}}=\sum_{\boldsymbol{r}\in
R^{k}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{k}$;

$\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}^{k}\leftarrow\mathcal{T}_{\boldsymbol{i}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{k}/\mathcal{P}_{\boldsymbol{i}}$
for all $k\in[K]$; // E-step

Update tensor $\mathcal{Q}^{k}$ for all $k\in[K]$ using Equations (1), (4),
and (5); // M-step

Update mixture ratio $\eta^{k}$ using Equation (12) for all $k\in[K]$; //
M-step

until  _Convergence_ ;

return _$\mathcal{P}$_

Algorithm 1 Non-negative Tensor Mixture Learning

####  3.2 Analysis of computational complexity

In the following, we discuss the computational complexity of the proposed EM
low-rank approximation without a mixture, specifically when $K=1$. It is
straightforward to see that the computational complexity for $K>1$ is simply
the sum of the complexity of each of the EM low-rank approximations used in
the mixture.

In our proposed framework, we compute the sum over the visible variables
$\boldsymbol{i}\in\Omega_{I}$ of the tensor $\mathcal{M}^{k}$ in the M-step,
using Equations (4), (5), and (12) for which the computational cost increases
rapidly as the size of the tensor increases. However, if the input tensor
$\mathcal{T}$ is sparse, which is a reasonable assumption for most density
estimation tasks, we can reduce the computational complexity of the summation
because the tensor $\mathcal{M}^{k}$ is also sparse for indices
$\boldsymbol{i}$ as the definition of
$\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}^{k}=\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}$.
More specifically, we replace the sum over the visible variables
$\sum_{\boldsymbol{i}\in\Omega_{I}}$ with
$\sum_{\boldsymbol{i}\in\Omega^{\mathrm{o}}_{I}}$ where
$\Omega^{\mathrm{o}}_{I}$ is the set of indices of nonzero values of the
tensor $\mathcal{T}$. The resulting time computational complexity is $O(\gamma
NDR)$ for EM-CP decomposition and $O(\gamma DNR^{D})$ for EM-Tucker and EM-
Train decomposition where $N$ is the number of nonzero values, $D$ is the
number of modes in the tensor, $\gamma$ is the number of iteration, assuming
ranks are $(R,\dots,R)$ for all low-rank models.

Furthermore, the computational cost of the tensor train decomposition can be
reduced to $O(\gamma DNR^{2})$ by computing the sum over the latent variables
$\boldsymbol{r}\in\Omega_{R}$ as follows. Firstly, we introduce the following
tensors

$\displaystyle\mathcal{G}^{\left(\rightarrow
d\right)}_{i_{1},\dots,i_{d},r_{d}}=\sum_{r_{d-1}}\mathcal{G}^{\left(\rightarrow
d-1\right)}_{i_{1},\dots,i_{d-1},r_{d-1}}\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}},\quad\mathcal{G}^{(d\leftarrow)}_{i_{d+1},\dots,i_{D},r_{d}}=\sum_{r_{d+1}}\mathcal{G}^{(d+1)}_{r_{d}i_{d+1}r_{d+1}}\mathcal{G}^{(d+1\leftarrow)}_{i_{d+2},\dots,i_{D},r_{d+1}}$
(13)

with $\mathcal{G}^{(\rightarrow 1)}=\mathcal{G}^{(1)}$,
$\mathcal{G}^{(D-1\leftarrow)}=\mathcal{G}^{(D)}$, and
$\mathcal{G}^{(\rightarrow 0)}=\mathcal{G}^{(D\leftarrow)}=1$. Each complexity
is $O(R_{d})$ to get $\mathcal{G}^{(\rightarrow d)}$ and
$\mathcal{G}^{(d\leftarrow)}$ when we compute them in the order of
$\mathcal{G}^{(\rightarrow 2)},\mathcal{G}^{(\rightarrow
3)},\dots,\mathcal{G}^{(\rightarrow D)}$, and
$\mathcal{G}^{(D-2\leftarrow)},\mathcal{G}^{(D-3\leftarrow)},\dots,\mathcal{G}^{(1\leftarrow)}$,
respectively. The low-rank tensor $\mathcal{P}$ also can be obtained with the
cost $O(R_{1})$ since it holds $\mathcal{P}=\mathcal{G}^{(\rightarrow
D)}=\mathcal{G}^{(0\leftarrow)}$. Using these tensors, the update rule can be
written as

$\displaystyle\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}^{\mathrm{o}\backslash
d}}\frac{\mathcal{T}_{\boldsymbol{i}}}{\mathcal{P}_{\boldsymbol{i}}}\mathcal{G}^{(\rightarrow
d-1)}_{i_{1},\dots,i_{d-1},r_{d-1}}\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}}\mathcal{G}^{(d\leftarrow)}_{i_{d+1},\dots,i_{D},r_{d}}}{\sum_{\boldsymbol{i}\in\Omega_{I}^{\mathrm{o}}}\frac{\mathcal{T}_{\boldsymbol{i}}}{\mathcal{P}_{\boldsymbol{i}}}\mathcal{G}^{(\rightarrow
d)}_{i_{1},\dots,i_{d},r_{d}}\mathcal{G}^{(d\leftarrow)}_{i_{d+1},\dots,i_{D},r_{d}}}$
(14)

for $\Omega_{I}^{\mathrm{o}\backslash
d}=\Omega_{I}^{\mathrm{o}}\cap\Omega_{I}^{\backslash d}$. We used the relation
$\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}=\mathcal{T}_{\boldsymbol{i}}{\Phi}_{\boldsymbol{i}\boldsymbol{r}}$
and
${\Phi}_{\boldsymbol{i}\boldsymbol{r}}=\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}/\mathcal{P}_{\boldsymbol{i}}$
to get the above update rule. Not all elements in Equation (13) are necessary
for updating tensors by Equation (14). We provide the EM-Train factorization
in Algorithm 2 in the supplementary material.

Figure 2: The M-step for the typical tensor tree structure is decomposed into
two exactly solvable many-body approximations, which are enclosed by dotted
lines.

Reordering tensor modes for EM-Train  The tensor train decomposition results
are influenced by the order of the modes in the tensor. To quantify the
influence on and potentially enhance the performance of these tensor
decompositions, we reorder the tensor modes based on normalized mutual
information (NMI) between pairwise features in the data. We describe the
definition of NMI in Section D.2 in the supplementary material. To illustrate
this, we consider a dataset with five features. First, we select the two
modes, $j_{1}$ and $j_{2}$, with the highest NMI. These become the middle
modes in the rearranged order. Next, we choose the mode $j_{3}$ with the
second-highest NMI with $j_{1}$, placing it to the left of $j_{1}$. Then, we
select the mode $j_{4}$ with the highest NMI with $j_{2}$ among the remaining
unselected features (i.e., modes), placing it to the right of $j_{2}$. This
process is repeated until all features (modes) are selected, resulting in the
tensor modes being rearranged from ($1,2,3,4,5$) to
($j_{5},j_{3},j_{1},j_{2},j_{4}$). The effectiveness of the reordering is
examined in the supplementary material.

####  3.3 EM algorithm for more general low-rank structures

We now discuss how to find the solution for the many-body approximation
required in the M-step when a more complex low-rank structure is assumed in
the model. By the basic property of the logarithm function, the function
$L_{\mathrm{MBA}}$ to be optimized in the M-step can be decoupled into
independent solvable parts with closed form. As an example, we here see the
tensor network state described in Figure 2, which is known as a typical tensor
tree structure [18]. The objective function of many-body approximation in the
M-step can be decoupled as follows:

$\displaystyle L_{\mathrm{MBA}}(\mathcal{M};\mathcal{Q})$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\mathcal{G}_{r_{1}r_{2}r_{5}}A_{i_{1}r_{1}}B_{i_{2}r_{2}}C_{r_{5}r_{6}}D_{i_{3}r_{3}}\mathcal{H}_{r_{3}r_{6}r_{4}}E_{r_{4}i_{4}}}$
$\displaystyle=L_{\mathrm{MBA}}(\mathcal{M}^{\mathrm{Tucker}};\mathcal{Q}^{\mathrm{Tucker}})+L_{\mathrm{MBA}}(\mathcal{M}^{\mathrm{Train}};\mathcal{Q}^{\mathrm{Train}})$
(15)

where we define
$\mathcal{M}_{i_{1}i_{2}r_{1}r_{2}r_{5}}^{\mathrm{Tucker}}=\sum_{i_{3}i_{4}r_{3}r_{4}r_{6}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}$
and
$\mathcal{M}_{i_{3}i_{4}r_{3}r_{4}r_{6}}^{\mathrm{Train}}=\sum_{i_{1}i_{2}r_{1}r_{2}r_{5}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}$.
We can optimize both terms in the final line by the closed-form solution
introduced in Section 2.1. We provide theoretical support for this procedure,
including normalization conditions, in Section D.1 in the supplementary
material. By flexibly combining CP, Tucker, and Train decompositions, our EM-
based framework can approximate tensors with a wide variety of low-rank
structures.

####  3.4 Adaptive noise learning

Our general framework enables us to not only deal with typical low-rank
structures and their mixture but also allow regularization and stabilization
terms to be added to the model. For example, we here define model
$\mathcal{P}$ as the convex combination of a low-rank tensor
$\mathcal{P}^{\text{low-rank}}$ and a normalized uniform tensor
$\mathcal{P}^{\mathrm{noise}}$ whose elements all are $1/|\Omega_{I}|$ as

$\displaystyle\mathcal{P}=(1-\eta^{\mathrm{noise}})\mathcal{P}^{\text{low-
rank}}+\eta^{\mathrm{noise}}\mathcal{P}^{\mathrm{noise}},$

where $|\Omega_{I}|$ is the number of elements of the input tensor
$\mathcal{T}$, that is, $|\Omega_{I}|=I_{1}I_{2}\dots I_{D}$. The value
$\eta^{\mathrm{noise}}$ indicates the magnitude of global shift or background
noise in the data. This is a learnable parameter from the data, not a
hyperparameter. It is obvious from the general theory of the EM algorithm
described in Section 3.1 that the convexity of the E and M steps remain even
in the presence of the noise term. We show in the next section by numerical
experiments that the noise term is practically important for stabilization.

###  4 Numerical Experiments

We empirically examined the effectiveness of our framework using eight real-
world categorical datasets. We downloaded these datasets from the repositories
described in Table 5 and divided samples into 70% training, 15% validation,
and 15% test samples to form train, validation, and test tensors,
respectively. Each method decomposes the training tensor with ranks tuned such
that the cross entropy between the reconstruction and the validation tensor is
minimized. Performance was evaluated by the cross entropy between the
reconstruction and the test tensor.

Table 1: Negative log-likelihood per sample obtained on the test dataset. | MPS | BM | LPS | CPTrainON
---|---|---|---|---
SolarFlare | 6.54(0.14) | 6.38(0.09) | 6.35(0.06) |  5.98(0.04)
Monk | 6.63(0.26) | 6.81(0.01) | 6.56(0.23) |  6.39(0.05)
SPECT | 14.83(2.05) | 12.04(0.44) | 12.01(0.22) |  11.60(0.27)
Lympho. | 13.03(0.04) | 12.34(0.20) | 12.37(0.17) |  12.13(0.35)
Votes | 11.80(0.26) | 10.64(0.38) | 10.43(0.12) |  10.37(0.04)
Hayesroth | 5.20(1.01) | 5.20(1.01) | 4.71(0.01) |  4.70(0.01)
Tumor | 9.54(0.19) | 9.28(0.17) | 9.18(0.04) |  9.11(0.04)
Chess2 |  10.38(0.01) | 10.75(0.22) | 10.40(0.00) | 10.44(0.01)

We compared our approach to the following baseline methods: Matrix Product
States (MPS), Born Machine (BM), and Locally Purified State (LPS), which are
also tensor-based methods [10]. The learning rate of the baseline methods was
tuned so that the reconstruction tensor minimizes the KL divergence from the
validation data. In contrast, the proposed methods do not require a learning
rate. We ran each of the procedures ten times with random initialization and
reported mean values and the standard error in Table 1. Dataset details,
experimental setup, and computational environment are described in detail in
Section C. While the proposed framework can explore a variety of low-rank
structures and their mixtures, we examine here the performance of CPTrainON, a
mixture of CP and Train decompositions with adaptive noise term and tensor
mode re-ordering. The generalization performance of each model is shown in
Table 1. We see that CPTrainON has the best generalization performance on all
datasets except Chess2. A more detailed comparison, including EMCP, can be
found in Section B in the supplementary material.

Figure 3: Validation error for each number of parameters for each method. The
symbol ‘O’ indicates that the model includes mode reordering, and ‘N’
indicates that the model includes the noise term.

Furthermore, we evaluated the validation error of the CP, Train, and Tucker
decompositions including the adaptive noise term by varying the number of
parameters in Figure 3. For the Tucker decomposition, we only considered three
relatively small datasets, SolarFlare, Monk, and, Hayesroth, due to its high
computational cost. No method alone shows superior performance on all
datasets, however, we observe that the mixture of CP and Train generally
performs well by combining the modeling capabilities of the two decomposition
approaches. As a result, for Chess2, SolarFlare, Monk, and SPECT datasets,
CPTrainON achieved the lowest validation error, while TrainN had the lowest
validation error for Lymphography datasets, and CPN had the lowest error for
Votes. Notably, our framework readily incorporates all these procedures
providing a framework that allows decompositions with various low-rank
structures and their combinations without tuning the learning rate.

Table 2: Generalization performance | CPTrain | CPTrainN
---|---|---
SolarFlare | 6.34(0.27) | 5.98(0.06)
Monk | 6.31(0.07) | 6.43(0.07)
SPECT | 11.22(0.11) | 11.49(0.23)
Lympho. | 22.13(11.02) | 11.83(0.41)
Votes | 11.13(0.65) | 10.39(0.12)
Hayesroth | 4.74(0.04) | 4.76(0.00)
Tumor | 9.28(0.14) | 9.24(0.20)
Chess2 | 10.33(0.01) | 10.44(0.01)

We also compared the results using CPTrain with and without the adaptive noise
term to verify its usefulness in Table 2. Although the noise term does not
change the generalization performance significantly, we see in Figure 4 of the
supplementary material that the models not including the adaptive noise term
had significantly larger errors when the models were specified with large
numbers of parameters and thus prone to overfitting. The full results of the
experiment are given in Table 6 in the supplementary material.

###  5 Conclusion

While it is a well-established principle that both real- and complex-valued
tensors can be approximated with various low-rank structures based on local
singular value decomposition (SVD), a unified framework for non-negative low-
rank approximation has not been developed. Consequently, nonnegative tensor
factorizations typically require piecemeal tailored implementations that use
gradient methods. This study introduces an EM-based unified framework that
decouples the non-negative low-rank structure into multiple solvable many-body
approximations and applies a closed-form solution locally, thereby eliminating
the need for gradient methods. Our framework not only uniformly handles
various low-rank structures but also supports adaptive noise terms and
mixtures of low-rank tensors without losing the convergence guarantee. The
flexibility of our modeling stems from the generality of the EM algorithm and
a formulation that can take advantage of the logarithmic function properties
that appear in the KL divergence. Empirical results demonstrate that our
framework achieves superior generalization performance compared to baseline
methods.

Limitations and broader impacts  The proposed framework only works on non-
negative tensors. The theoretical analysis supporting the generalization
performance of the mixture models remains in future work. While this study
discussed low-rank structures with tree structures, such as CP, Tucker, and
Train decomposition, and their combinations, we did not discuss tensor
networks with loops. Our work aims to establish a fundamental methodology for
machine learning, and there is no direct risk of misuse or ethical issues.
Please also see Section D.4 in the supplementary material.

### References

  * mis, [1987] (1987).  Congressional Voting Records.  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C5C01P.
  * mis, [1989] (1989).  Solar Flare.  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C5530G.
  * Amari, [2016] Amari, S.-i. (2016).  Information geometry and its applications, volume 194.  Springer.
  * Bain and Hoff, [1994] Bain, M. and Hoff, A. (1994).  Chess (King-Rook vs. King).  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C57W2S.
  * Cheng et al., [2019] Cheng, S., Wang, L., Xiang, T., and Zhang, P. (2019).  Tree tensor networks for generative modeling.  Phys. Rev. B, 99:155131.
  * Chi and Kolda, [2012] Chi, E. C. and Kolda, T. G. (2012).  On tensors, sparsity, and nonnegative factorizations.  SIAM Journal on Matrix Analysis and Applications, 33(4):1272–1299.
  * Dempster et al., [1977] Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society, 39(1):1–22.
  * Ghalamkari and Sugiyama, [2021] Ghalamkari, K. and Sugiyama, M. (2021).  Fast Tucker rank reduction for non-negative tensors using mean-field approximation.  In Advances in Neural Information Processing Systems, volume 34, pages 443–454, Virtual Event.
  * Ghalamkari et al., [2023] Ghalamkari, K., Sugiyama, M., and Kawahara, Y. (2023).  Many-body approximation for non-negative tensors.  In Advances in Neural Information Processing Systems, volume 36, pages 257–292, New Orleans, US.
  * Glasser et al., [2019] Glasser, I., Sweke, R., Pancotti, N., Eisert, J., and Cirac, I. (2019).  Expressive power of tensor-network factorizations for probabilistic modeling.  Advances in neural information processing systems, 32.
  * Hayes-Roth and Hayes-Roth, [1989] Hayes-Roth, B. and Hayes-Roth, F. (1989).  Hayes-Roth.  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C5501T.
  * Huang and Sidiropoulos, [2017] Huang, K. and Sidiropoulos, N. D. (2017).  Kullback-Leibler principal component for tensors is not NP-hard.  In 2017 51st Asilomar Conference on Signals, Systems, and Computers, pages 693–697. IEEE.
  * Iblisdir et al., [2007] Iblisdir, S., Orus, R., and Latorre, J. (2007).  Matrix product states algorithms and continuous systems.  Physical Review B, 75(10):104305.
  * Jensen, [1906] Jensen, J. L. W. V. (1906).  Sur les fonctions convexes et les inégalités entre les valeurs moyennes.  Acta mathematica, 30(1):175–193.
  * Kim et al., [2008] Kim, Y.-D., Cichocki, A., and Choi, S. (2008).  Nonnegative Tucker decomposition with alpha-divergence.  In 2008 IEEE International Conference on Acoustics, Speech and Signal Processing, pages 1829–1832.
  * Kurgan et al., [2001] Kurgan, L. A., Cios, K. J., Tadeusiewicz, R., Ogiela, M., and Goodenday, L. S. (2001).  Knowledge discovery approach to automated cardiac spect diagnosis.  Artificial intelligence in medicine, 23(2):149–169.
  * Kırbız and Günsel, [2014] Kırbız, S. and Günsel, B. (2014).  A multiresolution non-negative tensor factorization approach for single channel sound source separation.  Signal Processing, 105:56–69.
  * Liu et al., [2018] Liu, Y., Long, Z., and Zhu, C. (2018).  Image completion using low tensor tree rank and total variation minimization.  IEEE Transactions on Multimedia, 21(2):338–350.
  * Novikov et al., [2021] Novikov, G. S., Panov, M. E., and Oseledets, I. V. (2021).  Tensor-train density estimation.  In Uncertainty in artificial intelligence, pages 1321–1331. PMLR.
  * Orús, [2014] Orús, R. (2014).  A practical introduction to tensor networks: Matrix product states and projected entangled pair states.  Annals of Physics, 349:117–158.
  * Phan and Cichocki, [2008] Phan, A. H. and Cichocki, A. (2008).  Fast and efficient algorithms for nonnegative Tucker decomposition.  In Advances in Neural Networks-ISNN 2008: 5th International Symposium on Neural Networks, ISNN 2008, Beijing, China, September 24-28, 2008, Proceedings, Part II 5, pages 772–782. Springer.
  * Silverman, [2018] Silverman, B. W. (2018).  Density estimation for statistics and data analysis.  Routledge.
  * Strehl and Ghosh, [2002] Strehl, A. and Ghosh, J. (2002).  Cluster ensembles—a knowledge reuse framework for combining multiple partitions.  Journal of machine learning research, 3(Dec):583–617.
  * Takeuchi et al., [2013] Takeuchi, K., Tomioka, R., Ishiguro, K., Kimura, A., and Sawada, H. (2013).  Non-negative multiple tensor factorization.  In 2013 IEEE 13th International Conference on Data Mining, pages 1199–1204. IEEE.
  * Thrun, [1991] Thrun, S. (1991).  The monk’s problems: A performance comparison of different learning algorithems.  Technical Report of Carnegie Mellon University.
  * Tomczak, [2021] Tomczak, J. M. (2021).  Deep Generative Modeling.  Cham: Springer International Publishing.
  * Wu, [1983] Wu, C. J. (1983).  On the convergence properties of the em algorithm.  The Annals of statistics, pages 95–103.
  * Yang et al., [2011] Yang, Z., Zhang, H., Yuan, Z., and Oja, E. (2011).  Kullback-leibler divergence for nonnegative matrix factorization.  In International Conference on Artificial Neural Networks, pages 250–257. Springer.
  * [29] Zwitter, M. and Soklic, M. (1988a).  Lymphography.  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C54598.
  * [30] Zwitter, M. and Soklic, M. (1988b).  Primary Tumor.  UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C5WK5Q.

## Supplementary Material

###  Appendix A Proofs

####  A.1 Proofs for exact solutions of many-body approximation

First, we show the known solution formulas for the best CP rank-$1$
approximation that globally minimizes the KL divergence from the given tensor.

######  Theorem 1 (Optimal M-step in CP decomposition [12]).

For a given non-negative tensor
$\mathcal{M}\in\mathbb{R}^{I_{1}\times\dots\times I_{D}\times R}_{\geq 0}$,
its many body approximation with interactions as described in Figure 1(a) is
given as

$\displaystyle
A^{(d)}_{i_{d}r}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}r}}{\mu^{1/D}\left(\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}r}\right)^{1-{1/D}}},\quad\mu=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{r\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}r}$

In the following, we provide proofs of the closed-form solutions of many-body
approximation in Figure 1(b) and (c). When a factor in a low-body tensor is
multiplied by $\nu$, the value of the objective function of many-body
approximation remains the same if another factor is multiplied by $1/\nu$,
which we call the _scaling redundancy_. The key idea in the following proofs
is reducing the scaling redundancy and absorbing the normalizing conditions of
the entire tensor into a single factor. This enables the decoupling of the
normalized condition of the entire tensor into independent conditions about
factors. This trick is also used in Section D.1 to decouple more complicated
low-rank structures into a combination of CP, Tucker, and Train
decompositions.

######  Theorem 2 (The closed form of the optimal M-step in Tucker
decomposition).

For a given tensor $\mathcal{M}\in\mathbb{R}_{\geq 0}^{I_{1}\times\dots
I_{D}\times R_{1}\times\dots\times R_{D}}$, its many-body approximation with
interaction described in Figure 1(b) is given as

$\displaystyle\mathcal{G}_{\boldsymbol{r}}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}},\quad
A^{(d)}_{i_{d}r_{d}}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\sum_{\boldsymbol{r}\in\Omega_{\boldsymbol{r}}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}.$

######  Proof 1.

The objective function of the many-body approximation is

$\displaystyle
L_{\mathrm{MBA}}(\mathcal{M};\mathcal{Q}^{\mathrm{Tucker}})=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{\mathrm{Tucker}}}$
(16)

where

$\displaystyle\mathcal{Q}^{\mathrm{Tucker}}_{i_{1}\dots i_{D}r_{1}\dots
r_{D}}=\mathcal{G}_{r_{1}\dots r_{D}}A^{(1)}_{i_{1}r_{1}}\dots
A^{(D)}_{i_{D}r_{D}}.$

Since many-body approximation parameterizes tensors as discrete probability
distributions, we optimize the above objective function with normalizing
condition
$\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{\mathrm{Tucker}}=1$.
Then, we consider the following Lagrange function:

$\displaystyle\mathcal{L}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\mathcal{G}_{\boldsymbol{r}}A^{(1)}_{i_{1}r_{1}}\dots
A^{(D)}_{i_{D}r_{D}}}-\lambda\left(\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{G}_{\boldsymbol{r}}A^{(1)}_{i_{1}r_{1}}\dots
A^{(D)}_{i_{D}r_{D}}-1\right)$

To reduce the scaling redundancy and decouple the normalizing condition, we
introduce scaled factor matrices $\tilde{A}^{(d)}$ as

$\displaystyle\tilde{A}^{(d)}_{i_{d}r_{d}}=\frac{A^{(d)}_{i_{d}r_{d}}}{a^{(d)}_{r_{d}}},\quad\mathrm{where}\
\ a^{(d)}_{r_{d}}=\sum_{i_{d}}A_{i_{d}r_{d}},$ (17)

and the scaled core tensor,

$\displaystyle\tilde{\mathcal{G}}_{\boldsymbol{r}}=\mathcal{G}_{\boldsymbol{r}}a^{(1)}_{r_{1}}\dots
a^{(D)}_{r_{D}}.$

The normalizing condition
$\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{\mathrm{Tucker}}=1$
guarantees the normalization of the core tensor $\tilde{\mathcal{G}}$ as

$\displaystyle\sum_{\boldsymbol{r}\in\Omega_{R}}\tilde{\mathcal{G}}_{\boldsymbol{r}}=1.$
(18)

The tensor $\mathcal{Q}^{\mathrm{Tucker}}$ can be represented with the above
introduced tensors as

$\displaystyle\mathcal{Q}^{\mathrm{Tucker}}_{\boldsymbol{i}\boldsymbol{r}}=\mathcal{G}_{\boldsymbol{r}}A^{(1)}_{i_{1}r_{1}}\dots
A^{(D)}_{i_{D}r_{D}}=\tilde{\mathcal{G}}_{\boldsymbol{r}}\tilde{A}^{(1)}_{i_{1}r_{1}}\dots\tilde{A}^{(D)}_{i_{D}r_{D}}.$

We optimize $\tilde{\mathcal{G}}$ and $\tilde{A}^{(d)}_{i_{d}r_{d}}$ instead
of $\mathcal{G}$ and $A^{(d)}_{i_{d}r_{d}}$. Thus the Lagrange function can be
written as

$\displaystyle\mathcal{L}$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\tilde{\mathcal{G}}_{\boldsymbol{r}}\tilde{A}^{(1)}_{i_{1}r_{1}}\dots\tilde{A}^{(D)}_{i_{D}r_{D}}}-\lambda\left(\sum_{\boldsymbol{r}}\tilde{\mathcal{G}_{\boldsymbol{r}}}-1\right)-\sum_{d=1}^{D}\sum_{r_{d}}\lambda^{(d)}_{r_{d}}\left(\sum_{i_{d}}\tilde{A}^{(d)}_{i_{d}r_{d}}-1\right).$
(19)

The condition

$\displaystyle\frac{\partial\mathcal{L}}{\partial\mathcal{G}_{\boldsymbol{r}}}=\frac{\partial\mathcal{L}}{\partial
A^{(d)}_{i_{d}r_{d}}}=0$

leads equations

$\displaystyle\tilde{\mathcal{G}}_{\boldsymbol{r}}=\frac{1}{\lambda}\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}},\quad\tilde{A}^{(d)}_{i_{d}r_{d}}=\frac{1}{\lambda_{d}}\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}.$

The values of Lagrange multipliers are identified by the normalizing
conditions (17) and (18) as

$\displaystyle\lambda=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}},\quad\lambda^{(d)}_{r_{d}}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}.$

######  Theorem 3 (The closed form of the optimal M-step in Train
decomposition).

For a given tensor $\mathcal{M}\in\mathbb{R}_{\geq 0}^{I_{1}\times\dots
I_{D}\times R_{1}\times\dots\times R_{D-1}}$, its many-body approximation with
interactions described in Figure 1(c) is given as

$\displaystyle\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d,d-1}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}}$

for $d=1,\dots,D$, assuming $r_{0}=r_{D}=1$.

######  Proof 2.

The objective function of the many-body approximation is

$\displaystyle
L_{\mathrm{MBA}}(\mathcal{M};\mathcal{Q}^{\mathrm{Train}})=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{\mathrm{Train}}}$

where

$\displaystyle\mathcal{Q}^{\mathrm{Train}}_{i_{1}\dots i_{D}r_{1}\dots
r_{D}}=\mathcal{G}^{(1)}_{i_{1}r_{1}}\mathcal{G}^{(2)}_{r_{1}i_{2}r_{2}}\dots\mathcal{G}^{(D)}_{r_{D-1}i_{D}}$

Since many-body approximation parameterizes tensors as discrete probability
distributions, we optimize the above objective function with normalizing
condition
$\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{\mathrm{Train}}=1$.
Then, we consider the following Lagrange function:

$\displaystyle\mathcal{L}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\mathcal{G}^{(1)}_{i_{1}r_{1}}\mathcal{G}^{(2)}_{r_{1}i_{2}r_{2}}\dots\mathcal{G}^{(D)}_{r_{D-1}i_{D}}}-\lambda\left(\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{G}^{(1)}_{i_{1}r_{1}}\mathcal{G}^{(2)}_{r_{1}i_{2}r_{2}}\dots\mathcal{G}^{(D)}_{r_{D-1}i_{D}}-1\right)$

To decouple the normalizing condition and make the problem easier, we
introduce scaled core tensors
$\tilde{\mathcal{G}}^{(1)},\dots,\tilde{\mathcal{G}}^{(D-1)}$ that are
normalized over $r_{d-1}$ and $i_{d}$ as

$\displaystyle\tilde{\mathcal{G}}^{(d)}_{r_{d-1}i_{d}r_{d}}=\frac{g^{(d-1)}_{r_{d-1}}}{g^{(d)}_{r_{d}}}\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}},$

where we define

$\displaystyle
g^{(d)}_{r_{d}}=\sum_{r_{d-1}}\sum_{i_{d}}\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}}g^{(d-1)}_{r_{d-1}},$

with $g_{r_{0}}=1$. We assume $r_{0}=r_{D}=1$. Using the scaled core tensors,
the tensor $\mathcal{Q}^{\mathrm{Train}}$ can be written as

$\displaystyle\mathcal{Q}^{\mathrm{Train}}_{i_{1}\dots i_{D}r_{1}\dots
r_{D}}=\mathcal{G}^{(1)}_{i_{1}r_{1}}\mathcal{G}^{(2)}_{r_{1}i_{2}r_{2}}\dots\mathcal{G}^{(D)}_{r_{D-1}i_{D}}=\tilde{\mathcal{G}}^{(1)}_{i_{1}r_{1}}\tilde{\mathcal{G}}^{(2)}_{r_{1}i_{2}r_{2}}\dots\tilde{\mathcal{G}}^{(D)}_{r_{D-1}i_{D}}$

with

$\displaystyle\tilde{\mathcal{G}}^{(D)}_{r_{D-1}i_{D}}=\frac{1}{g^{(D-1)}_{r_{D-1}}}\mathcal{G}^{(D)}_{r_{D-1}i_{D}}.$
(20)

The matrix $\tilde{\mathcal{G}}^{(D)}$ is normalized, satisfying
$\sum_{r_{D-1}}\sum_{i_{D}}\tilde{\mathcal{G}}_{r_{D-1}i_{D}}^{(D)}=1$. Thus,
the Lagrange function can be written as

$\displaystyle\mathcal{L}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\tilde{\mathcal{G}}^{(1)}_{i_{1}r_{1}}\tilde{\mathcal{G}}^{(2)}_{r_{1}i_{2}r_{2}}\dots\tilde{\mathcal{G}}^{(D)}_{r_{D-1}i_{D}}}$
$\displaystyle-\sum_{d=1}^{D-1}\lambda^{(d)}_{r_{d}}\left(\sum_{r_{d-1}}\sum_{i_{d}}\tilde{\mathcal{G}}^{(d)}_{r_{d-1}i_{d}r_{d}}-1\right)$
$\displaystyle-\lambda^{(D)}\left(\sum_{r_{D-1}}\sum_{i_{d}}\tilde{\mathcal{G}}^{(D)}_{r_{D-1}i_{D}}-1\right).$
(21)

The critical condition

$\displaystyle\frac{\partial\mathcal{L}}{\partial\mathcal{G}^{(d)}_{r_{d-1}i_{d}r_{d}}}=0$

leads the equation

$\displaystyle\tilde{\mathcal{G}}^{(d)}_{r_{d-1}i_{d}r_{d}}=\frac{1}{\lambda^{(d)}_{r_{d}}}\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d,d-1}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}},$

where the values of multipliers are identified by the normalizing conditions
in Equation (20) as

$\displaystyle\lambda^{(d)}_{r_{d}}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}^{\backslash
d}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}.$

####  A.2 Proofs for EM-algorithm

We prove the propositions used to derive the EM algorithm for non-negative
tensor mixture learning in Section 3. Furthermore, we show the convergence of
the proposed method in Theorem 4.

######  Proposition 1.

For any tensors $\Phi^{1},\dots,\Phi^{K}$ that satisfies
$\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}=1$,
the cross entropy $L(\hat{\mathcal{Q}})$ in Equation (7) can be bounded as
follows:

$\displaystyle
L(\hat{\mathcal{Q}})\geq\overline{L}(\hat{\mathcal{Q}},\Phi)=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}\log{\frac{\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}}}.$

######  Proof 3.

For any tensors $\Phi^{1},\dots,\Phi^{K}$ that satisfies
$\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}=1$,
we can transform the cross entropy $L(\hat{\mathcal{Q}})$ as follows:

$\displaystyle L(\hat{\mathcal{Q}})$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{T}_{\boldsymbol{i}}\log{\mathcal{P}_{\boldsymbol{i}}}$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{T}_{\boldsymbol{i}}\log{\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{T}_{\boldsymbol{i}}\log{\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\frac{\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}}}$
$\displaystyle\geq\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}\log{\frac{\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}}}=\overline{L}(\mathcal{Q},\Phi)$

where the following relation, called the Jensen inequality [14], is used:

$\displaystyle
f\left(\sum_{m=1}^{M}\lambda_{m}x_{m}\right)\leq\sum_{m=1}^{M}\lambda_{m}f(x_{m})$
(22)

for any convex function $f:\mathbb{R}\rightarrow\mathbb{R}$ and real numbers
$\lambda_{1},\dots,\lambda_{M}$ that satisfies $\sum_{m=1}^{M}\lambda_{m}=1$.
The inequality is adaptable because the logarithm function is a convex
function.

######  Proposition 2.

In E-step, the optimal update for $\Phi$ is given as

$\displaystyle\tilde{\Phi}^{k}_{\boldsymbol{i}\boldsymbol{r}}=\frac{\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}.$

######  Proof 4.

We put Equation (10) into the lower bound $\overline{L}$ in Equation (8),

$\displaystyle\overline{L}(\hat{\mathcal{Q}},\tilde{\Phi})$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\tilde{\Phi}^{k}_{\boldsymbol{i}\boldsymbol{r}}\log{\frac{\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\tilde{\Phi}^{k}_{\boldsymbol{i}\boldsymbol{r}}}}$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\frac{\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}\log{\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}$
$\displaystyle=\sum_{\boldsymbol{i}\in\Omega_{I}}\mathcal{T}_{\boldsymbol{i}}\log{\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\hat{\mathcal{Q}}^{k}_{\boldsymbol{i}\boldsymbol{r}}}$
$\displaystyle=L(\hat{\mathcal{Q}}).$

Jensen’s inequality in Equation (22) shows that

$\displaystyle
L(\hat{\mathcal{Q}})=\overline{L}(\hat{\mathcal{Q}},\tilde{\Phi})\geq\overline{L}(\hat{\mathcal{Q}},\Phi)$
(23)

for any tensors $\Phi\in\boldsymbol{\mathcal{D}}$ where
$\boldsymbol{\mathcal{D}}=\set{\left(\Phi^{1},\dots,\Phi^{K}\right)\mid\sum_{k}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}=1}$.
Thus, the tensors $\tilde{\Phi}=(\tilde{\Phi}^{1},\dots,\tilde{\Phi}^{K})$ are
optimal.

######  Proposition 3.

In M-step, the optimal update for $\eta$ that optimizes $J(\eta)$ in Equation
(9) is given as

$\displaystyle\eta^{k}=\frac{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}}{\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}}$

######  Proof 5.

We optimize the decoupled objective function

$\displaystyle
J(\eta)=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}\log{\eta^{k}}$

with conditions $\sum_{k=1}^{K}\eta^{k}=1$ and $\eta^{k}\geq 0$. Thus we
consider the following Lagrange function

$\displaystyle\mathcal{L}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}\log{\eta^{k}}-\lambda\left(\sum_{k=1}^{K}\eta^{k}-1\right)$

The condition $\partial\mathcal{L}/\partial\eta^{k}=0$ leads to

$\displaystyle\eta^{k}=\frac{1}{\lambda}\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}$

where the normalization identifies the multiplier $\lambda$ as

$\displaystyle\lambda=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{k=1}^{K}\sum_{\boldsymbol{r}\in\Omega_{R^{k}}}\mathcal{T}_{\boldsymbol{i}}\Phi^{k}_{\boldsymbol{i}\boldsymbol{r}}$

######  Theorem 4.

Mixture EM-tensor factorization always converges regardless of the choice of
low-rank structure and mixtures.

######  Proof 6.

We prove the convergence of the EM algorithm from the fact that the objective
function $L$ is bounded and that the E-step and M-step maximize the lower
bound with respect to $\Phi$ and $\hat{\mathcal{Q}}$, respectively. The E-step
in iteration $t$ updates $\Phi_{t-1}$ by optimal $\tilde{\Phi}_{t}$ to
maximize the lower bound for $\Phi$ such as

$\displaystyle
L(\hat{\mathcal{Q}}_{t})=\overline{L}(\hat{\mathcal{Q}}_{t},\tilde{\Phi}_{t})\geq\overline{L}(\hat{\mathcal{Q}}_{t},\Phi_{t-1}).$

The M-step in iteration $t$ updates $\hat{\mathcal{Q}}_{t}$ by
$\hat{\mathcal{Q}}_{t+1}$ to maximize lower bound for $\hat{\mathcal{Q}}$ such
as

$\displaystyle\overline{L}(\hat{\mathcal{Q}}_{t+1},\tilde{\Phi}_{t})\geq\overline{L}(\hat{\mathcal{Q}}_{t},\tilde{\Phi}_{t}).$

Again, the E-step in iteration $t+1$ updates $\Phi_{t}$ by
$\tilde{\Phi}_{t+1}$ to maximize lower bound for $\Phi$ such as

$\displaystyle
L(\hat{\mathcal{Q}}_{t+1})=\overline{L}(\hat{\mathcal{Q}}_{t+1},\tilde{\Phi}_{t+1})\geq\overline{L}(\hat{\mathcal{Q}}_{t+1},\Phi_{t}).$

Combining the above three relations, we get

$\displaystyle
L(\hat{\mathcal{Q}}_{t+1})=\overline{L}(\hat{\mathcal{Q}}_{t+1},\tilde{\Phi}_{t+1})\geq\overline{L}(\hat{\mathcal{Q}}_{t+1},\tilde{\Phi}_{t})\geq\overline{L}(\hat{\mathcal{Q}}_{t},\tilde{\Phi}_{t})=L(\hat{\mathcal{Q}}_{t}).$

Thus, it holds that $L(\hat{\mathcal{Q}}_{t+1})\geq L(\hat{\mathcal{Q}}_{t})$.
The algorithm converges because the cost function is bounded and monotonically
increasing.

###  Appendix B Additional experimental results

In this section, we discuss additional experimental results that could not be
included in the main text due to page limitations.

##### Effect of the adaptive noise term

We compare the validation error of the proposed model with and without the
adaptive noise term in Figure 4. Focusing on the vertical axis, we can see
that the validation errors for models without the adaptive noise term are huge
where the number of parameters is large which makes the learning process
unstable. Although the EM algorithms converge in theory, as shown in Theorem
4, the objective function often becomes NaN when the model has no adaptive
noise term. This is because of numerical instability due to extremely small
values in the logarithm function. The noise term, which adds values to all
elements of the tensor, eliminates this problem. While adding a small uniform
constant to values in a logarithmic function is often used as a heuristic to
stabilize learning, our method provides a principled approach to learning such
constant from the data.

Figure 4: Validation error for each number of parameters for each dataset with
adaptive noise learning (solid line) and without adaptive noise learning
(dashed line).

##### Detailed comparison of EM-based methods

We provide the comparison between Tucker and TuckerN in Table 3. We
additionally provide the full results of the experiments in Table 6 to compare
EM-based methods, including the exiting EMCP decomposition [12].
Interestingly, for some datasets, the simple EMCP shows the best
generalization performance. However, as shown in Figure 4, EMCP without the
noise term is prone to overfitting and, as a result, exhibits instability such
that the value of the validation error increases rapidly when the number of
parameters is large, i.e., the CP rank is large.

Table 3: Generalization performance | Tucker | TuckerN
---|---|---
SolarFlare | 6.78(0.15) |  6.60(0.04)
Monk | 6.80(0.00) |  6.70(0.04)
Hayesroth |  4.52(0.00) | 4.62(0.00)

##### Direct comparison between MPS and EM-Train

Both the proposed EMTrain and the baseline MPS are the same model, with
different optimization methods. Although the ranks in both models are vector
values, the official implementation of the MPS 2 assumes that each element of
the rank is the same while the rank of proposed methods has been tuned
following the way described in Section C.2. Therefore, for a direct and more
fair comparison, we have conducted our comparison experiments again fixing the
range of train ranks for EMTrainN and EMTrainON to be $(r,...,r)$ for
$r=1,\dots,8$. The results are provided in Table 4 where we again observe that
EMTrainON outperforms MPS for all datasets except Chess2.

Table 4: Generalization performance | MPS | TrainN | TrainON
---|---|---|---
SolarFlare | 6.54(0.14) |  6.02(0.04) | 6.07(0.06)
Monk | 6.63(0.26) | 6.66(0.10) |  6.43(0.07)
SPECT | 14.83(2.05) | 12.08(0.00) |  11.22(0.11)
Lymphography | 13.03(0.04) | 12.69(0.19) |  12.10(0.12)
Votes | 11.80(0.26) |  10.31(0.04) | 10.56(0.16)
Hayesroth | 5.20(1.01) | 4.94(0.06) |  4.85(0.06)
Tumor | 9.54(0.19) |  9.23(0.11) | 9.52(0.00)
Chess2 |  10.38(0.01) | 10.54(0.01) | 10.52(0.01)

###  Appendix C Experimental details

####  C.1 Experimental setup

We download six categorical tabular datasets SolarFlare, Monk, SPECT,
Lymphography, Hayesroth, and Chess2 from the UCI
database,111https://archive.ics.uci.edu/,  none of which contain any missing
values, and two preprocessed categorical tabular datasets Votes and Tumor from
the official repository of baselines [10]. All these datasets, except
Hayesroth and Chess2, are used in the paper of the baseline methods [10]. Each
dataset contains $N$ categorical samples $x_{n}=(i_{1},...,i_{D})$. Both the
number of samples, $N$, and the sample dimension, $D$, vary across datasets as
seen in Table 5. In the original datasets, each $i_{d}$ represents a
categorical quantity such as color, location, gender, etc. By mapping these to
natural numbers, each feature $i_{d}$ is converted to a natural number from
$1$ to $I_{d}$, where $I_{d}$ is the degree of freedom in the $d$-th feature.
We randomly select 70% of the $N$ samples to create the training index set
$\Omega^{\mathrm{train}}$, 15% of samples to create the validation index set
$\Omega^{\mathrm{valid}}$, and the final 15% of the samples form the test
index set $\Omega^{\mathrm{test}}$. Some datasets may contain exactly the same
samples. To deal with such datasets, we suppose that these indices sets may
contain multiple identical elements.

We create empirical tensors $\mathcal{T}^{\mathrm{train}}$,
$\mathcal{T}^{\mathrm{valid}}$, and $\mathcal{T}^{\mathrm{test}}$, where each
value $\mathcal{T}^{\ell}_{\boldsymbol{i}}$ is defined as the number of
$\boldsymbol{i}$ in the set $\Omega^{\ell}$ for
$\ell\in\set{\mathrm{train},\mathrm{valid},\mathrm{test}}$. They are typically
very sparse tensors. The above procedure to create empirical tensors is
consistent with the discussion at the beginning of Section 3. We normalize
each tensor by dividing all the elements by the sum of the tensor to map them
to a discrete probability distribution.

During the training phase, by optimizing the log-likelihood

$\displaystyle
D(\mathcal{T}^{\mathrm{train}},\mathcal{P})=\sum_{\boldsymbol{i}\in\Omega^{\mathrm{train}}}\mathcal{T}_{\boldsymbol{i}}^{\mathrm{train}}\log{\mathcal{P}_{\boldsymbol{i}}},$

we decompose the tensor $\mathcal{T}^{\mathrm{train}}$ to obtain the
reconstructed tensor $\mathcal{P}$, which approximates
$\mathcal{T}^{\mathrm{train}}$. We adjust hyper-parameters such as tensor
ranks, bounds, and learning rates to minimize the distance
$D(\mathcal{T}^{\mathrm{valid}},\mathcal{P})$. Finally, we evaluate the
generalization error $D(\mathcal{T}^{\mathrm{test}},\mathcal{P})$, where
$\mathcal{P}$ is the reconstruction approximating
$\mathcal{T}^{\mathrm{train}}$ with tuned rank. The proposed method and the
baseline method have initial value dependence. Hence, all calculations were
repeated five times to evaluate the mean and standard deviation of the cross-
entropy.

####  C.2 Implementation detail

We describe the implementation details of methods in the following. All
methods are implemented in Python 3.12.3. For reproducibility, we provide our
source code for all experiments in the supplementary material.

##### Proposed method

The pseudocodes are described in Algorithms 1 and 2. All tensors and mixture
ratios were initialized with a uniform distribution between $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 iteration was below 10e-6. The source code used
in the experiments is available in the supplementary material. We manually
determined the values of the ranks to be searched for each dataset so ensure
that we observe underfitting, better fitting, and overfitting regimes for each
validation dataset. In the EM tensor-train model, ranks of the central core
tensors are adjusted to be equal to or larger than the ranks of the core
tensors at the edges. This is because modes with large mutual information are
gathered in the center of the train-structure due to reordering. The searched
rank ranges are available in exp_config.py in the supplementary material. To
reorder tensor modes, we use the greedy method, which can be found in MI.py in
the supplementary material.

##### Baselines

We downloaded the source code for Positive Matrix Product State (MPS), Real
Born Machine (BM), and Real Locally Purified State (LPS) from the official
repository.
222https://github.com/glivan/tensor_networks_for_probabilistic_modeling,  The
license of the code is described in the repository as MIT License. We varied
the learning rates from 1.0e-4, 1.0e-3, … to 1.0e+4 for training data. We then
used the learning rate that yielded 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 20,000. We performed 5 experiments at each
bond and evaluated the mean and standard deviation. We varied the bond of the
model from $1$, $2$, … to $8$, and evaluated the test data with the bond value
that best fit the validation data.

####  C.3 Additional information for reproducibility

##### Environment

Experiments were conducted on 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 learning rate of the baselines, was less than 240 hours, using 88 threads
of parallel computing.

##### Dataset detail, license, and availability

We downloaded real-world datasets described in Table 5 through the Python
package ucimlrepo. All of these datasets are licensed under a Creative Commons
Attribution 4.0 International (CC BY 4.0) license, as seen on each web page in
the UCI database. 1 The imported data were directly converted to tensors by
the procedure described in Section C.1. For the Monk dataset, we concatenate
targets and features to make the data sparser, which is shown in prepro.py in
the supplementary material.

Table 5: Datasets used in experiments. The sparsity of dataset is defined as
$N/|\Omega_{I}|$.
| # Feature | # Non-zero values | Tensor size | Sparsity |
---|---|---|---|---|---
| $D$ | $N$ | $|\Omega_{I}|$ | $N/|\Omega_{I}|$ |
Solarflare [2] | 9 | 1067 | 41472 | 0.0257 | Link
Monks [25] | 7 | 432 | 864 | 0.5000 | Link
SPECT [16] | 23 | 267 | 4194304 | 1.7e-04 | Link
Lymphography [29] | 18 | 148 | 113246208 | 1.3e-06 | Link
Votes [1] | 17 | 376 | 86093442 | 4.4e-6 | Link
Hayesroth [11] | 5 | 160 | 256 | 0.6250 | Link
Tumor [30] | 17 | 301 | 2654208 | 1.1e-4 | Link
Chess2 [4] | 6 | 25057 | 65536 | 0.382 | Link

###  Appendix D Additonal remarks

####  D.1 Technical detail for complicated low-rank structures

We decoupled the tensor many-body decomposition into independent problems in
Section 3.3 using the basic property of logarithmic functions. We show here
that the normalization for the model is satisfied when we use the solution
formula of the many-body approximation for each decoupled problem.

In the following, we discuss the decomposition as described in Figure 2 as an
example, while the generalization to arbitrary tree low-rank structures is
straightforward. When we decouple the many-body approximation in Equation
(3.3) into the many-body approximation corresponding to the Tucker and Train
decompositions, we need to guarantee that the normalizing condition

$\displaystyle\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}=1$
(24)

is satisfied where we define

$\displaystyle\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}=\mathcal{G}_{r_{1}r_{2}r_{5}}A_{i_{1}r_{1}}B_{i_{2}r_{2}}C_{r_{5}r_{6}}D_{i_{3}r_{3}}\mathcal{H}_{r_{3}r_{6}r_{4}}E_{r_{4}i_{4}}.$
(25)

As explained at the beginning of Section A.1, we reduce scaling redundancy by
normalizing each factor and decomposing the Lagrange function into independent
parts. More specifically, we define a single root tensor and introduce
normalized factors that sums over the edges that lie below from the root.
Although the choice of the root tensor is not unique, we let tensor
$\mathcal{G}$ be the root tensor and introduce

$\displaystyle\tilde{A}_{i_{1}r_{1}}$
$\displaystyle=\frac{1}{a_{r_{1}}}A_{i_{1}r_{1}},\quad\tilde{B}_{i_{2}r_{2}}=\frac{1}{b_{r_{2}}}B_{i_{2}r_{2}},\quad\tilde{C}_{r_{5}r_{6}}=\frac{h_{r_{6}}}{c_{r_{5}}}C_{r_{5}r_{6}},$
(26) $\displaystyle\tilde{D}_{i_{3}r_{3}}$
$\displaystyle=\frac{1}{d_{r_{3}}}D_{i_{3}r_{3}},\quad\tilde{E}_{i_{4}r_{4}}=\frac{1}{e_{r_{4}}}E_{i_{4}r_{4}},\quad\tilde{\mathcal{H}}_{r_{3}r_{4}r_{6}}=\frac{d_{r_{3}}e_{r_{4}}}{h_{r_{6}}}\mathcal{H}_{r_{3}r_{4}r_{6}}$
(27)

where each normalizer is defined as

$\displaystyle a_{r_{1}}$ $\displaystyle=\sum_{i_{1}}A_{i_{1}r_{1}},\quad
b_{r_{2}}=\sum_{i_{2}}B_{i_{2}r_{2}},\quad
c_{r_{5}}=\sum_{r_{6}}C_{r_{5}r_{6}}h_{r_{6}},$ $\displaystyle d_{r_{3}}$
$\displaystyle=\sum_{i_{3}}D_{i_{3}r_{3}},\quad
e_{r_{4}}=\sum_{i_{4}}E_{i_{4}r_{4}},\quad
h_{r_{6}}=\sum_{r_{3}r_{4}}d_{r_{3}}e_{r_{4}}\mathcal{H}_{r_{3}r_{4}r_{6}},$

then it holds that

$\displaystyle\sum_{i_{1}}\tilde{A}_{i_{1}r_{1}}=\sum_{i_{2}}\tilde{B}_{i_{2}r_{2}}=\sum_{r_{6}}\tilde{C}_{r_{5}r_{6}}=\sum_{i_{3}}\tilde{D}_{i_{3}r_{3}}=\sum_{i_{4}}\tilde{E}_{i_{4}r_{4}}=\sum_{r_{3}r_{4}}\tilde{\mathcal{H}}_{r_{3}r_{4}r_{6}}=1.$

We define the tensor $\tilde{\mathcal{G}}$ as
$\tilde{\mathcal{G}}_{r_{1}r_{2}r_{5}}=a_{r_{1}}b_{r_{2}}c_{r_{5}}\mathcal{G}_{r_{1}r_{2}r_{5}}$
and putting Equations (26) and (27) into Equations (25) and (24), we obtain
the condition for the root tensor $\mathcal{G}$ as

$\displaystyle\sum_{r_{1}r_{2}r_{5}}\tilde{\mathcal{G}}_{r_{1}r_{2}r_{5}}=1.$

Then, the tensor $\mathcal{Q}$ can be written as

$\displaystyle\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}=\tilde{\mathcal{G}}_{r_{1}r_{2}r_{5}}\tilde{A}_{i_{1}r_{1}}\tilde{B}_{i_{2}r_{2}}\tilde{C}_{r_{5}r_{6}}\tilde{D}_{i_{3}r_{3}}\tilde{\mathcal{H}}_{r_{3}r_{6}r_{4}}\tilde{E}_{r_{4}i_{4}}.$
(28)

The above approach to reduce scaling redundancy is illustrated in Figure 5.
Finally, the original optimization problem with the Lagrange function

$\displaystyle\mathcal{L}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}}-\lambda\left(1-\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}\right)$

is equivalent to the problem with the Lagrange function

$\displaystyle\mathcal{L}=\mathcal{L}^{\mathrm{Tucker}}+\mathcal{L}^{\mathrm{Train}}$

where

$\displaystyle\mathcal{L}^{\mathrm{Tucker}}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\tilde{G}_{r_{1}r_{2}r_{5}}\tilde{A}_{i_{1}r_{1}}\tilde{B}_{i_{2}r_{2}}\tilde{C}_{r_{5}r_{6}}}+\lambda^{\mathcal{G}}\left(\sum_{r_{1}r_{2}r_{5}}\tilde{\mathcal{G}}_{r_{1}r_{2}r_{5}}-1\right)$
$\displaystyle+\sum_{r_{1}}\lambda^{A}_{r_{1}}\left(\sum_{i_{1}}\tilde{A}_{i_{1}r_{1}}-1\right)+\sum_{r_{2}}\lambda^{B}_{r_{2}}\left(\sum_{i_{2}}\tilde{B}_{i_{2}r_{2}}-1\right)+\sum_{r_{5}}\lambda^{C}_{r_{5}}\left(\sum_{r_{6}}\tilde{C}_{r_{5}r_{6}}-1\right),$

which is equivalent to the Lagrange function for the Tucker decomposition
given in Equation (19) and

$\displaystyle\mathcal{L}^{\mathrm{Train}}=\sum_{\boldsymbol{i}\in\Omega_{I}}\sum_{\boldsymbol{r}\in\Omega_{R}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}\log{\tilde{\mathcal{H}}_{r_{3}r_{4}r_{6}}\tilde{D}_{i_{3}r_{3}}\tilde{E}_{i_{4}r_{4}}}$
$\displaystyle\ \
+\sum_{r_{3}}\lambda^{D}_{r_{3}}\left(\sum_{i_{3}}\tilde{D}_{i_{3}r_{3}}-1\right)+\sum_{r_{4}}\lambda^{E}_{r_{4}}\left(\sum_{i_{4}}\tilde{E}_{i_{4}r_{4}}-1\right)+\sum_{r_{6}}\lambda^{\mathcal{H}}_{r_{6}}\left(\sum_{r_{3}r_{4}}\tilde{\mathcal{H}}_{r_{3}r_{4}r_{6}}-1\right),$

which is also equivalent to the Lagrange function for the Train decomposition
given in Equation (2) assuming $G^{(D)}$ is a normalized uniform tensor. For
simplicity, we define tensors

$\displaystyle\mathcal{M}_{i_{1}i_{2}r_{1}r_{2}r_{5}}^{\mathrm{Tucker}}=\sum_{i_{3}i_{4}}\sum_{r_{3}r_{4}r_{6}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}},\quad\mathcal{M}_{i_{3}i_{4}r_{3}r_{4}r_{6}}^{\mathrm{Train}}=\sum_{i_{1}i_{2}}\sum_{r_{1}r_{2}r_{5}}\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}},$

then, solve these independent many-body approximations by the closed-form
solution by Equations (4) and (5) for given tensors
$\mathcal{M}^{\mathrm{Tucker}}$ and $\mathcal{M}^{\mathrm{Train}}$,
respectively, and multiply solutions to get optimal tensor $\mathcal{Q}$ as
Equation (28), which satisfied the normalizing condition in Equation (24).

Figure 5: We normalize all tensors except for the root tensor, which is
enclosed in a bold line. We then push the normalizer of each tensor,
$a,b,c,d,e$, and $h$ on the root tensor. The root tensor absorbs scaling
redundancy. This procedure decouples the Lagrangian into two independent
problems.

####  D.2 Definition of normalized mutual information

This section describes the definition of the mutual information [23] used in
the paper. For a given normalized non-negative tensor
$\mathcal{T}\in\mathbb{R}_{\geq 0}^{I_{1}\times\dots\times I_{D}}$, we define
its normalized mutual information between two modes $d,l\in[D]$ as

$\displaystyle{\rm H}_{d,l}=\frac{{\rm I}_{d,l}}{\sqrt{{\rm H}_{d}{\rm
H}_{l}}}$

where mutual information ${\rm I}\in\mathbb{R}^{D\times D}$ and individual
entropy ${\rm H}\in\mathbb{R}^{D}$ are defined as

$\displaystyle{\rm
I}_{d,l}=\sum_{i_{d},i_{l}}\mathcal{T}^{(d,l)}_{i_{d}i_{l}}\log{\frac{\mathcal{T}^{(d,l)}_{i_{d}i_{l}}}{\mathcal{T}_{i_{d}}^{(d)}\mathcal{T}_{i_{l}}^{(l)}}},\quad{\rm
H}_{d}=\sum_{i_{d}}\mathcal{T}_{i_{d}}^{(d)}\log{\mathcal{T}^{(d)}_{i_{d}}}$

for marginalized tensors

$\displaystyle\mathcal{T}^{(d,l)}_{i_{d}i_{l}}=\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d,l}}\mathcal{T}_{\boldsymbol{i}},\quad\mathcal{T}^{(d)}_{i_{d}}=\sum_{\boldsymbol{i}\in\Omega_{I}^{\backslash
d}}\mathcal{T}_{\boldsymbol{i}}.$

input : Non-negative tensor $\mathcal{T}$ and train rank
$R^{\textrm{Train}}=(R_{1},\dots,R_{D-1})$.

Initialize $\mathcal{Q}^{\mathrm{Train}}$;

repeat

$\mathcal{M}_{\boldsymbol{i}\boldsymbol{r}}^{\mathrm{Train}}\leftarrow\mathcal{T}_{\boldsymbol{i}}\mathcal{Q}_{\boldsymbol{i}\boldsymbol{r}}^{\mathrm{Train}}/\mathcal{P}_{\boldsymbol{i}}^{\mathrm{Train}}$
for $\boldsymbol{i}\in\Omega_{I}^{\mathrm{o}}$; // E-step

for  _$d\leftarrow 1$ to $D$_ do

Obtain $\mathcal{G}^{\left(\rightarrow d\right)}$ and
$\mathcal{G}^{\left(D-d\leftarrow\right)}$ using Equation (13);

Update $\mathcal{G}^{(1)},\dots,\mathcal{G}^{(D)}$ using Equation (14);

$\mathcal{P}^{\mathrm{Train}}_{\boldsymbol{i}}\leftarrow\mathcal{G}^{(\rightarrow
D)}_{\boldsymbol{i}}$ for $\boldsymbol{i}\in\Omega_{I}^{\mathrm{o}}$;

until  _Convergence_ ;

return _$\mathcal{G}^{(1)},\dots,\mathcal{G}^{D}$_

Algorithm 2 EM Train decomposition

####  D.3 Notation, terminology, and evaluation metric

Although all symbols and technical terms are properly introduced in the main
text, we provide our notation here for readability.

The symbol $\mathbb{R}_{\geq 0}$ denotes the set of non-negative real numbers.
The set of all natural numbers less than or equal to a natural number $K$ is
denoted by $[K]$. We use the Landau symbol $O$ for computational time
complexity.

##### Tensors

We refer to a tensor whose elements all are nonnegative as a nonnegative
tensor. The number of axes of a tensor is called its modes. The number of
modes is called the order. For example, a vector is a first-order tensor, and
a matrix is a second-order tensor. Tensors are denoted by calligraphic capital
letters, as $\mathcal{T},\mathcal{P},\mathcal{Q},\mathcal{M},\dots$. A non-
negative tensor whose sum over all indices is $1$ is called a normalized
tensor. Tensors and matrices whose sums over some indices equal $1$ are marked
with a tilde as $\tilde{\mathcal{Q}}$. Although the beginning and last core
tensors of the tensor-train format are matrices rather than tensors, they are
denoted by the calligraphic form $\mathcal{G}^{(1)}$ and $\mathcal{G}^{(D)}$
for notational convenience. The element-wise product of a normalized non-
negative tensor $\mathcal{Q}^{k}$ and a weight $\eta^{k}$ is written with a
hat, e.g., $\hat{\mathcal{Q}}^{k}=\eta^{k}\mathcal{Q}^{k}$.

##### EM algorithm

The solution spaces for the E and M steps are represented by
$\boldsymbol{\mathcal{D}}$ and $\boldsymbol{\mathcal{B}}$, respectively. Since
the solution space of the M-step is equivalent to that of the many-body
approximation for tensors, we use $\boldsymbol{\mathcal{B}}$ to denote both.
The objective function of the EM algorithm is $L$, and its lower bound is
written using overlined $L$ as $\overline{L}$.

##### Indices

The visible variables $\boldsymbol{i}=(i_{1},\dots,i_{D})$ and hidden
variables $\boldsymbol{r}=(r_{1},\dots,r_{V})$ are denoted as lower
subscripts. The index set of visible variables and hidden variables are
denoted as $\Omega_{I}$ and $\Omega_{R}$, respectively. Multiple tensors are
represented using superscripts with parentheses, such as
$\mathcal{G}^{(1)},\dots,\mathcal{G}^{(D)}$. Indices for mixtures are
expressed as superscripts without parentheses. For example, the mixing ratio
is represented by $\eta=(\eta^{1},\dots,\eta^{K})$. The symbol
$\Omega^{\mathrm{o}}_{I}$ is for the set of indices of nonzero values of the
tensor $\mathcal{T}$, that is,
$\Omega^{\mathrm{o}}_{I}=\set{\boldsymbol{i}\mid\mathcal{T}_{\boldsymbol{i}}\neq
0}\subseteq\Omega_{I}$.

##### Metric

The proposed framework and baselines optimize the KL divergence between
tensors [28], which is defined as

$\displaystyle
D_{KL}(\mathcal{T},\mathcal{P})=\sum_{\boldsymbol{i}\in\Omega_{I}}\left\\{\mathcal{T}_{\boldsymbol{i}}\log{\frac{\mathcal{T}_{\boldsymbol{i}}}{\mathcal{P}_{\boldsymbol{i}}}}-\mathcal{T}_{\boldsymbol{i}}+\mathcal{P}_{\boldsymbol{i}}\right\\}$

where $\mathcal{T}$ and $\mathcal{P}$ can be non-normalized or normalized
tensors. For a given tensor $\mathcal{T}$, the optimization of the KL
divergence for the tensor $\mathcal{P}$ is equivalent to the maximization of
the cross-entropy, which is also equivalent to minimizing the negative log-
likelihood per sample.

####  D.4 Limitations

In addition to the description in Section 5, the limitations are further
discussed below. We have empirically examined the effectiveness of the
proposed method only for discrete density estimation. The method cannot be
used directly in situations where there are missing values in the data. The
number of ranks that need to be tuned is larger in mixture models than in the
non-mixture low-rank model. The method to decouple a complicated low-rank
structure into solvable CP, Tucker, and Train decompositions is not unique.
Thus, establishing an efficient decoupling method is also a subject for future
work.

Table 6: Negative log-likelihood per sample obtained on the test dataset. The symbol O means that the model includes mode reordering and N means the model includes the adaptive noise term. | CP | CPN | Train | TrainN | TrainO | TrainON | CPTrain | CPTrainN | CPTrainO | CPTrainON
---|---|---|---|---|---|---|---|---|---|---
SolarFlare | 6.04(0.19) | 5.95(0.02) | 6.07(0.05) | 6.02(0.06) | 6.13(0.02) | 6.06(0.05) | 6.34(0.27) | 5.98(0.06) | 6.54(0.73) | 5.98(0.04)
Monk | 6.29(0.10) | 6.49(0.08) | 6.40(0.35) | 6.72(0.10) | 6.57(0.32) | 6.39(0.07) | 6.31(0.07) | 6.43(0.07) | 6.36(0.11) | 6.39(0.05)
SPECT | 11.28(0.56) | 11.46(0.19) | 12.67(1.15) | 12.08(0.00) | 11.60(0.08) | 11.67(0.12) | 11.22(0.11) | 11.49(0.23) | 11.35(0.00) | 11.60(0.27)
Lympho. | 13.02(0.00) | 12.58(0.00) | 17.34(3.20) | 12.86(0.11) | 13.53(1.24) | 12.21(0.13) | 22.13(11.02) | 11.83(0.41) | inf(nan) | 12.13(0.35)
Votes | 10.83(0.43) | 10.34(0.07) | 11.98(2.01) | 10.30(0.02) | 12.86(0.90) | 10.49(0.01) | 11.13(0.65) | 10.39(0.12) | 17.26(14.42) | 10.37(0.04)
Hayesroth | 4.52(0.00) | 4.84(0.03) | 4.56(0.02) | 4.86(0.05) | 4.48(0.06) | 4.74(0.06) | 4.74(0.04) | 4.76(0.00) | 4.65(0.01) | 4.70(0.01)
Tumor | 9.11(0.02) | 9.21(0.06) | 9.39(0.00) | 9.23(0.12) | 9.38(0.00) | 9.51(0.01) | 9.28(0.14) | 9.24(0.20) | 9.13(0.08) | 9.11(0.04)
Chess2 | 10.33(0.01) | 10.53(0.01) | 10.35(0.01) | 10.54(0.01) | 10.33(0.01) | 10.52(0.01) | 10.33(0.00) | 10.44(0.01) | 10.33(0.00) | 10.44(0.01)