Many-Body Approximation for Tensors
Authors: Kazu Ghalamkari
Mahito Sugiyama
Venue: Submitted to ICLR 2023 (2023)
Keywords: Tensor decomposition
Energy based model
Tensor networks

Abstract: We propose a nonnegative tensor decomposition with focusing on the relationship between the modes of tensors. Traditional decomposition methods assume low-rankness in the representation, resulting in difficulties in global optimization and target rank selection. To address these problems, we present an alternative way to decompose tensors, a many-body approximation for tensors, based on an information geometric formulation. A tensor is treated via an energy-based model, where the tensor and its mode correspond to a probability distribution and a random variable, respectively, and many-body approximation is performed on it by taking the interaction between variables into account. Our model can be globally optimized in polynomial time in terms of the KL divergence minimization, which is empirically faster than low-rank approximations keeping comparable reconstruction error. Furthermore, we visualize interactions between modes as tensor networks and reveal a nontrivial relationship between many-body approximation and low-rank approximation.

Under review as a conference paper at ICLR 2023
MANY-BODY APPROXIMATION FOR NON-NEGATIVE TENSORS
Anonymous authors Paper under double-blind review
ABSTRACT
We propose nonnegative tensor decomposition with focusing the relation- ship between the modes of tensors. Traditional decomposition methods assume low-rankness in the representation, resulting in difficulties in global optimization and target rank selection. To address these problems, we present an alternative way to decompose tensors, a many-body approximation for tensors, based on an information geometric formulation. A tensor is treated via an energy -base model, where the tensor and its mode correspond to probability distribution and ran- dom variable, respectively, and many-body approximation is performed on it by taking the interaction between variables into account. Our model can be globally optimized in polynomial time in terms of the KL divergence minimization, which is empirically faster than low-rank approximations keeping comparable recon- struction error. Furthermore, we visualize interactions between modes as tensor networks and reveal nontrivial relationship between many-body approximation and low-rank approximation.
INTRODUCTION
Tensors are generalization of vectors and matrices. Data in various fields such as neuroscience [rro & Hunyadi 2022), bioinformatics (Luo et al., 2017), signal processing (Cichocki et al. 2015). and computer vision (Panagakis et al., 2021) are often stored in the form of tensors, and features are extracted from them. Tensor decomposition and its non-negative version (Shashua & Hazan 2005) are popular methods that extract features by approximating tensors by the sum of products of smaller tensors. These smaller tensors are often called factors. It usually tries to minimize the difference between the tensor reconstructed from obtained factors and an original tensor, called the reconstruction error.
In most of tensor decomposition approaches, low-rank structure is typically assumed, where a given tensor approximated by linear combination of a small number of bases. Such decomposition requires the following two information. First, it requires the structure, which specifies the type of decomposition such as CP decomposition (Hitchcock) 1927) and Tucker decomposition (Tucker) 1966). In recent years, fensor networks (Cichocki et al. (016) have been introduced, which can intuitively and flexibly design the structure including tensor train decomposition (Oseledets) 2011) tensor ring decomposition (Zhao et al. 2016). and tensor tree decomposition (Murg al. 2010).
Second, it requires the rank value, the number of bases used in the decomposition. Since larger ranks increase the capability of the model while increasing the computational cost, the user required to find the appropriate rank in this tradeoff problem. Since the above tensor decomposition via minimization of the reconstruction error is non-convex, which causes initial value dependence (Kolda & Bader] 2009. Chapter 3), the problem of finding an appropriate setting of the low-rank structure is highly nontrivial in practice as it is hard to locate the cause the decomposition does not perform well. As a result, to find proper structure and rank, the user often needs to perform decomposition multiple times with various settings, which is time and memory consuming.
Instead of the low-rank structure that has been the focus of attention in the past, in this paper, we propose a novel formulation of tensor decomposition, called many-body approximation, that focuses on the relationship among modes of tensors. We determine the structure of decomposition based on the existence of the interactions between modes. The proposed method requires only the decomposition structure naturally determined by the interactions between the modes and does not

---

Under review as a conference paper at ICLR 2023
require the rank value, which traditional decomposition methods also require 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 energy-based model in machine learning, which has been used in Legendre decomposition (Sugiyama et al. 2018 2016) that decomposes tensors via convex optimization. Technically, finds factors of tensor by treating tensor as probability distribution and enforcing some of its natural parameters to be zero. We point out that interactions in the energy function H-)- can be represented using natural parameters of distribution, and we can successfully formulate many-body approximation as special case of Legendre decomposition by setting some of natural parameters to be zero. The advantage of this approach is that many-body approximation can be also achieved by convex optimization that minimizes the Kullback-Leibh (KL. divergence (Kullback & Leibler 195T). Our approach, describing interactions between modes using energy functions, is different from existing methods that focus on interactions between mode matrices Vasilescu Terzopoulos, 2002 Vasilescu 201T) or block tensors (Vasilescu et al. 2021)
Furthermore, we introduce way of representing tensor interactions, which visualizes the presence or absence of interactions between modes. We discuss the correspondence between our representation and the tensor network and point out that operation called coarse- sousse-rained transformation (Levin & Nave 2007). 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 decomposition.
We summarize our contribution as follows:
By focusing on the interaction between modes of tensors, introduce alternative rank- free tensor decomposition, many-body approximation. This decomposition is realized by convex optimization.
We present way of describing tensor many-body approximation, interaction representation, a diagram that shows interactions within a tensor. This diagram can be transformed into tensor networks, which tells us the relationship between many-body approximation and existing low-rank approximation.
We empirically show that many-body approximation is faster than low-rank approximation with competitive reconstruction errors.
2 TENSOR MANY-BODY APPROXIMATION
Our proposal, tensor many-body approximation, based on the formulation of Legendre decomposi- tion for tensors. We first review Legendre decomposition and its optimization in Section 2.1 We introduce interactions between modes and its visual representation to prepare for many-body approxi- mation Section 2.2 Using interactions between modes, we define many-body approximation in Section 2.3 Finally, we transform the interaction representation into a tensor network and point out the connection between many-body approximation and existing low-rank decomposition methods in Section [.4]
In the following discussion, we consider D-order non negaaiiee tensors whose size is (II...., Ip). We assume the sum of all elements in P is for simplicity, while this assumption can be eliminated using the general property of Kullback-Leiih (KI.) divergence, AD (P, P,P, Q) = DRI((P, XQ), for any real number X.
2.1 REMINDER TO LEGENDRE DECOMPOSITION AND ITS OPTIMIZATION
Legendre decomposition is a method to decompose a non-negative tensor by regarding the tensor discrete distribution and representing it with limited number of parameters. We describe non-negative tensor P using natural parameters f = H
(g....))) and its energy function
P Pii....i exp (Hi,,.,,io)) H Hipm..
0
(I)
2

---

Under review as a conference paper at ICLR 2023
where P1..... has role of normalization. Here it is clear that a tensor corresponds to distribution whose sample space is its index set; that is, the value of each element is regarded as the probability of realizing the corresponding index (Sugiyama et al. 2017).
As we can see in equation] 1 we can uniquely identify tensors from natural parameters 0. We can compute the natural parameter from a given tensor
Ip
GooD log Pej.....b
using the Möbius function μ 8:88 (-,, +1}, where S is the of indices, defined inductively as follows:
ia i'2 for all d (I....,D). Zzz=ia P else if is i' for all d {I..... D}.
P
otherwise.
The above modelling for non-negative tensors is an instance of the log-linear model on posets (Sugiyama et al. (017)) Since distribution described by equation I belongs to the exponential family, we can also identify each tensor by expectation parameters 7) i M7, ppp using the Möbius inversion formula as
LD Σ D=iD
I
Io
Th
Pi....b'
P Pin...de
bionnd Th Bi
h==1
See Supplemental Materials for examples of the above calculation. Since distribution is determined by specifying either 6 p-parameters or 7-parameters, they form two coordinate systems called the 8-coordinate system and the cooooddiittte system, respectively. By using the dual flatness, Legendre decomposition achieves convex optimization as shown in the following.
2.1.1 OPTIMIZATION
Legendre decomposition approximates a tensor by setting some 9 values to be zero, which corresponds to dropping some parameters for regularization. Let B be the set of indices of 0 parameters that are not imposed to be 0. Then Legendre decomposition coincides with a projection of a given nonnegative if (6, ( BD)
tensor P onto the subspace B {0 O.......h>
Let us consider projection of a given tensor P onto B. The space probability distributions is not a Euclidean space. Therefore, it is necessary consider geometry of probability distributions, which is studied in information geometry. It is known that a subspace with linear constraints on natural parameters 0 is flat, called e-flat (Amari 2016 Chapter 2). The subspace B is e-flat, meaning that the logarithmic combination, or called e-geodesic, R E {(1 t) log Qi + t log Q2 o(t) 0 < t < 1} of any two points 21, Q2 E B is included in the subspace B, where @(t) a normalizer. There is always unique point P on the e-flat subspace that minimizes the KL divergence from any point P.
P = arg min DKz(P, 2:055
(4)
This projection is called the m-projection. The m-projection onto e-flat subspace is convex optimization. We define two vectors 8 = (fg)se and 7 = (7)6 (ee B- We write as (B] the number of elements in these vectors since it is equal to the cardinality of B. The derivative of the KL divergence and the Hessian matrix G E IR [|| are given as
2 d08 DKL = B, G = 7hnax(1 Ji max(in Tii.....ppp Tj.....D (5)
where 8 and rj are the expectation parameters of Q and P. respectively, and u (ii,..., in) (J.....,ju jD) E B. This matrix is also known as the negative Fisher information matrix. Using update 8 in each iteration t
gradient descent with second-order derivative,
2+1 = 0.6 -G-1 - B
(6)
3

---

Under review as a conference paper at ICLR 2023
P
Minimizing
D (P, P)
P
m- projection
11
2/2
23
4
geodesic
e-flat space
B
b.
d.
mm
mm
H (k.)
H (k.m.p)
21
2/2
23
24
1k
1p
ik
Figure 1: (a) An illustration of optimization of Legendre decomposition. Interaction representations
corresponding to (c) equation 10 and (d) equation I1
The distribution Q++1 is calculated from the updated natural parameters 0 This step finds point Qr+1 that is closer to the destination along with the e-geodesic from Q, to P. We can also calculate the expected value parameters 7/++1 from the distribution. By repeating this process until convergence, we can always find the globally optimal solution satisfying equation [ This procedure illustrated in Figure. [aa. See Supplemental Materials for more detail on the optimization.
2.2 INTERACTION AND ITS REPRESENTATION OF TENSORS
In this subsection, we introduce interactions between modes and its visual representation to prepare for many-body approximation. The following discussion enables intuitively describe relationships between modes and formulate our novel rank-free tensor decomposition.
First introduce n-body parameters, which is generalized concept of one-body and two-body parameters in (Ghalamkari & Sugiyama, 2022). Let n of a n-body parameter be the number of non-one indices; for example, 0,,2,1, one-body parameter, 04,3,1, two-body parameter and 01,2,4.3 is a three-body parameter. We also use the following notation for n-body parameters:
0.h 01,,,, Ii....i
(k..)) .1..... Li Lim. ie
(kom.p) 01..... Pestmi
4
for n = 1, and 3, respectively. We write the energy function H with n-body parameters as m-
D
Hag---i .i Ho +
ΣΣ
(k.m) Hngt
ΣΣΣ
H (k.mp)) Sastmip
+...
H (......
(8)
m=lk=1
where the n-th order energy is introduced as
(9)
For simplicity, we suppose that S 1 < 1/ < ... Ln S D holds. We set Ho = O1.... We say that n-body interaction exists between modes (,,.., Ln there indices ... it, satisfying H(ll....da) +0. O. Aiy....
The first term Hy in equation 8 is called the normalized factor or the partition function. The terms (k) are called bias in machine learning and magnetic field or self-energy in statistical physics. The terms m)) are called the weight of the Boltzmann machine machine learning and two-body interaction electron-electron interaction in physics. To visualize the existence interactions within a tensor, we newly introduce a diagram called interaction representation, which is inspired by factor graphs in graphical modelling (Bishop &
4

---

Under review as a conference paper at ICLR 2023
Nasrabadi 2006 Chapter 8). The graphical representation of the product of tensors is widely known as tensor networks. However, displaying the relations between the modes of a tensor as a factor graph is novel approach. We represent the n-body interaction black square, connected with modes. We describe examples of the two-hody interaction between modes (k; m) and the three-body interaction among modes (k, m,p) in Figure [b). Combining these interactions, the energy function including all two-body interactions is shown in Figure c). and the energy function including all two-body and three-body interactions is shown in Figure (d) for D 4.
This visualization allows us to 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 the following section, we reduce some of n-body interactions, that is, H Hlgu..nt 0, by fixing pll.....a)) = 0 for all indices ( it,)) E (2...., 41 x ... x f
each parameter " foutti
2.3 MANY-BODY APPROXIMATION
Our proposed method, tensor many-body approximation, approximate a given tensor with assuming the existence of dominant interactions between the modes of the tensor and ignoring other interactions. Since this operation can be understood as setting some natural parameters of the distribution to be zero, can be achieved by convex optimization through the theory of Legendre decomposition. As we see below, approximated tensors are represented without the summation symbol Σ This property different from existing low-rank approximations except for rank-1 approximation.
As example, we consider two types of approximations of nonnegative tensor P by tensors represented in Figure [c). (d). If all energies greater than 2nd-order or those than 3rd-order
(((t...do) = 0 for n in equation 8 are ignored, that is, ".....ii
2 or n 3, P is approximated as follows:
Piniiaota P22
Na = 11.22 is $1.24 12 12.34 ig.4
(1.3) (1.4)
(1.4) (2.3) (2.4) (3.4)
(10)
Plsiziiaii P33
Xi1.22 X21 Xig Aig.
(1,2.3) (1.2.4) (1,3.4)
(11)
where each factor on the right-hand side is represented as
X
(k)
exp
1
H!hh
(k)
H (k.m)
H
2
3
( (k)
H (m)
H (p)
(kmmpp)
V cxp
+
km
3
tp
+ H (k))
2 ie
1
+
2
1
H (m.p)
+ H (k.p)
2
+
(K (k.m.p)
The partition function, or the normalization factor, is given Z exp (-01,11.1), which do not depend on indices (i1, i2, ig. 14). When the tensor P is approximated by pmmm the set B contains only indices of n(s mi)-oody parameters.
In the above discussion, we consider many-body approximation with all n-body parameters, while our formulation allows us to use only a part of n-body interactions as shown in the following. We consider the situation where only one-body interaction and two-body interaction between (k, k + 1) exist for all E (I....,D) D} (D + 1 implies 1 for simplicity). Figure 2(a) shows the interaction
representation of the approximated tensor. As can confirm by substituting 0 for (k.) if + k 1, we can describe the approximated tensor as
Pii...i paye Pire.ii
((1)
X
(12)
where
x
i2.f441
2/2
cxp
H. H((+t1) it ie.
H(()) (((++1))
iet1
(13)
with normalization factor = exp ( (-O)....)) When the tensor approximated by paye the set
B contains only all one-body parameters and two-body parameters 11 for d E {1, D) We call this approximation cyclic two-body approximation since the order of indices in equation 12 is
cyclic. We show the connection between cyclic two-body approximation and existing tensor ring decomposition in the following subsection.
5

---

Under review as a conference paper at ICLR 2023
b.
(
i1
(1)
1
14
12
X
X (2)
P
14
X
(4)
X
(2
12
X (3)
23
ceeraatioonReepreeenttttoo
Tensor Network
Tensor Ring Decomposition
Figure 2: (a) Interaction representation of example of cyclic two-body approximation and its transformed tensor network for 1) 4. (b) Tensor network of tensor ring decomposition.
2.4 CONNECTION TO TENSOR NETWORK
Tensor interaction representation is a diagram that focuses on the relationship between modes. Tensor networks, which are well known as diagrams that focus on factors after decomposition, represent tensor undirected graph, whose nodes correspond matrices tensors and edges the modes of summation tensor products (Cichocki et al., (016).
Our tensor interaction representation has tight connection to tensor networks, and we convert tensor interaction representation to a tensor network, For the conversion, we use a hyper-diagonal tensor Q2 that is defined as Sijj d1zōzu.ks where 6;j = 1 if j and O otherwise. The tensor is often represented by e in tensor networks. the community of tensor network, the tensor Ω appears in the CNOT gate and special case Z spider (Nielsen & Chuang 2010). The tensor network in Figure 2 represents the following formula
D πI
ΣV 4
Σ
((d)
Ris 1soosdara)
(14)
lada+s
ja
where jD+1 h1, iD+1 i1, 1D+1 4. Substituting the definition of Ω in equation 14 realize that the tensor network corresponds to equation 12
We point out that the tensor network representation of cyclic two-body approximation is similar to the tensor network of the tensor ring decomposition. The tensor ring decomposition is an extension of the tensor train decomposition, and its representation is shown in Figure (bbb using tensor network. In fact, if we consider the region enclosed by the dotted line in the tensor network as new tensor, the tensor network of the cyclic two-body approximation coincides with the tensor network of the tensor ring decomposition(See (See more details in the Supplemental Materials). This operation, in which multiple tensors are regarded as a new tensor tensor network, called renormalization or coarse- coarse-graining transformation (Evenbly & Vidal 2015).
Comparing the number parameters The number elements an input tensor Ii x I2 x x Ip. After the cyclic two-body approximation, the number of parameters is given
B = 1+ Σ la - - d=1
1)(1+11 -1
(15)
d-1
where we assume Ip+1 = Ij. The first term is for normalizer, the second is the number of one-body parameters, and the final term the number of two-body parameters. In contrast, in the tensor ring decomposition with the target rank (Rı,..., Rp), the number of parameters is given | = R2 I RK+11 The ratio of the number of parameters of these two methods [ proportional to R/ 2 if we assume R₂ R and Ia for all d {1...., D} for simplicity. Therefore, when the target rank is small and the size the imput tensor is large, the proposed method has more parameters than the tensor ring decomposition.

---

Under review as a conference paper at ICLR 2023
b.
A
B
B
$1
9
16
g
3g
Interaction Interaction-Represeniii Representation
Tensor Network
Tensor Tree Decomposition
Figure 3; (a) Interaction representation corresponding to equation 16 and its transformed tensor network for D (b) Tensor network of variant of tensor tree decomposition.
2.4.1 OTHER EXAMPLE OF MANY-BODY APPROXIMATION AND ITS TENSOR NETWORK
In the same way, we can find correspondence between another example of many-body approximation and the existing low-rank approximation. For D 9, we consider three-body and two-body interactions among (i1, i2; is), (i4, 4s, is), and (i7, is, ig) and three-body approximation among (ig; ie; ig). We provide the interaction representation of the target energy function Figure 3a). In this approximation, the decomposed tensor can be described as
Pii..... = Ai is Ciy. 9i S9.99
(16)
In the same way in the case of the cyclic two-body approximation, we can convert the interaction representation to tensor network, as described in Figure 3a). A tensor network of tensor tree decomposition in Figure b) emerges when the region enclosed by the dotted line is replaced with a new tensor (shown with tilde) in Figure 3 a). Such tensor tree decomposition is used in generative
modeling (Cheng et al. 2019), computational chemistry (Murg et al. 2015) and quantum many-body
physics (Shi et al. 2006).
As we have seen above, by transforming tensor interaction representation to tensor networks and applying coarse- goanss-griinng.... we reveal the relationship between tensor many-body approximations and low-rank approximations.
2.5 MANY-BODY APPROXIMATION AS GENERALIZATION OF MEAN-FIELD APPROXIMATION
has been already pointed out that any tensor P can be represented by vectors (d) E R for d E {L...., D} as P Pip
(D) Fin if and only if all n(> 2)-body 0-parameters are 0 (Ghalamkari & Sugiyama, 2021. The right-hand side is equal to the Kronecker product of D the rank-1 approximation
vectors (x(1) ,..., (D) and therefore this approximation is equivalent since the rank of the tensor that can be represented by the Kronecker product is always 1, which 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.
is also known to correspond
2.6 COMPUTATIONAL COMPLEXITY
We analyze the computational complexity of many-body approximation. In many-body approxima- tion, the overall complexity is dominated by the update of 8, which includes matrix inversion of G. The complexity of computing the inverse of an n x n matrix is O(n²). Therefore. the computational complexity of many-boot approximation is ()) where y is the number of iterations.
This complexity can be reduced if we reshape tensors so that the size of each mode becomes small. For example, let consider 3-order tensor whose size is (J2,22 J2, J2) and its cyclic two-
7

---

Under review as a conference paper at ICLR 2023
a.
330,30,30,00,30
R15
b.
20, R10
C.
(30,30,30,30,30
d.
(20,20,20,20,00220)
0.020
0.50008
Relative Error
0.0100
Relative Error
Relative Emor
0.50006
Relative Error
0.50020
0.015
0.50015
0.0075
0.50004
0.50010
0.0050
0.010
0.50002
0.50005
0.0025
0.005
0.50000
0.50000
10a
10*
10a
10
10a
10
10a
10
parameters
# parameters
8 parameters
# parameters
Running time (sec)
10a
Bunning time (sec)
Running time (sec)
10a
10a
10a
Running time
10a
10
10a
10a
10
10
10
10
10
10
10
10a
10
parameters
parameters
parameters
parameters
APG
MM
HALS
IraMM
proposal
Figure 4: (a)(b) Results for low ring rank tensor. (c)(d) Results for tensors sampled from uniform distribution. The vertical red dotted line is [B] (See equation 15
body approximation. In this case, the time complexity is ((yJ since it holds that [B] J (See equation [1).. In contrast, if we reshape the input tensor to a 6-order tensor whose size is JJ, J, J), the time complexity becomes O(7:9) since it holds that [B] J
This technique of reshaping a tensor into a larger-order tensor is used practically not only in the proposed method but also in various methods based on tensor networks, such tensor ring decompo- sition (Malik& & Becker 2021
3 EXPERIMENTS
As seen in Section 2.4 many-body approximation has a close connection to low-rank approxima- tion. For example, tensor ring decomposition, if we impose that decomposed factors can be represented as products with hyper- hyper-diagonal tensors Q. this decomposition is equivalent to a cyclic two-body approximation (see Figure 2. Therefore, to examine our conjecture that cyclic two-body approximation is as capable of approximating as tensor ring decomposition, we empirically examine the efficiency and effectiveness of cyclic two-body approximation compared with tensor ring de- composition. As baselines, NTR-APG, NAR-HALS, NTR-MU, NTRMMM and NTR-IraMM (Yu et al., [021]] (022). These methods minimize the reconstruction error defined with the Frobenius norm by the gradient method. See Supplemental Materials for implementation detail.
five existing methods of non-negative tensor ring decomposition.
We evaluate the approximation performance by the relative error - ||| for input tensor T and a reconstructed tensor T with the Frobenius norms || ||... Since all the existing methods are based on nonconvex optimization, we plot the best score (minimum relative error) among restarts with random initialization. In contrast, the score of our method is obtained by single run it is convex optimization and such restarts are fundamentally unnecessary. We compare the total running time of them.
Synthetic data We performed experiments on four synthetic datasets. The first two are synthetic data with low tensor ring rank. This setting is often used in evaluation tensor ring decomposition. We create D core tensors of size Rt x I R by sampling from uniform distribution. Then a tensor with the size 1 and the tensor ring rank (R.... R) is obtained by the product of these D tensors. Results for R 15, D 5, I 30 shown in Figure [a... and those for R 10, D 6. I = 20 in Figure (b). Relative error and computation time are plotted with gradually increasing the target rank of the tensor ring decomposition, which is compared to the score of our method, plotted as the
8

---

Under review as a conference paper at ICLR 2023
4DLFD
TT_ TT_ChartRes
TT. Origami
TT. Paint
0.242
0.32
Relative Error
Relative Error
0.3200
Relative Emror
0.33
Relative Error
0.32
0.236
0.3150
0.30
0.31
0.3100
0.30
0.29
0.230
400
600
1200
1800
500
Serueu * parameters
1000
2000
para meters
parameters
# parameters
10a
Running time (sec)
Running (sec)
Running time (sec)
10
10a
10a
104
Running
10a
10a
10
300
400
500
600
1200
1800
500
1000
1500 2000
500
1000 1500 2000 parameters
parameters
meram
# parameters
APG
MU
MM
-
HALS
IraMM
proposal
Figure 5: Experimental results for real datasets. The vertical red dotted line is |B] (See equation 15p
cross point of horizontal and vertical red dotted lines. Please note that our method does not have the concept of the rank, thus the score of our method invariant to changes of the target rank unlike other methods. If the cross point of red dotted lines is lower than other lines, the proposed method is better than other methods.
In addition to the above case in which assumed the low-rankness, also generated synthetic datasets without such an assumption. We created a tensor of size 30 and tensor of size 20 by sampling from uniform distribution and performed the same experiment. Results are shown in Figure cc) and Figure d). In all experiments, the proposed method is superior to comparison partners both efficiency and effectiveness. should be noted that the relative error of the proposed method is smaller even when the target rank of the tensor ring decomposition is large and the number of parameters is several times larger than the proposed method.
Real data Next, we evaluate our method on real data. produced from 4D Light Field Dataset (Honauer al., 2016 1996). STT_ChartRees TT Or irmmimi and TT Pa Painn nt are 7-order tensors, which is produced from TocyoTech Hyperspectral Image Dataset (Monno al. 2015 2017). Each tensor has been reshaped to reduce the computational complexity. See the dataset details in the Supplemental Materials. The proposed method is always faster than baselines with keeping the competitive relative errors. In baseline methods, slight change of the target rank induce significant increase of the reconstruction error due to the nonconvex nature of them.
4DLFD is a 9-order tensor, which is Gorder et al. 1996 Levoy & Hanrahan
CONCLUSION
We propose many-body approximation for tensors, which decomposes tensors with focusing on the relationship between modes represented by an energy-based model. It approximates tensors by ignoring the energy corresponding some interactions, which can be viewed as generalization of mean-field approximation that considers only one-body 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. Furthermore, introduce way of visualize interactions between modes, called interaction representation, to see activated interactions between modes. We have established transformation between representation and tensor networks, which reveals the nontrivial connection between many-body approximation and the classical tensor low-rank tensor decomposition.

---

Under review as a conference paper at ICLR 2023
REPRODUCIBILITY STATEMENT
We provide implementation and dataset details in Supplemental Material. We provide both of code for the proposed method and code for comparison methods in Supplemental Material.
ETHICS STATEMENT
This study theoretical analysis tensors and we believe that our theoretical discussion would not have any negative societal impacts.
REFERENCES
S. Amari. Information Geometry and Its Applications. Springer, 2016.
Christopher Bishop and Nasser Nasrabadi. Puttern recognition and machine learning, volume Springer, 2006.
Song Cheng, Lei Wang. Tao Xiang, and Pan Zhang. Tree tensor networks for generative modeling. Physical Review B, 191151115111, 2019.
Andrzej Cichocki, Danilo Mandic, Lieven De Lathauwer, Guoxu Zhou, Qibin Zhao, Cesar Caiafa, and Huy Anh Phan. Tensor decompositions for signal processing applications: From two-way to multiway component analysis. IEEE signal processing magazine, 32(2):145-163, 2015.
Andrzej Cichocki, Namgil Lee, Ivan Oseledets, Anh-Huy Phan, Qibin Zhao, Danilo Mandic, et al. Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions. Foundations and Trends@ in Machine Learning. 94-55):249-44 429, 2016.
Aybüke Erol and Borbála Hunyadi. Tensors for neuroimaging: A review applications of tensors to unravel the mysteries of the brain. Tensors for Data Processing, pp. 427-482, 2022.
Glen Evenbly and Guifre Vidal. Tensor network renormalization. Physical review letters, 115(18): 180405, 2015.
K. Ghalamkari and M. Sugiyama. Fast tucker rank reduction for non-negative tensors using mean- field approximation. In Advances in Neural Information Processing Systems, volume 34, pp. 443 454, Virtual Event, December 2021.
K. Ghalamkari and M. Sugiyama. Fast rank-1 NMF for missing data with KL divergence. In Proceedings of the 25th International Conference on Artificial Intelligence and Statistics, pp. 2927-2940, Virtual Event, March 2022.
Steven Gortler, Radek Grzeszezuk, Richard Szeliski, and Michael Cohen. The lumigraph. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pp- 13-54, 1996.
Frank L Hitchcock. The expression of a tensor or polyadic as sum of products. Journal of Mathematics and Physics, 6(1-4):164-189, 1927.
Katrin Honauer, Ole Johannsen, Daniel Kondermann, and Bastian Goldluccke. A dataset and evaluation methodology for depth estimation on 4d light fields. In Asian Conference Computer Vision, pp. 19-34. Springer, 2016.
Tamara G Kolda and Brett W Bader. Tensor decompositions and applications. SIAM review, 51(3): 455-500, 2009.
Solomon Kullback and Richard A Leibler. On information and sufficiency. The annals of mathemati- cal statistics, 22(1):79-86, 1951.
Michael Levin and Cody P Nave. Tensor renormalization group approach to two-dimensional classical lattice models. Physical review letters, 10010, 2007.
10

---

Under review as a conference paper at ICLR 2023
Mare Levoy and Pat Hanrahan. Light field rendering. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques, pp. 31 42, 1996.
Yuan Luo, Fei Wang, and Peter Szolovits. Tensor factorization toward precision medicine. Briefings in bioinformatics, 18(3):511-55) 2017.
Osman Asif Malik and Stephen Becker. A sampling- based method for tensor ring decomposition. In International Conference on Machine Learning, pp- 7400-7411. PMLR, 2021.
Yusuke Monno, Daisuke Kiku, Masayuki Tanaka, and Masatoshi Okutomi. Adaptive residual interpolation for color and multispectral image demosaicking. Sensors, 17(12):2787, 2017.
Yusukex Monno, Sunao Kikuchi, Masayuki Tanaka, and Masatoshi Okutomi. A practical one-shot multispectral imaging system using a single image sensor. IEEE Transactions on Image Processing, 24(10):3048-33) 3059, 2015.
Valentin Murg, Frank Verstracte, Örs Legeza, and Reinhard M Noack. Simulating strongly correlated quantum systems with tree tensor networks. Physical Review B, 82(20):205105, 2010.
Valentin Murg, Frank Verstraete, Reinhold Schneider, Peter R Nagy, and Legeza. Tree tensor network state with variable tensor order: An efficient multireference method for strongly correlated systems. Journal of Chemical Theory and Computation, 1 1133):027--0036, 2015.
Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quantum Information: 10th Anniversary Edition. Cambridge University Press, 2010. doi: 111111//C809885119776666
Ivan V Oseledets. Tensor-train decomposition. SIAM Journal on Scientific Computing, 33(5):
2295-2317, 2011.
Yannis Panagakis, Jean Kossaifi, Grigorios Chrysos, James Oldfield, Mihalis Nicolaou, Anima Anandkumar, and Stefanos Zafeiriou. Tensor methods in computer vision and deep learning. Proceedings of the IEEE, 109-5),863-890, 2021.
Amnon Shashua and Tamir Hazan. Non-negative tensor factorization with applications to statistics and computer vision. In Proceedings of the 22nd international conference Machine learning, pp. 792-799, 2005.
Y-Y Shi, L-M Duan, and Guifre Vidal. Classical simulation of quantum many-body systems with a tree tensor network. Physical review 74(2):022320, 2006.
M. Sugiyama. H. Nakahara, and K. Tsuda. Information decomposition structured space. In 2016 IEEE International Symposium on Information Theory (ISIT), pp. 575-579, Barcelona, Spain, July 2016.
M. Sugiyama, H. Nakahara, and K. Tsuda. Tensor balancing statistical manifold. In Proceedings of the 34th International Conference on Machine Learning (ICML), volume 70 of Proceedings of Machine Learning Research, pp. 2270-3279, Sydney, Australia, August 2017.
M. Sugiyama, H. Nakahara, and K. Tsuda. Legendre decomposition for tensors. In Advances in Neural Information Processing Systems 31, pp. 8825 8825-8835 8835,, Montréal, Canada, December 2018.
Ledyard R Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3): 279-311, 1966.
M. A. O. Vasilescu. Multilinear projection for face recognition via canonical decomposition. In Proc. IEEE Inter. Conf. on Automatic Face Gesture Recognition (FG 2011), pp. 476-483, Mar 2011. doi: 10.1 10..///..111.57711111
M. A. O. Vasileseu and D. Terzopoulos. Multilinear analysis of image ensembles: TensorFaces. In Proc. European Conf. on Computer Vision (ECCV 2002), pp. 447-460, Copenhagen, Denmark, May 2002.
M Alex O Vasilescu, Eric Kim, and Xiao S Zeng. Causalx: Causal explanations and block multilinear factor analysis. In 2020 25th International Conference on Pattern Recognition (ICPR), pp. 10736- 10743. IEFE, 2021.
11

---

Under review as a conference paper at ICLR 2023
YuYuan Yu, Kan Xie, Jinshi Yu, Qi Jiang, and Shengli Xie. Fast nonnegative tensor ring decomposition based on the modulus method and low-rank approximation. Science China Technological Sciences, 64(9):1843-1833, 2021.
Yuyuan Yu, Guoxu Zhou, Ning Zheng, Yuning Qiu, Shengli Xie, and Qibin Zhao. Graph-regularized non-negative tensor-ring decomposition for multiway representation learning, IEEE Transactions on Cybernetics, 2022.
Qibin Zhao, Guoxu Zhou, Shengli Xie, Liqing Zhang, and Andrzej Cichocki. Tensor ring decomposi- tion. arXiv preprint wrXiv:1606.0055000000 2016.
A CYCLIC TWO-BODY APPROXIMATION AND RING DECOMPOSITION
We can interpret cyclic two-body approximation as tensor ring decomposition with constraints as de- scribed below. Non-negative tensor ring decomposition approximates given tensor P E IR 12
1xxxxx1g xID
with D core tensors x( x())x
(D) with
(())
R220
R Ra for each d {1.2. ..... D)
R
R₂
RD
Pe Phi...i P P
Σ
Σ
Xrp.h 2r1.
(D) Xrp 1.DDD
(18)
1==1r==1
rn=1
where (RI..... RD) is called tensor ring rank. The decomposition is described in Figure 2 c). The eyclic two-body approximation also approximates the tensor P in the form of equation 18 imposing
an additional constraint that each core tensor (d) is decomposed as
Ia
Xra aafara
(d)
Σ
x
ma Amaiiars
(19)
m,=1
for each d E f1. D}. where Piyk 8eydguOka. We assume To = TD for simplicity. We obtain equation] 12 by substituting equation 19 into equation [8]
This constraint enables us to perform convex optimization, From Kronecker's ' delta 8, ra = is holds in equation 19 thus X (d)) a tensor with the size Ta-1 x Ta x la- Tensor ring rank after the cyclic two-hody approximation is (J,,...: Tp) since the size of core tensors coincides with tensor ring rank.
B IMPLEMENTATION DETAIL
We describe the implementation details of methods in the following.
Proposed method Our method is implemented in Julia 1.8. We use natural gradient method for cyclic two-body approximation. The natural gradient method uses the inverse of the Fisher infor- mation matrix to perform second-order optimization in non-euclidean space. For non-normalized tensors, we conduct the following procedure. First, we compute the total sum of elements of an input tensor. Then, normalize the tensor. After that, we conduct Legendre decomposition for the normalized tensor. Finally, we get the product of the result of the previous step and the total sum we compute initially. The termination criterion the same as the original implementation of Legendre Decomposition by /Sugiyama et al. (2018), that is, it terminates if ||m! 7 || < 10-5, where B the expectation parameters on f-th step and r) is the expectation parameters of an input tensor, which are defined in Section 1. The overall procedure is described in Algorithm [ Note that this algorithm is based on Legendre decomposition by /Sugyyama al. (2018).
Baseline methods We implemented baseline methods by translating MATLAB code provided by the authors into Julia code for fair comparison. As we can see from their original papers, -TR-APG. NTR-HALS, NTR-MU, NTR-MM and NTR-IraMM have an inner and outer loop to find local solution. We repeat the inner loop 100 times. We stop the outer loop when the difference between the relative error of the previous and the current iteration is less than 10e-4. NTR-MM and NTR-IraMM require diagonal parameters matrix E. We define E = wI where is an identical matrix and + = 0.1. The NTR-IraMM method performs low-rank approximation the matrix obtained by mode expansion of an input tensor. The target rank is 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).
12

---

Under review as a conference paper at ICLR 2023
Algorithm 1: Many-body approximation
MANYBODYAPROXIMATTON
B)
f- Total sum of
Obtain normalized input tensor P T/8
|| "1"""" denotes element-wise division
Compute 7 of P using equation
Initialize 0
e.g. 0, for all E
t 1+1
repeat
Compute 2g using the current parameter 0.7h with equation 1
Compute Tt. from Q, using equation 3
Compute the inverse of the Fisher information matrix using equation 5
0
0
G )
until |
B - f <E
// We set 10 in our implementation;
T Q: 8
|/ " denotes element-win multiplication
*
return T
Environment Experiments were conducted on Ubuntu 20.04.1 with a single core of 2.I Intel Xeon CPU Gold 5218 and 128GB memory.
C DATASET DETAIL
We describe the details of each dataset in the following.
Synthetic Datasets For all experiments synthetic datasets, we change the target ring-rank as 2....... for baseline methods.
(r.....r) for
Real Datasets 4DLFD is originally a (9,9, 9, 512, 512,3) 3) tensor, which is produced from 4D Light Field Dataset described in Honauer al. (2016). Its license is Creative Commons Attribution- NonCommercial-Shaee
4.0 International License. We use dino images and their depth and disparity map in training scenes. We concatenate them to produce tensor. We re- shaped the tensor as (6, 8, 6, 8, 6, 8, 6, 12). For baseline methods, we chose the target ring- 2, 2, 5, 2, (2, 5, 2, 2,2,2
rank as (2, 3, 2, 2, 2, 2, 2, 2, (2, 2,3, 2, 2). (2,2, 2, 2,, TT TT_ VT_ChartRes is origi-
(2,
(2, 2, 2, 2, 2, 2, 2, 2, 7). (2, 2, nally a (736, 736, 31) tensor, which is produced from TokyoTech 31-band Hyperspectral Image Dataset. We use ChartRes mat. We reshaped the tensor as (23, 8, 23, 4, For baseline methods, we chose the target ring-rank as 2, 2,
(2,2,2, 5), 2,,
2, (3, 23,2, 2,2,, 2,2,)) (3, 2, 2, 3, 2, 2, 10), 23 2, 2, 6). (4,,,2,2,2.2.6), VT_Origami and TT_Paint TokyoTech 59-band Hyperspectral Image Dataset. We use Origami. mat and Paint mat. In [TT_Origami. 0.0016% of elements were negative, hence we preprocessed all elements of TT rT_Origami by subtracting -0.000764, the smallest value in TT TT_Origami: to make all elements non-negative. We reshaped the tensor as (8, 8,
2,,,2,32,,,33
2,, 2, 2, 12), 2,,2,2,3, 15), (3, (3,2, 2, 2, 16). originally (512, 512, 59) tensors, which produced from
8, 59). For baseline methods, we chose the target ring-rank as (2, 2, 2, 2, 2, r) for r 2, 3, ... 15. These reshaping reduces the computational complexity as described Section 2.6 to complete the experiments in a reasonable time.
D PROJECTION THEORY IN INFORMATION GEOMETRY
We explain concepts of information geometry used in this study, including natural parameters, expectation parameters, model flatness, and convexity optimization. In the following discussion, consider only discrete probability distributions.
13

---

Under review as a conference paper at ICLR 2023
(8, n) coorrdinat and geodesics In this study, we map a normalized D-order negative tensor x1x to discrete probability distribution with D random variables. Let lI be the set of discrete probability distributions with D random variables. The entire space U 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 two points is a straight line. In non- Euclidean space, such a 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, P2 U, e- and m-geodesics be defined as
P E R 22
{ Re log R₂ = (1 t) log Pi + tlo P }. { R. = (1 - t)P1 + tP2 }
respectively, where 0 S S 1 and (tt) a normalization factor to keep R, to be distribution.
We can parameterize distributions p E U by a parameter called the natural parameter. We have described the relationship between distribution and natural parameter 0 (01.......... ............. JI) in equation [] The natural parameter 0 serves as a coordinate system U4 since any distribution in 2 specified by determining 0. Furthermore, we can also specify a distribution by its expectation parameter 17 (f)............................. ..... ........ ........))) which corresponds to expected values of the distribution and an alternative coordinate system of U. The definition of the expectation parameter rj is described in equation 3 (I-cordiinates and n-coordinates are orthogonal with each other, which means that the Fisher information matrix G has the following property, Ga = an /00, and (G (G-1) (G-")a, = 00,/ /0 and m- geodesics also be described using these parameters as follows.
{ O' 0 (1 1)0 1 + 10 }, {n' 17 (1 t)n² 2 in P3
where P and 17" are f- and cocorddnatee of distribution P E U.
Flattness and projections A subspace is called e-flat when any e-geodesic connecting two points in a subspace is included in the subspace. The vertical descent of an m m-geodesic from a point P U onto e-flat subspace B, is called m-projection. Similarly, e-projection is obtained when we replace all e with yn and m with e. The flatness of subspaces guarantees the uniqueness of the projection destination. The projection destination P or P obtained by m- or e-projection onto B, or Bm minimizes the following KL divergence.
P arg min (PI,P, 2), P = arg min DKI(O, QEB.
QEBm.
The KL divergence from discrete distributions P E to Q E u is given as
Ip
(P,
Σ 11=
P. Pp..... log Ap=1
Σ
(20)
Ci....op
It is known that a subspace with linear constraints on natural parameters 0 is e-flat (Amari 2016 Chapter 2). The proposed many-body approximation performs m-projection onto the subspace B C. 21 with some natural parameters fixed to be O. From this lincar constraint, e-flat. Therefore, the optimal solution of the many-body approximation is always unique. When a space is e-flat and m-flat at the same time, say that the space is dually-flat. ti is dually-flat.
know that B is
Natural gradient method e(m)-flatnes: guarantees that cost functions to be optimized equa- tion 20 are convex. Therefore, (n(e)-projection onto e(m))-lat subspace can be implemented by gradient method using a second-order gradient. We call this gradient method the natural gradient method. The Fisher information matrix G appears by second-order differentiation of the KI. diver- gence (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 R4xx--Io with 0 and v using the Möbius function, defined in Section 2.1 We provide examples here. In equation 2 The Möbius function is used to find the natural parameter 0 from distribution P. For example, if D 2, 3, it holds that
0t.,62 log Pi,,.2 log P P,, -1,2 logg P1, log P,,-1,b2-1,
0, = log Pii, - log P,, P - log P P - loo P21
+ log P, -1,21,22-1,8 1,40
log P - log P, -2 log Ps 11,6 -1,
14

---

Under review as a conference paper at ICLR 2023
where we assume Pa,t P,,0 and P11.2.00 01.00 Po.12.53 = Note that, to identify
the value of fes....ke. we need only Piç Piç.....% with (f1.....*k) k2) E {i - 1, x {12 - 1, X
(ia - 1, (d)- In the same way, using equation 3 we can find the distribution P by the expectaion parameter v7. For example, if D 2, holds that
Pi 111,62 1,62 01 + t1igh,
Pindaiis = 111i221 - +.i2 - 121. - d2:
+ +11+1,22 +1.+1 11ist+ +hs+
where we assume +1,42 77,,1++1 = 0 and +1.64 1,, = 1h, 0.
15


Reviews:
[{"review_id": "j64XUirAJG", "reviewer": "Reviewer_QF9P", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors' starting point is the \"Legendre Decomposition for Tensors\" paper. The Legendre decomposition is modified to approximate a data tensor by \"dividing\" the original data tensor N-ways and approximating the N-blocks which \"interact\" with each other through a subset of the tensor modes. The authors assume that tensor blocks have dominant interactions while others can be deemed to be irrelevant and set to zero. The authors refer to this modified approximation as a many-body approximation.  The approach was tested on simulation data and 4D light field data etc. The authors report on some of the experiments lower reconstruction errors relative to prior research.\n\n\n", "summary_of_the_review": "The paper introduces a modified Legendre decomposition that subdivides a data tensor into multiple tensor blocks to approximate a tensor whole.  It also introduces an interaction tensor network diagram for better visualization.  However, the paper is missing an algorithm, a concrete example, some important discussions, and comparisons to prior work.\n\nI look forward to reading the final version of this paper when all outstanding issues discussed above are fixed.\n\n______\n______\n______\n______\n______\n---  UPDATED REVIEW  given the authors' response:  ---\n---\nThe paper's title, text, fig 3, and experiment description originally led me to believe that a tensor was approximated by multiple bodies  that interact through a subset of modes, such as the parts that make up a body. I assumed that the division into multiple bodies was achieved by setting some of the natural parameters of a modified Legendre decomposition to zero which was then followed by factorization of each body with a Legendre decomposition. The authors informed me that this was not the case.  \n\nThe paper does state that it approximates a data tensor by formulating a \"many-body approximation as a special case of Legendre decomposition by setting some of the natural parameters to be zero\" (pg.2, paragraph 2).  However, the multi-body language is unrelated to the multi-bodies/particles one encounters in physics and unrelated to the multi-bodies/parts that make up an object.  Apparently, the \"bodies\" and \"interactions\" are artifacts of a sequential implementation that do not have a  physical meaning and do not hold true if one employs a parallel implementation.  Despite repeated requests, the authors have refused to define the meaning of a 2-body or N-body.  Why?\n\nThe name \"interaction diagram\" is a misnomer. A more appropriate name is an operation order diagram.   All the mode interactions are active in a Tucker decomposition which is self-evident in an M-mode SVD (parallel Tucker), but it is not immediately apparent if one looks at the operation order of a tensor ring (sequential Tucker). In Fig. 2, \"the mode interaction diagram\" for a Tucker decomposition depicts the operation order of a sequential implementation, ie a tensor ring. \n\nThe algorithm in the supplement is the natural gradient descent algorithm that does not integrate equations 1,3,5.  The equations are in the main body of the paper and contain undefined variables that would have had to be specified if they were part of the algorithm.  It would have been nice if the authors would have explained the meaning of the superscripts and subscripts which can vary from author to author.\n\nThe experiments are inconclusive. The authors have not solved any AI/ML/CV/CG problems. Hence, it is unclear if their approach is useful. The authors have asserted that their approximated images are still usable, but have not displayed any synthesized or reconstructed images.  Why?\n\n\nConclusion:\n---\nThe authors' claim that they are performing a multi-body approximation is unsubstantiated.  The terms \"multi-bodies\" and \"interactions\"  are unrelated to the multi-bodies/particles in physics and these terms seem to be artifacts of a sequential optimization.  Beyond the inappropriate nomenclature that is misleading, the paper talks in generalities without being grounded in a concrete AI/ML/CV/CG example and without providing any intuition.  The algorithm is missing many details. \n\nCurrently, this paper is not ready for publication at ICLR.\n\n\nOther comments:\n---\n1. The title is misleading and a type of flag planting.  The authors are not performing a multibody approximation where the bodies are physics particles or the parts of an object.  Hence, a more appropriate title might be \"Legendre Tensor Decomposition envisioned as an M-body Approximation,\" and even this is a stretch.\n\n2. > Our approach, describing interactions between modes using energy functions, is different from existing methods that focus on interactions between mode matrices (Vasilescu & Terzopoulos 2002; Vasilescu 2011) or block tensors (Vasilescu et al. 2007)\n\nThe above is inaccurate.  Object parts may not interact with each other at all, or may interact with each other through a subset of modes.  The following is a more accurate description:\n\nTensor factor analysis was introduced in computer vision with the TensorFaces papers (Vasilescu and Terzopoulos 2002; Vasilescu 2011) and in computer graphics with the TensorTexture paper (Vasilescu and Terzopoulos 2004). Vasilescu etal employ energy functions and focus on interactions between the many bodies that make up a tensor whole, such as the parts that make up an object (Vasilescu etal 2021; Vasilescu ICPR 2022). \n\nBy comparison, we define multi-body as  ...  (please provide a definition).\n\n___\n___\n\nGiven the other glowing review ratings, it is unlikely that there will be a reply forthcoming to any of my inquiries. Hence, I have gone ahead and just downgraded the paper. Best of luck and congratulations on your ICLR 2023 paper.", "strength_and_weaknesses": "Strengths:\n1. This modified non-negative tensor approximation subdivides a data tensor into multiple tensor blocks that approximate the tensor whole.  The tensor blocks have modes that are entirely their own, and modes that they share with other tensor blocks.\n\n2. Mode interaction is visualized through a new network diagram.\n\nWeakness:\n1. There is no algorithmic description.  From Appendix B, I see that a tensor is reshaped and split into several tensor blocks each approximated by a tensor ring. Presumably, some of the irrelevant mode interactions between the tensor rings (tensor blocks) are set to zero.  Which mode interactions were set to zero? Would any reshaping and tensor block subdivision work equally well?\n\n2. The authors have not clearly articulated \n\n+ (i) the criteria for splitting a data tensor into $N$ data blocks,\n+ (ii) the criteria for determining which mode interactions are deemed to be dominant or irrelevant and set to zero, or \n+ (iii) the choice of tensor ring ranks.\n \n4. It is unclear what the experimental results are demonstrating.  What is the AI, ML, CV, or CG problem being solved?   The goal is to find meaningful representations that increase recognition rates or result in high-quality reconstructed or synthesized images, etc.  Reconstruction errors are not always good proxy indicators for the previously mentioned tasks. For CG, the authors may wish to display reconstructed or newly synthesized images versus reconstruction errors versus decomposition speed versus reconstruction speed. \n\nIf the goal is compression, this should be stated.  The authors may wish to compare their work with jpeg or other state-of-the-art compression techniques.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "_vT54uid-YL", "reviewer": "Reviewer_GSHo", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes a new kind of tensor decomposition (TD) framework for non-negative tensors which is inspired by many-body interactions in physics. Unlike most traditional TDs which are hard to fit since they correspond to non-convex optimization problems, the proposed framework corresponds to a *convex* optimization problem. The authors propose a gradient descent-based scheme for solving it.\nAnother notable feature of the proposed framework is that it *does not* have a notion of rank---a parameter which is difficult to choose in more standard TDs. \n\n### Update after rebuttal ###\n\nThe authors have done a good job of addressing all the issues I brought up in my original review. The optimization scheme makes sense now. What the various n-body constraints imply in terms on constraints on the $\\theta$ parameters has also been clarified. The connections to other decompositions is clearer. The explanation of physics terms has also been improved which makes the paper more readable for non-physicists. \n\nWith these improvements, I think the paper is good enough for publishing in ICLR. I have increased my score from 3 to 8.\n\nA few minor things I noticed when reading through the paper again after the rebuttal that the authors can fix in the camera ready version as appropriate:\n- When mentioning natural parameters and dual flatness for the first time, it would be great if you provided a reference to the supplement where these things are discussed more in depth.\n- In Eq (8), the order of the two summations should be changed, right? The first one should be $\\sum_{k=1}^D \\sum_{m=1}^{k-1}$ instead of what it is now, and the second one should be $\\sum_{k=1}^D \\sum_{m=1}^{k-1} \\sum_{p=1}^{m-1}$.\n- In sentence above Eq (10), there's a word missing (indicated in bold): \"If all energies greater than 2nd-order or those *greater* than 3rd-order...\". \n- In the various places where you refer to the supplementary material, it would be great if you indicated which section in the supplementary material.\n- Section \"Synthetic data\" at the bottom of page 8, there's a stray period that should be removed. (\"...are shown in Figure 4 (a)*.*, and those...\")\n- Section \"Real data\" on page 9: You should refer to Fig 5 at some point, since it contains the results for this experiment.\n- On page 14, in one of the paragraph headings: \"Flattness\" -> \"Flatness\"", "summary_of_the_review": "I think the proposal itself is novel and interesting. It has many benefits like a convex optimization formulation and a lack of rank parameters to tune. Additionally, the experiment results (speed-ups in particular) are promising. These considerable strengths are unfortunately brought down by lack of clarity in the text. For example, the gradient descent optimization scheme is hard to follow, the connection between various parameters is not explained properly, and the connection between the presented framework and more traditional decompositions is not explained well. Also, the use of various physics terms that are non-standard in ML and numerical linear algebra add to the confusion. \n\nI think once these clarity issues are sorted out, this will be a strong and interesting paper. I encourage the authors to work on improving the clarity, and I hope that my comments will be helpful for doing that.", "strength_and_weaknesses": "**--- Strengths ---**\n\n**S1.** The novelty of the proposal is a strong strength. I haven't seen anything like this proposal before. I think Sec 2.5 explains the novelty in a nice way: The framework does rank-1 decompositions, where the notion of \"rank-1\" is now generalized beyond the traditional tensor rank.\n\n**S2.** The fact that the decomposition corresponds to a *convex* problem is another considerable strength, since most other tensor decompositions correspond to notoriously difficult non-convex optimization problems.\n\n**S3.** The fact that no rank parameters need to be chosen is another important strength. Choosing an appropriate rank is typically also difficult for more traditional tensor approaches. Intuitively, it seems like it might be easier to choose appropriate levels of interaction between modes in the proposed framework than it is to choose rank in more traditional approaches.\n\n**S4.** Although the experiments are not very comprehensive, the speed-up results in the few experiments that are provided are impressive.\n\n**--- Weaknesses ---**\n \nThe proposal in the paper is very interesting. However, unfortunately there are points that aren't explained properly and that are difficult to follow.\n\n**W1.** The optimization procedure---arguably a key advantage of the method---is difficult to understand. \n- Concepts like \"natural parameters\", \"flats\" and \"natural gradient\" are never properly defined. Since most of the ICLR audience are unlikely to be familiar with these concepts, they should be explained. If there's not enough room in the main paper, this could be put in the appendix.\n- The vectors $\\theta^B$ and $\\eta^B$ are, as far as I can tell, length-$|\\mathcal{B}|$ vectors that contain the entries from $\\theta$ and $\\eta$ corresponding to indices in $B$, i.e., indices for which $\\theta$ is zero: $\\theta_b = 0$ for all $b \\in B$. With this in mind, it seems like the stepping scheme in Eq (6) only updates entries in $\\theta$ that are already zero. For a distribution $Q_1$ to lie in the e-flat as shown in Fig 1 (a), doesn't all entries $\\theta_b$ for $b \\in B$ have to be fixed to zero? In that case, why are those entries updated in Eq (6)? It seems like all this optimization procedure does is update parameters we know should be zero, while leaving the other parameters unchanged...\n- The complexity analysis in Sec 2.6 is hard to follow since it doesn't refer back to the particular computational steps discussed earlier.\n- What termination criteria is used in the gradient descent method? This is important when interpreting the experiment results.\n\n**W2.** The practical connection between the $\\theta$ parameters (e.g., in Eq (7)) and the $H$ parameters (e.g., Eq (9)) isn't explained well. The optimization is done over the parameters in $\\theta$, but the conditions (e.g., set $n$-body interactions to zero for $n \\geq 3$) are imposed on the $H$ parameters, which are sums of $\\theta$ parameters. It is not clear to me how such conditions on the $H$ parameters translate to conditions on the $\\theta$ parameters. For example, if $H_{i_k, i_m}^{(k,m)}$ is zero, then what indices should be added to the zero-index set $B$? The parameter count in Eq (15) counts the number of parameters in the $H$ factors---how does this relate to the number of $\\theta$ parameters that are actually being optimized?\n\n**W3.** The comparison between specific $n$-body approximations and conventional tensor decompositions is vague. For example, in Sec 3 you say \"if we impose that decomposed factors can be represented as products with hyper-diagonal tensors $\\Omega$, this decomposition is equivalent to a cyclic two-body approximation.\" It's not clear to me what this means. How should the symbol $\\simeq$ in the figures be interpreted? Since fitting a tensor ring is a non-convex problem, it's clear that they can't be equivalent.\n\n**W4.** Examples of things throughout the text that aren't well-defined/hard to follow:\n- Sec 2.1: \"dual flatness and orthogonality\" of coordinate systems. What does this mean?\n- Sec 2.3: \"element-wise cyclic product of matrices\". \n- While terms like \"partition function\" may be standard in physics, they will be unfamiliar to most members of the ICLR audience.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "5gBc01d2rL", "reviewer": "Reviewer_w947", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In this paper, a new tensor decomposition for nonnegative tensors is proposed interpreting the tensor as a probability measure over a multidimensional discrete space given by its indices commonly known as a random field. \nThe authors used tools from Physics in which an Energy function is defined containing terms that corresponds to potentials associated to first order, second order, etc. interactions between variables in the random field. \nIn this model, instead of using the concept of rank of the decomposition, the complexity is regulated by the order and map of interactions between variables (indices). \nFor the optimization of the model to fit it to some data, they applied Information Theory tools so that the optimization problem is reduced to a convex one and a natural gradient-based algorithm is derived to minimize a Kulback-Liebler distance.\nThe paper demonstrates that, for example, by considering a cyclic 2-order interaction the obtained model is equivalent to a Tensor Ring (TR) network.\nExperimental results using only a cyclic 2-order interaction model are included demonstrating the good properties in terms of accuracy and computation time of the proposed learning algorithm.\n", "summary_of_the_review": "A novel tensor decomposition for nonnegative tensors is proposed and an efficient algorithm supported by information theory ideas is provided. Although the experimental evaluation could be highly improved, I think the community working on tensor decompositions will appreciate this work because it provides new models that are rank-free and with provable convergent algorithms.", "strength_and_weaknesses": "Strengths:\n-\tOriginal model with nice properties such as its rank-free feature and the convexity of the associated optimization problem.\n-\tStrong theoretical ground based on Information Theory is provided.\n-\tVery well written and clear presentation of the model and algorithm\n\nWeaknesses:\n-\tThe model only applies to nonnegative tensors, which limits its applicability to many tensor data problems. Current title: “MANY-BODY APPROXIMATION FOR TENSORS” suggests a general model. A more appropriate title would be “MANY-BODY APPROXIMATION FOR NONNEGATIVE TENSORS”\n-\tThe equivalence between the proposed model (cyclic Many Body) and Tensor Ring (TR) should be explored more deeply. For example, given a cyclic Many Body model, it is shown that it is equivalent to a TR network. However, it is not clear what is the rank of the obtained TR model. Also, is the obtained TR-rank the minimal rank? Is there a way to check the actual rank of the proposed model? \n-\tThe experimental results are rather limited:\n1)\tThey only consider a cyclic 2-order interactions model in the experiments, which is equivalent to a TR network. I think more experiments considering higher order interactions would provide rich insights about the expression power of the new model compared to classical tensor decompositions.\n2)\tI found the experiments with totally random tensors not very useful because in that case there is no structure on the data tensor itself. The value of this experiment, if any, should be better explained in the paper.\n3)\tThe results on real datasets are not properly discussed. The behavior of the proposed algorithm is not well understood. For example, the error on the 4DLFD dataset is larger than other algorithms without any clear explanation about this result. In general, it seems that HALS always provides the best results in terms of error. Is there any explanation about it?\n\nMinor issues:\n-\tThere is a typo or a problem of visualization in Figure 1. The KL distance should be between P and \\bar{P}.\n-\tIn section 2.5, I think the following sentence needs to be fixed: “any tensor P can be represented by …” -> “any tensor P can be approximated by … ”, because in general not any tensor can exactly represented as a rank-1 tensor.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject