Non-negative low-rank approximations for multi-dimensional arrays on statistical manifold

Although low-rank approximation of multi-dimensional arrays has been widely discussed in linear algebra, its statistical properties remain unclear. In this paper, we use information geometry to uncover a statistical picture of non-negative low-rank approximations. First, we treat each input array as a probability distribution using a log-linear model on a poset, where a structure of an input array is realized as a partial order. We then describe the low-rank condition of arrays as constraints on parameters of the model and formulate the low-rank approximation as a projection onto a subspace that satisfies such constraints, where parameters correspond to coordinate systems of a statistical manifold. Second, based on information-geometric analysis of low-rank approximation, we point out the unexpected relationship between the rank-1 non-negative low-rank approximation and mean-field approximation, a well-established method in physics that uses a one-body problem to approximate a many-body problem. Third, our theoretical discussion leads to a novel optimization method of non-negative low-rank approximation, called Legendre Tucker rank reduction. Because the proposed method does not use the gradient method, it does not require tuning parameters such as initial position, learning rate, and stopping criteria. In addition, the flexibility of the log-linear model enables us to treat the problem of non-negative multiple matrix factorization (NMMF), a variant of low-rank approximation with shared factors. We find the best rank-1 NMMF formula as a closed form and develop a rapid rank-1 NMF method for arrays with missing entries based on the closed form, called A1GM.


Introduction
Tabular data in the real world is stored in computers as multidimensional arrays such as tensors or matrices. The low-rank approximation [1] is a popular task that approximates a multidimensional array by a low-rank array with a linear combination of fewer bases or principal components to extract features, find patterns, and store data efficiently with small memory storage [2]. A lot of applications are known in computer-vision [3], recommender-systems [4], denoising [5,6], and machine learning [7]. The low-rank representation is typically obtained by minimizing the reconstruction error from the input array based on singular value decomposition (SVD) or its extension to tensors, high-order singular value decomposition (HOSVD).
Non-negative low-rank approximation, a variant of low-rank approximations imposing non-negativity on both input and output, has been widely studied and applied to many practical problems because of the usefulness of non-negative representation in image processing, counting data, and audio analysis [8]. For instance, non-negative Tucker rank decomposition (NTD) approximates a tensor with a mode-product of a non-negative core tensor and matrices [9]. It is applied in various fields of data mining [10,11], bio-informatics [12], and mechanical engineering [13]. While nonnegative low-rank approximation has many applications, the non-negativity causes SVD-based optimization to be infeasible. As a result, it is implemented by gradient methods that iteratively reduce reconstruction errors [14][15][16]. However, the gradient methods require appropriate tuning for initial values, stopping criterion, and learning rates, which is often nontrivial in practice.
The present study proposes non-gradient-based methods for non-negative low-rank approximation for multidimensional arrays based on information geometric analysis. We treat each input array as a probability distribution using a log-linear model on a poset [17], where a structure of an input array is realized as a partial order. We describe the conditions for a tensor to be rank-reduced with θ -and η-parameters of the model, which are canonically used as coordinate systems of a dually flat manifold in information geometry. We capture the problem of non-negative low-rank approximation using such parameterization and realize it as a projection onto the subspace satisfying the constraints (low-rank space). In this novel formulation, we show that the rank-r approximation can be realized by combining the rank-1 approximation on each mode. By using the known formula of the tensor best rank-1 approximation [18], we developed a fast low-Tucker rank approximation method-Legendre Tucker rank reduction (LTR)-which is not based on the gradient method.
In addition, if the target rank is 1, we show that the low-rank space consists of the product of independent distributions. Approximating a joint distribution by the product of independent distributions is called mean-field approximation, which is a well-established method in physics [19]. The information-geometric formulation of low-rank approximation gives us a new perspective with which to view rank-1 approximation as a mean-field approximation. The rank-1 approximation is used when we want to quickly extract only the most significant factor from a large-scale input.
The flexibility of the log-linear model allows us to extend our formulation to related tasks, such as non-negative multiple matrix factorization (NMMF), that conduct simultaneous factor-sharing decomposition of multiple matrices; this task appears in purchase forecast systems [24] and recommender systems [25]. We find the best rank-1 NMMF formula as a closed form by describing the rank-1 condition with shared factors in a dually flat coordinate system. As an application of the formula, we develop an efficient method for rank-1 non-negative matrix factorization (NMF) for matrices with missing entries. The key idea is to transform the cost function of NMF for missing data into NMMF by increasing the missing values, enabling us to use the closed form.
Finally, we explain the geometric relationship between tensor rank-1 approximation and tensor balancing based on an information-geometric analysis. Tensor balancingthe operation of rescaling the axial sums of a tensor-has already been formulated using the log-linear model on posets [17]. We suppose that both the balancing and rank-1 conditions are imposed simultaneously on a tensor. In that case, all parameters of the distribution are determined, so the distribution is uniquely determined. This implies that balancing of a rank-1 tensor is unique.
A preliminary version of this paper was presented at the NeurIPS2021 [26] and AISTATS2022 [27]. In [27], we have discussed NMMF for three input matrices, but in this study we extend it for four input matrices and find the more general rank-1 approximation formula in a closed form. This enables us to exactly solve the missing NMF with a rank-2 weight matrix, as shown in Sect. 3.2.3. Moreover, we firstly point out the information-geometric relationship between rank-1 approximation and tensor balancing in Sect. 3.3.8.

Preliminaries
We provide the notations in Sect. 2.1 and briefly explain our modeling framework, the log-linear model on a partially ordered set (poset), which was introduced in [17], in Sect. 2.2.

Notation
Tensors are denoted by calligraphic capital letters such as T and P. Matrices are denoted by bold capital letters like X and Y. Each element is denoted as T i 1 ,...,i d and X i j . Mode-k expansion of a tensor T is denoted by T (k) . Vectors are denoted by lowercase bold alphabets like a and b. We regard vectors and matrices as particular cases of tensors. The total sum of a vector, a matrix, or a tensor is represented as S(·). The i-th component of a vector a is written in its non-bold letter as a i . The Kronecker product of two vectors a and b is denoted by (a⊗b), which is a rank-1 matrix, and each element is defined as (a ⊗ b) i j = a i b j . The I × J all-one and all-zero matrices are denoted by 1 I ×J and 0 I ×J , respectively. The identity matrix is denoted by I. The transpose of a matrix X is denoted by X T . The element-wise product of the two matrices A and B is denoted by A • B. We denote a pair of natural numbers n and m (≥ n) with [n, m] = {n, n + 1, . . . , m − 1, m} and abbreviate [1, m] as [m]. We use P a (k) :b (k) to denote the subtensor obtained by fixing the range of kth index to only from a to b. The set difference of B and A is denoted by B\A. When we use the Kullback-Leibler (KL) divergence D(P, T ) for non-negative tensors T and P, it is defined as follows: [9] D(P, T ) = We define Rank(X) as the rank of a matrix X.

Modeling
This paper discusses non-negative low-rank approximations of tensors, matrices, and multiple matrices. To treat various structures of such inputs, we use the log-linear model on posets, which is a flexible statistical model to treat various structures by designing appropriate partial orders [17]. Furthermore, we describe the low-rank conditions using the model's parameters. Consequently, we formulate the low-rank approximation as a projection onto a subspace satisfying the conditions that some natural parameters become zero. We provide an overview of the log-linear model on posets in Sect. 2.2.1 and the general projection theory in information geometry in Sect. 2.2.2.

Reminder of the log-linear model on poset
The log-linear model on a poset [17] is a generalization of Boltzmann machines [28], where we can flexibly design interaction between variables using partial orders. A poset ( , ≤) is a set of elements associated with a partial order ≤ on , where the relation "≤" satisfies the following three properties: For all x, y, z ∈ , (1) x ≤ x, We consider a discrete probability distribution p on a poset ( , ≤), which is treated as a mapping p : → (0, 1) such that x∈ p(x) = 1. Each element p(x) is assumed to be strictly larger than zero. We assume that the structured domain has the least element ⊥; that is, ⊥≤ x for all x ∈ . The log-linear model for a distribution p on ( , ≤) is defined as where θ(⊥) corresponds to the normalizing factor (partition function). The convex quantity defined as the sign inverse ψ(θ) = −θ(⊥) is called the Helmholtz free energy of p. This model belongs to the exponential family and θ corresponds to natural parameters except for θ(⊥).
The log-linear model's natural parameter θ uniquely identifies the distribution p. Using θ as a coordinate system in the set of distributions, which is a typical approach in information geometry [29], we can draw the following geometric picture: Each point in the θ -coordinate system corresponds to a distribution. Moreover, because the log-linear model belongs to the exponential family, we can also identify a distribution by expectation parameters defined as Using the Möbius function [30], inductively defined as otherwise, each distribution can be described as We write p η if p is determined by the expectation parameter η. We can also identify each point using the η-coordinate system. Each expectation parameter η(x) is consistent with the expected value E[F x (s)] for the function F x (s), such that F x (s) = 1 if x ≤ s and 0 otherwise [31]. Technically, η(⊥) is not an expectation parameter as it is always 1 from the definition. In addition, the θ -coordinate and the η-coordinate are orthogonal with each other, which guarantees that we can combine these coordinates together as a mixture coordinate and a point specified by the mixture coordinate also identifies a distribution uniquely [29].

Projection theory in information geometry
We consider the set of distributions S. We achieve non-negative low-rank approximations by projection onto a subspace Q ⊆ S. In a subspace parameterized by θand η-coordinate systems, two types of geodesics-m-and e-geodesics from a point q 1 ∈ Q to q 2 ∈ Q-can be defined as respectively, where 0 ≤ t ≤ 1 and φ(t) is a normalization factor to keep r t to be a distribution. 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-geodesic from a point p ∈ S to q in an e-flat subspace Q e 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 destination r m or r e obtained by m-or e-projection onto Q e or Q m minimizes the following KL divergence: We assume that distributions in S are parameterized by N natural parameters. Let Q be the set of distributions satisfying the condition f (θ 1:n ) = 0 for a linear function f (·) on a part of the natural parameters θ 1:n = (θ 1 , . . . , θ n ), where we assume that this part is from 1 to n (≤ N ) without loss of generality. Since a subspace defined by linear constrains in natural parameters is e-flat [29,Chapter 2.4], the m-projection from p onto Q is unique. This m-projection does not change the rest of the part of expectation parameters η n+1:N = (η n+1 , . . . , η N ) [32,Chapter 11.3]. We call this property expectation conservation law in m-projection, which is a key idea for getting analytical solutions to non-negative rank-1 approximations.

The proposed methods
We propose our low-rank approximation methods with a focus on their relationship to information geometry. In the following, we regard normalized inputs array as a distribution using the log-linear model on posets. The sample space of the distribution is the index set of the array. The probability of realizing an index is the corresponding value of the multidimensional array on the index. For example, for a positive matrix X ∈ R 2×3 >0 , the sample space of the model is the index set = [2] × [3] and the probability is given as p(i, j) = X i j / i j X i j for random variables (i, j) ∈ .
The distribution is parameterized by θ -parameter and η-parameters. We can describe low-rank conditions using θ and η to formulate the low-rank approximation as a projection onto the subspace satisfying the condition.
The advantage of our approach is that low-rank approximation via our formulation always becomes a convex optimization due to the flatness of its subspace. The limitation of our approach is that since our proposed methods are based on the loglinear model, they only treat positive arrays and cannot handle zero or negative values. Although experimental results show the usefulness of our methods, even for inputs array including zeros, we always assume that every element of the input is strictly positive in our theoretical discussion.
Since our approach is based on information geometric analysis, we naturally adopt the KL divergence as a cost function. It is known that optimizing the KL divergence prevents overfitting for sound data, noisy data, and purchasing data generated from Poisson distribution [33]. This is why many studies of robust low-rank approximations minimize the KL divergence [34,35].
Before discussing the low-rank approximations of tensors, we begin with a discussion of a particular case of tensors: matrices. We introduce the best rank-1 NMMF formula as a closed form in Sect. 3.1. Then, as an application of the formula, we introduce faster rank-1 NMF method for missing data in Sect. 3.2. By extending these discussions to tensors, we introduce a fast non-gradient-based low-Tucker-rank approximation method, LTR, in Sect. 3.3.

Best rank-1 NMMF
In this section, to find an example of how information geometric discussion enables us to formulate non-negative low-rank approximation without gradient methods, we focus on NMMF, a variant of matrix factorization. First, we define the rank-1 NMMF in Sect. 3.1.1. We then provide the best rank-1 approximation formula in Sect. 3.1.2. In Sects. 3.1.3-3.1.4, we introduce information geometric formulation of NMMF using the log-linear model and derive the closed form solution. Although they can be viewed as a proof of the closed form solution, we provide it in the main text because it includes various interesting properties between NMMF and information geometry, particularly the characterization of rank-1 approximation via parameters of the exponential family, which we believe is valuable for further development of this line of research.

Rank-1 NMMF simultaneously decomposes four matrices
into four rank-1 matrices w ⊗ h, a ⊗ h, w ⊗ b and c ⊗ b, respectively, using non-negative vectors w ∈ R I The cost function of NMMF is defined as and the task of rank-1 NMMF is to find vectors w, h, a, b and c that minimize the above cost. We assume that the scaling parameters α, β and γ are non-negative real numbers. We provide a sketch of the task in Fig. 1.

A closed formula of the best rank-1 NMMF
We give the following closed form of the best rank-1 NMMF that exactly minimizes the cost function in Eq. (4), which is one of our main theoretical contributions. This formula efficiently extracts only the most dominant shared factors from four input matrices.

Theorem 1 (the best rank-1 NMMF) For any four positive matrices
and three parameters α, β, γ ≥ 0, five non- , and c ∈ R L ≥0 that minimize the cost function in Eq. (4) is given as The time complexity to obtain the best rank-1 NMMF is O(I J + N J + I M + L M) because all we need is to take the summation of each column and row of the matrices X, Y, Z and U. If N = M = 0, our result in Theorem 1 coincides with the best rank-1 NMF minimizing the KL divergence from an input matrix X shown in [36]. There is a known closed formula for the best rank-1 approximation that minimizes the KL divergence from a tensor [18]. The above case for N = M = 0 corresponds to the formula with d = 2, where d is the order of the tensor.

Posets for NMMF
The input of NMMF is a tuple (X, Y, Z, U), where X ∈ R I ×J >0 , Y ∈ R N ×J >0 , Z ∈ R I ×M >0 , and U ∈ R L×M >0 . For simplicity, we normalize them beforehand so that their sum is 1; that is, S(X) + S(Y) + S(Z) + S(U) = 1. It is straightforward to eliminate this assumption using the property of the KL divergence, μD (X, Y) = D (μX, μY), for any non-negative number μ. We model these four matrices using a single discrete distribution on a partial ordered sample space.
To make a one-to-one mapping from (X, Y, Z, U) to a probability mass function p, we prepare the index set as where the subspace X corresponds to the index of X, Y to Y, Z to Z, and U to U. Then, we define the following partial order "≤" between each element (s, t) in the where (s, t) ∈ . The θ -parameters {θ 21 , . . . , θ N +I +L,J +M } are identified so that they satisfy , and l ∈ [L], and the normalizing factor θ 11 is uniquely determined from them. Figure 2 illustrates the partial order for the input triple (X, Y, Z, U).
There are other possible ways to model (X, Y, Z, U) as a probability distribution using different partial order structure. However, the solution formula that we obtain in Theorem 1 does not depend on the modeling.
Using results provided in [17] based on the incidence algebra between θ -and ηparameters, we can also obtain the η-parameters using the following formula: To make the following discussion clear, for all i ∈ [I ], j ∈ [J ], n ∈ [N ], m ∈ [M] and l ∈ [L], we define

Derivation of the exact solution of rank-1 NMMF
We define the terms one-body and many-body parameters to describe low-rank conditions for multiple matrices and tensors. Let a one-body parameter be a parameter of which at least d − 1 indices are 1; for example, θ 1,3 and η 5,1 are one-body parameters for d = 2, where d is the order of the tensor. Parameters other than one-body parameters are called many-body parameters. Gray-colored nodes in Fig. 2a correspond to one-body parameters. We treat the matrix case in this section-that is, d is always 2-and we will discuss the general case in Sect. 3.3.5. This categorization of parameters is inspired by the study of the Boltzmann machine [28], which is a special case of the log-linear model [17], where a one-body parameter corresponds to a bias and a many-body parameter corresponds to a weight. In the Boltzmann machine, if all weights are zero, the distribution modeled by the Boltzmann machine is represented by the product of independent distributions. In the same way, if all the many-body parameters are 0, the distribution is represented by the product of independent distributions. Approximating a joint distribution by a product of independent distributions is called mean-field approximation, and this approach is frequently used in physics. Further, as the following proposition shows, the rank of a multidimensional array is 1 if all many-body parameters are 0. We will discuss the relationship between the mean-field approximation and the rank-1 approximation in more detail in Sect. 3.3.5.

Proposition 1 (simultaneous rank-1 θ -condition) A tuple (X, Y, Z, U) is simultaneously rank-1 decomposable if and only if its all many-body θ -parameters are 0.
See Supplement for proof. We call a subspace that satisfies simultaneous rank-1 condition simultaneous rank-1 subspace. From the viewpoint of information geometry, we can understand the best rank-1 NMMF as follows (shown in Fig. 2b). The input of NMMF (X, Y, Z, U) corresponds to a point in the space described by the θ -coordinate system. The best rank-1 NMMF is an m-projection onto the simultaneous rank-1 subspace from the input point.
Since the m-projection is a convex optimization, we can get the projection destination by a gradient method. However, it requires appropriate settings for initial values, stopping criterion, and learning rates.
Our closed analytical formula of the projection destination in Theorem 1 addresses all the drawbacks of the gradient-based optimization. According to the expectation conservation law in this m-projection onto simultaneous rank-1 subspace, one-body η-parameters do not change in the m-projection. That is, for any where η is the expectation parameter of input, and η is the expectation parameter after the m-projection onto simultaneous rank-1 subspace. By the definition of expectation parameters, we obtain The expectation conservation law in Eq. (6) guarantees that the values of the lefthand sides do not change before the m-projection, and they do not change after the m-projection either. Since the sum of each matrix X, Y, Z, and U is represented by one-body η-parameters, the sum of each matrix does not change in the m-projection. Using these facts and multiplying these equations together, we can derive Theorem 1. A complete proof of Theorem 1 is provided in the Supplementary Material.

Rank-1 missing NMF
As an application of the closed form in Theorem 1, we develop an efficient method to solve rank-1 NMF for missing data. Rank-1 NMF for a given matrix X ∈ R I ×J ≥0 with missing values (rank-1 missing NMF) is the task of finding two non-negative vectors w ∈ R I ≥0 and h ∈ R J ≥0 that minimize a weighted cost function D (X, w ⊗ h) defined as for a binary weight matrix ∈ {0, 1} I ×J . The weight matrix indicates the position of missing values; that is, i j = 0 if the entry X i j is missing, i j = 1 otherwise. If the binary weight matrix satisfies Rank ( ) ≤ 2, we can find the exact solution for rank-1 missing NMF. After we mention the relationship between NMMF and missing NMF in Sect. 3.2.1, we demonstrate a way to find the best rank-1 missing NMF when it holds Rank ( ) ≤ 2 in Sect. 3.2.2-3.2.3. We also develop an efficient method for the general cases to find an approximate solution based on the closed formula. The proposed method is described in Sect. 3.2.4. 1

Connection between NMMF and missing NMF
Our discussion is based on the two following fundamental facts. First, we can regard NMMF as a special case of missing NMF. We assume that a binary weight matrix ∈ {0, 1} N +I +L,J +M and an input matrix K ∈ R N +I +L,J +M ≥0 are given in the form of where All of the elements of E and F are missing. We consider the rank-1 approximation of K as In this situation, the cost function of missing NMF is equivalent to that of NMMF [37]: The second fundamental fact is the homogeneity of rank-1 missing NMF, which ensures a factorization after row or column permutations can be reproduced by permutations after the factorization. See Supplement for proof.
Proposition 2 (Homogeneity of rank-1 missing NMF) Let NMF 1 ( , X) be the best rank-1 matrix w ⊗ h, which minimizes the cost function in Eq. (7). For any permutation matrices G and H, it holds that Therefore, using the closed formula of the best rank-1 NMMF in Theorem 1, we can solve the rank-1 missing NMF when we can relocate the position of missing values to the form Eq. (8) by row and column permutations.

Rank-1 missing NMF for grid-like missing
We introduce the term grid-like, which is defined as follows. As we describe in this section, we can regard missing NMF as NMMF that requires only three matrices X, Y, Z as input, which corresponds to the case L = 0 in Theorem 1.  (1) and j ∈ S (2) , 1 otherwise, then is said to be grid-like.
It holds that rank ( ) = 2 if is grid-like, but the converse is not true. We discuss the case rank ( ) = 2 but not grid-like in Sect. 3.2.3. Real-world tabular datasets tend to have missing values only on certain rows or columns. Therefore, the binary weight matrix often becomes grid-like in practice (we show example datasets in Sect. 4). Figure 3 illustrates examples of matrices with grid-like missing values.
When ∈ {0, 1} I +N ,J +M is grid-like, we can transform it in the form given in Eq. (8) with L = 0 using row and column permutations. Let S (1) (1) and | S (2) |= C (2) be the row and column index sets for zero entries in . For the block at the upper right of whose row and column indices are specified as where G and H are corresponding permutation matrices to G and H, respectively. We can obtain G and H as follows. First, we focus on row permutation G. We want to include each row j ∈ S (1) (1) , in B (1) by row permutation G, which can be achieved by any one-to-one mapping from S (1) The corresponding permutation matrix is given as where R k↔l is a permutation matrix, which switches the k-th row and the l-th row; that is, In the same way, any one-to-one mapping from S (2) (2) and B (2) (2) . The corresponding permutation matrix is given as which is also a symmetric matrix.
The above discussion leads to the following procedure of the best rank-1 missing NMF for an input matrix K if a binary weight matrix is grid-like. The first step is to find proper permutations G and H to collect the missing values in the upper right corner. In the next step, we obtained NMF 1 (G H, GKH) using the closed formula of the best rank-1 NMMF. In the final step, we operated the inverse permutations of G and H to the result of the previous step; that is, Note that G −1 = G T = G and H −1 = H T = H always holds since these permeation matrices are orthogonal and symmetrical.

Rank-1 missing NMF with Rank (8) ≤ 2
In this subsection, we show that we can relocate missing values into the form of Eq. (8) by column and row permutations if the binary weight satisfies rank ( ) = 2 and solve the rank-1 missing NMF as a rank-1 NMMF.
As we can confirm immediately, there are only two cases when the rank of a binary matrix is 1. In the first case, all of the elements i j are 1. This does not happen in our case because the number of missing values is assumed to be strictly larger than 0, resulting in = 1. In the second case, if there are rows or columns with all zero elements in , the matrix rank of can be 1. However, since rows and columns with all zero elements do not contribute to the cost function (7), we ignore such rows and columns. As a result, we discuss only the case of Rank( ) = 2.
We consider a rank-2 weight matrix ∈ {0, 1} N +I +L,J +M . There are two linear independent column bases if and only if Rank ( ) = 2 since row-rank and columnrank are always the same. We define them as a ∈ {0, 1} N +I +L and b ∈ {0, 1} N +I +L . We also assume that a = 0 and b = 0 since a zero-vector cannot be a basis. Then, for the binary matrix = [c (1) , . . . , c (J +M) ], it should be possible to write any column c (i) as c (i) = α i a + β i b using two bases a and b. Since c (i) is a binary vector, possible domains of α i and β i are limited and we analyze the domains in the following by separating them into three cases. In all three cases, we can rearrange into the form of Eq. (8) be a, b, or a + b. Since a and . , a, 1, . . . , 1, b, . . . , b], where H is a permutation matrix corresponding to the column permutation. After the column permutation, we conduct row permutation as follows. First, we define There is a one-to-one mapping G from S ∩ B c to S c ∩ B. The corresponding permutation matrix is given as where B c = [J + M]\B. By operating the permutation G on bases, we obtain, Using the fact G1 = 1, finally, we obtain which means G H is in the form of Eq. (8). Case 2: The bases are not disjoint, but one of them is one vector If the bases a and b are not disjoint but a = 1, it holds that  . , a, b, . . . , b].
In the same way as with Case 1, we obtain G H = [ã, . . . ,ã,b, . . . ,b], which corresponds to the form of Eq. (8) with I = 0, corresponding to NMMF for only two matrices Z and U.
The above discussion leads to the following procedure of the best rank-1 missing NMF for an input matrix K if the rank of the binary weight matrix is 2. The first step is to find the bases a and b of and find proper permutations G and H as described above to collect the missing values in the corners. The remaining step is the same as described in the final paragraph in Sect. 3.2.2.

Rank-1 missing NMF for the general case
If the rank of the binary weight matrix is strictly larger than 2, the above procedure cannot be directly applied. To treat any matrices with missing values, our idea is increase missing values so that the rank of becomes 2. However, the optimal way to increase missing values is not obvious to make rank-2. Then, we increase missing values so that becomes grid-like since the optimal way to increase missing values is clear. Although this strategy is counter-intuitive because we lose some information, which may cause a larger reconstruction error, we gain the efficiency instead of using our closed form solution in Theorem 1 with L = 0, and, as we empirically show in Sect. 4.1.2, the error increase is not significant in many datasets. Examples of this step are demonstrated in Fig. 3.
In the worst case, the number of missing values after this step becomes k 2 for k missing values. If every row or column has at least one missing value, all indices are missing after this step, for which our algorithm does not work. Thus, our method is not suitable if there are too many missing values in a matrix.
We illustrate an example of the overall procedure of A1GM in Fig. 4 and show its algorithm in Supplement. Since the time complexity of each process of A1GM is at most linear with respect to the number of entries of an input matrix, the time complexity is O((I + N )(J + M)) for input matrix X ∈ R (I +N )×(J +M) .

Legendre Tucker rank reduction
In this section, by designing posets appropriately, we extend the discussion so far to tensors and derive Legendre Tucker rank Reduction (LTR), which is a non-gradient method for non-negative low-Tucker-rank approximation. First, we define the task in Sect. 3.3.1. Then, we overview our fundamental ideas for LTR in Sect. 3.3.2, following the introduction to the algorithm of LTR in Sect. 3.3.3 and derive the LTR algorithm in Sects. 3.3.4-3.3.6, pointing out the relationship between the rank-1 approximation and mean field approximation. 2 Finally, we mention the relationship between the proposed LTR and related works in Sects. 3.3.7 and 3.3.8.

Low-Tucker-rank approximation for tensors
First we define the Tucker rank of tensors and formulate the problem of non-negative low-Tucker-rank approximation. The Tucker rank of a dth-order tensor P ∈ R I 1 ×···×I d is defined as a tuple (Rank(P (1) ), . . . , Rank(P (d) )), where each P (k) ∈ R I k × m =k I m is the mode-k expansion of the tensor P [38][39][40]. We describe the definition of mode-k expansion at the end of this subsection. If the Tucker rank of a tensor P is (r 1 , . . . , r d ), it can always be decomposed as with a tensor G ∈ R r 1 ×···×r d , called the core tensor of P, and vectors a (k) [41]. The core tensor and these vectors are often called factors.
In this paper, we say that a tensor is rank-1 if its Tucker-rank is (1, . . . , 1). Nonnegative low-Tucker-rank approximation approximates a given non-negative tensor P by a non-negative lower-Tucker-rank tensor T that optimizes the cost function D(P, T ). The task is not a non-negative factorization that imposes nonnegativity on factors, but a low-rank approximation that allows negative factors [42,43].

Mode-k expansion
The mode-k expansion [41] of a tensor P ∈ R I 1 ×···×I d is an operation to convert P into a matrix P (k) ∈ R I k × d m=1(m =k) I m focusing on kth mode. The relation between tensor P and its mode-k expansion P (k) is formally given as

Idea of LTR
In this subsection, we overview our core idea for LTR. LTR reduces Tucker rank of an input tensor P by known closed-formula of the best rank-1 approximation. As an example, here we reduce the Tucker rank of a positive tensor P ∈ R I 1 ×I 2 ×I 3 >0 to at most (r 1 , r 2 , I 3 ) as shown in Fig. 5. In the space of positive tensors, there exists a subspace F 1 consisting of positive tensors of the Tucker rank at most (r 1 , I 2 , I 3 ) and a subspace F 2 consisting of positive tensors of Tucker rank at most (I 1 , r 2 , I 3 ). We want to find a low-Tucker-rank tensor in F 1 ∩ F 2 that approximates P as close as possible.
First, we map a tensor to a probability distribution. Then, using natural parameters of the distribution, we describe sufficient conditions for reducing the Tucker rank of tensors, called bingo rule. For m = 1 and 2, we define e-flat subspace B (m) ⊂ F m , called bingo space, that satisfies the bingo rule. The projection from P onto B (1) can be conducted by using known rank-1 approximation formula onto the subtensor of P. Also, the projection from the point on B (1) onto B (1) ∩ B (2) can be conducted in the same way.
Bingo spaces only cover tensors generated by applying rank-1 approximations to subtensors along with each mode of a tensor. Therefore, the search space is smaller to at most (r 1 , r 2 , I 3 ) by the proposed method LTR. F 1 is the set of positive tensors with the Tucker rank at most (r 1 , I 2 , I 3 ) and F 2 with the Tucker rank at most (I 1 , r 2 , I 3 ). The best approximation tensor exists in F 1 ∩ F 2 , enclosed by the dotted lines. For m = 1, 2, there exist e-flat bingo spaces B (m) ⊂ F m . The projection onto B m can be performed by dividing P into subtensors along with mode-m direction and replacing each subtensor with its rank-1 approximation. The choice of bingo space is not unique than the traditional low-rank approximation, which approximates the tensor with an appropriately chosen basis and its coefficients, and there is no guarantee that LTR finds the best approximation; however, we can guarantee that LTR finds a tensor in the selected bingo spaces that minimizes the KL divergence from an input tensor. Such a smaller search space derived by the bingo rule makes our algorithm efficient without a gradient method. We discuss this point in more detail in Sect. 3.3.6.

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 [18]. The optimal rank-1 tensor P of P is given by where each s (k) = (s Now we introduce LTR, which iteratively applies the above rank-1 approximation to subtensors of a tensor P ∈ R I 1 ×···×I d . When we reduce the Tucker rank of P to (r 1 , . . . , r d ), LTR performs the following two steps for each k ∈ [d]: Step 1: We construct C = {c 1 , . . . , c r k } ⊆ [I k ] by random sampling from [I k ] without replacement, where we always assume that c 1 = 1 and c l < c l+1 for every l ∈ [r k − 1].
Step 2: For each l ∈ [r k ], if c l = c l+1 −1 holds, we replace the subtensor P c (k) of P by its rank-1 approximation obtained by Eq. (11).
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 Supplement.

Computational complexity of LTR
Step 1 of LTR requires O(r 1 + r 2 + · · · + r d ) since we only need to sample integers from 1, 2, . . . , I k for each k ∈ [d] using the Fisher-Yates method. Since the above procedure repeats the best rank-1 approximation at most r 1 r 2 . . . r d times, the worst computational complexity of LTR is O(r 1 r 2 . . . r d I 1 I 2 . . . I d ).

Posets for LTR
We derive the LTR algorithm by information geometric formulation of low-Tuckerrank approximation. The discussion is based on the log-linear model on poset. For simplicity, we normalize input tensor beforehand so that the sum is 1. To regard any positive tensor P ∈ R I 1 ×···×I d >0 as a distribution, we introduce the following partial order "≤" between each elements (i 1 , . . . , i d ) in the index set d = [I 1 ] × · · · × [I d ] of the tensor P, See Fig. 6 as an example for the poset d with I 1 = I 2 = I 3 = 3 and d = 3. We regard P as a discrete probability distribution whose sample space is the index set of P by log-linear model on ( d , ≤). Any positive normalized tensor can be described by natural parameters as A parameter vector (θ ) i 1 ,...,i d = (θ 2,...,1 , . . . , θ I 1 ,...,I d ) uniquely identifies the normalized positive tensor P. Therefore, (θ ) i 1 ,...,i d can be used as an alternative representation of P. Note that θ 1,...,1 is the normalizing factor and not independent from natural parameters. In our modeling in Eq. (13), which clearly belongs to the exponential family, each value of the vector of η-parameters (η) i 1 ,...,i d is written as follows and uniquely identifies a normalized positive tensor P: The normalization condition is realized as η 1,...,1 = 1. As shown in [17], by using the Möbius function [30] inductively defined as (a) (b) Fig. 6 a A poset ( 3 , ≤) corresponding to 3 × 3 × 3 tensor. The parameters on the gray nodes are one-body parameters. b The non-negative rank-1 approximation is formulated as a m-projection from input tensor to rank-1 space, minimizing KL divergence otherwise. each distribution P can be described as using the η-coordinate system. Note that, to identify the value of P i 1 ,...,i d , we need only

Information geometric view of rank-1 approximation
Before we dive into the general case of non-negative low-Tucker-rank approximation, we focus on the problem of rank-1 approximation for positive tensors and show the fundamental relationship with the mean-field theory. We describe the necessary and sufficient condition for the rank of a tensor to be 1 using θ -and η-parameters. We formulate tensor rank-1 approximation as a projection onto a subspace consisting of positive rank-1 tensors, which is called 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.
In the following, we use the one-body parameters and many-body parameters defined in Sect. 3.1.4 to describe the low-rank condition of tensors. For clarity, we also use the following notation for one-body parameters of a dth-order tensor, The rank-1 condition for positive tensors is described as follows using many-body θ parameters, and we have also succeeded in describing the rank-1 condition using the η-parameter. See Supplement for proof.

Proposition 4 (rank-1 condition as η-form) For any positive dth-order tensor
, rank(P) = 1 if and only if all of its many-body η-parameters are factorizable as Since the rank-1 space is described by linear constraints in natural parameters, the rank-1 space is e-flat [29,Chapter 2.4]. Using Propositions 3 and 4, we can derive the projected destination without the gradient method, which reproduces the closed formula of the best rank-1 approximation minimizing KL divergence [18] from the view of information geometry. We also use the expectation conservation low in the m-projection to prove the following proposition (See Appendix D in Supplement for proof). In this case, one-body η-parameters do not change during the m-projection as follows:

Proposition 5 (m-projection onto factorizable subspace) For any positive tensor
, its m-projection onto the rank-1 space is given as Since the m-projection minimizes the KL divergence, it is guaranteed that P obtained by Eq. (18) minimizes the KL divergence from P within the set of rank-1 tensors. If a given tensor is not normalized, we need to divide the right-hand side of Eq. (18) by the (d − 1)-th power sum of all entries of the tensor in order to match the scales of the input and the output tensors, which is consistent with Eq. (11).

Rank-1 approximation as a mean-field approximation
We consider a rank-1 positive tensor P ∈ R I 1 ×···×I d >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, by substituting 0 for all many-body parameters in our model in Eq. (13), we obtain where s (k) ∈ R I k is a positive first-order tensor normalized as I k j k =1 s (k) j k = 1; we can then regard s (k) as a probability distribution with a single random variable j k ∈ [I k ].
The above discussion means that any positive a 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 [19]. Nowadays, it appears in a wide range of domains such as statistics [44], game theory [45,46], and information theory [47]. From the viewpoint of information geometry, [48] developed a theory of mean-field approximation for Boltzmann machines [28], which is defined as p(x) = exp( i b i x i + i j w i j x i x j ) for a binary random variable vector x ∈ {0, 1} n with a bias parameter b = (b) i ∈ R n and an interaction parameter W = (w i j ) ∈ R n×n . To illustrate that a rank-1 approximation can be regarded as a meanfield 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 w i j = 0 for all i and j, which is called a factorizable subspace. Since the distribution with the constraint w i j = 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 us to know the expectation value η i ≡ E[x i ] of an input distributions and requires O(2 n ) computational cost [49], 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 e-projection finds the distribution P e that minimizes D(P e ; P) for a given distribution P and the projection is conduced by solving the mean-field equations η i = σ (b i + j w i j η j ) numerically, where σ (·) is a sigmoid function. There is no theoretical guarantee that the e-projection destination is uniquely determined, since the factorizable subspace is e-flat but not m-flat. The factorizable subspace has a special property such that we can calculate the expectation value η i ≡ E[x i ] from a distribution as η i = tanh −1 b i and can also compute the distribution from the expectation value as b i = 1 2 log 1−η i 1−η i . The analogy of rank-1 approximation 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, m-projection can be performed directly in our formulation with O(I 1 . . . I d ), which is computationally infeasible in the case of Boltzmann machines due to O(2 n ) cost. Moreover, as we prove in Proposition 7 in in the supplementary material, the rank-1 space has the same property as the factorizable subspace of Boltzmann machines; that is, parameters can be easily computed from the dual parameter in a closed form:  rank (8, 5, 2), we firstly define three bingos on mode-2 as shown in a since 8 − 5 is 3, and approximate contiguous blocks filled in green and blue in b by rank-1 tensors using Eq. (11). In the same way, c shows the case where the target rank is (7,5,2). We also define single bingo on mode-1 since 8 − 7 is 1. A subtensor approximated by Eq. (11) is filled in yellow. We assume that we project a tensor onto B (1) , followed by projecting it onto B (2) . After the second m-projection, the θ -parameters on red panels seem to be overwritten. However, these values remain zero after the second m-projection

Bingo rule for general low-Tucker-rank approximation
In this subsection, we extend the above discussion to arbitrary Tucker rank. First, we relax the θ -representation of the rank-1 condition, which was described in Proposition 3. We introduce the bingo rule as a sufficient condition for the tensor to be rank-reduced. Next, we formulate the low-Tucker-rank approximation as a projection onto the subspace that satisfies this bingo rule. Finally, we show that the projection can be achieved by rank-1 approximations for subtensors of input tensor without using the gradient method.
See Supplement for proof. Therefore, for any tensor P ∈ R I 1 ×···×I d >0 such that it has b k bingos on each mode-k, we can always guarantee that its Tucker rank is at most (I 1 − b 1 , . . . , I d − b d ). We define bingo space as a subspace consisting of tensors with bingos. We denote a bingo space by B. The bingo space is always e-flat since a subspace defined by linear constrains in natural parameters is e-flat [29,Chapter 2.4]. For a given bingo space B, the set of indices (i 1 , . . . , i d ) that imposes bingos θ i 1 ,...,i d = 0 is called the bingo-index set and denoted by B .
Finally, we prove that LTR successfully reduces the Tucker rank by extending the above discussion. We formulate low-Tucker-rank approximation as an m-projection onto a specific bingo space. This bingo space is constructed in Step 1 of LTR. Then, in Step 2, we 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, and then discuss 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 1 , . . . , I k−1 , r k , I k+1 , . . . , I d ) for an input positive tensor P ∈ R I 1 ×···×I d >0 and r k < I k . Let B (k) be the set of tensors having I k − r k bingos on mode-k and B (k) be the set of the bingo indices for mode-k constructed in Step 1 of LTR: Note that P ∈ B (k) implies that the Tucker rank of P is at most (I 1 , . . . , I k−1 , r k , I k+1 , . . . , I d ). Let P (k) be the destination of the m-projection from P onto B (k) andθ ,η be its corresponding parameters of θ -and η-coordinates. From the definition of m-projection and the conservation low of η, the parameters of tensor P (k) satisfỹ As described in Sect. 3.3.4, we need η-parameters on only (i 1 , . . . , Therefore, all we have to do to reduce the Tucker rank is to change the elements P i 1 ,··· ,i d for only (i 1 , . . . , i d ) / ∈ˆ B (k) . Such adjustable parts of P can be divided into some contiguous blocks; we call each of these a subtensor of P on mode-k. In Fig. 7a, b, for example, we can find two subtensors P 3 (1) :5 (1) and P 7 (1) :8 (1) . By conducting the rank-1 approximation introduced in Sect. 3.3.3 onto each subtensor, we obtain P (k) satisfying Eq. (21), since Proposition 3 and Eq. (17) hold in these rank-1 approximations. When the rank of only two modes are reduced Let us assume that the target Tucker rank of mode-k is r k < I k and that of mode-l is r l < I l . In this case, we need to consider two bingo spaces B (k) and B (l) associated with bingo-index set B (k) and B (l) . Let P (k) be the resulting tensor of m-projection of P onto B (k) and P (k,l) be the resulting tensor of m-projection of P (k) onto B (l) . To get P (k,l) , let us consider m-projection from P (k) ∈ B (k) to the bingo space B (l) . In this projection, the natural parameters that are set to be 0 in the previous m-projection onto B (k) from P seem to be overwritten (see red panels in Fig. 7c). However, as shown in the Proposition 8 in the supplementary material, after the rank-1 approximation of a tensor where some one-body θ -parameters are already zero, these parameters remain zero. As a result, the θ -and η-parameters of P (k,l) satisfỹ where η is expectation parameter of P. That means P (k,l) is resulting tensor of m- In conclusion, we can obtain the projected tensor onto B by 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 Fig. 8b. For general case Based on the above discussion, we can derive Step 2 for the general case of low-Tucker-rank approximation. We formulate non-negative low-Tucker-rank approximation as a m-projection onto the intersection of bingo spaces on each mode The m-projection destination is given as an iterative application of m-projection d times, starting from P onto subspace B (1) , and from there onto subspace B (2) , . . . , and finally onto subspace B (d) . Note that each mprojection needs only rank-1 approximation for subtensors on each mode. The result of LTR does not depend on the projection order. Since the m-projection minimizes the KL divergence from the input onto the bingo space, LTR always provides the best lowrank approximation in the specified bingo space B; that is, for a given non-negative tensor P, the output T * of LTR satisfies that 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.

Relationship to Legendre decomposition
Our theoretical analysis is closely related to Legendre decomposition [50], which also uses information geometric parameterization of tensors and solves the problem of tensor decomposition by a projection onto a subspace. However, their concept differs from ours in the following regard. 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. Using the conservation law for η-coordinates, we obtain an exact analytical representation of the coordinates of the projection destination without using a gradient method. Figure 8 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. In addition, although we derived the algorithm based on the discussion using (a) (b) Fig. 8 a The relationship among rank-1 approximation, Legendre decomposition [50], and mean-field equation, where we assume that the same bingo space is used in Legendre decomposition. A 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. O is an initial point of Legendre decomposition, which is usually a uniform distribution. P t is a tensor of the t-th step of gradient descent in Legendre Decomposition. b The m-projection of a common space of two different bingo spaces from P can be achieved by m-projection into one bingo space and then m-projecting into the other bingo space the dually flat manifold on the (θ, η)-coordinate, we do not have to know the value of (θ, η) during the algorithm, which also makes a difference from related works in [17,31].

Connection between rank-1 approximation and balancing
So far, we have seen that the rank of tensors can be reduced by describing the lowrank condition with the many-body θ -parameters. Related to this task, it has been reported that characterizing tensors by one-body η-parameters enables an operation called tensor balancing [17], which is often solved by the Sinkhorn algorithm [51], and its quantum information geometry is also studied [52]. Therefore, we summarize the relationship between the rank-1 approximation, which constrains many-body θparameters, and tensor balancing, which constrains one-body η-parameters. First, we introduce tensor balancing for a non-negative tensor P ∈ R I 1 ×···×I d ≥0 . There are two kinds of balancing: fiber balancing and slice balancing. In this section, we discuss only the slice balancing. See the supplementary material for discussion of fiber balancing. Given d vectors c (k) ∈ R I k ≥0 for k ∈ [d], the task of c-slice balancing is to rescale a tensor so that the sum of each slice satisfies (25) for i k ∈ [I k ]. Note that there is no solution when i c i , so we always assume that i c (k) i = 1 for any k ∈ [d] without loss of generality. The information geometric formulation of tensor balancing has already been performed in [17]. Recall the definition of η-parameters, the above condition for c-slice balancing can be expressed by one-body η-parameters of a tensor as η (k) Let us define c-slice balancing space Q c as the set of c-slice balanced tensor, yielding Q c = {p η | η satisfies the condition (26) }. By considering c-slice balancing and rank-1 approximation simultaneously in the framework of information geometry, we can derive the following property: the c-balanced rank-1 tensor always uniquely exists. As discussed in Sect. 2.2, we can identify a distribution using the mixture coordinate system (θ, η) that combines θ -and η-coordinates. On the intersection Q c ∩ B 1 , all one-body parameters are identified by c-balancing condition in Eq. (26) and other parameters are identified by rank-1 condition in Proposition 3, where B 1 is the rank-1 space. Now, balancing conditions and rank-1 condition specify all parameters; therefore, the mixture coordinate (θ, η) uniquely identifies the rank-1 c-balanced tensor. See Supplement for proof.

Theorem 2 The intersection Q c ∩ B 1 is a singleton.
To get the intuition of geometric structure across conditions on rank-1 approximation and balancing, we see a simple case of I 1 = I 2 = 2 and d = 2. We illustrate a simple case of d = 2 as 3D plots in Fig. 9. Let us consider the c-balanced matrix P ∈ R 2×2 ≥0 with n = 2. We obtain θ -and η-coordinate of P ∈ B c using P 22 : where log is element-wise logarithm. Remember that θ 11 corresponds to the normalizing factor and η 11 = 1. The subspace consisting of c-balanced matrices can be drawn as a convex curve in a three-dimensional space by regarding P 22 as a mediator variable. The curve becomes a straight line in the θ -coordinate only when c (1) = c (2) = (0.5, 0.5). In contrast, the set of rank-1 matrices is identified as a plane (θ 21 , θ 12 , 0) in the θ -coordinate since θ 22 = 0 ensures rank(P) = 1 and on the plane (η 21 , η 12 , η 21 η 12 ) in the η space (See Propositions 3 and 4). We can observe that cslice balancing space Q c and rank-1 space B 1 cross a point, which is shown in Fig. 9. It is coherent with Theorem 2. The cross point changes dynamically by c (1) and c (2) . In addition, Fig. 9 shows that we obtain the same matrix by rank-1 approximation of any matrix on Q c since rank-1 approximation does not change the values of one-body η-parameters.

Relation between A1GM and LTR
We summarize the relationship between the two proposed algorithms A1GM and LTR, which are based on the log-linear model on posets and its convex optimization via information geometry. The difference between A1GM and LTR is the structure of posets behind algorithms. After designing proper posets, these two algorithms perform m-projection in common, where some natural parameters become zero. Interestingly, θ -and η-parameters are not computed explicitly during the procedure of both algorithms. This is because the trick of describing the low-rank condition in a dual-flat coordinate system and using a conservation law for the parameters allows us to know the projection destination in a closed form.
It is also known that the constraints of tensor balancing can be described in terms of expectation parameters, as we see in Sect. 3.3.8. In our framework, the task to be solved, such as low-rank approximation or balancing, is described as a constraint in a dual-flat coordinate system. We expect that higher-rank approximations for multiple matrices could also be possible by defining bingos as well as LTR.
In summary, our approach formulates tasks as convex optimizations by taking the input data structure as proper poset and describing the constraints of the task in a dual-flat coordinate system.

Numerical experiments
We empirically examined the efficiency and the effectiveness of the proposed methods A1GM and LTR using synthetic and real-world datasets. See Supplement for implementation and dataset details.

Experiments for A1GM
We used three types of data to empirically investigate the efficiency and effectiveness of A1GM: (i) synthetic data with missing values at the upper right corner, (ii) synthetic data with random grid-like missing values, and (iii) real data with grid-like and nongrid-like missing values. It is guaranteed that A1GM always finds the best solution for any data of (i) and (ii). Thus, we only investigate efficiency in our experiments for (i) and (ii). By contrast, for data (iii), the reconstruction error can be worse than the existing methods due to increased missing values in A1GM. Therefore, in the experiment for data (iii), we investigate both efficiency (running time) and effectiveness (reconstruction error). We use KL-WNMF as a comparison method [53]. KL-WNMF is a commonly used gradient method that reduces the KL-based cost in Eq. (7) by multiplicative updates. Although faster NMF methods, such as ALS [54] and ADMM [55], have been developed, they are just as fast as a multiplicative update when the target rank is small [56]. Moreover, as we will show in Sects. 4.1.1 and 4.1.2, KL-WNMF converges within only 2-4 iterations in our experiments. Thus, these techniques are considered ineffective for speeding up rank-1 KL-WNMF. This is why we only compared A1GM with simple KL-WNMF. We implemented KL-WNMF by referring to the original paper [53]. The stopping criterion of KL-WNMF follows the implementation of the standard NMF in scikit-learn [57]. The initial values of KL-WNMF are determined by sampling from a uniform continuous distribution from 0 to 1.

Synthetic datasets
Missing values in the top right corner We prepared synthetic matrices X ∈ R N ×N and their weights ∈ {0, 1} N ×N . We assumed that each input weight is in the form of Eq. (8) with L = 0. We measured the running time to obtain rank-1 decomposition of X with varying the matrix size N . Figure 10a shows that A1GM is an order of magnitude faster than the existing gradient method. The number of iterations of the existing method until convergence was between 2 and 4. A1GM just applies the closed formula in Theorem 1 to parts of input matrices.

Random grid-like missing values
We also prepared synthetic matrices and its binary weight matrices ∈ {0, 1} N ×N . We assumed that every input weight matrix is gridlike, and we set the ratio of missing values to be 5 percent. We measured the running time of A1GM to complete the best rank-1 missing NMF compared with KL-WNMF by varying the matrix size N . Figure 10b shows that our method is always faster than the gradient method. The number of iterations of the existing method required for convergence was between 2 and 4. Note that, in these datasets, A1GM does not need to increase missing values.

Real datasets
We used 20 real datasets. We downloaded tabular datasets that have missing values from the Kaggle databank 3 or UCI dataset. 4 If a dataset contains negative values, we converted them to their absolute values. Zero values in a matrix were replaced with the average value of the matrix to make them all positive. We evaluated the relative error as D (X, A1GM(X)) D (X, WNMF(X)), where WNMF(X) and A1GM(X) are the rank-1 reconstructed matrices by KL-WNMF and A1GM, respectively, and the binary weight matrix indicates the locations of the missing values of X. We also compared the relative running time of A1GM to KL-WNMF.
The results are summarized in Table 1. In the table, the column increase rate means the ratio of the number of missing values after addition in A1GM to the original number of missing values. If increase rate is 1, it means that the location of missing values of the dataset is originally grid-like. For such datasets, it is theoretically guaranteed that our method A1GM always provides the best rank-1 missing NMF, which minimizes the KL divergence in Eq. (7). It is reasonable that the matrices reconstructed by KL-WNMF and by A1GM are the same since the cost function (7) is convex in w and h. The number of iterations of the existing method required for convergence was between 2 and 4 for real datasets.
We can see that A1GM is much faster than KL-WNMF for all the datasets. Moreover, the relative error remains low even if missing values of datasets are not grid-like for most of datasets. In some real data, missing values are likely to be biased towards a particular row or column. As a result, they become grid-like by just adding a small number of missing values. In these cases, our proposed method can conduct rank-1 missing NMF rapidly with competitive errors to KL-WNMF. By contrast, a large amount of information is lost after the increasing missing value step for some datasets (large increase rates). As a result, our method is not suitable for obtaining an accurately reconstructed rank-1 matrix, even though it is much faster than the existing method.

Experiments for LTR
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 non-negative tensor decomposition method [58] whose cost function is either the Least Squares (LS) error (NTD_LS) or the KL divergence (NTD_KL). The second method is sequential non-negative Tucker decomposition (lraSNTD), which is known as the faster of the two methods [59]. Its cost function is the LS error.
Results on synthetic data We created tensors with d = 3 or d = 5, where every I k = n. We changed n to generate various sizes of tensors. Each element was sampled from the uniform continuous distribution from 0 to 1. To evaluate the efficiency, we measured the running time of each method. To evaluate the accuracy, we measured the  The horizontal axis is n 3 for input (n, n, n) tensor. c, d The horizontal axis is the number of elements of the core tensor LS reconstruction error, defined as the Frobenius norm between input and output tensors. Figure 11a shows the running time and the LS reconstruction error for randomly generated tensors with d = 3 and n = 30 varying the target Tucker tensor rank. Figure 11b 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 but still retains the competitive approximation accuracy.
As described in Sect. 3.3.6, 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.

Conclusion
We have discussed non-negative low-rank approximations from an information geometric viewpoint. Specifically, using a log-linear model with a partial order structure in the sample space, we have treated a tensor as a probability distribution and successfully characterized the rank of the tensor in a dually flat coordinate system. We have pointed out that rank-1 approximation of a non-negative tensor can be geometrically captured as a mean-field approximation that reduces the many-body problem to the one-body problem, which is frequently used in statistical physics. We also have constructed a new algorithm LTR that can perform low-Tucker-rank approximation by appropriately designing constraints for the dually flat coordinate system.
We have also analyzed the problem of non-negative multiple matrix factorization (NMMF) in the same way and found the closed formula of globally optimal rank-1 NMMF. By focusing on the equivalence of the NMMF and the decomposition of matrices with missing values, we have constructed a new algorithm A1GM that can quickly compute approximate solutions to rank-1 approximations for matrices with missing values without using the gradient method. A1GM can extract the largest principal components from matrices with missing values.
Our work will form the basis for further investigation between linear algebraic matrix operations, statistics, and machine learning via information geometry.