Many-body Approximation for Non-negative Tensors
Authors: Kazu Ghalamkari
Mahito Sugiyama
Yoshinobu Kawahara
Venue: NeurIPS 2023 poster (2023)
Keywords: Tensor decomposition
Energy based model
Tensor networks
TL;DR: We formulate rank-free tensor decomposition via energy-based modeling of interactions between tensor modes. We analyze the relationship between our model and existing low-rank approximation models using tensor networks.

Abstract: We present an alternative approach to decompose non-negative tensors, called many-body approximation. Traditional decomposition methods assume low-rankness in the representation, resulting in difficulties in global optimization and target rank selection. We avoid these problems by energy-based modeling of tensors, where a tensor and its mode correspond to a probability distribution and a random variable, respectively. Our model can be globally optimized in terms of the KL divergence minimization by taking the interaction between variables (that is, modes), into account that can be tuned more intuitively than ranks. Furthermore, we visualize interactions between modes as tensor networks and reveal a nontrivial relationship between many-body approximation and low-rank approximation. We demonstrate the effectiveness of our approach in tensor completion and approximation.

Many-body Approximation for Non-negative Tensors
Kazu Ghalamkari
RIKEN AIP
kazu. ghalamkari@rr.
Mahito Sugiyama National Institute Informatics SOKENDAI mahito@nii.com ac.
Yoshinobu Kawahara
Osaka University
RIKEN AIP
kawaharaßist. osaka-aakkaooo u.. ac.
Abstract
We present an alternative approach to decompose non-negative tensors, called many- body approximation. Traditional decomposition methods assume low-rankness in the representation, resulting in difficulties in global optimization and target rank selection. We avoid these problems by energy-based modeling of tensors, where a tensor and its mode correspond to probability distribution and a random variable, respectively. Our model can be globally optimized in terms of the KL divergence minimization by taking the interaction between variables (that is, modes), into account that can tuned more intuitively than ranks. Furthermore, we visualize interactions between modes tensor networks and reveal a nontrivial relationship between many-body approximation and low-rank approximation. We demonstrate the effectiveness of our approach in tensor completion and approximation.
Introduction
Tensors are generalizations of vectors and matrices. Data in fields such as neuroscience [8] bioinfor- maties [2] signal processing (7), and computer vision (30 are often stored in the form of tensors, and features are extracted from them. Tensor decomposition and its non- non-negan negative version (33 are popular methods to extract features by approximating tensors by the sum of products of smaller tensors. These methods usually seek to minimize the reconstruction error, the difference between an original tensor and tensor reconstructed from obtained smaller tensors.
In tensor decomposition approaches, low-rank structure is typically assumed. where given tensor essentially represented by linear combination of small number of bases. Such decomposition requires two kinds of information: structure and the number of bases used in the decomposition. The structure specifies the type of decomposition such CP decomposition [[5] and Tucker decomposi- tion 37 Tensor networks, which were initially introduced in physics, have become popular in the
machine learning community recently because they help intuitive and flexible model design, including tensor train decomposition [29]. tensor ring decomposition 41, and tensor tree decomposition [25].
Traditional tensor structures physics, such as MERA 33 and PEPS [4) are also used machine learning. The number of bases is often referred to the rank. Larger ranks increase the capability of
the model while increasing the computational cost, resulting the tradeoff problem rank tuning that requires careful treatment by the user. Since these low-rank structure-basse decompositions via reconstruction error minimization are non-convex, which causes initial value dependence [17 Chapter 3), the problem of finding an appropriate setting of the low-rank structure is highly nontrivial
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

a.
b.
(km))
H kk
Figure I: Interaction representations corresponding to (a) second- and third-order energy (b) two- and (c) three-body approximation.
in practice as it is hard to locate the cause if the decomposition does not perform well. As result, in order to find proper structure and rank, the user must often perform decomposition multiple times with various settings, which is time- and memory-consuningg Although there is an attempt avoid rank tuning by approximating it by trace norms (2) it also requires another hyperparameter: the weights of the unfolding tensor. Hence, such an approach does not fundamentally solve the problem.
Instead the low-rank structure that has been the focus of attention in the past, in the present paper we propose a novel formulation of non-negative tensor decomposition, called many-body approximation, that focuses on the relationship among modes tensors. The structure of decomposition can be naturally determined based on the existence of the interactions between modes. The proposed method requires only the decomposition structure and does not require the rank value, which traditional decomposition methods also require as a hyperparameter and often suffer to determine.
To describe interactions between modes, we follow the standard strategy in statistical mechanics that uses an energy function H(-)) to treat interactions and considers the corresponding distribution exp (-H(-)). This model is known to be an mnrgy-bbsedd model in machine learning [99] and is exploited in tensor decomposition a tensor as a discrete probability distribution and reduces the number of parameters by enforcing some of them to be zero in optimization. We explore this energy-base approach further and discover the family of parameter sets that represent interactions between modes in the energy function K(-)) How to choose non-zero parameters in Legendre decomposition has been open problem; by addressing this problem and propose many-body approximation as a special case of Legendre decomposition. Moreover, although Legendre decomposition is not factorization of tensors in general, our proposal always offers factorization, which can reveal patterns in tensors. Since the advantage of Legendre decomposition is inherited our proposal, many-body approximation can be achieved by convex optimization that globally minimizes the Kullback-Leibler (KL) divergence (18).
Legendre decomposition [36]. Technically, it parameterizes
start
Furthermore, we introduce way of representing mode interactions, that visualizes the presence or absence of interactions between modes as a diagram. We discuss the relation to the tensor network and point out that operation called coarse coarse-grained transformation 20 in which multiple tensors are viewed as a new tensor reveals unexpected relationship between the proposed method and existing methods such as tensor ring and tensor tree decomposition.
We summarize our contribution as follows:
By focusing on the interaction between modes of tensors, we introduce rank-free tensor decompo- sition, called many-body approximation. This decomposition is achieved by convex optimization. We present way of describing tensor many- bany-body approximation called interaction representation, which is diagram that shows interactions within a tensor and can be transformed into tensor networks.
Many-body approximation can perform tensor completion via the em-algorithm, which empirically shows better performance than low-rank based existing methods.
2 Many-body Approximation for Tensors
Our proposal, many-body approximation, is based on the formulation of Legendre decomposition for tensors, which we first review in Section 2.1 Then we introduce interactions between modes and its visual representation in Section 2.2 many-body approximation in Section 2.3 and transformation of interaction representation into tensor networks in Section 2.4 In the following discussion, we consider
D-order non-negative tensors with the size (J,,..., ., JD), and we denote by [K] = (1,2,.... K} for positive integer K. We assume the sum of all elements in given tensor to be for simplicity,

---

while this assumption can be eliminated using the general property of the Bullback-Leibbe (KL) divergence, ID Q) (xI(AP, XQ), for any R20-
2.1 Reminder Legendre Decomposition and its optimization
Legendre decomposition is method for decomposing non-negative tensor by regarding the tensor as a discrete distribution and representing with a limited number of parameters. We describe
non-negative tensor P using parameters 0 (02............. (g.... and its energy function H. as
Pip = exp (- -H Hig...dD = Σ
where 0 has role of normalization. sample space of distribution and the value of each entry of the tensor is regarded as the probability of realizing the corresponding index (35]
clear that the index set of tensor corresponds to
As we see in Equation ( we can uniquely identify tensors from the parameters 0. Inversely, we can compute 0 from given tensor as
4
In
A...ko log Pi Pip...bb
(2)
=1=1
==
using the Möbius function μ USS 1 1-1, 1, +1}, where S is the set of indices, defined inductively as follows:
if tha Ad [D], else if ia S d'd. Vi E [D),
E IP Σ Dizia Pi...
otherwise.
Since distribution described by Equation I belongs to the exponential family, O corresponds to natural parameters of the exponential family, and we can also identify each tensor by expectation parameters (17........... M......n) using Möbius inversion as
Ti Tp
P Pip....i
===1
where 7 1 because of normalization. See Section A in Supplement for examples of the above calculation. Since distribution is determined by specitying either 0 e-parameters or 77parameters, they form two coordinate systems, called the d-coordinate system and the n-coordinate system, respectively.
2.1.1 Optimization
Legendre decomposition approxinates tensor by setting some 9 values be zero, which corresponds to dropping some parameters for regularization. It achieves convex optimization using the dual flatness of (- and 77-ooodinaaee systems. Let be the set of indices of parameters that are not imposed be 0. Then Legendre decomposition coincides with a projection of given nonnegative tensor P onto the subspace B {Ω O.,.....p = O (i1..... iD) & B for d-parameters of 2).
Let us consider projection of a given tensor P onto B. The space of probability distributions, which not Euclidean space, is studied information geometry. By taking geometry of probability distributions into account, we can introduce the concept of flatness for a set of tensors. The subspace B e-flat if the logarithmic combination, or called e-geodesic, R E {(1 - () log Q + log Q2 s(())
0 < t < 1} of any two points Q1, Q2 E B is included in the subspace B, where o(t) is a normalizer. There is always a unique point P on the e-flat subspace that minimizes the KL divergence from any point P.
= arg min min (PL(P,))
(4)
2.0EB
3

---

This projection is called the m-projection. The m-projection onto an e-flat subspace is convex optimization. We define two vectors 0 = (f))BEB and n78 = (17. Be and write the dimensionality of these vectors as []] since is equal to the cardinality of B. The derivative of the KL divergence and the Hessian matrix G IR [||[ are given as
00 DKL (K Q) = 878 -)" Gu = 7max (11 d. max((oo.)) ".i.......p Jor
(5)
where o7 and r) B are the expectation parameters of Q and P. respectively, and u (ii,.... iD). ): (J1....,jnn jD) E B. This matrix G is also known as the negative Fisher information matrix. Using gradient descent with second-order derivative, we can update 0 in each iteration t as
- = 0 - - (n 1)
(6)
The distribution Q++1 is calculated from the updated natural parameters 0.+11 This step finds point Q++1 E B that is closer the destination P along with the e-geodesic from Q P. We can also calculate the expected value parameters 7/+11 from the distribution. By repeating this process until convergence,
always find the globally optimal solution satisfying Equation (4). See Section A in Supplement for more detail on the optimization. The problem of how to design and interpret B has not been explored before, and firstly give the reasoning of B the presence or absence of particular interactions between modes, which naturally leads to the design guideline of B. The computational complexity of Legendre decomposition is (()) where y is the number of iterations.
2.2 Interaction and its representation of tensors
In this subsection, introduce interactions between modes and its visual representation to prepare for many-body approximation. The following discussion enables us to intuitively describe relationships between modes and formulate our novel rank-free tensor decomposition.
First, we introduce n-body parameters, a generalized concept of one-body and two-body parameters in M. In any n-body parameter, there are non-one indices; for example, 1,,2,1,1 is one- body parameter, 11,3,1,1 is a two-body parameter, and 01,2,1,3 is a three-body parameter. We also use the following notation for n-body parameters:
(k.m) Ae
(k.mp.)) Sesmm
1 = 01....1.........11
f1.....14
for n = 1,2 and 3, respectively. We write the energy function H with n-body parameters as
D
D
m-1
D m- k-
Hip io = Ho + Σ H(m)
Σ
(k.m) fe
+
Σ 2
H
(o.k.m) fp
m=
m
k=1
m=1
k=1 pm1
where the n-th order energy is introduced as
(8)
8 1/
For simplicity, we suppose that 15 144- .... that an n-body interaction exists between modes l,,..., L if there are indices it..... it, satisfying 1(((....d)) it.. # The first term, Ho in Equation { is called the normalized factor or the partition function. The terms (k) are called bias in machine learning and magnetic field self-energy in statistical physics. The terms H (k.m) are called the weight of the Boltzmann machine in machine learning and two-body interaction or dlectrnn-electrog interaction in physics.
D holds. We set Ho - We say
To visualize the existence of interactions within tensor, we introduce a diagram called interaction representation, which inspired by factor graphs graphical modeling Chapter 8). The graphical representation of the product of tensors is widely known tensor networks. However, displaying the relations between modes of a tensor as factor graph is our novel approach. We represent the n- bo-o interaction black square,
connected with modes. We deseribe examples of the

---

a.
X
b.
C.
8
x
14
X
X
P
44 (4
82
$g
Figure 2: (a)(eeft) Interaction representation of an example of cyclic two-body approximation and (a)(right) its transformed tensor network for D 4. Each tensor is enclosed by square and each mode is enclosed by a circle. black circle between modes mean interaction existence. (b) Tensor network of tensor ring decomposition. (c) The relationship among solution spaces of tensor-ring decomposition Bring. two-body approximation B2, and cyclic two-body approximation B Beye-
a hyper diagonal tensor. Edges through
two-body interaction between modes (k: m) and the three-body interaction among modes (k; m. p) in Figure a). Combining these interactions, the diagram for the energy function including all two-body interactions is shown in Figure bb), and all two-body and three-body interactions are included in Figure [e) for D 4. This visualization allows us intuitively understand the relationship between modes of tensors. For simplicity, we abbreviate one-body interactions in the diagrams, while we always assume them. Once interaction representation is given, we can determine the corresponding decomposition of tensors.
In Boltzmann machines, usually consider binary (two-level) variables and their second-orde interactions. In our proposal, we consider multi-level D variables, each of which can take a natural number from to La for d E [D]. Moreover, higher-order interactions among them are allowed. Therefore, our proposal is form of a multi-level extension of Boltzmann machines with higher-order interaction, where each node Boltzmann machines corresponds the tensor mode. the following section, we reduce some of n-bodd interactions, that is, ((lnn.dda)) = 0 - by fixing each parameter
(a) 0 for all indices (i, (4,,).., it,) {2..... 4. x .... {2..... ........ I.n
2.3 Many-body approximation
We approximate given tensor by assuming the existence of dominant interactions between the modes of the tensor and ignoring other interactions. Since this operation can be understood as enforcing some natural parameters of the distribution be zero, it can be achieved by convex optimization through the theory of Legendre decomposition. We discuss how we can choose dominant interactions in Section 3.1
As an example, we consider two types of approximations of nonnegative tensor P by tensors represented in Figure [b, e). If all energies greater than second-orde or those greater than third- order in Equation ( are ignored that is, H (l(l..dd)) 0 for n > or > 3 the tensor Piy.iai8ah4 approximated as follows:
Phiiziiais 222
x (1.2) (1.3) (1.4) (2.3) (2.4) (3,4)
i hise 11.24 ig.23
12. is.4
Phrizidsis 333
(1.2.3)
(12244) (1,3,4) (2.3,4)
Xi
Xiy.a.d
Xizi
where each of matrices and tensors on the right-hand side is represented as
(kum)
1
exp
((())
(k)
+
(k.))
+
1
H
(m))
77
3
3
(k.mp))
V
1
exp
H
(k)
H(mm (m)
(())
3
H. (k)m)
sm km
H (m.p)
+ Hr (R.p) + H (k.m.))
1
H.
The partition function, or the normalization factor, given as Z exp (-00,,,)) which does not depend on indices (6, f2, ig, i4). Each (k.m) (resp. X (k.m..) is factorized representation for the relationship between k-th and m-th (resp. k-th, m-th and p-th) modes. Although our model can be
5

---

transformed into linear model by taking the logarithm, our convex formulation enables us to find the optimal solution more stable than traditional linear low-rank-base based nonconvex approaches. Since we do not impose any low-rankness, factorized representations, such as X (k,r) and x (k.mp)) can be full-rank matrices or tensors.
For a given tensor P, we define its m-body approximation pSm as the optimal solution of Equation 4 that is, pmm = arg min D. DKL(P, 2) where the solution space B, is given as
B = { Q ii.....i O if Bs,....ip is n(> m--body parameters for 0 p-parameters of )}
Note that PSD P always holds. For any natural number m < D, it holds that Bm C mm+1- Interestingly, the two-body approximation for a non-negative tensor with I ... ID = 2 is equivalent to approximating the empirical distribution with the fully connected Boltzmann machine.
In the above discussion, we consider many-body approximation with all the n-body parameters, while we can relax this constraint and allow the use of only a part of n-body parameters in approximation. Let us consider the situation where only one-body interaction and two-body interaction between modes (d, + 1) exist for all d E [D](D + 1 implies for simplicity). Figure 2(a) shows the interaction representation of the approximated tensor. As we can confirm by substituting 0 for H (R.1) if 1 +
1, we can describe the approximated tensor as
Pii = pere Phi
x! x 1112 121d
X
where
(1) f4.4441
(1 H!(h io
1 + (((b1))
=
Vz
exp
2
+
2
with the normalization factor Z = exp (-01,..,p) When the tensor P is approximated by peye, the set B contains only all one-body parameters and two-body parameters ig-daja o(d,411) for d [D]. We call this approximation cyclic two-body approximation since the order of indices in Equation ( is cyclic.
It holds that B eye C B2 for the solution space of cyclie two-body approximation B Beyc (Figure ()).
2.4 Connection to tensor networks
Our tensor interaction representation is an undirected graph that focuses on the relationship between modes. By contrast, tensor networks, which well known as diagrams that focus on smaller tensors after dccomposition, represent a tensor as an undirected graph, whose nodes correspond to matrices or tensors and edges to summation over a mode in tensor products (6).
We provide examples in our representation that can be converted to tensor networks, which implies tight connection between our representation and tensor networks. For the conversion, we use hyper-diagonal tensor Z defined as Liyk ó1,ō,1ōm , where d, = if i = j and 0 otherwise. The tensor Z is often represented by in tensor networks. In the community of tensor networks, the tensor appears in the CNOT gate and special case of the Z spider 28 The tensor network Figure 2[a) represents Equation for D 4.
remarkable finding that the converted tensor network representation cyclic two-body approx- imation and the tensor network of tensor ring decomposition, whose tensor network is shown in Figure bb, have similar structure, despite their different modeling. If we consider the region enclosed by the dotted line in Figure a) as new tensor, the tensor network of the cyclic two-body approximation coincides with the tensor network of the tensor ring decomposition shown in Fig- ure bb). This operation, in which multiple tensors are regarded as a new tensor in tensor network, is called coarse-graining transformation [9]
Formally, cyclic two-body approximation coincides with tensor ring decomposition with a specific constraint as described below. Non- Non-negative tensor ring decomposition approximates a given tensor
P E R p with D core tensors (),..,, x..,) E RR Ra-1 Ra
20
20
as
R
Rp
Phi...i
Σ n==1
Xppa Xrp idpap
(1)
(D)
(10)
r==1

---

where (Ri...., RD) is called the tensor ring rank. The cyclic two-body approximation also approxi- mates the tensor P in the form of Equation (T)p, imposing an additional constraint that each core tensor x (d) decomposed
Ia Σ m==1
X ra
Imasaara
(11)
1.m
for each E [D]. We assume ro = FD for simplicity. We obtain Equation Equation (11) into Equation (1O).
by substituting
This constraint enables us to perform convex optimization. This means that we find a subclass eye that can be solved by convex optimization in tensor ring decomposition, which has suffered from the difficulty of non-convex optimization. In addition, this is simultaneously a subclass of two-body approximation, shown in Figure[[(e..
From Kronecker's delta 8, ra = ig holds in Equation 11 thus, ( is a tensor with the size (La-1, La: La). Tensor ring rank after the cyclic two-body approximation (LI...., ID) since the size of core tensors coincides with tensor ring rank. This result firstly reveals the relationship between Legendre decomposition and low-rank approximation via tensor networks.
Here we compare the number of parameters of cyclic two-body approximation and that of tensor ring decomposition. The number of elements of an input tensor is 11 .. Ip. After cyclic two-body approximation, the number [B] of parameters is given as
D = ++ Σ - (la- d=1
(12)
d=1
where we assume JD+1 = I1 The first term is for the normalizer, the second for the number of one-body parameters, and the final term for the number of two-body parameters. In contrast, in the tensor ring decomposition with the target rank (Rı,-. RD), the number of parameters is given as E RalaBa+1. The ratio of the number of parameters of these two methods is proportional to I/R² if we assume R₂ = R and Ia = I for all d E [D] for simplicity. Therefore, when the target rank is small and the size of the input tensor large, the proposed method has more parameters than the tensor ring decomposition. The correspondence of many-body approximation to existing low-rank approximation not limited tensor ring decomposition. We also provide another example of the correspondence for tensor tree decomposition in Section C.2 in Supplement.
2.5 Many-body approximation as generalization of mean-field approximation
Any tensor P can be represented by vectors (1) x (D)
(2) R la as P.,.....i 0. x/ J.
(D)
if and only if all n² 2)-body f-parameters are 0 [10 12 The right-hand side is equal to the Kronecker product of D vectors x ,.... (P); therefore, this approximation is equivalent to the rank-1 1 approximation since the rank of the tensor that can be represented by the Kronecker product is always 1, which is also known correspond to mean-field approximation. In this study, propose many-body approximation by relaxing the condition for the mean-field approximation that ignores n² 2)-body interactions. Therefore many-body approximation is generalization of rank-1 approximation and mean-field approximation.
The discussion above argues that the one-body approximation - that is, the rank-1 approximation that minimizes the KL divergence from a given tensor - a convex problem. Although finding the rank-1 tensor that minimizes the Frobenius norm from given tensor is known to be Na-hard problem [14], it has been reported that finding the rank-1 1 tensor that minimizes the KL divergence from given tensor is convex problem [16]. Therefore, our claim is not inconsistent with those existing studies. It should be also noted that, except for the rank-1 case, there is no guarantee that traditional low-rank decomposition that optimizes the KL divergence is convex problem 5
2.6 Low-body tensor completion
Since our proposal is an alternative to the existing low-rank approximation, we can replace it with our many-body approximation in practical applications. As representative application, we demonstrate
7

---

Algorithm 1: Low-body tensor completion
LBTC(T. Q1)
I/ Q is observed indices
t+1
Initialize pt
I| c.g. random tensor
Pá Tn
repeat
pe+1 MAY (P', B)
|| m-step, See Algorithm 2
Ph1 -- Tn
|| e-step
rest f |pt+1 - P'11 |/|1//11'/|
1 ( 1-1 1+1
until res" rest-1 and 1
// We set 10-5 in implementation;
return P
Table 1: Recovery fit score for tensor completion.
Scenario 1
Scenario 2
Scenario 3
LBTC
0.917 0.00120.00152
0.874:0.00153
HaLRTC
0.864
0.842
0.825
SiLRTC
0.863
0.841
0.820
SiLRTCTT
0.948
0.915
0.868
PTRCRW
0.0400000000
0.901
0.844=0.0000
missing value estimation by the em-algorithm [26]. where we replace the m-step of low-rank approximation with many-body approximation. We call this approach low-body tensor completion (LBTC). The em-algorithh with tensor low-rank approximation does not guarantee the uniqueness of the m-step because the model space is not flat, while LBTC ensures that the model space is flat, and thus the m-step of LBTC is unique. See more discussion information geometry for LBTC in Section C.4 We provide its pseudo-code in Algorithm [ Instead of hyper-parameten related to ranks that traditional tensor completion methods require, LBTC requires only information about interactions to be used, but it can be intuitively tuned, also seen in Section 3.1 We examine the performance of LBTC in Section 3.2
3 Experiments
We conducted three experiments see the usefulness and effectiveness of many-body approximation. Datasets and implementation details are available in Supplement.
3.1 An example of tuning interaction
We extracted 10 color images from Columbia Object Image Library (COIL-100) 27 and constructed tensor P I2-00 Tx40x00x3xxx where each mode represents the width, height, color, image index, respectively. Then decomposed P with varying chosen interactions and showed the reconstructed images, Figure 3 a) shows ps4, which coincides with the input tensor P itself since its order is four. Reconstruction quality gradually dropped when we reduced interactions from Figure 3 to Figure3jd) We can explicitly control the trend of reconstruction by manually choosing interaction, which one of characteristics of the proposal. In Figure 3 B(e,, each reconstructed image has monotone taste because the color is not directly interacted with the pixel position (width and height). In Figure although the whole tensor is monotone because the color independent of any other modes, the reconstruction keeps the objects' shapes since the interaction among width, height, and image index are retained. See Figure 7 in Supplement for the role of the interaction. Interactions between modes are interpretable, in contrast to the rank of tensors, which is often difficult to interpret. Therefore, tuning of interactions is more intuitive than that of ranks.
We also visualize the interaction between color and image a heatmap of matrix X (c.) E Bax 10 in Figure 3 el),
e), which is obtained in the case of Figure (c.. We normalize each column of X (c.)) with the

---

10
5
82
0.79
Figure 3: (a.b) Four-body and three-body approximation for the COIL dataset. (c.d) Appropriate selection of interactions controls the richness of color. Visualization of interaction between (e) color and image index, and (f) width, height, and image index.
L2 norm for better visualization. As expected, indices corresponding to red objects (that is, 1, 4, and have large values the third red-color row. We also visualize the interaction among width, height, and image index, which corresponds
tensor ( E 40xx40x 10 without any information of colors, as a series of heatmaps in Figure 3(). They provide the gray-scale reconstruction of images. The above discussion implies that many-body approximation captures patterns features of inputs.
3.2 Tensor completion by mony-body approximation
To examine the performance of LBTC, we conducted experiments on a real-world datasct: the traffic speed dataset in District 7, Los Angeles County, collected from PeMS (3]] This data was obtained from fixed-point observations of vehicle speeds every five minutes for 28 days. We made 28 24 12 x tensor T, whose modes correspond date, hour, minute, and lane. As practical situations, we assumeed that a disorder of a lane's speedometer causes random and continuous deficiencies. We prepared three patterss of missing cases:
with disorder in the speedometer in Lane (Scenario 1), a case with disorders in the speedometers in Lanes 1, 2, and 3 (Scenario 2), and case with disorders in all speedometers (Scenario 3). The missing rates for scenarios 1, and 3 were 9 percent, 27 percent, and 34 percent, respectively. We evaluate the recovery fit score 1 1- T, - T- || for the reconstruction and missing indices т.
In the proposal, we need to choose interactions to use. Intuitively, the interaction between date and minute and that between minute and lane seem irrelevant. Therefore, we used all the two-body interactions except for the above two interactions. In addition, the interaction among date, hour, and minute also included because it intuitively relevant.
As baseline methods, we use RIT, HaLRTC (1], SiLRTC-TT 2, and PTRCRW (38 SiLRTC- TT and PTRCRW are based on tensor train and ring decompositions, respectively, and PTRCRW is known as the state of the art. The reconstructions are available in Figure iin Supplement with ground truth. We tuned hyperparameters of baseline methods and compared their best results with LBTC. The resulting scores (see Table show that our method is superior to other methods. Moreover, we can intuitively tune interactions the proposal, while other baselines need typical hyperparameter tuning as seen in Figure in Supplement. In our method, interactions among ansossiccatee modes could worsen the performance as seen in Table 2 in Supplement. Since the proposal and PTRCRW are non-convex optimizations, we ran each of them 10 times with random initialization and reported mean values and the standard error. Many-body approximation itself is always convex, while our tensor completion algorithm LBTC is non-convex due to the use of the em-algorithm.

---

3.3 Comparison with ring decomposition
As seen in Section 2.4 many-body approximation has close connection low-rank approximation. For example, in tensor ring decomposition, if we impose the constraint that decomposed smaller tensors can be represented as products with hyper-diagonal tensors I, this decomposition equivalent to a eyclic two-body approximation. Therefore, to examine our conjecture that cyclic two-body approximation as capable of approximating as tensor ring decomposition, we empirically examined the effectiveness of cyclic two-body approximation compared with tensor ring decomposition. As baselines, we used five existing methods of non-negative tensor ring decomposition: NTR-APG, We also describe the efficiency in
NTR-HALS, NTR-MU, NTR-MM, and NTR-IraMM [39 40 Supplement.
We evaluated the approximation performance by the relative error I| -| for an input tensor ' and reconstructed tensor T. Since all the existing methods based nonconvex optimization, we plotted the best score (minimum relative error) among five restarts with random initialization. In contrast, the score of our method is obtained by single run as it is convex optimization and such restarts are fundamentally unnecessary.
Synthetic data We performed experiments two syn- thetic datasets with low tensor ring rank. We create sixth- and seventhth-order tensors with ring ranks (15, 15)
30,,30200 30,300 RIS
120,00,20 20,00,0
0.020
Reletive Ereo
0.0100
Reistvv Etrar
0.015
HALS
and (10, ..., 10), respectively. In Figure 4a), relative er- rors are plotted with gradually increasing the target rank of the tensor ring decomposition, which is compared to the score of method, plotted as the cross point of horizontal and vertical red dotted lines. Please note that method does not have the concept of the rank, the score of our method is invariant to changes of the target rank unlike other methods. If the cross point of red dot- ted lines is lower than other lines, the proposed method
0.0075
8.0050
0.010
0.0025
0.005
10
10'
Perameters
Perameters
ChartBass
Tr.rggami
0.3200
Refetive Enror
0.3150
Reletive
better than other methods. In both experiments, the proposed method is superior to comparison partners in terms of the reconstruction error.
0.3100
1200 1800
500 1500 2000
Parameters
* Paremeters
Real data Next, we evaluated our method on real data. TT_ChartRe and TT_ FT__rrigat 0riiamiiiii are seventh-order ten- sors that are produced from TokyoTech Hyperspectral Image Dataset 22 224. Each tensor has been reshaped to reduce the computational complexity. As seen in Fig-
Figure 4: Experimental results for (a) synthetic data and (b) real datasets. The vertical red dotted line is [B] (see Equa- tion (12).
b)), the proposed method keeps the competitive relative errors. In baselines, slight change of the target rank can induce a significant increase of the reconstruction error due to the nonconvexity. These results mean that we climinate the instability of non-negative tensor ring decomposition by our convex formulation.
4 Conclusion
We propose many-body approximation for tensors, which decomposes tensors by focusing on the relationship between modes represented by an energy-based model. This method approximates and regularizes tensors by ignoring the energy corresponding to certain interactions, which can be viewed as generalized formulation of mean-field approximation that considers only one- interactions. Our novel formulation enables us to achieve convex optimization of the model, while the existing approaches based on the low-rank structure are non-convex and empirically unstable with respect to the rank selection. Furthermore, we introduce a way of visualizing activated interactions between modes, called interaction representation. We have demonstrated transformation between our representation and tensor networks, which reveals the nontrivial connection between many-body approximation and the classical tensor low-rank tensor decomposition. We have empirically showed the intuitive model's 's designability and the better usefulness in tensor completion and approximation tasks compared baseline methods.
Limitation Proposal only works on non-negative tensors. We have examined the effectiveness of LBTC on only traffic datasets.
10

---

Acknowledgments and Disclosure of Funding
This work was supported by RIKEN, Special Postdoctoral Researcher Program (KG), JST, CREST Grant Number JPMJCR22DS, Japan, and ISPS KAKENHI Grant Number JP21H03503 (MS), and IST, CREST Grant Number JPMJCR1913, Japan (YK).
References
[1] Amari, S. (2016). Information Geometry and Applications. Springer.
[2] Bengua, J. A., Phien, H. N., Tuan, H. D., and Do, M. N. (2017). Efficient tensor completion for color image and video recovery: Low-rank tensor train. IEEE Transactions Image Processing, 26(5):2466-2477 2479..
[3] Bishop, C. M. and Nasrabadi, N. M. (2006). Pattern recognition and machine learning, volume 4. Springer.
[4] Cheng, S., Wang, L.,, and Zhang, (2021). Supervised learning with projected entangled pair states. Physical Review B, 103(12):125117.
[5] 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
[6] Cichocki, A., Lee, N., Oseledets, I., Phan, A.-H., Zhao, Q. Mandic, D. P., et al. (2016). Tensor networks for dimensionality reduction and large-scale optimization: Part low-rank tensor decompositions. Foundations and Trends@ in Machine Learning, 9445))249-429,
[7] Cichocki, A., Mandic, D., De Lathauwer, L., Zhou, G., Zhao, Q., Caiafa, C., and Phan, H. A. (2015). Tensor decompositions for signal processing applications: From two-way to multiway component analysis. IEEE signal processing magazine, 3252)11
[8] Erol, A. and Hunyadi, B. (2022). Tensors for neuroimaging: review on applications of tensors unravel the mysteries of the brain. Tensors for Data Processing, pages 427-482.
[9] Evenbly, G. and Vidal, G. (2015). Tensor network renormalization. Physical review letters, 1 118(15)180445.
[10] Ghalamkari, K. and Sugiyama, M. (2021). Fast tucker rank reduction for non non-negative ten- sors using mean-field field approximation. In Advances Neural Information Processing Systems, volume 34, pages 443445, Virtual Event.
[11] Ghalamkari, K. and Sugiyama, M. (2022). Fast rank-1 NMF for missing data with KL diver- gence. In Proceedings of the 25th International Conference on Artificial Intelligence and Statistics, pages 2927-2940, Virtual Event.
[12] Ghalamkari, K. and Sugiyama, M. (2023). Non-negative low-rank approximations for multi- dimensional arrays on statistical manifold. Information Geometry.
[13] Hansen, S., Plantenga, T., and Kolda, T. G. (2015). Newton-based optimization for kullback- leibler nonnegative tensor factorizations. Optimization Methods and Software, 1005):1002-1029.
[14] Hillar, C. J. and Lim, L.-H. (2013). Most tensor problems are np-hard. Journal of the ACM (JACM), 60(6):1-39,
[15] Hitchcock, F. L. (1927). The expression a tensor or a polyadic as a sum of products. Journal of Mathematics and Physics, 6 6(1-4):164-189.
[16] 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.
[17] Kolda, T. G. and Bader, B. (2009). Tensor decompositions and applications. SIAM review, 5153):555-500.

---

[18] Kullback, S. and Leibler, R. A. (1951). On information and sufficiency, The annals of mathematical statistics, 22(1):79-86.
[19] LeCun, Y., Chopra, S., Hadsell, R., Ranzato, M., and Huang, F. (2006). tutorial on energy- based learning. Predicting struciured data, 1(0).
[20] Levin, M. and Nave, C. P. (2007). Tensor renormalization group approach to two-dimensional classical lattice models. Physical Review Letters, 19(12):120001
[21] Liu. J., Musialski, P., Wonka, P., and Ye, J. (2013). Tensor completion for estimating missing values in visual data. IEEE Transactions Pattern Analysis and Machine Intelligence, 35(1):208- 220.
[22] Luo, Y., Wang. F., and Szolovits, P. (2017). Tensor factorization toward precision medicine. Briefings in Bioinformatics, 18(3):511-514.
[23] Monno, Y. Kiku, D., Tanaka, M., and Okutomi, M. (2017). Adaptive residual interpolation for color and multispectral image demosaicking. Sensors, 17(1))27787.
[24] Monno, Y., Kikuchi, S., Tanaka, M., and Okutomi, M. (2015). A practical one-shot multispectral imaging system using single image sensor. IEEE Transactions on Image Processing. 24(10):3048- 3059.
[25] Murg, V., Verstracte, F., Legeza, Ö., and Noack, R. M. (2010). Simulating strongly correlated quantum systems with tree tensor networks. Physical Review B, 82(20):205105.
[26] Narita, A., Hayashi, K., Tomioka. R., and Kashima, H. (2012). Tensor factorization using auxiliary information. Data Mining and Knowledge Discovery, 25:298-3244
[27] Nayar and Murase, H. (1996). Columbia object image library: Co0l-100. Technical Report CUCS-006 006-96, Department of Computer Science, Columbia University.
[28] Nielsen, M. A. and Chuang, I. L.. (2010). Quantum Computation and Quantum Information: 10th Anniversary Edition. Cambridge University Press.
[29] Oseledets, I. V. (2011). Tensor-train decomposition. SIAM Journal on Scientific Computing, 3315)
[30] Panagakis, Y., Kossaifi, J., Chrysos, G. G., Oldfield, J., Nicolaou, M. A., Anandkumar, A., and Zafeiriou, S. (2021). Tensor methods in computer vision and deep learning. Proceedings of the IEEE, 109(5):863-8)) 890.
[31] Ran, B., Tan, H., Feng, J., Liu, Y., and Wang, W. (2015). Traffic speed data imputation method based on tensor completion. Computational Intelligence and Neuroscience, 2015:22-22 22.
[32] Reyes, J. A. and Stoudenmire, F. M. (2021). Multi-scale tensor network architecture for machine learning. Machine Learning: Science and Technology, 2(3):035036.
[33] Shashua, A. and Hazan, T. (2005). Non-negative tensor factorization with applications to statistics and computer vision. In Proceedings of the 22nd International Conference on Machine learning (ICML), pages 792-799, Bonn, Germany.
[34] Shcherbakova, E. and Tyrtyshnikov, E. (2020). Nonnegative tensor train factorizations and some applications. In Large-Scale Scientific Computing: 12th International Conference. LSSC 2019, Sozopol, Bulgaria, June 10-14, 2019, Revised Selected Papers 12, pages 156-164 Springer.
[35] Sugiyama, M., Nakahara, H., and Tsuda, K. (2017). Tensor balancing on statistical manifold. Proceedings of the 34th International Conference on Machine Learning (ICML), volume 70 of Proceedings of Machine Learning Research, pages
Sydney, Australia.
[36] Sugiyama, M., Nakahara, H., and Tsuda, K. (2018). Legendre decomposition for tensors. In Advances Neural Information Processing Systems 31, pages 8825-8833 Montréal, Canada.
[37] Tucker, L. R. (1966). Some mathematical notes on three-mode factor analysis. Psychometrika, 3113):279-311.
12

---

[38] Yu, J., Zhou, G., Li, C., Zhao, Q. and Xie. S. (2021a). Low tensor-ring rank completion by parallel matrix factorization. IEEE Transactions on Neural Networks and Learning Systems, 32(7):300-33333.
[39] Yu, Y. Xie, K., Yu, J., Jiang, Q. and Xie, S. (2021b). Fast nonnegative tensor ring decomposition based on the modulus method and low-rank approximation. Science China Technological Sciences, 64(9): 1843-18848
[40] Yu, Y., Zhou, G., Zheng, N., Qiu, Y., Xie, S., and Zhao, Q. (2022). Graph-regularizeed non- negative tensor-ring decomposition for multiway representation learning. IEEE Transactions on Cybernetics.
[41] Zhao, Q. Zhou, G., Xie, S., Zhang, L., and Cichocki, A. (2016). Tensor ring decomposition. arXiv preprint 1606.055335
13

---

Appendix
Table of Contents
A Projection theory in information geometry
14
B Formal definitions
16
Theoretical remarks
Convexity and uniqueness many- approximation
Another example: Tensor decomposition
17
Reducing computational complexity by reshaping tensor
17
Information geometric view of LBTC
17
D
Related work
18
E
Additional numerical results
E.
Additional results for COIL dalasel
18
E.2
Additional results for LBTCC
18
E.3
Additional results for comparison with ring decomposition
20
E,4
Comparison with traditional low-rank approximation
20
Dataset details
22
F.1
For experiment with RGB images
22
F2
For experiments with LRTC
22
F33
For comparison with non-negative tensor ring decomposition
22
G
Implementation detail
23
Proposed method
23
Baselines
23
References for Supplementary Materials
25
A Projection theory in information geometry
We explain concepts of information geometry used in this study, including natural parameters, expectation parameters, model flatness, and convexity of optimization. In the following discussion, we consider only discrete probability distributions.
(0, cooordiiatt and geodesies In this study, we treat a normalized D-order non-negative tensor P E R. xT discrete probability distribution with D random variables. Let ie be the set of discrete probability distributions with D random variables. The entire space 24 is non-Euclidean space with the Fisher information matrix G as the metric. This metric measures the distance between two points. In Euclidean space, the shortest path between
points is a straight line. In non- Euclidean space, such shortest path is called a geodesic. In the space U, two kinds of geodesics can be introduced: e-geodesics and m-geodesics. For two points PI, Pa U. e- and m-geodesics are defined as
R2 log Rr = (1- - log Pi + log P: }(1), R,
R, = - (1- + tP2
respectively, where 0 < < 1 and ((1) normalization factor to keep R, to be distribution.
We can parameterize distributions E 2 by parameters known as natural parameters. In Equation have described the relationship between a distribution P and natural parameter vector 0
14

---

(f2.......... ,0 Ip). The natural parameter 0 serves a coordinate system of U. since any distribution in U specified by determining 0. Furthermore, we can also specify a distribution P by its expectation parameter vector 7 (,2.....,........ "...... ). which corresponds to expected values of the distribution and an alternative coordinate system Uf. The definition of the expectation parameter is described in Equation (3). f-coordinates and c-coordinates are orthogonal with each other, which means that the Fisher information matrix G has the following property, Gu,v = / 100, and (G "),,) = 08w/ /mm. e- and m-geodesics can also be described using these parameters as follows.
-0'10'= =1 (1 - 1)0- + 10"² }.
(n'lnn (1-1)=) + +
where 0 and n are 0- and n-coordinaae of distribution P E U.
subspace called e-flat when any e-geodesic connecting two points in subspace is included in the subspace. The vertical descent of an m-geodesic from point P E u onto e-flat subspace B. is called m-projection. Similarly, e-projection is obtained when we replace all with m and m with e. The flatness of subspaces guarantees the uniqueness of the projection destination. The projection destination P or P obtained by minimizes the following KL divergence,
Flatness and projections
e-projection onto B. or B..
P = arg min (K 2), p = arg min DKL QEB..
QEB
The KI. divergence from discrete distributions P E U to Q E is given as
In
DKL Q) =
Σ
Σ ip=1
P
ii..in in log
(13)
Qi.....
A subspace with some of its natural parameters fixed at is e-flat [ Chapter 2], which is obvious from the definition of e-flatness. The proposed many-body approximation performs m-projection onto the subspace B C 24 with some natural parameters fixed to be 0. From this linear constraint, we know that B is e-flat. Therefore, the optimal solution of the many-body approximation is always unique. When space is e-flat and m-flat at the same time, we say that the space is dually-flat. The set of discrete probability distributions 24 is dually-flat.
Natural gradient method c(m)-flatness guarantees that cost functions to be optimized in Equa- tion (3) are convex. Therefore. m(e)-projection onto an e(m))-fat subspace can be implemented by gradient method using second-order gradient. This second-order gradient method is known as the natural gradient method [42]. The Fisher information matrix G appears by second-orde differentiation of the KL. divergence (see Equation (5)). We can perform fast optimization using the update formula in Equation (6). using the inverse of the Fisher information matrix.
Examples for Möbius function In the proposed method, we need to transform the distribution P E R4xx--xxi with and n using the Möbius function, defined Section 2.1 We provide examples here. In Equation 2 The Möbius function is used to find the natural parameter @ from distribution P. For example, if D = 2, it holds that
0 0.r log BigPi,, Pi, _ log -,, log P P1 log P, -1,
0 = log Pin log p, -iig - log P,n.i--1.is log P1n-i2
+ +logP,, 1 1,is +log Pir,iz P + log P1i.i2-1,is-
log P: P -,,2,,,-1,,-11
where we assume Po.iz Pi,.0 and Punig Pi, = Po = 1. To identify the value of Oit....ar we need only P....... with (i (i.....4)) 4 (6 In the same way, using Equation 3 For example, if = 2, it holds that
x {iz 1, x {ia (2)
find the distribution P by the expectation parameter n.
P Mii 1,t1is B1. + 7/ 41,it
Pantaids = Mh 12.55 - -Th 41.12. - f21.ss Th 12.651
+ M,, +., + Th, 41.22att + 12, 1hit1,i2tistt
where we assume 7/1+1 7/1 0 and 71 Mal 7h,.2..441 0.
15

---

a.
b.
A
B
g
88
p
bg
17
26
$7
47 8
Interaction- Interaction-RRpprsenttiie Representation
Tensor Network
Tensor Tree Decomposition
Figure 5: (a) Interaction representation corresponding to Equation (14 and its transformed tensor network for D = 9. (b) Tensor network of variant of tensor tree decomposition.
B Formal definitions
We provide give formal definitions of Legendre decomposition (35
follows:
Definition 1. For a normalized non-negative tensor P and a set of indices B, Legendre decomposition P is Q defined as
arg arg DKI(P, 0.08
where B 12 01.... if (i,, (p))...,ip)
ip) $ B for 0-parameters of 2}.
As in the following definition, particular case of the Legendre decomposition is the many-body decomposition.
Definition 2. For a normalized non-negative tensor P and a natural number m, we define its m-body approximation psm as
Q arg DKz
2:2EB
where B used in B contains all the indices of mm)body parameters.
C Theoretical remarks
C.1 Convexity and uniqueness of many-body approximation
is widely known that the maximum likelihood estimation of ordinary Boltzmann machines without hidden variables is a convex problem. Since we optimize the KI. divergence, the proposed many- body approximation
also a maximum likelihood estimation that approximates non-negative tensor, which is regarded as empirical distribution, by extended Boltzmann machine without hidden variables, as described in Section 2.2 As with ordinary Boltzmann machines, the maximum likelihood estimation of extended Boltzmann machines can be globally optimized. Theorem I is general proof of this using information geometry.
The solution of many-body approximation is always unique, and the objective function of many- body approximation is convex.
Theorem 1.
Proof. As Definition shows, the objective function of the proposed many-body approximation is the KI. divergence from an empirical distribution (a given normalized tensor) to subspace B. We can immediately prove that subspace B is e-flat from the definition of the flatness [ Chapter 2] for any B. The KL divergence minimization from a given distribution to e-flat space is called m-projection (see the second paragraph in Section A). The m-projection onto e-flat subspace known to be convex and the destination is unique (see the third paragraph in Section A). Thus, the optimal solution of the many- body approximation always unique, and the objective function convex.
16

---

m-flat
D
m-flat
D
Not unique
e-flat
Low-rank manifold B
Low-body manifold B
Figure 6: The em-algorithm. Dashed(Solid) arrows are e(m)-projection.
C.2 Another example: Tensor tree decomposition
As seen in Section 2.4, cyclic two-body approximation coincides with tensor ring decomposition with constraint. In this section, we additionally provide another example of correspondence between many-body approximation and the existing low-rank approximation. For = 9, let us consider three-body and two-body interactions among (i1, i2, is), (i1, 15, dG), and (ł7, ig, ig) and threc-body approximation among (1g, ig). ig). We provide the interaction representation of the target energy function in Figure a). In this approximation, the decomposed tensor be described as
Pi.....o A2, Be Ag Csy Sig.5.
(14)
Its converted tensor network is described in Figure D a). It becomes tensor network of tensor tree decomposition in Figure (bbb when the region enclosed by the dotted line is replaced with new tensor (shown with tilde) by coarse-graine graining Such tensor tree decomposition is used in generative modeling [43] computational chemistry 48 and quantum many-body body physics [50].
C.3 Reducing computational complexity by reshaping tensor
As described in Section 2.1.1 the computational complexity of many-body approximation is B))) where is the number of iterations, because the overall complexity is dominated by the update of 6, which includes matrix inversion of and the complexity of computing the inverse of an n matrix is (n²)) This complexity can be reduced if we reshape tensors so that the size of each mode becomes small. For example, let us consider third-order tensor whose size is (J2, J2, J2) and its cyclic two-body approximation. In this case, the time complexity is O()J12) since it holds that B| = (JJ)) (see Equation (12)). By contrast, we reshape the input tensor to sixth-order tensor whose size is JJ, J, J, J), the time complexity becomes (()Jo) since it holds that B (J2).
This technique of reshaping a tensor into a larger-order tensor is widely used practically in various methods based tensor networks, such as tensor ring decomposition [47). We use this technique in the experiment in Section 3.3 to handle large tensors in the both of proposal and baselines.
C.4 Information geometric view of LBTC
Traditional methods for tensor completion are based on the following iteration, called em- algorithm (51.
e-step: Pa Ta:
m-step: P +- low ---.oorrakkppproxiiaiiioooo
with an initialized tensor P. where Ω is the set of observed indices of given tensor T with missing values. From the information geometric view, m-step does m-projection onto the model manifold B consisting of low-rank tensors, and e-step does e-projection onto the data man- ifold D {PI Pa- To It is straightforward to prove that e-step finds point p arg millpe (xI(P, (P, (P)P) for the obtained tensor P by the previous m-step 52 The iteration of e- and m-step finds the closest point between two manifolds B and D. Since we can immediately
prove that tQ (1 t) 22 belongs to D for any Q1, Q2 D and t < 1, the data manifold D
17

---

Figure 7: (a) Many-body approximation with three-body interaction among width, height, and index image can reconstruct the shapes of each image. (b) Without the interaction, many-body approximation cannot reconstruct the shapes of each image.
m-flat, which guarantees the uniqueness of the c-step. See the definition of m-flat in Section However, since low-rank subspace is not e-flat, that m-projection is not unique. LBTC conducts tensor completion by replacing low-rank approximation in m-step with many-bodk approximation. Since the model manifold of many-body approximation is always e-flat and data manifold is m--fat, each step in the algorithm unique.
D Related work
As discussed in Section I we assume low-rankness in the structure of the tensor in conventional decompositions. However, since tensors have many modes, there are too many kinds of low-rank structures. Examples include, the CP decomposition [15], which decomposes a tensor by the sum of Kronecker products of vectors; the Tucker decomposition (37]. which decomposes tensor by core tensor and multiple matrices; the tensor train decomposition [29], which decomposes a tensor by multiple core tensors; and tensor-ring decomposition [41], which assume periodic structure on tensor train decomposition are well known and popularly used in machine learning.
However, in general, it is difficult find out the low-rank structure of given data set. Therefore, it difficult to say how to choose low-rank decomposition model from those options. To alleviate this difficulty in model selection and provide intuitive modeling, we have introduced the concept of body and formulated a decomposition that does not rely on ranks.
E
Additional numerical results
E.1
Additional results for COIL dataset
In Section 3.1 we see that interactions related to color control color-richness the reconstruction. We provide additional experimental results in Figure 7 to the role of the interaction among width, height, and image index. The figure shows the influence the interaction for reconstruction. While Figure b), without the interaction, could not capture the shapes of images, Figure 7 () reconstruct each shape properly. Based on these results, is reasonable that the interaction is responsible for these shapes.
E.2 Additional results for LBTC
We provide the reconstructions by the experiments in Section 3.2 in Figure with ground truth.
In traditional low-rank approximation, tuning hyper-parameters based on prior knowledge difficult because its semantics are unclear. By contrast, the interactions in many-body approximation be tuned easily by considering the expected correlation among modes of given tensor. To support this claim, we also confirm that LBTC scores get better when we activate interactions between modes that are expected to be correlated based on prior knowledge. In this experiment, we use the same datasets Section 3.2 We prepare two set of interactions, Interactions A and B, shown on the left in Figure 9 In Interaction A, we activate interaction between modes such as (date, hour) and (hour, minute, lane), which we expect be correlated. In Interaction B, we activate interaction between modes such as (date, minute) and (minute, lane), which we do not expect to be correlated. For fair comparison, we chose these interactions so that the number parameters used in Interactions A and B are as close as possible, where the number of parameters for Interaction A was 1847, and the number of parameters for Interaction was 1619. seen in Table [ LBTC with Interaction always has better score than Interaction B, which implies that many-body approximation has better
18

---

performance when we choose interactions among correlated modes. We conducted the experiment 10 times and reported mean values with standard error. We also visualize the reconstruction Figure y to show that LBTC with Interaction A can complete tensors more accurately than that with Interaction B.
As seen in Figure [4 and Section G.2 existing methods require ppper-prraetters that are difficult to tune since we usually cannot find any correspondence between these values and prior knowledge about the given tensor. In contrast, as demonstrated in this section, our energy-based model is intuitively tunable that can skip time-consuming nyper-parameter tuning.
Table 2: Dependency of recovery fit score of LBTC on choice of interaction.
Scenario 1
Scenario 2
Scenario 3
Interaction A
0.900 0.00125
0.889 + 0.000632 0.819 + 0.000699
0.864 + 0.001264 0.805 + 0. 0.000527
Interaction B
0.820 = 0.00143
a.
Scenario
Proposal |
90
Lane
90-
Lane 2
Sood
30
30
80
30
30
b.
Scenario
Sped
22
29
Scenario
90
Lame
Lanc
Nee trqu
60
4
29
90
8
1 1 oo M
1pp1 10000
/ w 1 MMY
8
22
Figure& 8: Ground truth and reconstruction by LBTC in the experiment in Section Filled areas are missing parts, and blue lines are estimated values by the proposal. Green lines are ground truth.
19

---

Scenario
Sconaria
Scenario
A
B
Figure Interaction dependency of LBTC. Filled areas are missing parts, and blue lines are estimated values by the proposal with Interactions A and B, respectively. Green lines are ground truth. Due to space limitations, we show the reconstruction Lane from day to 20, where many missing parts exists.
E.3 Additional results for comparison with ring decomposition
We additionally provide Figure TO which shows running time compassion the experiment in Section As described in the main text, repeat five times each baseline method with random initialization while the proposal runs only once. We compare the total running time of them. As seen in Figure [00] for both synthetic and real datasets, the proposed method always faster than baselines, keeping the competitive relative errors. The flatness of solution space enables many-body approximation be solved with the Newton method. In contrast, NTR is solved by first-order derivative. The quadratic convergence of the Newton method is well known, and it requires only few iterations find the best solution in the flat solution space B Beye. That is the reason why the proposed method is fast and stable.
E.4
Comparison with traditional low-rank approximation
(30,30,30,30,30) R15
(20,20,20,20,20.0)) R10
b.
ChartRess
TT. TT.Orrgaai
0.020
0.3200
Relative Erro
0.0100
Relative Ener
Relelive Error
Relative Error
IraM Proposal
0.015
0.0075
0.32
0.3150
0.0050
0.31
0.0025
0.005
0.3100
10
10a
10
10
600
1200
1800
500
1000
1500
Parameters
Parameters
Parameters
Parameters
Running Time (sec.)
10
Running Time (sec.)
Running Time (sec.)
Running Time (sec.)
1033
10a
10a
102
10a
10a
10a
10a
10a
600
1200
1800
500
1000 1500
2000
Parameters
Parameters
Parameters
Parameters
Figure 10: Comparison of reconstruction capabilities (top) and computational speed (bottom) of cyclic two-body approximation and tensor ring decomposition for synthetic low-tensor ring datasets (a) and real datasets (b). The vertical red dotted line is [B]] (See Equation (22
low-rank approximation, we often gradually increase the rank enlarge model capability. In the same way, if we have no prior or domain knowledge about a tensor, we can gradually increase of m-body approximation and get more accurate reconstruction. To demonstrate that, we conducted We compared the reconstruction ability
additional experiments with the same dataset in Section 3.2
20

---

a.
b.
1.00
1.00
NTTF Poppol
Proposal
Prsporal
0.95
0.95
0.95
0.95
E 0.90
E: 0.00
0.90
0.80
10
0.80
10" Perameters
0.80
10a
107
10"
0.80 10a
107
10*
Parameters
Parameters
Parameters
Figure 11: We compared the fit score of proposal pDl, paye ps? ps3 with baselines. (a) The ranks for baselines assumed as (r, (r,)..,?) r) are tuned taking the tensor into size account.
with changing the value of r. (b) The ranks for baselines
One body
Cyclic Approximation
Body Approximation
Approximation
Origial
Spuad (mipk)
14
21
21
11 Hour
17
23
0
11
17
23
11 Hour
17
23
0
11 Hour
17
23
5
11
17
23
Hour
Hour
Figure 12: Heatmaps for reconstructed tensor p paye ps pS3 and input tensor P. We visualize the traffic on only the first lane at 00 minutes every hour, which corresponds to 1,
of the proposed m-body approximation with non-neg non-negative tucker decomposition (NTD) and non-negative tensor train decomposition (NTTF) 34
In many applications based matrix and tensor factorization, their objective functions often the fit score (RMSE) and they try to find good representations to solve the task by optimizing the objective. Therefore, for fair evaluation without depending on specific applications, we examine the performance by fit score.
For m-body approximation, we gradually increased the order of the interaction by changing to 1, 2.33 We ran the proposal only once since the proposal does not depend on initial values. For baselines, gradually enlarged model capacity by increasing the target rank. We repeated each baseline 20 times with random initialization. We defined the Tucker and train ranks for the baseline method in ways. First, assumed that all values in ranks are the increase the value of r. However, since the size an input tensor is (28, 24, 12, 4) and each mode size is not unique, it can be unnatural rank setting, which may also adjusted target rank by taking the tensor size into account. In particular, the rank is set to be 11,1,1,1,, (6,2,2.1), 1), (6,5,2,1), (8,5,4.2), (10,7,4,2), (12,9,42),, (16, 12,6.2), (18.15,8,2) (20,18,11.3) for NTD and (1,1,1), (1,2,1), (1,2.2), (2,2,2), (4,2,2), (6,2,2), (6,3,2), (8,3,2), (10,3,2), (10,4.2), (15,4,2), (20,5.2), (20,8 (30,10.2), (30,15,8) for NTTF.
that is, (r, (r,...,,r) and
worse RMSE. To avoid that,
We evaluated the approximation performance with the fit score 1 |P P and the reconstructed tensor P. We show these results in Figure 11 While baselines are often trapped in worse local minimum, the proposal reconstructs the original tensor better. We also show the reconstruction by the proposed method as heatmaps in Figure 12 Even one-body approximation able to reconstruct the rush hour around 17:00. When the number of interactions increases, observe the trend that the time the busy hours on weekends are different from those on weekdays more clearly.
P| for the input
To see the effectiveness of proposal for data mining, we also plotted two-body interaction between day and lane extracted by the two-body approximation of the above traffic data. This is matrix
21

---

( E 2 As seen in Figure 13 the value in the first column is always the largest and the value the fourth column always the smallest. From this result, we can interpret that the data corresponding to :, ;, were taken in the overtaking lane and the data for Pl were taken the overtaken lane.
F
Dataset details
We describe the details of each dataset below.
F1 For experiment with RGB images
We downloaded the COIL-100 dataset from the official website We used 10 images with object numbers 4, 5, 7, 10, 17, 23, 28, 38, 42, and 78. This selection aimed to include both monochromatic and multi-colored objects so that we can see how the proposed method changes the reconstruction we deactivate the interactions related to color. We reshaped the size of each image to 40 x 40.
F.2 For experiments with LRTC
We downloaded traffic speed records in District 7, Los Angeles County, from PeMS ? The data are in the public domain. The whole period of the data lasts for 28 days from February 4 to March 2013. By measuring the speed vehicles four lanes at a certain point on the highway 12 times five-minute intervals, 12 4 matrix is generated per hour. Twenty-four of these matrices are generated per day, so the dataset is a 24 x 12 x 4 tensor. Further, since we used these tensors for 28 days, the size of the resulting tensor is 28 x 24 12 We artificially generated defects in this tensor. As practical situations, we assumed that a disorder of lane's 's speedometer causes random and continuous deficiencies. We prepared three patterns of missing cases: case with a disorder in the speedometer in Lane 1 (Scenario 1), a case with disorders in the speedometers in Lanes 1, 2, and 3 (Scenario 2), and case with disorders all speedometers (Scenario 3). The missing rates for Scenarios 1, 2, and 3 were 9 percent, 27 percent, and 34 percent, respectively. The source code for imparting deficits are available in the Supplementary Materials.
F.3 For comparison with non-negative tensor ring decomposition
Synthetic datasets We created D core tensors with the size (R, I, R,) by sampling from uniform distribution. A tensor with the size Io and the tensor ring rank ......) R) was then obtained by the product of these D tensors. Results for R 15, D 5, 30 are shown in the left column in Figure[4 and those for R = 10, D = 6, / = 20 in the right column in Figure For all experiments on synthetic datasets, change the target ring-rank as (r: for 2.3,.... ,9 for baseline methods.
Real datasets TT_ChartRes originally a (736, 736,31 31) tensor, produced from TokyoTech 31- band Hyperspectral Image Dataset. We downloaded ChartRes mat from the official website and reshaped the tensor as (23, 8, 4, 23, 8, 4, 31)... For baseline methods, we chose the target ring- rank as (2, 2, 2, 2, (2.2, (2.2,,,, 2,8), (3. 2, 2, 3, 2, 2, 5). (2,
2, 9),
(3,2, 2, 2, 2,2,9). (3, 2, 2, 3, 2, 2, 10), 2,, 2, 15), or (3, 2, 2, 16). (512, 512, 59) tensor, which is produced from Tokyo'Tech 59-band Hyperspectral Image Dataset. We downloaded Origami. mat from the official website In TT IT-_rrigam, 0.0016 percent of ele- ments were negative, hence we preprocessed all elements of TT_ "T_Origaai by subtracting - -.0000664 the smallest value in TT_Origami, make all elements non-negative. We reshaped the tensor as (8, 8, 8, 8, 8, 8, 59). For baseline methods, we chose the target ring-rank as (2, 2, 2, 2, 2, 2, for r 2, ... 15. These reshaping reduces the computational complexity as described Section to complete the experiment in reasonable time.
(3, 2, 2, 3, 2, 6). 2,2, (3, 2, 2, 2,
(3, 2,,,,,4,2,,, 2,, 23,2,2,4,2,,,3
TT rT_Origami is originally
http:://pem.doo.c.. gov/
22

---

(d-1) interaction in two body approximation
Interaction intensity
1.0
0.9
0.8
0.7
Lane
Lane
0.6
Lane
Lane
0.5
10
15
20
25
Day
Figure 13: The obtained interaction between day and lane in two-body approximation.
G Implementation detail
We describe the implementation details of methods in the following. All methods are implemented in Julia 1.8.
Environment Experiments were conducted on Ubuntu 20.04.1 with a single core of 2.1GHz Intel Xeon CPU Gold 5218 and I 128GB of memory.
G.1 Proposed method
Many-body approximation We a natural gradient method for many-body approximation. The natural gradient method uses the inverse of the Fisher information matrix to perform second-order conducted the following
optimization in non-Euclidean space. For non-normalized tensors, procedure. First, we computed the total sum of elements of an input tensor. Then, we normalized the tensor. After that, we conducted Legendre decomposition for the normalized tensor. Finally, we obtained the product of the result of the previous step and the total sum we computed initially. That pro- cedure was justified by the general property of the KI. divergence AD (P. Q) = D (XP, XQ), for any À x- The termination criterion is the same as the original implementation of Legendre decomposition by [36]; that is, it terminates if |n ; I| < 10-5, where n,," is the expectation parameters t-th step and n is the expectation parameters of an input tensor, which are defined in Section 2.1. The overall procedure is described in Algorithm 2 Note that this algorithm is based on Legendre decomposition by [36].
LBTC
The overall procedure described in Algorithm with our tolerance threshold. We randomly initialized each missing element in p(t=1) by Gaussian distribution with mean 50 and variance 10.
G.2 Baselines
Baseline methods for tensor completion We implemented baseline methods for tensor completion by translating pseudo codes in original papers into Julia code. As baseline methods, used HaLRTC, SiLRTC 21). SiLRTC-TT (2). and PTRCRW 38). SiLRTC and HaLRTC are well-known classical methods based HOSVD. SiLRTC-TT is based tensor train decomposition. PTRCRW is based on tensor-ring decomposition and is recognized as the state of the art. HaLRTTC requires real value p and a weight vector a. We used the default setting on Q described in the original paper and tuned the value as p 10-5, 5 110- for d = 1, 2, 3, This is the threshold truncated SVD for the mode-d matricization of the input tensor. Following the original paper, we assumed =
102. SiLRTC requires a hyperparameter na
and we tune 7 in the range 10*). SiLRTC-TT requires a weight vector a and real number f. We set the value of a as the default setting described in the original paper. We tune the value of in the range of [0.001, 0.01, 0.025, 0.1, 0.25, 1 J. PTRC-RW requires two weight vectors, a and B, step length and tensor ring rank as hyper parameters. Following the original paper, we used the default setting for a, B, and d. Then we tuned the tensor ring rank R, R, R, R) within the range E 2,,2,... 6). assuming that it has the same values in each element, which is also assumed in the original paper.
of (10-2 2 10
23

---

Algorithm 2: Many-body approximation
MANTBODYAPPROXMMTTTOTT B) + Total sum of
Obtain normalized input tensor P - T/ss
|| "1""/" denotes element-wise division
Compute '7 of P using Equation 3
Initialize 0
I/ e.g. O. for all E B
6-1
repeat
Compute Q, using the current parameter 0 with Equation
Compute 7, from Q using Equation (3)
Compute the inverse of the Fisher information matrix G using Equation
0
G- - -n)") ))")
(+1 1.11 1.+11
until || B -1p1 <E
// We set 10 in implementation:
T Q: 8
",""* denotes element-wise multiplication
return T
Since PTRC-RW is not convex method, we repeated it 10 times for each with randomly sampled initial values by Gaussian distribution with mean 50 and variance 10. We provide Figure 14] which shows the sensitivity of the tensor completion performance with respect to changes of the above syper-parameters. All of the above baseline methods include iterative procedures. These iterations were stopped when the relative error between the previous and the current iteration is less than 10-5 or when the iteration count exceeds 1500 iterations.
Baseline methods for tensor ring decomposition We implemented baseline methods for tensor ring decomposition by translating MATLABB code provided by the authors into Julia code for fair comparison. we can see from their original papers, NTR-APG, NTR-HALS, NTR-MU, NTR-MM, and NTR-IraMM have an inner and outer loop to find a local solution. We repeated the inner loop 100 times. We stopped the outer loop when the difference between the relative error the previous and the current iteration is less than 10-4 NTR-MM and NTR-IraMM require diagonal parameter matrix E. We define = w/ where I an identical matrix and w 0.1. The NTR-IraMM method performs low-rank approximation to the matrix obtained by mode expansion of input tensor. The target rank was set to be 20. This setting is the default setting in the provided code. The initial positions of baseline methods were sampled from uniform distribution (0, 1).
Baseline methods for low-rank approximation We implemented NNTF by translating the pseudo- code in the original paper (4] into Julia code. The algorithm requires NMF in each iteration. For that, we also implement the most standard method, multiplicative update rules, following the paper [46]. We use default values of hyper parameters of sklearn 49] for NMF. For NTD, translate TensorLy implementation [55 to Julia code, whose cost function is the Frobenius norm
between input and output.
24

---

HaLRTC
SILRTC
SILRTC-TT
PTRC-RW
1.0
1.0
1.0
Scenario
0.5
0.5
0.6
0.4
0.4
0.4
0.4
0.2
0.2
0.2
0.0
10
10
10"h
10
10"
10"
10
10
10
10
10
10
Rank
1.0
Scenario
0.0
0.0
0.0
0.6
0.4
0.4
0.4
0.4
0.2
022
0.2
0.2
0.0
10
10
Rank
1.0
Scenario 3
0.8
0.8
0.6
0.6
0.6
0.8
04
0.4
0.4
0.4
0.4
0.0
0.0
Rank
Figure 14: Hyper-parameter dependencies for the performance of baselines in the experiment in
Section 3.2
References for Supplementary Materials
[42] Amari, S. 1998. "Natural gradient works efficiently in learning." Neural computation 10(2):251- 276.
[43] Cheng, Song, Lei Wang, Tao Xiang and Pan Zhang. 2019. "Tree tensor networks for generative modeling." Physical Review B 1911511151131.
[44] Kim, Yong-Deok and Seungjin Choi. 2007. Nonnegative tucker decomposition. In 2007 IEEE Conference on Computer Vision and Pattern Recognition. IEEE pp. 1-8.
[45] Kossaifi, Jean, Yannis Panagakis, Anima Anandkumar and Maja Pantic. 2019. "TensorLy: Tensor Learning Python." Journal of Machine Learning Research 20(26):1-6.
URL: http:///iiil.orrgppperrs/////// 2ttp://///////22...... 8-277.22mmm
[46] Lee, Daniel and H Sebastian Scung. 2000. Algorithms for non-negative matrix factorization. In Advances Neural Information Processing Systems 13, Denver, United States: pp. 556-562.
[47] Malik, Osman Asif and Stephen Becker. 2021. A sampling-hased method for tensor ring decomposition. In International Conference Machine Learning. PMLR pp. 7400-7411.
[48] Murg, Valentin, Frank Verstraete, Reinhold Schneider, Peter R Nagy and O Legeza. 2015. "Tree tensor network state with variable tensor order: An efficient multireference method for strongly correlated systems." Journal of Chemical Theory and Computation 113):-027-1036.
[49] Pedregosa, F., G. Varoquaux, A. Gramfort, Michel, Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot and E. Duchesnay. 2011. "Scikit-learn: Machine Learning Python." Journal of Machine Learning Research 12:2825-8833
25

---

[50] Shi, Y-Y, L-M Duan and Guifre Vidal. 2006. "Classical simulation of quantum many-body systems with tensor network." " Physical review 74(2):022320.
[51] Walczak, B. and D.L. Massart. 2001. "Dealing with missing data: Part 1." Chemometrics and Intelligent Laboratory Systems 58(1):15-27.
URL: htt.://wwww..ttteooo.o..o./tttta
[52] Zhang, Sheng, Weihong Wang. James Ford and Fillia Makedon. 2006. Learning from incomplete ratings using non- non-negative matrix factorization. In Proceedings of the 2006 SIAM International Conference on Data Mining, SIAM Bethesda, United States: pp. 549-553.
26


Reviews:
[{"review_id": "eldMMpP0wG", "reviewer": "Reviewer_xQBv", "summary": "This paper proposes a novel many-body decomposition, based on Legendre decomposition, for non-negative tensor. Instead of specifying a rank parameter, the proposed method use energy-based model to the interactions among modes.", "questions": "See weakness section", "limitations": "1. As the authors already stated in the paper, the proposed method can only be applied to non-negative tensors.\n2. Compared to rank-based tensor decomposition, the proposed method seems to require more expert knowledge about the application itself.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The paper is overall well-written and easy to follow. Especially the figures that explain the interactions.\n2. The work is technically sound and aims at an interesting angle of tensor decomposition which is worth more attention from the research community.\n3. The proposed method doesn't require the rank parameter which can save a headache of parameter tuning in many existing tensor analyses.", "weaknesses": "1. The novelty of the proposed method is limited given the existing Legendre decomposition.\n2. The empirical advantages of the proposed methods are not sufficiently demonstrated. \n3. The rank parameter is not necessarily required in traditional tensor methods, for example, tensor nuclear norm based approach.\n4. The details of the proposed algorithm, LBTC, should be more emphasized in the main text. \n5. How is the computational cost of the proposed approach, compared to the traditional tensor decompositions? There were many fast methods have been proposed for rank based decomposition, such as tensor CUR decompositions.\n6. If I understand right, the proposed method is indeed a generalized rank-1 approximation.\n7. \"interactions between modes\" -> \"interactions among modes\"\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your positive feedback. We provide our point-by-point response to each of your comments in the following. \n\n### In the \"Weaknesses\" section\n\n\n> 1. The novelty of the proposed method is limited given the existing Legendre decomposition. \n\n\nAlthough our proposal is based on Legendre decomposition, we have the following nontrivial novel contributions:\n\n- As described in lines 48-49, Legendre decomposition does not provide factorization, while the proposed many-body approximation factorizes tensors into a product form that enables us to observe components of the data. For example, as seen in Figure 3, our proposal factorizes the COIL image dataset into two components, the color tendency of each object (e) and the shape of each object (f). We also provide another example with the traffic dataset in Figure 13 in Appendix.\n\n- As described in lines 44-45, we have revealed that the $n$-body parameters control interaction among modes. Such a tight connection between natural parameters and interaction between modes was not known in literature of Legendre decomposition.\n\n- We have proved a theoretical connection between Legendre decomposition and low-rank approximation (Section 2.4).\n\n\n> 2. The rank parameter is not necessarily required in traditional tensor methods, for example, tensor nuclear norm based approach.\n\nThank you for pointing this out. You are right that nuclear norm-based methods do not require the rank parameter. We will add the following text in line 32:  \nAlthough there is an attempt to avoid rank tuning by approximating it by trace norms [2], it also requires another hyperparameter, the weights of the unfolding tensor, hence such an approach does not fundamentally solve the problem.\n\n--\n\nPlease note that nuclear norm-based methods require other hyperparameters instead of ranks. In our experiments, the baseline methods SiLRTC and SiLRTC-TT in Section 3.2 are typical nuclear norm-based methods, and their performance depends on hyperparameters, as explained in Section E.2 and Figure 14. These dependencies were also noted in their original papers (for example please refer to Section VI. A in [2]).  \n\n[2] Bengua, J. A., Phien, H. N., Tuan, H. D., and Do, M. N. (2017). Efficient tensor completion for color image and video recovery: Low-rank tensor train. IEEE Transactions on Image Processing, 26(5):2466–2479\n\n\n> 3. The empirical advantages of the proposed methods are not sufficiently demonstrated.\n\nWe have demonstrated the empirical advantages of the proposed method in terms of hyperparameter tuning and model selection in Sections 3.2 and C.2, initial value dependency in Figure 11 in Appendix, and tensor completion in Section 3.2, which can be significant merits in practical situations.\n\n> 4. The details of the proposed algorithm, LBTC, should be more emphasized in the main text.\n\nThank you for your suggestion. In the camera-ready version, we will add a new section 2.6 to describe LBTC, to which we put lines 269-280 and Algorithm 2.\n\n> 5. How is the computational cost of the proposed approach, compared to the traditional tensor decompositions? There were many fast methods have been proposed for rank based decomposition, such as tensor CUR decompositions.\n\nAs described in line 117, our method is based on the natural gradient, and its computational complexity is cubic with respect to the number of parameters to be optimized, that is, the total complexity is $O(\\gamma |B|^3)$, where $\\gamma$ is the number of iterations. Due to the second-order convergence property of Newton's method, the iteration number $\\gamma$ is usually small (Section C.3 in Appendix). For many low-rank decomposition methods, the computational complexity is cubic with respect to the rank. For example, the computational cost of tensor train decomposition of a $d$-order tensor of size $n×n× \\dots ×n$ is $O(dnr^3)$, where the rank is $(r,r,…,r)$ [a].\n\nPlease note that your suggested acceleration techniques based on SVD, such as CUR decomposition, cannot directly incorporate the non-negative condition in decomposition. As seen in Figure 10 in Appendix, we have also compared runtime of cyclic two-body approximation methods, NTR-MM and NTR-lraMM, that are known as faster methods than traditional tensor ring decomposition [32].  \n\n[a] Oseledets, Ivan V. \"Tensor-train decomposition.\" SIAM Journal on Scientific Computing 33.5 (2011): 2295-2317.  \n[32] Yu, Y., Xie, K., Yu, J., Jiang, Q., and Xie, S. (2021b). Fast nonnegative tensor ring decomposition based on the modulus method and low-rank approximation. Science China Technological Sciences, 64(9):1843–1853.\n\n\n> 6. If I understand right, the proposed method is indeed a generalized rank-1 approximation.\n\nYour understanding is correct. We have discussed that point in Section 2.5. We have shown in this paper the significant potential of a generalization of rank-1 approximation by formulating it as a many-body approximation, which has comparable results with traditional non-convex low-rank approximation.\n\n> 7. \"interactions between modes\" -> \"interactions among modes\"\n\nThank you for pointing out that. We will revise that in the camera-ready version. \n\n### In the \"Limitations\" section\n\n> Compared to rank-based tensor decomposition, the proposed method seems to require more expert knowledge about the application itself\n\nIn low-rank approximation, we often gradually increase the rank to enlarge model capability. In the same way, if we have no prior or domain knowledge about a tensor, we can gradually increase $m$ of $m$-body approximation and get more accurate results. Hence, even when the user knows nothing about the data, the proposed method is also easy to use. "}, {"review_id": "mcOYEO99wN", "reviewer": "Reviewer_59t8", "summary": "The authors introduce a new, energy-model based approach to the decomposition of non-negative tensors. They compare their method to mainstay techniques.", "questions": "My questions all relate to the similarities and structures of the other techniques.", "limitations": "This is not an issue. The work presented by the authors is a general tool.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The technique seems solid and reasonable, but how solid and reasonable is hard to assess (see below).", "weaknesses": "As presented, it's hard for me to determine how novel and distinct the presented method is. I work extensively with tensors, but as practical tools, so I am not fully versed in the various techniques and advantages for the various decomposition technique. For me to determine this, I would need a much more extensive literature review (perhaps adding a table) and to compare the new technique much more thoroughly and extensively to existing techniques (to highly the similarities, differences, advantages, and disadvantages).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments about our novelty and impact. \n\nAccording to your suggestion, we have prepared a table comparing the proposed method with major tensor decomposition methods, the CP decomposition, Tucker decomposition, tensor train decomposition, and tensor ring decompositions. If there is any other information to be added, please let us know.\n\n\n| Models | Global Optimal | Uniqueness  of solution | Required parameter | Remarks |\n| ------------- | ------------- | ------------- | ------------- | ------------- |\n| CP decomp. | No | No  | CP-rank $r$  | NP-hardness, ill-posed |\n| Tucker decomp. | No | No  | Tucker rank $(r_1, \\dots,r_D)$  | |\n| Tensor train decomp. | No | No  | Train rank  $(r_1, \\dots,r_{D-1})$   | Optimal rank is imbalanced |\n| Tensor ring decomp. | No | No  | Ring rank $(r_1, \\dots, r_D)$  | Cyclic structure |\n| Many-body Approx. | __Yes__ | __Yes__  | Interactions, or $m$ of $m$-body | Only for non-negative tensors|"}, {"review_id": "2clnbZvNQz", "reviewer": "Reviewer_oVis", "summary": "This work introduces a novel approach to non-negative tensor decomposition, termed \"many-body approximation,\" which specifically addresses the relationship among modes of tensors. It is formulated as a variant of Legendre decomposition and realized through globally minimizes the Kullback-Leibler divergence. Additionally, this work explores the connections between the proposed method and existing methods, using interaction representation. A series of experimental results substantiate the superiority of the proposed method over other existing ones in dealing with tensor completion applications.", "questions": "Is it possible to derive some theoretical results to support the merits of the proposed \"many-body approximation\" method in dealing with tensor completion and approximation problems?", "limitations": "See the Weakness\\Questions above.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper is comprehensive and includes all the essential components. The structure is logical, and the design principle is effectively and adequately illustrated.\n2. The proposed method is a novel rank-free method, which is particularly advantageous as identifying an appropriate low-rank structure can be highly challenging in some practical applications.\n3. This paper visualizes the presence or absence of interactions between modes and further demonstrates the interpretability of interactions among all the modes.", "weaknesses": "1. The proposed decomposition is only applicable to non-negative tensors. This may restrict the method's effectiveness and generalizability when dealing with tensors that contain negative values.\n2. Some theoretical analysis about the proposed \"many-body approximation\"method is missed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your positive feedback. \n\n> Some theoretical analysis about the proposed \"many-body approximation\" method is missed. Is it possible to derive some theoretical results to support the merits of the proposed \"many-body approximation\" in dealing with tensor completion and approximation problems?\n\nAs described in lines 112-114 in the main text and the second paragraph titled “Flatness and projections”  in Section A, we have derived theoretical results of many-body approximation using information geometry, the uniqueness of the solution and the global convergence of the natural gradient method. These theoretical results directly lead to practical merits, that is, we do not have to run the algorithm multiple times with random initial values to find a better approximation.\n"}, {"review_id": "MERgg49Ya1", "reviewer": "Reviewer_gp7Q", "summary": "The paper proposes a new type of tensor decomposition, based on considering an interaction network among the different modes of a tensor. Nonnegative tensors are considered and factorized based on a Legendre decomposition, which may be viewed as representing the tensor elements by an underlying probability distribution. The paper then proposes to model the tensor in terms of interactions, each of which couples two or more modes. These interactions may be represented by a tensor network, which contains no rank modes, e.g., T_{ijkl}=U_{ij}V_{jk}W_{kl}. By contrast, typical tensor decompositions, include auxiliary/contracted modes. Using the theory of probability distributions, the authors argue that the resulting optimization problem to factorize the tensor according to interaction network, is convex. A gradient based method is used for optimization. Connections are made to other types of tensor networks, demonstrating that the interaction network arises by imposing additional structure on such networks. Experimental results show that, on the datasets considered, the interaction/many-body approximation method can be competitive in terms of accuracy.\n", "questions": "Please provide a counterargument to my assertion above that the optimization problem is in fact nonconvex. If I have somehow overlooked something on this, I would be much more supportive of the paper. If its indeed nonconvex, the paper needs to be corrected and re-evaluated.\n", "limitations": "None.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": " * The proposed tensor decomposition type appears to be new and is interesting / has connections to probability theory and physics.\n * The use of probability theory to demonstrate convexity of the optimization problem is innovative, however, it appears to me that the assertion of convexity is incorrect, see weaknesses. [After reviewing the rebuttal, I follow at least the high-level motivation for why the approach is convex and revised my evaluationp]\n\n", "weaknesses": " * The Legendre decomposition is defined vaguely in the paper and supplementary material, and as far as I can tell this is not a broadly known concept (discussed only in a recent line of ML papers). This makes the paper not self-contained, and makes it difficult to ensure correctness of the arguments. The optimization method and proof of convexity of the optimization problem are also proven hastily / with lack of good definition/theorem/lemma/proof structure, in particular restate the theorems being used from prior work (e.g., when discussing results based on linear dependence of parameters in the distribution).\n * My main concern is that from a tensor point of view, it seems impossible that decomposing a tensor into a many-body approximation for an arbitrary interaction network leads to a convex optimization problem. Here is a counterexample. Consider a sixth order tensor with 3 interaction terms, each acting on two independent modes, T_{ijklmn}=U_{ij}V_{kl}W_{mn}. If the dimensions of modes j, l, and n are 1, we've obtained a rank-1 approximation of the tensor T'_{ikm}=T_{i1k1m1}=U_{i1}V_{k1}W_{m1}. Its well known that rank-1 approximation, which also corresponds to the tensor eigenvalue problem, is nonconvex (often has many local minima). [After reviewing the rebuttal, I see that KL-divergence enables convexity based on prior results.]\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate your careful reading and constructive comments. As described below, our proposal is ensured to be convex.  \n\n> My main concern is that from a tensor point of view, it seems impossible that decomposing a tensor into a many-body approximation for an arbitrary interaction network leads to a convex optimization problem. Its well known that rank-1 approximation, which also corresponds to the tensor eigenvalue problem, is nonconvex.\n\nOur argument is correct because we use not the Frobenius norm but the KL divergence. As you demonstrate and we describe in Section 2.5, our proposal includes rank-1 approximation. However, instead of the Frobenius norm, our proposal optimizes the KL divergence from an empirical distribution (a given normalized non-negative tensor) onto a model manifold (a set of interaction-reduced tensors). When the cost function is defined with the Frobenius norm, the best rank-1 approximation is NP-hard, as seen in Equation (28) in [a], which is non-convex. In contrast, when the cost function is the KL divergence, it becomes a convex problem, as shown in [b]. To avoid the misleading, we will add the following sentence in line 243 in the camera-ready version:  \nFinding the rank-1 tensor that minimizes the Frobenius norm from a given tensor is known to be an NP-hard problem [a]. However, it has been reported that finding the rank-1 tensor that minimizes the KL divergence from a given tensor is a convex problem [b].\n\n[a] DOI: 10.1145/2512329  \n[b] DOI: 10.1109/ACSSC.2017.8335432\n\n\n> The Legendre decomposition is defined vaguely in the paper and supplementary material, and as far as I can tell this is not a broadly known concept. This makes the paper not self-contained, and makes it difficult to ensure correctness of the arguments. \n\nWe agree that Legendre decomposition is not a broadly known concept. Thus, we will make a new section, “Formal definitions” in Appendix and add the following formal definition of Legendre decomposition in the camera-ready version.\n\n__Definition 1__: Legendre Decomposition [29]\n\nFor a normalized non-negative tensor $\\mathcal{P}$ and a set of indices $B$, many-body approximation of $\\mathcal{P}$ is $\\mathcal{Q}$ defined as\n$ \\mathcal{Q} = argmin_{\\mathcal{Q} \\in \\boldsymbol{\\mathcal{B}}} D_{KL} (\\mathcal{P},\\mathcal{Q}) $ \nwhere $\\boldsymbol{\\mathcal{B}} = \\set{ \\mathcal{Q} \\mid \\theta_{i_1,\\dots,i_D} = 0 \\ \\ \\mathrm{ if }\\ \\ (i_1,\\dots,i_D)\\notin B \\text{ for } \\theta\\text{-parameters of } \\mathcal{Q}\\}$.\n\nIn the same way, we will give a formal definition of many-body decomposition as __Definition 2__  in the camera-ready version.\n\n> The optimization method and proof of convexity of the optimization problem are also proven hastily / with lack of good definition/theorem/lemma/proof structure, in particular restate the theorems being used from prior work.\n\nWe appreciate your constructive comments for readability.\n\n- For the optimization method :\nThe optimization method is thoroughly described in Section 2.1.1, Section A, and Algorithm 1. This is the Newton method that optimizes KL divergence from a given distribution. We are happy to address that further if you could provide unclear points specifically.\n\n- For convexity of the proposal:\nWe have discussed in Section A to carefully introduce the necessary definitions to prove the argument, i.e., geodesic, flatness, and projection. We will add the following theorem in Appendix B with proof to further clarify our claim.\n\n__Theorem 1__: The solution of many-body approximation is always unique, and the objective function of many-body approximation is convex.\n\n__Proof__:\nAs we see in Definition 2, the objective function of the proposed many-body approximation is the KL divergence from an empirical distribution (a given normalized tensor) to a subspace $\\boldsymbol{\\mathcal{B}}$. We can immediately prove that a subspace $\\boldsymbol{\\mathcal{B}}$ is $e$-flat from the definition of the flatness [1, Chapter 2] for any $B$. The KL divergence minimization from a given distribution to $e$-flat space is called $m$-projection (see the second paragraph in Section A). The $m$-projection onto $e$-flat subspace is known to be convex and the destination is unique (see the third paragraph in Section A). Thus, the optimal solution of the many-body approximation is always unique, and the objective function is convex. \n\n--\n\nIt is widely known that maximum likelihood estimation of ordinary Boltzmann machines without hidden variables is a convex problem. Since we optimize the KL divergence, the proposed many-body approximation is also maximum likelihood estimation that approximates a non-negative tensor, which is regarded as an empirical distribution, by an extended Boltzmann machine without hidden variables, as described in lines 147-153. As with ordinary Boltzmann machines, the maximum likelihood estimation of extended Boltzmann machines can be globally optimized. __Theorem 1__ is a general proof of this using Information Geometry. We will add this description in Appendix in the camera-ready version.\n\n\n> when discussing results based on linear dependence of parameters in the distribution.\n\nIn the second paragraph titled “Flatness and projections” in Appendix A, the sentence “It is known that a subspace with linear constraints on natural parameters $\\theta$ is $e$-flat” describes the general condition for $e$-flatness. However, to prove the convexity of the proposal, we need only a particular condition “A subspace with some of its natural parameters fixed at 0 is $e$-flat” and it is obvious from the definition of $e$-flatness. To simplify, we replace that sentence in Appendix A with this sentence. In the same way, in line 98, “it is known that a subspace with linear constraints on natural parameters $\\theta$ is flat, called $e$-flat [1, Chapter 2]” will be replaced as “we can introduce a concept of flatness for a set of tensors” to increase the readability.\n"}, {"review_id": "emRVCT120W", "reviewer": "Reviewer_wHyr", "summary": "Authors propose a tensor decomposition approach based on the Legendre decomposition and convex optimization with natural gradient-based algorithm for non-negative tensors which is inspired by many-body interactions in physics. The tensor is interpreted as a probability measure over a multidimensional discrete space given by its indices, and the complexity is regulated by the order and map of interactions between variables (this makes it possible to avoid the difficult problem of choosing the optimal rank of the decomposition). Experimental results demonstrating the good properties in terms of accuracy and computation time of the proposed method in comparison with a number of alternative methods of low-rank approximations.", "questions": "1. [Line 23] \"More recently, tensor networks [4] have been introduced...\". I do not fully agree with this position. Tensor networks have been used for a very long time in the scientific community, including quantum applications (e.g., DMRG), and also the noted decompositions (CP, Tucker) represent a case of the tensor network. It seems to me that this paragraph should be substantially reformulated. It also seems appropriate to note such decompositions as MERA, PEPS, etc.\n\n2. [Line 268] It might be worth refining the positioning of the graphics so that the section title (3.2) isn't at the end of the page.\n\n3. [Line 330] \"error[s] are...\"\n\n4. [Lines 349-350] \"visualizing interactions between modes ...  to see activated interactions between modes\". Perhaps this sentence should be reformulated to avoid repeating the phrase.", "limitations": "The method can only be applied to tensors with non-negative elements. There are also likely to be problems associated with high computational complexity when trying to apply this approach to significantly multidimensional data arrays.", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "The proposed method looks novel and interesting, and it has benefits like a convex optimization formulation and a lack of rank parameters to tune. The work is well structured, written in a good style and easy to read.", "weaknesses": "1. From the presented numerical experiments, the actual advantages of the proposed approach over baselines are not quite clear to me. For the case of a dataset with graphic images, what is the actual quality of the low-rank approximation (it would be useful to give examples of reconstructed images)?\n\n2. As far as I know, there are effective rank-adaptive methods of low-rank approximations (rank-adaptive TT-cross method, rank-adaptive tensor ring, etc.), see, e.g., \"Adaptive rank selection for tensor ring decomposition\" (2021). Will the proposed approach be superior to such modern methods?\n\n3. I did not see estimates of the computational complexity of the proposed method. What dimensions limit the applicability of the approach, i.e., can we use it, for example, for the 100-dim tensors?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "I appreciate your comments. We provide our point-by-point response to each of your comments in the following. \n### In the \"Weaknesses\" section\n\n> 1. From the presented numerical experiments, the actual advantages of the proposed approach over baselines are not quite clear to me.\n\nMany-body approximation has two significant practical advantages over baselines, the global convergence and intuitive model design: \n- Many-body approximation always converges to the globally optimal solution as described in lines 112-114 and Section A. Thus we do not have to repeat algorithms to find better initial values. In contrast, presented results in Section 3-3 for baselines are best results after multiple trials with random initial values (see line 326).\n- Many-body approximation has an advantage in terms of model selection. For example, it does not require hyperparameter tuning as described in Figure 14 and Section E.2, while it still has comparable recovery fit scores compared to baselines. This can be a significant advantage of the proposed method in practice. In traditional low-rank approximation, hyper-parameter tuning based on prior knowledge is difficult because its semantics is unclear. In contrast, the interactions in many-body approximation can be tuned easily by considering the expected association among modes of a given tensor (see more discussion in Sections 3.1 and C.2 in Appendix).\n\n> 1. For the case of a dataset with graphic images, what is the actual quality of the low-rank approximation?\n\nPlease refer to the global response for the reconstruction images. Also, Figure 11 in Appendix shows that the proposal has the better reconstruction quality for the traffic dataset compared to two traditional low-rank approximations, non-negative Tucker decomposition and non-negative tensor-train decomposition.\n\n> 2. Will the proposed approach be superior to rank-adaptive methods?\n\nWe would like to argue that our proposal has two significant advantages over rank-adaptive methods: convexity and intuitive interaction tuning.\n\nRank-adaptive methods typically repeat low-rank approximation so that the reconstruction error becomes smaller than a threshold. This strategy requires a threshold as an additional hyperparameter and more computational cost, as shown in Table III in [a]. In addition, since it is non-convex optimization, each iteration for low-rank approximation cannot find the global optimal due to the initial value dependency, as described in Section IV in [a]. Moreover, interpretation or understanding of the meaning of the obtained rank is difficult as it is not on the original feature space. In contrast, since our proposal is based on convex optimization, there is no initial value dependency; thus, we do not have to run the algorithm multiple times with random initial values to find a better approximation. Furthermore, as described in Sections 3.1, 3.2, and C.1, the interaction is more interpretable than the traditional rank.\n\n[a] Sedighin, Farnaz, Andrzej Cichocki, and Anh-Huy Phan. \"Adaptive rank selection for tensor ring decomposition.\" IEEE Journal of Selected Topics in Signal Processing 15.3 (2021): 454-463.\n\n> 3. I did not see estimates of the computational complexity of the proposed method. \n\nAs described in line 117, our method is based on the natural gradient, and its computational complexity is cubic with respect to the number of parameters to be optimized. For example, the number of parameters for cyclic two-body approximation is described in equation (12). Also, we provide a computational speed comparison in Figure 10 in Appendix. \n\n> 3. What dimensions limit the applicability of the approach, i.e., can we use it, for example, for the 100-dim tensors?\n\nIt depends on not the dimensionality but the number of non-zero natural parameters. Please also refer to the discussion about computational complexity in line 117. For huge-size tensors and much higher-order interactions, we need to consider efficient variations of the natural gradient, which can be our future work.\n\n### In the \"Questions\" section\n\n> 1. Tensor networks have been used for a very long time in the scientific community, including quantum applications (e.g., DMRG), and also the noted decompositions (CP, Tucker) represent a case of the tensor network. It seems to me that this paragraph should be substantially reformulated. It also seems appropriate to note such decompositions as MERA, PEPS, etc.\n\nThank you for pointing it out. We agree that tensor networks have a long history and have been used to evaluate wave functions and partition functions for physics. What we wanted to say in that paragraph is that tensor networks, introduced in physics, have recently been used in the machine learning community as well. To avoid misleads, we will revise that paragraph as follows:  \nMore recently, tensor networks, initially introduced in physics, have become popular in the machine learning community because it helps intuitive and flexible model design, including tensor train decomposition [24], tensor ring decomposition [34], and tensor tree decomposition [2]. Nowadays, traditional tensor structures for physics, such as MERA [b] and PEPS [c], are also used in machine learning.\n\n[b] 10.1088/2632-2153/abffe8 [c] 10.1103/PhysRevB.103.125117\n\nWe will also addressee issues 2, 3, and 4 in the camera-ready version. We appreciate your careful reading. \n\n### In the \"Limitations\" section\n\n> There are also likely to be problems associated with high computational complexity when trying to apply this approach to significantly multidimensional data arrays.\n\nAlthough this is a general problem of tensor factorizations, we would like to clarify that, as you can see in Figure 10 in Appendix, our methods work on tensors whose sizes are ten million scales, 20×20×20×20×20×20 = 64,000,000, with Intel Xeon CPU Gold 5218 and 128GB of memory (See Section E). Moreover, we have provided a way to handle larger tensors in Section B-2 in Appendix.\n"}]

Author rebuttal:
We thank all the reviewers for their insightful feedback and constructive suggestions. We responded to individual reviewers and look forward to further discussions.

Also, to address the concern raised by Reviewer wHyr, we submit a PDF file of reconstructed images of the COIL dataset to compare proposal and low-rank approximations, non-negative Tucker decomposition (NTD) and non-negative tensor train decomposition (NTT).

Meta-review: Overall the reviewers all appreciate the idea and a decomposition but all flagged issues in clarity. Given the reviewers all vote for acceptance and agree this will be of interest to the community, I recommend acceptance. I strongly suggest the authors update the papers according to the discussion with the reviewer and revisit the paper for clarity.
Decision: Accept (poster)