Fast Tucker Rank Reduction for Non-Negative Tensors Using Mean-Field Approximation
Year: 2022

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.

heading: 1 Introduction
text: A multidimensional array, or tensor, is a 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 [30]. In 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 [9, 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, non-negative low-rank approximation has also been developed, not only for matrices such as NMF [26], but also for tensors [27]. In particular, non-negative Tucker decomposition (NTD) [21] and its efficient variant lraSNTD [43] approximate a 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 classification [23], recommendation [37], and denoising [12], efficient low-rank approximation remains fundamentally challenging. Even the simplest case, the rank-1 approximation in terms of minimizing the Least Squares (LS) error between a given tensor and a low-rank tensor, is known to be NP-hard [16]. Various methods have been developed to efficiently find approximate solutions in polynomial runtime [7, 8, 10, 25, 42]. If we use the Kullback–Leibler (KL) divergence instead of the LS error as a cost function, we can alleviate the problem as the best rank-1 approximation can be obtained in the closed formula [19]. 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-Tucker-rank approximation method for 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 (NeurIPS 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 a 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 a sufficient condition for such parameters, called the bingo rule, to reduce Tucker ranks of tensors. We then show that low-rank approximation is achieved by m-projection; this is one of the two canonical projections in information geometry [2], where a 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-negative tensors can be exactly solved by a mean-field approximation, a well-established method in physics that approximates a joint distribution by independent distributions [40], as we can represent any non-negative rank-1 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 bym-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 a projection onto a 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 in our theoretical discussion.

heading: 2 The Proposed Rank Reduction Algorithm
text: We propose a Tucker rank reduction algorithm, called the Legendre Tucker rank reduction (LTR), which transforms a given non-negative tensor P ∈ RI1×···×Id≥0 into a tensor that approximates P with arbitrary reduced Tucker rank (r1, . . . , rd) specified by the user. The term “Legendre” 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.

heading: 2.1 Problem setup and notation
text: First we define the Tucker rank of tensors and formulate the problem of non-negative Tucker rank reduction. The Tucker rank of a dth-order tensor P ∈ RI1×···×Id is defined as a tuple (Rank(P(1)), . . . , Rank(P(d))), where each P(k) ∈ RIk× ∏ m6=k Im is the mode-k expansion of the tensor P and Rank(P(k)) denotes the matrix rank of P(k) [11, 15, 39]. See Supplement for definition of the mode-k expansion. If the Tucker rank of a tensor P is (r1, . . . , rd), it can be always decomposed as
P = ∑r1 i1=1 · · · ∑rd id=1 Gi1,...,ida (1) i1 ⊗ a(2)i2 ⊗ · · · ⊗ a (d) id
with a tensor G ∈ Rr1×···×rd , called the core tensor of P , and vectors a(k)ik ∈ RIk , ik ∈ {1, . . . , rk}, for each k ∈ {1, . . . , d}, where ⊗ denotes the Kronecker product [24]. If every element of the Tucker rank is the same as 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 a given non-negative tensor by a non-negative lower Tucker rank tensor. We denote by [d] = {1, 2, . . . , d} for a positive integer d and denote by Pa(k):b(k) the subtensor obtained by fixing the range of kth index to only from a to b.

heading: 2.2 The LTR Algorithm
text: 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 P is given by
P = λs(1) ⊗ s(2) ⊗ · · · ⊗ s(d), (1)
Algorithm 1: input :Tensor P , target Tucker rank r = (r1, . . . , rd) output :Rank reduced tensor Q LTR(P ,r)
(I1, . . . , Id)← the size of the input tensor P foreach k = 1, . . . , d do
Construct {c1, . . . , crk} ⊆ [Ik] by random sampling from [Ik] without replacement, where we always assume that c1 = 1 and ci < ci+1. foreach l = 1, . . . , rk do if cl 6= cl+1 − 1 then
Replace the subtensor P c (k) l :c (k) l+1−1 of P by its rank-1 approximation as P c (k) l :c (k) l+1−1 ← BESTRANK1(P c (k) l :c (k) l+1−1 )
Q ← P return Q
BESTRANK1(P) (I1, . . . , Id)← the size of the input tensor P foreach k = 1, . . . , d do
foreach ik = 1, . . . , Ik do s (k) ik ←∑I1i1=1 · · ·∑Ik−1ik−1=1∑Ik+1ik+1=1 · · ·∑Idid=1 Pi1,...ik−1,ik,ik+1,...,id
λ← sum of all elements of P P ← λs(1) ⊗ s(2) ⊗ · · · ⊗ s(d) return P
where each s(k) = (s(k)1 , . . . , s (k) Ik ) with k ∈ [d] is defined as
s (k) ik = I1∑ i1=1 · · · Ik−1∑ ik−1=1 Ik+1∑ ik+1=1 · · · Id∑ id=1 Pi1,...,id
and λ is the (1− d)-th power of the sum of all elements of P . Now we introduce LTR, which iteratively applies the above rank-1 approximation to subtensors of a tensor P ∈ RI1×···×Id . When we reduce the Tucker rank of P to (r1, . . . , rd), LTR performs the following two steps for each k ∈ [d]: Step 1: We construct C = {c1, . . . , crk} ⊆ [Ik] by random sampling from [Ik] without replacement, where we always assume that c1 = 1 and cl < cl+1 for every l ∈ [rk − 1]. Step 2: For each l ∈ [rk], if cl 6= cl+1 − 1 holds, we replace the subtensor Pc(k)l :c(k)l+1−1 of P by its rank-1 approximation obtained by Equation (1).
Step 1 of LTR requires O(r1 + r2 + · · ·+ rd) since we only need to sample integers from 1, 2, . . . , Ik for each k ∈ [d] using the Fisher-Yates method. Since the above procedure repeats the best rank-1 approximation at most r1r2 . . . rd times, the worst computational complexity of LTR is O(r1r2 . . . rdI1I2 . . . Id). The choice of C in Step 1 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, we can directly choose these indices for C instead of random sampling to obtain a more accurate reconstructed tensor.
We provide the algorithm of LTR in algorithmic format in Algorithm 1.

heading: 3 Theoretical Analysis of LTR
text: To theoretically guarantee that LTR always reduces the Tucker rank to (r1, . . . , rd), 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.

heading: 3.1 Modeling tensors as probability distributions
text: 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-linear model on partially ordered sets (poset) introduced by Sugiyama et al. [35] for our probabilistic formulation as they have studied tensors on a statistical manifold [36] using the formulation. We see that the low Tucker rank tensor approximation for positive tensors in terms of the KL divergence is formulated as a convex optimization problem.
In the following, we assume that an input positive tensor P ∈ RI1×···×Id>0 is always normalized so that the sum of all elements is 1, and we regard P as a discrete probability distribution whose sample space is the index set of tensors Ωd = [I1] × · · · × [Id]. Any positive normalized tensor can be described by canonical parameters (θ)i1,...,id = (θ1,...,1, . . . , θI1,...,Id) as
Pi1,...,id = exp (∑i1 i′1=1 · · · ∑id i′d=1 θi′1,...,i′d ) . (2)
The condition of normalization is exposed on θ1,...,1 with Ω+d = Ωd\{(1, . . . , 1)} as
θ1,...,1 = − log  ∑ (i1,...,id)∈Ω+d exp (∑i1 i′1=1 · · · ∑id i′d=1 θi′1,...,i′d ) . (3) This belongs to the log-linear model, and a parameter vector (θ)i1,...,id in Equation (2) uniquely identifies the normalized positive tensor P . Therefore, (θ)i1,...,id can be used as an alternative representation of P . If we see each probabilistic distribution – a normalized tensor in our case – as a point in a statistical manifold, (θ)i1,...,id corresponds to a coordinate of the manifold, which is a typical approach in information geometry and known as the θ-coordinate system [2].
It is known that any distribution in an exponential family log pθ(x) = C(x) + ∑k i=1 θiFi(x)− ψ(θ) for the normalization factor ψ(θ) and canonical parameters (θ)i can be also uniquely identified by expectation parameters ηi = Epθ [Fi] [4]. In our modeling in Equation (2), which clearly belongs to the exponential family, each value of the vector of η-parameters (η)i1,...,id is written as follows and uniquely identifies a normalized positive tensor P:
ηi1,...,id = ∑I1 i′1=i1 · · · ∑Id i′d=id Pi′1,...,i′d . (4)
Hence we can also use (η)i1,...,id as a coordinate system of the set of distributions, known as the η-coordinate system in information geometry [2]. As shown in Sugiyama et al. [35], by using the Möbius function [32, 34] inductively defined as
µ i′1,...,i ′ d
i1,...,id =  1 if ik = i′k for all k ∈ [d], −∏dk=1∑i′k−1jk=ik µj1,...jdi1,...,id if ik < i′k for all k ∈ [d], 0 otherwise,
each distribution P can be described as Pi1,...,id = ∑ (i′1,...,i ′ d)∈Ωd µ i′1,...,i ′ d i1,...,id ηi′1,...,i′d . (5) using the η-coordinate system. The normalization condition is realized as η1,...,1 = 1.
In information geometry, it is known that θ- and η-coordinates are connected via Legendre transformation [2]. The remarkable property of this pair of coordinate systems is that they are orthogonal with each other and we can combine them to define the (θ, η)-coordinate as a mixture coordinate. If we specify either θi1,...,id or ηi1,...,id for every (i1, . . . , id) ∈ Ωd, we can always uniquely identify a normalized positive tensor P .

heading: 3.2 Bingo rule: θ representation of Tucker rank condition
text: We theoretically derive the relationship between the Tucker rank of tensors and conditions on the θ-coordinate system. We achieve Tucker rank reduction using such conditions instead of directly imposing constraints on elements of tensors and solving it using classical methods such as the Lagrange multipliers. The advantage of our approach is that we can formulate low-rank approximation as a projection of a distribution P on a subspace of the (θ, η)-space, which always becomes a convex optimization due to the flatness of its subspace [2].
Definition 1 (Bingo). Let (θ)(k)ij = (θ (k) 11 , . . . , θ (k) IkK ) with K = ∏ m 6=k Im be the θ-coordinate representation of the mode-k expansion of a tensor P ∈ RI1×···×Id>0 . If there exists an integer i ∈ [Ik] \ {1} such that θ(k)ij = 0 for all j ∈ [K] \ {1}, we say that P has a bingo on mode-k. Proposition 1 (Bingo and Tucker rank). If there are bk bingos on mode-k, it holds that
Rank(P(k)) ≤ Ik − bk.
Therefore, for any tensor P ∈ RI1×···×Id>0 such that it has bk bingos for each kth index, we can always guarantee that its Tucker rank is at most (I1 − b1, . . . , Id − bd).

heading: 3.3 Projection theory for LTR
text: 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 θ- and η-coordinate systems, two types of geodesics, m- and e-geodesics from a point P to Q in a manifold can be defined as
{ Rt | Rt = (1− t)P + tQ } , { Rt | logRt = (1− t) logP + t logQ− φ(t) } , respectively, where 0 ≤ t ≤ 1 and φ(t) is a normalization factor to keep Rt to be a distribution. In the above definition, we regard P and Q as its corresponding coordinates points (θ)i1,...,id or (η)i1,...,id . A subspace is called e-flat when any e-geodesic connecting any two points in a subspace is included in the subspace. The vertical descent of an m-geodesic from a point P to Q in an e-flat subspace Me is called m-projection. Similarly, e-projection is obtained when we replace all e with m and m with e. The flatness of subspaces guarantees the uniqueness of the projection destination. The projection destinationRm orRe obtained by m- or e-projection onto Me or Mm minimizes the following KL divergence,
Rm = argmin Q∈Me DKL(P;Q), Re = argmin Q∈Mm DKL(Q;P).
LTR conducts an m-projection from an input positive tensor P onto the low-rank space B, which is the set of tensors, each of which satisfies a given set of bingos. In practical situations, a Tucker rank constraint (r1, . . . , rd) 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 (r1, . . . , rd). We call the low-rank subspace B bingo space. Since the low-rank space B is e-flat, the m-projection that minimizes the KL divergence is uniquely determined.
In m-projection into a space where some θ-parameters are constrained to be 0, every η-parameter with the same index with some restricted θ-parameter does not change [2]. Therefore, we always have
ηi1,...,id = η ′ i1,...,id if θ′i1,...,id 6= 0, (6)
where ηi1,...,id is the η-coordinate of an input positive tensor P and η′i1,...,id is that of the destination of m-projection onto B from P . Using the conservation law for the η-coordinate is our key insight to conduct the efficient m-projection onto the bingo space B.
Relationship to Legendre decomposition Legendre decomposition [36] is a tensor decomposition 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 a subspace that has some constraints on the θ-coordinate is taken as the initial state and moves by gradient descent inside the subspace to minimize the KL divergence from an input tensor. This operation is an e-projection, where the constrained θ-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 η-coordinates, as shown in Equation (6). Using this conservation law for η-coordinates, we can obtain an exact analytical representation of the coordinates of the projection destination without using a gradient method. Figure 1(a) illustrates the relationship between our approach and Legendre decomposition. Moreover, the Tucker rank is not discussed in the Legendre decomposition, so it is not guaranteed that Legendre decomposition reduces the Tucker rank, which is in contrast to our method.

heading: 3.4 Rank-1 approximation as mean-field approximation
text: Here we focus on the problem of rank-1 approximation for positive tensors and show the fundamental relationship with the mean-field theory. In the following, we 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 θ, η are corresponding parameters of θ- and η-coordinates. Let a one-body parameter be a parameter of which at least d− 1 indices are 1; for example, θ1,1,3,1 and η1,5,1,1 are one-body parameters for d = 4. Parameters other than one-body parameters are called many-body parameters. These namings come from the Boltzmann machine [1], which is a special case of the log-linear model [35], where a one-body parameter corresponds to a bias and a many-body parameter to a weight. We also use the following notation for one-body parameters of a dth-order tensor,
θ (k) j ≡ θ1, . . . , 1︸ ︷︷ ︸
k−1 ,j,1, . . . , 1︸ ︷︷ ︸ d−k , η (k) j ≡ η1, . . . , 1︸ ︷︷ ︸ k−1 ,j,1, . . . , 1︸ ︷︷ ︸ d−k for each k ∈ [d].
The rank-1 condition for positive tensors is described as follows using many-body θ parameters as a special case of Proposition 1.
Proposition 2 (rank-1 condition on θ). For any positive tensor P , rank(P) = 1 if and only if all of its many-body θ parameters are 0.
Note that the bingo constraints on θ-coordinates become unique if the target Tucker rank is (1, . . . , 1), hence the projection always finds the best rank-1 tensor that minimizes the KL divergence from an input tensor. This means that our information geometric formulation of the rank-1 approximation leads to the same result with the closed formula of the best rank-1 approximation given by Huang and Sidiropoulos [19]. Indeed, using the factorizability of η-parameters for rank-1 tensors, we can reproduce the closed formula in [19] (see Supplement for its proof).
We consider a rank-1 positive tensor P ∈ RI1×···×Id>0 and show that it is represented as a product of independent distributions, which leads to an analogy with the mean-field theory. In the rank-1 space, the normalization condition for dth-order tensors imposed on θ parameters is given as
θ1,...,1 = − log d∏ k=1 1 + Ik∑ ik=2 exp  ik∑ i′k=2 θ (k) i′k  (7)
by assigning 0 to every many-body θ 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 (2), we obtain
Pj1,...,jd = exp ( θ1,...,1 ) d∏ k=1 exp  jk∑ j′k=2 θ (k) j′k (7)= d∏ k=1
exp (∑jk
j′k=2 θ
(k) j′k ) 1 + ∑Ik ik=2 exp (∑ik i′k=2 θ (k) i′k ) ≡ d∏ k=1 s (k) jk ,
where s(k) ∈ RIk is a positive first-order tensor normalized as∑Ikjk=1 s(k)jk = 1. Since exactly one element is used in each s(k) to determine Pj1,...,jd , we can regard s(k) as a probability distribution with a single random variable jk ∈ [Ik]. The above discussion means that any positive rank-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 a wide range of domains such as statistics [31], game theory [6, 28], and information theory [5]. From the viewpoint of information geometry, Tanaka [38] developed a theory of mean-field approximation for Boltzmann machines [1], which is defined as p(x) = exp( ∑ i bixi + ∑ ij wijxixj) for a binary random variable vector x ∈ { 0, 1 }n with a bias parameter b = (b)i ∈ Rn and an interaction parameter W = (wij) ∈ Rn×n. 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-field approximation of Boltzmann machines is formulated as the projection from a given distribution onto the e-flat subspace consisting of distributions whose interaction parameters wij = 0 for all i and j, which is called a factorizable subspace. Since the distribution with the constraint wij = 0 for all i and j can be decomposed into a product of independent distributions, we can approximate a given distribution as a product of independent distribution by the projection onto a factorizable subspace. The m-projection onto the factorizable subspace requires knowing the expectation value ηi ≡ E[xi] of an input distributions and requires O(2n) computational cost [3], so we usually approximate it by replacing the m-projection with e-projection. The e-projection is usually conducted by a 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 s(k) by projecting P onto the subspace such that all many-body θ 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, mprojection can be directly performed in our formulation, which is computationally infeasible in the case of Boltzmann machines due to O(2n) cost. Moreover, the rank-1 space has the same property as the factorizable subspace of Boltzmann machines; that is, η 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 1(a).

heading: 3.5 Derivation of LTR
text: 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 a 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 a 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 (I1, . . . , Ik−1, rk, Ik+1, . . . , Id) for an input positive tensor P ∈ RI1×···×Id>0 . Let B(k) be the set of tensors satisfying bingos and ΩB(k) be the set of bingo indices for mode-k constructed in Step 1 of LTR:
B(k) = { P | θi1,...,id = 0 for (i1, . . . , id) ∈ ΩB(k) } . (8)
Note that P ∈ B(k) implies that the Tucker rank of P is at most (I1, . . . , Ik−1, rk, Ik+1, . . . , Id). Let P̃ be the destination of the m-projection and θ̃, η̃ be its corresponding parameters of θ- and
η-coordinates. From the definition of m-projection and the conservation low of η in Equation (6), the tensor P̃ should satisfy
θ̃i1,...,id = 0 for (i1, . . . , id) ∈ ΩB(k) , (9) η̃i1,...,id = ηi1,...,id for (i1, . . . , id) 6∈ ΩB(k) . (10)
As we can see from Equation (5), η-parameters identify a tensor P . In particular, to identify each element Pi1,...,id , we require η-parameters on only (i1, . . . , id) ∈ {i1, i1 + 1}× {i2, i2 + 1}× · · · × {id, id + 1}. For example, if d = 2, Pi1,i2 = ηi1+1,i2+1 − ηi1+1,i2 − ηi1,i2+1 + ηi1,i2 [35]. This leads to the fact that, for Ω̂B(k) such that
Ω̂B(k) = { (i1, . . . , id) | {i1, i1 + 1} × · · · × {id, id + 1} * ΩB(k) } , (11) it holds that Pi1,...,id = P̃i1,...,id for (i1, . . . , id) ∈ Ω̂B(k) . Therefore, all we have to do to reduce the Tucker rank is to change values for only (i1, . . . , id) 6∈ Ω̂B(k) . Such adjustable parts of P can be divided into some contiguous blocks, and we call each of them a subtensor of P on mode-k. In Figure 2, for example, we can find two subtensors P3(1):5(5) and P7(1):8(1) . By conducting the rank-1 approximation introduced in Section 2 onto each subtensor, the rank-1 bingo condition ensures the θ condition in Equation (9) and the conservation low of η in Equation (6) for the best rank-1 approximation ensures the η condition in Equation (10). Therefore, we obtain the tensor P̃ satisfying Equation (9) and (10) simultaneously by the rank-1 approximations on each subtensor on mode-k, which always belongs to B(k).
When the rank of only two modes are reduced Let us assume that the target Tucker rank of mode-k is rk and that of mode-l is rl. In this case, we need to consider two bingo spaces B(k) and B(l) associated with bingo index sets ΩB(k) and ΩB(l) . Let P(k) be the resulting tensor of mprojection of P onto B(k) and P(k,l) be the resulting tensor of m-projection of P(k) onto B(l). The parameters of P(k,l) satisfy
θ̃i1,...,id = 0 if (i1, . . . , id) ∈ ΩB(k) ∪ ΩB(l) , (12) η̃i1,...,id = ηi1,...,id otherwise. (13)
It is obvious that the m-projection from P(k) onto B(l) is equivalent to the m-projection from P onto B = B(l) ∩ B(k) The conservation low of η-parameters ensures Equation (12) after the rank-1
approximation of subtensors of P(k) in the mode-l. In this operation, the part of θ-parameters which are set to be 0 in the previous m-projection onto B(k) from P seems to be overwritten (see red panels in Figure 2(b)). However, as shown in Proposition 6 in Supplement, when the rank-1 approximation of a tensor where some one-body θ-parameters are already zero, the values of such one-body θ-parameters after the rank-1 approximation remains to be zero. As a result, Equation (12) holds; that is, we can obtain the tensor P(k,l) 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 P(l,k) = P(k,l); that is, the projection order does not matter. The projection sketch is shown in Figure 1(b).
Using the above discussion, we can derive Step 2 for the general case of Tucker rank reduction. Theorem 1 (LTR). For a positive tensor P ∈ RI1×···×Id>0 , them-projection destination onto the bingo space B = B(1) ∩ · · · ∩ B(d) is given as an iterative application of m-projection d times, starting from P onto subspace B(1), then from there onto subspace B(2), . . . , and finally onto subspace B(d).
The result of LTR does not depend on the projection order. The conservation law of η-parameters 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, in the specified bingo space B, that is, for a given non-negative tensor P , the output Q∗ of LTR satisfies that
Q∗ = argmin Q∈B DKL(P;Q).
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 a smaller search space for low-rank tensors. This search space allows us to derive an efficient algorithm without a gradient method, which always outputs the globally optimal solution in the space.

heading: 4 Numerical Experiments
text: We empirically examined the efficiency and the effectiveness of LTR using synthetic and real-world 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 [41] whose cost function is either the Least Squares (LS) error (NTD_LS) or the KL divergence (NTD_KL). The second method is sequential nonnegative Tucker decomposition (lraSNTD), which is known as the faster of the two methods [43]. Its cost function is the LS error; see Supplement for implementation details. We also obtained almost the same results as in Figure 3 with the KL reconstruction error (see Supplement).
Results on Synthetic Data. We created tensors with d = 3 or d = 5, where every Ik = n. We change n to generate various sizes of tensors. Each element is 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 is defined as the Frobenius norm between input and output tensors. Figure 3(a) shows the running time and the LS reconstruction error for randomly generated tensors with d = 3 and n = 30 with varying the target Tucker tensor rank. Figure 3(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, 512, 512, 3) tensor [18] and AttFace is a (92, 112, 400) tensor [33]. AttFace is commonly used in tensor decomposition experiments [21, 22, 43]. For the 4DLFD dataset, we chose the target Tucker rank as (1,1,1,1,1), (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,8,1), (3,3,16,16,1), (3,3,20,20,1), (3,3,40,40,1), (3,3,60,60,1), and (3,3,80,80,1). For the AttFace dataset, we chose (1,1,1), (3,3,3), (5,5,5), (10,10,10), (15,15,15), (20,20,20), (30,30,30), (40,40,40), (50,50,50), (60,60,60), (70,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 lraSNTD. Nevertheless, our experiments show that the approximation accuracy of LTR is competitive with other methods. This means that NTD and lraSNTD do not effectively treat linear combinations of bases.

heading: 5 Conclusion
text: We have derived a new probabilistic perspective to rank-1 approximation for tensors using information geometry and shown that it can be viewed as mean-field approximation. Our new geometric understanding leads to a 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 a theoretical analysis of tensors and we believe that our theoretical discussion will not have negative societal impacts.

heading: Acknowledgement
text: This work was supported by JSPS KAKENHI Grant Number JP20J23179 (KG), JP21H03503 (MS), and JST, PRESTO Grant Number JPMJPR1855, Japan (MS).

Review summaries:
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.

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).

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.