Fast Tucker Rank Reduction for Non-Negative Tensors Using Mean-Field Approximation
Authors: Kazu Ghalamkari
Mahito Sugiyama
Venue: NeurIPS 2021 Poster (2021)
Keywords: Low-rank approximation
Tucker rank
tensor
mean-field theory
Information geometry

Abstract: We present an efficient low-rank approximation algorithm for non-negative tensors. The algorithm is derived from our two findings: First,  we show that rank-1 approximation for tensors can be viewed as a mean-field approximation by treating each tensor as a probability distribution. Second, we theoretically provide a sufficient condition for distribution parameters to reduce Tucker ranks of tensors; interestingly, this sufficient condition can be achieved by iterative application of the mean-field approximation. Since the mean-field approximation is always given as a closed formula, our findings lead to a fast low-rank approximation algorithm without using a gradient method. We empirically demonstrate that our algorithm is faster than the existing non-negative Tucker rank reduction methods and achieves competitive or better approximation of given tensors.

Fast Tucker Rank Reduction for Non-Negative Tensors Using Mean-Field Approximation
Kazu Ghalamkari!2 Mahito Sugiyama!* 1.2
National Institute of Informatics
The Graduate University for Advanced Studies, SOKENDAL {gkazu, aemii.ac..j.. jp.
Abstract
We present an efficient low-rank approximation algorithm for non-negale non-negative ten- sors. The algorithm is derived from our two findings: First, we show that rank- I approximation for tensors can viewed as mean-field approximation by treat- ing each tensor as probability distribution. Second, we theoretically provide a sufficient condition for distribution parameters to reduce Tucker ranks of tensors; interestingly. this sufficient condition can be achieved by iterative application of the mean-field approximation. Since the mean-field approximation always given closed formula, our findings lead fast low-rank approximation algorithm without using a gradient method. We empirically demonstrate that our algorithm is faster than the existing non-negative Tucker rank reduction methods and achieves competitive or better approximation given tensors.
Introduction
multidimensional array, or tensor, fundamental data structure in machine learning and statistical data analysis, and extraction of the essential information contained in tensors has been studied extensively 14 20
For second-order tensors that is, matrices low-rank approximation by singular value decomposition (SVD) is well established [13]. SVD always provides the best low-rank approximation in the sense of arbitrary unitarily invariant norms B0] contrast, the problem of low-rank approximation becomes much more challenging for tensors higher than the second order, where the question of how to define the rank of tensors is even nontrivial. To date, various types of ranks the CP-rank [17 24,, the Tucker rank g 39 and the tubal rank 29 have been proposed, and low-rank approximation of tensors in terms of one of the above two ranks has been widely studied. Furthermore, nog-negative low-rank approximation has also been developed, not only for matrices such as NMF (26), but also for tensors [27). In particular, non non-negan -og-negaevve Tucker decomposition (NTD) [21 and its efficient variant IraSNTD 43 approximate given non-negative tensor by a tensor with the lower Tucker rank.
While these approximations have been widely used in various domains such as image classifi- cation 233, recommendation (37), and denoising [12], efficient low-rank approximation remains fundamentally challenging. Even the simplest case, the rank-1 1 approximation in terms of mini- mizing the Least Squares (I.S) error between a given tensor and low-rank tensor, known to be NP-hard [16]. Various methods have been developed to efficiently find approximate solutions in polynomial runtime 7 10 25 (42. we use the Kullback-Leibh
(KI.) divergence instead of the LS error cost function, we can alleviate the problem as the best rank-1 approximation can be obtained in the closed formula [1)] However, the general case of low-rank approximation in terms of the KL divergence is also still under development.
In this paper, we present a fast low- Now-Tucker-ran approximation method for non-neen non-negative tensors. To date, the majority of low-rank approximation methods are based on gradient decent using the
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

derivative of the cost function, which often requires careful tuning of initialization and/or a tolerance threshold. In contrast, our method is not based on a gradient method; the solution is directly obtained based on closed formula, which we derive from information geometric treatment of tensors. Through an alternative parameterization of tensors by treating them as probability distributions in a statistical manifold, we theoretically provide sufficient condition for such parameters, called the bingo rule, to reduce Tucker ranks of tensors. We then show that low-rank approximation achieved by m-projection; this is one of the two canonical projections in information geometry (2), where distribution (corresponding to a given non-negative tensor) is projected onto the subspace restricted by the bingo rule (corresponding to the set of non-negative low-rank tensors).
The key insight is that rank-1 approximation for non- non-negative tensors can be exactly solved by a mean-field approximation, well-established method in physics that approximates a joint distribution by independent distributions [40] as we can represent any non-negative rank-1 I tensor by a product of independent distributions. Moreover, we show that the bingo rule, our sufficient condition for tensor Tucker rank reduction, can be achieved by iterative applications of the mean-field approximation. This, combined with the fact that mean-field approximation is computed by m-projection in the closed form, enables us to derive our fast low-Tucker-rank approximation method without using a gradient method. Our theoretical analysis has a close relationship to [36], whose proposal, called Legendre decomposition, also uses information geometric parameterization of tensors and solves the problem of tensor decomposition by projection onto subspace. Although we use the same information geometric formulation of tensors, they did not provide any connection to the Tucker ranks, and Tucker rank reduction is not guaranteed by their approach. A limitation of our method is that it only treats non-negative tensors and cannot handle tensors with zero or negative values. Although experimental results show the usefulness of our method, even for tensors including zeros, we always assume that every element of an input tensor is strictly positive our theoretical discussion.
2 The Proposed Rank Reduction Algorithm
We propose Tucker rank reduction algorithm, called the Legendre Tucker rank reduction (LTR), which transforms a given non-negative tensor P IR/20 1xxxx-l into tensor that approximates P with arbitrary reduced Tucker rank (ri,..., ra) specified by the user. The term "Legendree" of our algorithm comes from the fact that the theoretical support of our algorithm uses parameterization of tensors based on the Legendre transformation, which will be clarified in Section 3
2.1 Problem setup and notation
First we define the Tucker rank of tensors and formulate the problem of non-negative Tucker rank reduction. The Tucker rank of dth-order tensor P e R/1x Rgixx-11a is defined as a tuple (Rank(p(1)). (ank( where each p(k) E R Im Ims is the mode-k expansion of the tensor P and Rank(p(()) denotes the matrix rank of ((k) I 15 39]. See Supplement for definition of the mode-k expansion. If the Tucker rank of a tensor P is (T1..... .Ta) it can be always decomposed as
P Σ
9
(1)
a
a
(d)
with tensor G E R. called the core tensor of P, and vectors aip.
(k) E RRI, is E (1....,re) r
for each & E (1...., d), where @ denotes the Kronecker product 24
If every element of the Tucker rank is the same each other element, it coincides with the CP rank. In this paper, we say that a tensor is rank-1 if its Tucker rank is (1,..., 1). The problem with non-negative Tucker rank reduction is approximating given non-negative tensor by non-negative lower Tucker rank tensor. We denote by [d] = (1, (1.2..... d) for a positive integer and denote by P g000 the subtensor obtained by fixing the range of kth index to only from a to b.
2.2 The LTR Algorithm
In LTR, we use the rank-1 approximation method that always finds the rank-1 tensor that minimizes the KL divergence from an input tensor [19]. The optimal rank-1 tensor P of is given by
P = (1(
@ s
(d)
(1)

---

Algorithm 1:
input :Tensor P, target Tucker rank (T1, (fi,..., ra)
output: :Rank reduced tensor Q (TR(P,r)
(I1...., Ia) t- the size of the input tensor P foreach ,......d do
Construct {cI,..
Cra S [Tx] by random sampling from [Tg] without replacement,
where we always assume that 01 1 and o C+41-
foreach / - - do
cy * a++ then
Replace the subtensor P (0 -1 of P by its rank-1 approximation as 10 f BESTRANK (P (k) 10/11
return Q
BESTRANKI(P)
(I....., La) + the size of the input tensor P
foreach .....,, do
foreach is
IR do
-= Ei=i=1 Emtt=1 ... Eig=1 P 1a2dh
À - sum of all elements of P (1(1) x s
s(d)
return P
where each (4b) (5)1) (k)
...... $1
(k)
with & [d]
is defined as
1+41
(k)
Σ
Σ
Sip
Pa..dd
4=1
and À is the (1 d1-h power of the sum of all elements of P.
Now we introduce LTR, which iteratively applies the above rank-1 1 approximation to subtensors of tensor P E R X2, When we reduce the Tucker rank of P to (ri....,ra), LTR performs the following two steps for each
[d]:
Step 1: We construct C = (ct,..., Cr. C [Tx] by random sampling from [7x] without replacement, where we always assume that 6 = and c G-11 for every 7 [rx 1].
Step 2: For each E [rw], if c * c++ holds, we replace the subtensor P P0 of P by its rank-1 approximation obtained by Equation (I].
Step of LTR requires O(ri + + +... + ra) since we only need to sample integers from 1, 2, 2,..... Jg for each k [d] using the Fisher- Yates method. Since the above procedure repeats the best rank-1 approximation at most ri72 Ta times, the worst computational complexity of LTR is ((rir2...rall ... Iz). The choice of C' in Step is arbitrary, which means that another strategy can be used. For example, if we know that some parts of an input tensor are less important than others, directly choose these indices for C instead of random sampling obtain accurate reconstructed tensor.
We provide the algorithm of LTR in algorithmic format in Algorithm
3 Theoretical Analysis of LTR
To theoretically guarantee that LTR always reduces the Tucker rank to (ri,..., ra), in the following subsections, we introduce information geometric analysis of the low-rank tensor approximation by treating each tensor as a probability distribution. The proof of the propositions are given in the supplementary material.

---

3.1 Modeling tensors as probability distributions
Our key idea for the derivation of LTR is the information geometric treatment of positive tensors in a probability space. We use the special case of a log-lineat model on partially ordered sets (poset) introduced by Sugiyama al. (35] for our probabilistic formulation as they have studied tensors on statistical manifold (36] using the formulation. We see that the low Tucker rank tensor approximation for positive tensors in terms of the KI. divergence is formulated as a convex optimization problem.
In the following, we assume that an input positive tensor P RS0 R4xx is always normalized so that the sum of all elements is I, and we regard P as discrete probability distribution whose sample space is the index set of tensors Pnt [T1] described by canonical parameters (0) (O)
[Ta]. Any positive normalized tensor can be
(0I.... ......,,.... ....., 01 1a) as
fa
P Pir Ad
exp
0 0/
(2)
9==1
6==1
The condition of normalization is exposed on @I....! with 22 Sdd 1dd1(11...... 1)}
1 log Σ
exp Σ t
Og Osg....t'
(3)
This belongs to the log-linear model, and a parameter vector (0)s in Equation 2 uniquely identifies the normalized positive tensor P. Therefore, (0) can be used as an alternative representation of P. If we see each probabilistic distribution a normalized tensor in our case - as a point statistical manifold, (f)......44 corresponds to a coordinate of the manifold, which is typical approach in information geometry and known as the 8-coordinate system (2).
It is known that any distribution in an exponential family log (x()) = C(x) + Σ 0xF(())- (0) for the normalization factor (0) and canonical parameters (0); can be also uniquely identified by expectation parameters 7): Ep [F] 44. In our modeling in Equation (2). which clearly belongs to the exponential family, each value of the vector of 7)- 1-parameters (n..a.....f is written as follows and uniquely identifies a normalized positive tensor P:
T/c1 E R==
Pk...
Hence we can also use (T))..... as coordinate system of the set of distributions, known as the n-coordinate system in information geometry [2]. As shown in Sugiyama et al. (35), by using the Möbius function 32 34 inductively defined as
if i for all k [d],
Ihm Ej=1, p
if is < i/ for all & E [d], otherwise,
each distribution P can be described
Pii...de E Se
(5)
using the y-coordinate system. The normalization condition is realized as Th...., =1 1.
In information geometry, it is known that 0- and n-coordinates connected via Legendre transfor- mation (2). The remarkable property of this pair of coordinate systems is that they are orthogonal with each other and can combine them define the (0, 11) cooordinate mixture coordinate. If we specify either 0i..... or Tii,..., for every (i1,-.., ia) E Zd we can always uniquely identify a normalized positive tensor P.
3.2 Bingo rule: 0 representation of Tucker rank condition
We theoretically derive the relationship between the Tucker rank of tensors and conditions on the 9-coordinate system. We achieve Tucker rank reduction using such conditions instead of directly imposing constraints on elements of tensors and solving using classical methods such as the Lagrange multipliers. The advantage of our approach is that we can formulate low-rank approximation as projection of a distribution P on a subspace of the (0, space, which always becomes a convex optimization due to the flatness of its subspace (2)

---

Definition 1 (Bingo), Let (0))) = (0 ,....,opkkk with K = Immp Im be the O-coordinate representation of the mode-k expansion of tensor P E I8-00 xxxxx If there exists an integer i E [Ix] (1} such that (()) 0 for all j [K] [1), we say that P has bingo on mode-k. Proposition 1 (Bingo and Tucker rank). If there are by bingos on mode-k, it holds that
ij
Rank(p(() S bu.
Therefore, for any tensor P E 1850 guarantee that its Tucker rank is at most (I1 bi...... ......... Ia ba).
R-xxx-Xaa such that it has by bingos for each kth index, we can always
3.3 Projection theory for LTR
We achieve low Tucker-rank approximation by projection onto a submanifold; that is, our proposed algorithm TLR performs projection in the viewpoint of information geometry. In a statistical manifold parameterized by 6- and c-corddinate systems, two types of geodesics, m- and e-geodesics from a point P Q in manifold can be defined
{ R2 R2 (1 -1- 4)P +
R Ro (1 (1-t) t) log P tlo Q 6(t) }.
respectively, where 0 < t S 1 and @(t) is a normalization factor to keep R₁ to be distribution. In the above definition, we regard P and Q as its corresponding coordinates points (o)......s or (1))i.....ig. A subspace is called e-flat when any e-geodesic connecting any two points in subspace included in the subspace. The vertical descent an m-geodesic from a point to Q an e-flat subspace M. 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 R., or Re obtained by e-projection onto Me or Mm minimizes the following KL divergence,
Rm = argmin (RI. Q); R₂ = argmin (2, M
QeM,
LTR conducts an m-projection from an input positive tensor P onto the low-rank space B, which the set of tensors, each of which satisfies given set of bingos. In practical situations, Tucker rank constraint (f....., ..t.....ral ra) is given as an input parameter, and LTR first implicitly constructs B in Step 1 by translating the Tucker rank condition into bingos, where any tensor in B has the Tucker rank at most (ri,..., rT). We call the low-rank subspace B bingo space. Since the low-rank space is e-flat, the m-projection that minimizes the KL divergence is uniquely determined.
In m- mrrooeection into a space where some 6 d-parameters are constrained to be 0, every 77-parameter with the same index with some restricted f parameterr does not change (2) Therefore, we always have
Thi T if 1f ....... 0
where Mii....... is the 17-coordinate of an input positive tensor and v/ Wi.... is that of the destination of m-projection onto B from P. Using the conservation law for the -ooordinnte is key insight to conduct the efficient m-projection onto the bingo space B.
a tensor decomposition Relationship to Legendre decomposition Legendre decomposition 36
method based on an information geometric view. However, their concept differs from ours in the following aspect. In the Legendre decomposition, any single point in subspace that has some constraints on the (-coordinate is taken as the initial state and moves by gradient descent inside the subspace minimize the KL divergence from an input tensor. This operation is e-projection, where the constrained f-coordinates do not change from the initial state. In contrast, we employ the m-projection from the input tensor onto the low-rank space by fixing some 7-coordinates, as shown in Equation (6). Using this conservation law for c-corrdinates, we can obtain an exact analytical representation of the coordinates of the projection destination without using gradient method. Figure aa) illustrates the relationship between our approach and Legendre decomposition. Moreover, the Tucker rank is not discussed in the Legendre decomposition, so it not guaranteed that Legendre decomposition reduces the Tucker rank, which is in contrast to our method.
5

---

(a)
P
(b)
Mean field equation
PR.(P.Pa)
Rank-1
P
approximation
P(e)
(xPPPPP
60h Du8r
P.
Legendre decomposition
Pg)
Poxy
P
P
P
6
Rank-1 space
m-projection
e-projection
Figure 1: (a) The relationship among rank-1 approximation, Legendre decomposition (66, and mean-field equation, where assume that the same bingo space is used in Legendre decomposition. solid line illustrates m-projection with fixing one-body ') parameters. P is an input positive tensor and P is the rank-1 tensor that minimizes the KL divergence from P. an initial point Legendre decomposition, which usually uniform distribution. P. a tensor of the 1-th step gradient descent in Legendre Decomposition. (b) The m-projection of common space of two different bingo spaces from P can achieved by m-projection into one bingo space and then m-projecting into the other bingo space.
3.4 Rank-1 approximation mean-field approximation
Here we focus on the problem rank-1 approximation for positive tensors and show the fundamental relationship with the mean-field theory. In the following, refer to the subspace consisting of positive rank-1 tensors as a rank-1 space. We use the overline for rank-1 tensors; that is, P is a rank-1 tensor and B, 7) are corresponding parameters of 0- and n-coordinates. Let one-body parameter be a parameter of which at least parameters for d = 4. Parameters other than one-body parameters are called many-body parameters. These namings come from the Boltzmann machine (1) which
indices are 1; for example, f1,1,31 and 71.5.111 are one-body
special case of the log-linear model (35), where one-body parameter corresponds to a bias and many-body parameter to weight. We also use the following notation for one-body parameters of a dth-order tensor,
0/h) = 11,..,, 1
Mi
m M.... 1.4.i for each k: E [d].
The rank-1 condition for positive tensors is described as follows using many-body 0 parameters as special case of Proposition]] T
Proposition 2 (rank-1 condition on 0). For any positive tensor P, raak = if and only if all of its many-body 0 parameters are O.
Note that the bingo constraints on f-coordinates become unique if the target Tucker rank is (1,..., 1). hence the projection always finds the best rank-1 tensor that minimizes the KI. divergence from an input tensor. This means that our information geometric formulation of the rank-1 approximation leads to the same result with closed formula of the best rank-1 approximation given by Huang and Sidiropoulos [19]. Indeed, using the factorizability of 77-aarmeterss for rank-1 tensors, we can reproduce the closed formula in [9 (see Supplement for its proof).
We consider a rank-1 positive tensor P 1 and show that it represented product of independent distributions, which leads to an analogy with the mean-field field theory. In the rank-1 space, the normalization condition for dth-order tensors imposed on parameters is given as
0 o = log I
1
exp
=2

---

by assigning 0 to every many-body 9 parameter in Equation (3. Note that the empty sum is treated as zero. Next, by substituting 0 for all many-body parameters in our model in Equation we obtain
Juu 2
(K)
Pi exp (0)
II exp
II k=1
(k) Si
Σ
Iu in
=
exp
where ((h) E R is a positive first-order tensor normalized as Σ 2 jk 1. Since exactly one
1x (k)
element is used each ((x) to determine Pi.....aa. we can regard (A)) as a probability distribution with single random variable je E [4x]. The above discussion means that any positive rank-1 1 tensor can be represented as a product of normalized independent distributions.
The operation of approximating a joint distribution as a product of independent distributions is called mean-field approximation. The mean-field approximation was invented in physics for discussing phase transition in ferromagnets [40]. Nowadays, it appears in wide range of domains such as statistics [3. game theory 6 2 and information theory [5]. From the viewpoint of information geometry, Tanaka [38 developed theory of mean-field approximation for Boltzmann machines [. which is defined as p(x)) = exp(E; b,244 + Eij W1y04,r)) for a binary random variable vector æ E 1² with a bias parameter b = (b), E IR" and an interaction parameter W = (wyy) E RRnx,, To illustrate that a rank-1 approximation can be regarded as a mean-field approximation, we prepare the mean-field theory of Boltzmann machines, as follows.
The mean-fiell approximation of Boltzmann machines is formulated as the projection from given distribution onto the e-flat subspace consisting of distributions whose interaction parameters Wij = 0 for all i and j, which is called factorizable subspace. Since the distribution with the constraint wi = 0 for all i and j can be decomposed into a product of independent distributions, we can approximate given distribution as product of independent distribution by the projection onto a factorizable subspace. The m-projection onto the factorizable subspace requires knowing the expectation value тソ E[æ:]] of input distributions and requires O(2") computational cost (3), so we usually approximate it by replacing the m- a--rojecctoo with e-projection. The e-projection is usually conducted by self-consistent equation called mean-field equation.
The analogy of LTR and mean-field theory is clear. In our modeling, a joint distribution P is approximated by a product of independent distributions (kk) by projecting P onto the subspace such that all many-body 0 parameters are 0, leading to the rank-1 tensor P. Since we can compute expectation parameters by simply summing the input positive tensor in each axial direction, m- projection can be directly performed in our formulation, which is computationally infeasible in the case of Boltzmann machines due to O(2") cost. Moreover, the rank-1 space has the same property as the factorizable subspace of Boltzmann machines; that is, 1) can be easily computed from (see Supplement for more details). This is the first time that the relationship between mean-field approximation and the tensor rank reduction has been demonstrated. The sketch of the relationship among mean-field equations and the rank-1 approximation is illustrated in Figure a).
3.5 Derivation of LTR
Finally, we prove that LTR successfully reduces the Tucker rank by extending the above discussion. We formulate Tucker rank reduction as an m-projection onto specific bingo space. This bingo space is constructed in Step 1 of LTR. Then, in Step 2, we can perform the m-projection using the closed formula of the rank-1 approximation without gradient method. We first discuss the case in which the rank of only one mode is reduced, followed by discussing the case in which the ranks of two modes are reduced.
When the rank of only one mode is reduced Let us assume that the target Tucker rank is (I,,.... I--1. Tke Ikt1...., Ia) for an input positive tensor P E TR xxx Let B+a) be the set of tensors satisfying bingos and S2g00) be the set of bingo indices for mode-k constructed in Step of LTR:
B(k)
{ P O
0 for (f.....,ia) E &ggas }.
Note that P E B(O) implies that the Tucker rank of P is at most (I1,... I rk. Jeti..... Ia). Let P. be the destination of the m-projection and 0, ńj be its corresponding parameters of 0- and
7

---

(a)
(b) Target Tucker Rank (7,5,2)
Target Tucker Rank (8,5,)
S 22
C Ug
Figure 2: Examples of LTR for (8, 2) tensor. Gg is bingo index. Tensor values and their тH- parameters fgg do not change, and the 7-parameters rank is (8, 5, 2). Bingos are on mode-2. In this example, the number of contiguous blocks 2. (b) Target rank (7,5, 2). Bingos are on modes 1 and and the number of contiguous blocks is 1 and for modes and respectively. We assume that project tensor onto ( followed by projecting it onto (3)2) After the second m-projection, the 0 e-parameters on red pancls are overwritten. However, these values remain to be zero after the second m-projeccion in our LTR..
O% ( 122 also do not change. (a) Target
7--cordiiaatss From the definition of m-projection and the conservation low of rj in Equation 6 the tensor P should satisfy
0 0i 0 for (ii,..., ia) E Pgars
(9)
Th Tii for (f,,..., id) 9gg0s
(10)
As we can see from Equation (5). 17.17-17 p-pparmmeters identify a tensor P. particular, to identify each element P Pi.....ia> (ia, àd - 1}. For example, d 2, Pi = 1 7/6 41,62 1 + 76. [35] This leads to the fact that, for Ωggo; such that
require 7-parameters only (41,..., ia) E f8i 1 + 1) x f 1} x....>
ago = (f1,..., id) f + 13 x .... fla. La + 13 & 2gw }
(11)
holds that Pi.....ha = Pi.....d for (d1,...,ia) E Pgges. Therefore, all we have to do to reduce the Tucker rank is to change values for only (i...... ia) & Qgan Such adjustable parts of P can be divided into some contiguous blocks, and we call each of them a subtensor of on mode-k. In Figure 2 for example, we can find two subtensors P3 and P7 By conducting the rank-1 approximation introduced in Section[]oo the 0 condition in Equation (9) and the conservation low of 17 in Equation (6) for the best rank- I approximation ensures the condition in Equation (10). Therefore, we obtain the tensor Po satisfying Equation (9) and (TO) simultaneously by the rank- approximations on each subtensor on mode-k, which always belongs to (k).
each subtensor, the rank- bingo condition ensures
When the rank of only two modes are reduced Let us assume that the target Tucker rank of mode-k is тı and that of mode-l is rı. In this case, we need to consider two bingo spaces BBk) and (8() associated with bingo index sets Qg(s and Qgn. Let P(x) the resulting tensor of m-
projection of P onto ((k) and P(k.r) be the resulting tensor of m-projection of P(k) onto (1)), The parameters of P(u.. satisfy
0 = if (i1,..., .ag.....,,, ia) E Aggs Ufggn
(12)
= Tis.....s otherwise.
(13)
obvious that the m- m-projection from P(() onto B()) is equivalent to the m-projection from P onto B (()) n (3k)) The conservation low of p-parmmeterss ensures Equation (12 after the rank- 1

---

(a) Random (30,30,30,30,30) tensor
(b)
Target rank (10,10,10)
Running time (sec.)
101.
Running time (sec.)
10"
whagee
ore
8000
1011
InadRIDD
1914
L5 errer
1500
101
L5 error
6000
1475
10
10th
10
1450
10-3
1425
7--0-00000000 --w----------------00000 B-g
0
300
207
Tensor rank
Tensor rank
Tensor size
Tensor size
(c)
AttFace
(d) 4DLFD
Running time (sec.)
450
Running time (sec.)
Zgge
10
errar
LS error
1064
350
355
10
2010
10"
10'
30"
10
10
100
1000
10000
md
1000
Size of core tensor
Size of core tensor
Size of tensor
Size core tensor
Figure 3: Experimental results for synthetic (a. b) and real-world (c, d) datasets. Mean errors + standard error for 20 times iterations are plotted. (a) The horizontal axis r for target tensor rank (r, r,,rr,r, r, r). (b) The horizontal axis is n for input (n, n, n.) tensor. (c, The horizontal axis is the number elements of the core tensor.
approximation of subtensors of P(x) in the mode-l. In this operation, the part of o-parameters which set be 0 in the previous m-projection onto g(k) from P seems to be overwritten (see red panels in Figure 2)). However, as shown in Proposition 6 in Supplement, when the rank-1 approximation of tensor where some one-body 6-parameters are already zero, the values of such one-body 0 e-parameters after the rank-l approximation remains to be zero. As a result, Equation (12) holds; that is, we can obtain the tensor P(w.i) satisfying Equation (12) and (13) simultaneously by the rank-1 approximations on each subtensor of P(k) on mode-l. We can also immediately confirm Pu.A) P(,,)); that is, the projection order does not matter. The projection sketch is shown in Figure [b).
Using the above discussion, we can derive Step for the general case of Tucker rank reduction.
Theorem (LTR), For positive tensor 1 xIs, the m-projection destination onto the bingo space B (()) n... B(d) is given as an iterative application of m-projection d times, starting from P onto subspace B3)), then from there onto subspace ( ..... and finally onto subspace (dd).
The result of LTR does not depend on the projection order. The conservation law of 11-arameters during the m-projection also ensures that LTR conserves the sum in each axial direction of the input tensor (see Supplement for its proof).
Since the m-projection minimizes the KL divergence from the input onto the bingo space, LTR always provides the best low-rank approximation, the specified bingo space B, that is, for a given non-negative tensor P, the output Q" of LTT satisfies that
argmin DKL (P: QEB
The usual low-rank approximation without the bingo-space requirement approximates a tensor by a linear combination of appropriately chosen bases. In contrast, our method with the bingo-space requirement approximates a tensor by scaling of bases. Therefore, our method has smaller search space for low-rank tensors. This search space allows us to derive an efficient algorithm without gradient method, which always outputs the globally optimal solution in the space.
4 Numerical Experiments
We empirically examined the efficiency and the effectiveness of LTR using synthetic and real-worke datasets. We compared LTR with two existing non-negative low Tucker-rank approximation methods.

---

The first method is non-negative Tucker decomposition, which is the standard nonnegative tensor decomposition method [4 whose cost function is either the Least Squares (LS) error (NTD LS) or the KI. divergence (NTD_K)) KI.). The second method is sequential nonnegative Tucker decomposition (IraSNTD), which is known as the faster of the two methods 43]. Its cost function is the L. error; see Supplement for implementation details. We also obtained almost the same results as in Figure with the KL reconstruction error (see Supplement..
Results on Synthetic Data. We created tensors with d = 3 or d = 5, where every 4a = n. We change n to generate various sizes of tensors. Each element sampled from the uniform continuous distribution on [0, 1]. To evaluate the efficiency, we measured the running time of each method. To evaluate the accuracy, we measured the LS reconstruction error, which defined as the Frobenius norm between input and output tensors. FigureBjar aa) shows the running time and the LS reconstruction error for randomly generated tensors with d and n = 30 with varying the target Tucker tensor rank. Figure b)) shows the running time and the LS reconstruction error for the target Tucker rank (10, 10, 10) with varying the input tensor size n. These plots clearly show that our method is faster than other methods while keeping the competitive approximation accuracy.
Results on Real Data. We evaluated running time and the LS reconstruction error for two real-world datasets. 4DLFD is a (9, 9, 51,, 51,,512,33 3) tensor [88] and AttFace is a (92, 112,
400) tensor [33]. 43 For the 4DLFD
AttFace is commonly used in tensor decomposition experiments [2] 22 dataset, we chose the target Tucker rank as (1,1,1,1,I), (2,2,2,2,1), (3.3,4,4.1),, (3,3,5,5.1,, (3,3,6,6,1), (3,3,7,7.1). (3,3,8,811), (3.3,16,16.1) (3.3.20,20.I), (3,3,40,40.1), (3,3,60,60,1), and (3.3,80.80.1) For the AuFacee dataset, we chose (1,1,1), (3,3,3), (5,5.5), (10,10,10) 10), (15,15,15), (20,20,20), (30.30,30), (40,40,40), (50,50,50), (60,60,60), 770,70,70), and (80,80,80). See dataset details in the Supplement. In both datasets, LTR is always faster than the comparison methods, as shown in Figure 3 c, d), with competitive or better approximation accuracy in terms of the LS error.
As described above, the search space of LTR is smaller than that of NTD and IraSNTD. Nevertheless, our experiments show that the approximation accuracy of LTR is competitive with other methods. This means that NTD and IraSNTD do not effectively treat linear combinations of bases.
5 Conclusion
We have derived a new probabilistic perspective to rank-1 approximation for tensors using informa- tion geometry and shown that it can be viewed as mean-field approximation. Our new geometric understanding leads novel fast non-negative low Tucker-rank approximation method, called LTR, which does not use any gradient method. Our research will not only lead to applications of faster tensor decomposition, but can also be a milestone of the research of tensor decomposition to further development of interdisciplinary field around information geometry and the mean-field theory.
This study is theoretical analysis of tensors and we believe that our theoretical discussion will not have negative societal impacts.
Acknowledgement
This work was supported by JSPS KAKENHI Grant Number JP20J23179 (KG), JP21 (MS). and IST, PRESTO Grant Number JPMJPRI855, Japan (MS).
References
[1] Ackley, D. H., Hinton, G. F., and Sejnowski, T. (1985). A learning algorithm for Boltzmann machines. Cognitive Science, 911:147-116.
[2] Amari, S. (2016). Information Geometry and Its Applications. Springer.
[3] Anderson, R. and Peterson, C. (1987). A mean field theory learning algorithm for neural networks. Complex Systems, 1:995-100 1019.
[4] Arwini, K. A. and Dodson, C. T. (2008). Information geometry. Springer.
10

---

[5] Bhattacharyya, C. and Keerthi, S. S. (2000). Information geometry and Plefka's mean-fied field theory. Journal of Physics A: Mathematical and General, 33(7):1307.
[6] Caines, P. E., Huang, M., and Malhamé, R. (2006). Large population stochastic dynamic games: Closed-loop McKean- Vlasov systems and the Nash certainty equivalence principle. Communica- lions Information and Systems, 6(3):221-252.
McKean-Vlasov
[7] da Silva, A. P., Comon, P., and de Almeida, A. L. (2016), A finite algorithm to compute rank-1 tensor approximations. JEEE Signal Processing Letters, 23(7):959-963,
[8] da Silva, A. P., Comon, P., and de Almeida, A. L. F. (2015). Rank-1 tensor approximation methods and application to deflation. arXiv preprint arXiv:1508.05277
[9] De Lathauwer, L., De Moor, B., and Vandewalle. J. (2000a). A multilinear singular value decomposition. SIAM Journal Matrix Analysis and Applications,
[10] De Lathauwer, L., De Moor, B., and Vandewalle, J. (2000b). On the best rank-1 and rank-(r 1, 2..... m) approximation of higher-order tensors. SIAM Journal on Matrix Analysis and Applications, 2144):1324-1342.
[11] De Silva, V. and Lim, L.-H. (2008). Tensor rank and the ill-posedness of the best low-rank approximation problem. SIAM Journal on Matrix Analysis and Applications, 308)
[12] Dong, W., Li, G., Shi, G., Li, X., and Ma, Y. (2015). Low-rank tensor approximation with laplacian scale mixture modeling for multiframe image denoising. In Proceedings of the IEEE International Conference on Computer Vision, pages 442-449.
[13] Eckart, C. and Young. G. (1936). The approximation of one matrix by another of lower rank. Psychometrika, 133):211-218.
[14] Grasedyck, L., Kressner, D., and Tobler, C. (2013). A literature survey of low-rank tensor approximation techniques. GAMM- MMMM-MMitelllngenn, 36(1):53-78.
[15] Hackbusch, W. (2019). Tensor spaces and numerical tensor calculus. Number v. 56 in Springer series computational mathematics. Springer, 2nd edition.
[16] Hillar, C. and Lim, L.-H. (2013). Most tensor problems are np-hard. Journal of the ACM (JACM), 60(6):1-39.
[17] Hitchcock, F. L. (1928). Multiple invariants and generalized rank of p-way matrix or tensor. Journal of Mathematics and Physics, 7(1-4):39-79
[18] Honauer, K., Johannsen, O., Kondermann, D., and Goldluecke, B. (2016). A dataset and evaluation methodology for depth estimation 4d light fields. In Asian Conference Computer Vision, pages 19-34. Springer.
[19] 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
[20] Ji, Y., Wang, Q., L., X., and Liu, J. (2019). A survey on tensor techniques and applications in machine learning. IEEE Access, 7:162950-162990.
[21] Kim, Y.-D. and Choi, (2007). Nonnegative Tucker decomposition. In 2007 IEEE Conference on Computer Vision and Pattern Recognition, pages 1-8. IEEE.
[22] Kim, Y.-D., Cichocki, A., and Choi, S. (2008). Nonnegative Tucker decomposition with alpha- divergence. In 2008 IEEE International Conference on Acousties, Speech and Signal Processing. pages 1182--1832 IEEE.
[23] Klus, S. and Gelß, P. (2019). Tensor-based algorithms for image classification. Algorithms, 12(11):240.
[24] Kolda, T. G. and Bader, B. W. (2009). Tensor decompositions and applications. SIAM review, 51(3):555-500
11

---

[25] Lathauwer, L. D. (2000). A multilinear singular value decomposition. SIAM Journal on Matrix Analysis and Applications, 2154):1253
[26] Lee, and Seung, H. S. (2001a). Algorithms for non-negative matrix factorization. In Leen, T., Dietterich, T., and Tresp, V., editors, Advances in Neural Information Processing Systems, volume 13, pages 556 566-562.44 MIT Press.
[27] Lee, D. D. and Seung, H. S. (2001b). Algorithms for non-negative matrix factorization. In Advances in neural information processing systems, pages 556-562.
[28] Lions, P.-L. and Lasry, J.-M. (2007). Large investor trading impacts on volatility, Annales de I'Instiut Henri Poincare (C) Non Linear Analysis, 24(2):311-32)
[29] Lu, C., Feng, J.. Chen, Y., Liu, W., Lin, Z., and Yan, S. (2019). Tensor robust principal component analysis with new tensor nuclear norm. IEEE Transactions on Pattern Analysis and Machine Intelligence, 42(4):925-
[30] Mirsky, L. (1960). Symmetric gauge functions and unitarily invariant norms. The Quarterly Journal of Mathematics, 11(1):50-59.
[31] Peterson, C. (1987). A mean field theory learning algorithm for neural networks. Complex Systems, pages 995-1019.
[32] Rota, G.-C. (1964). On the foundations of combinatorial theory Theory of Möbius functions. Z Walhsseheinlichkeesstteeooi 2:340-368.
[33] Samaria, F. S. and Harter, A. C. (1994). Parameterisation of stochastic model for human face identification. In Proceedings of 1994 IEEE Workshop on Applications of Computer Vision, pages 142. IEEE.
[34] Spiegel, E. and O'Donnell, C. J. (1997). Incidence Algebras. Marcel Dekker.
[35] Sugiyama, M., Nakahara, H., and Tsuda, K. (2017). Tensor balancing statistical manifold. arXiv preprint arXiv:
[36] Sugiyama, M., Nakahara, H., and Tsuda, K. (2019). Legendre decomposition for tensors. Journal of Statistical Mechanics: Theory and Experiment, 2019(12:1220117.
[37] Symeonidis, P., Nanopoulos, A., and Manolopoulos, Y. (2008). Tag recommendations based on tensor dimensionality reduction. In Proceedings of the 2008 ACM conference on Recommender systems, pages 43-50.
[38] Tanaka, (1999). theory of mean field approximation. Advances in Neural Information Processing Systems, pages 351-360.
[39] Tucker, L. R. (1966). Some mathematical notes on three-mode factor analysis. Psychometrika, 31(3)377--31.
[40] Weiss, P. (1907). L'hypothèse du champ moléculaire et la propriété ferromagnétique. Journal de Physique Théorique et Appliquée, 6(1):661-690)
[41] Welling, M. and Weber, M. (2001). Positive tensor factorization. Pattern Recognition Letters, 22(12):1255-1261
[42] Zhang, T. and Golub, G. H. (2001). Rank-one approximation to high order tensors. SIAM Journal on Matrix Analysis and Applications, 23(2):534-550
[43] Zhou, G., Cichocki, A., and Xie, S. (2012). Fast nonnegative matrix/tensor factorization based low-rank approximation. IEEE Transactions on Signal Processing, 60(6):2928-2944
12


Reviews:
[{"review_id": "zbxuCgM356k", "reviewer": "Reviewer_F6qZ", "summary": "The paper presents a low-rank tensor approximation algorithm for non-negative tensors. The algorithm is based on two observations. First, rank-1 approximation for tensors can be viewed as a mean-field approximation by treating each tensor as a probability distribution. Second, the authors provide a sufficient condition for distribution parameters to reduce Tucker ranks of tensors and, interestingly, this sufficient condition can be achieved by iterative application of the mean-field approximation. Since the mean-field approximation is always given as a closed formula, a fast low-rank approximation algorithm can be formulated without using a gradient method. The authors empirically demonstrate on several datasets that the algorithm is faster than the existing non-negative Tucker rank reduction methods ( in terms of achieving comparable approximations of given tensors).", "questions": "", "limitations": "", "rating": 7, "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": "Pros\n* Strong theoretical and algorithmic contribution with an original connection between tensor factorizations and mean-field approximation theory.\n* The proposed algorithm is novel and interesting \n* I'd rate the significance of the work as high as the paper draws on ideas from information geometry and physics. Hence the paper will likely be of interest to diverse sets of researchers. \n* The paper is clearly and well written. \n\nCons \n* The paper is a bit sparse on the experimental side but that's not necessarily a con given the algorithmic/theoretical contribution as mentioned above. \n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The paper does not include any narrative of any negative societal impact. That said, I don't immediately see an immediate negative impact beyond the hypothetical possibility of tensor decompositions being used for unethical downstream tasks. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "KH01ClVMMA3", "reviewer": "Reviewer_8bMP", "summary": "This paper derives a tensor decomposition method using mean-field approximation and information geometry. The main claim of this paper is the speed of computation compared to some basic baseline methods. An interesting direction but there are several directions that the paper can improve further.", "questions": "", "limitations": "", "rating": 6, "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": "This paper is technically interesting and can be a useful addition to the tensor decomposition method. However, this new method does not show any additional benefits other than the gain in speed.\n\nExperiments are not convincing. First of all, I feel that the dimensions for synthetic data are small. I would like to see the speed for much larger tensors. It would be informative if the Tucker ranks of tensor have different values, rather than considering equal r. I believe that it will be challenging to identify the Tucker rank in cases where there is not equal r. \nFrom the given plots it is difficult to identify differences in accuracy between the proposed method and the NTD. Further, it is hard to understand what is meant by the size of the core tensor, can the authors explain it (and how it is decided)?\nCan authors include simulations with larger tensors that with different Tucker ranks on each mode?   Can the authors provide more details on the prediction accuracy of tensor decomposition? \n\nIt would be useful to know if the proposed methods mean-field approximation can be extended to other tensor decomposition methods such as tensor networks. Tensor networks such as tensor trains have shown to be more computationally robust for higher-order tensors. \n\nThe paper is written well. I am not highly familiar with all the methods used discussed in the paper, however, the paper seems to be technically correct. The plots in the experiments section should be large enough to see lines clearly. Can the authors give core-computational steps in algorithmic format for easy understanding? ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Only the advantage of the speed of decomposition is discussed. It is useful if reconstruction accuracy is also improved.  \n\nAnother limitation is the restriction to Tucker decomposition, since recently developed methods such as tensor trains can be more feasible to use with higher-order tensors.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1.5 hour"}, "extra_text": {}}, {"review_id": "GCEI0hgAGqq", "reviewer": "Reviewer_TkCd", "summary": "The paper considers the low Tucker-rank approximation problem, uses the mean-field theory to interpret rank-one approximation for tensors by treating tensors as probability distributions, and then develops a two-step Tucker-rank reduction algorithm which does not involve gradient descent. Theoretical analysis of the proposed algorithm is provided in the information geometry framework. In particular, a Bingo rule is introduced to relate Tucker rank and modal expansion of a tensor in another coordinate system, which paves the way to reformulate the Tucker rank reduction problem as an m-projection. ", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "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": "The proposed tensor-rank reduction method is interesting, and theoretical discussions from the perspective of mean-field theory has certain novelty. Numerical results on various types of data sets are supportive and convincing. This work is important for low-rank approximations in tensor recovery and could be insightful for other tensor rank reduction methods. Minor issues: In line 19-20, there are more tensor ranks other than CP-rank and Tucker-rank, e.g., tubal rank. In line 84, what does \"the total sum of $\\mathcal{P}$\" mean? In line 94, what are \"uninformative indices\"? Maybe some examples could be referred. In the caption of Figure 1, \"Bingo space\" is not defined. In lines 211-212, how can the distribution independence assumption be guaranteed for generic tensors in the original form? Does the performance of the algorithm rely on this assumption? In Figure 3, are there any reasons for choosing the tensors with the same dimension across the modes?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The proposed framework works only for nonnegative tensors due to the exponential involved coordinate transform, which may limit the use in practice. Moreover, if the number of modes or some dimension is very large, the first step of LTR involves a random sampling which may waste some iterations and computation of (1) could be slow. Stopping criteria for LTR could be mentioned as well. ", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1"}, "extra_text": {}}, {"review_id": "CnACPN2mSxa", "reviewer": "Reviewer_3L1X", "summary": "This paper considers the problem of low-rank approximation for non-negative tensors. The authors propose a Tucker rank reduction algorithm based on information geometry, and provide a probabilistic perspective through mean-field approximation.", "questions": "", "limitations": "", "rating": 4, "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": "As the paper is mainly focused on the theoretical analysis of low-rank tensor approximation, I think the main drawback is that there is no theoretical guarantee for both the approximation error and the computational complexity of the algorithm. The performance of the algorithm is only shown by numerical experiments. In addition, the proposed sufficient condition and the mean-field connection seem to be of only theoretical interest, and fail to provide insights and guidance for practical applications.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper gives a new way to do Tucker rank reduction for non-negative tensors. The new approach puts additional constraints on the form of the decomposition (characterized by "bingo-space" defined in the paper) and allows the algorithm to find the optimal solution under these constraints very efficiently. The approach is quite novel and has some reasonable guarantees. The authors should clarify these guarantees in detail in the revised version.
Decision: Accept (Poster)