Tensor Balancing on Statistical Manifold
Mahito Sugiyama Hiroyuki Nakahara Koji Tsuda
=1

roundnatbib




       
            
       
       
       
      


plainnat















mydefinition
  5.5pt 
  0pt 
   
   
   
  . 
  .5em 
   
mytheorem
  5.5pt 
  0pt 
   
   
   
  . 
  .5em 
   
mydefinition
definitionDefinition
exampleExample
mytheorem
theoremTheorem
propositionProposition
lemmaLemma


itemize
•  
 
 
 
 
 
 
 
 
 


enum2
enumerate

 enum2.  
 enum2
 
 
 
 
 
 
 
 
 
 






































Tensor Balancing on Statistical Manifold
    
  Mahito Sugiyama 

  National Institute of Informatics 

  JST, PRESTO 

   

  Hiroyuki Nakahara 

  RIKEN Brain Science Institute 

   

  Koji Tsuda 

  The University of Tokyo 

  RIKEN AIP; NIMS 

   


    December 30, 2023
=====================================================================================================================================================================================================================




 We solve tensor balancing, rescaling an Nth order nonnegative tensor by multiplying N tensors of order N - 1 so that every fiber sums to one.
 This generalizes a fundamental process of matrix balancing used to compare matrices in a wide range of applications from biology to economics.
 We present an efficient balancing algorithm with quadratic convergence using Newton's method and show in numerical experiments that the proposed algorithm is several orders of magnitude faster than existing ones.
 To theoretically prove the correctness of the algorithm, we model tensors as probability distributions in a statistical manifold and realize tensor balancing as projection onto a submanifold.
 The key to our algorithm is that the gradient of the manifold, used as a Jacobian matrix in Newton's method, can be analytically obtained using the Möbius inversion formula, the essential of combinatorial mathematics.
 Our model is not limited to tensor balancing, but has a wide applicability as it includes various statistical and machine learning models such as weighted DAGs and Boltzmann machines.




§ INTRODUCTION

Matrix balancing is the problem of rescaling a given square nonnegative matrix A ∈^n × n_≥ 0 to a doubly stochastic matrix RAS, where every row and column sums to one, by multiplying two diagonal matrices R and S.
This is a fundamental process for analyzing and comparing matrices in a wide range of applications, including input-output analysis in economics, called the RAS approach <cit.>, seat assignments in elections <cit.>, Hi-C data analysis <cit.>, the Sudoku puzzle <cit.>, and the optimal transportation problem <cit.>.
An excellent review of this theory and its applications is given by <cit.>.

The standard matrix balancing algorithm is the Sinkhorn-Knopp algorithm <cit.>, a special case of Bregman's balancing method <cit.> that iterates rescaling of each row and column until convergence.
The algorithm is widely used in the above applications due to its simple implementation and theoretically guaranteed convergence.
However, the algorithm converges linearly <cit.>, which is prohibitively slow for recently emerging large and sparse matrices.
Although <cit.> and <cit.> tried to achieve faster convergence by approximating each step of Newton's method,
the exact Newton's method with quadratic convergence has not been intensively studied yet.

Another open problem is tensor balancing, which is a generalization of balancing from matrices to higher-order multidimentional arrays, or tensors.
The task is to rescale an Nth order nonnegative tensor to a multistochastic tensor, in which every fiber sums to one, by multiplying (N - 1)th order N tensors.
There are some results about mathematical properties of multistochastic tensors <cit.>.
However, there is no result for tensor balancing algorithms with guaranteed convergence that transforms a given tensor to a multistochastic tensor until now.

Here we show that Newton's method with quadratic convergence can be applied to tensor balancing while avoiding solving a linear system on the full tensor.
Our strategy is to realize matrix and tensor balancing as projection onto a dually flat Riemmanian submanifold (Figure <ref>), which is a statistical manifold and known to be the essential structure for probability distributions in information geometry <cit.>.
Using a partially ordered outcome space, we generalize the log-linear model <cit.> used to model the higher-order combinations of binary variables <cit.>, which allows us to model tensors as probability distributions in the statistical manifold.
The remarkable property of our model is that the gradient of the manifold can be analytically computed using the Möbius inversion formula <cit.>, the heart of combinatorial mathematics <cit.>, which enables us to directly obtain the Jacobian matrix in Newton's method.
Moreover, we show that (n - 1)^N entries for the size n^N of a tensor are invariant with respect to one of the two coordinate systems of the statistical manifold.
Thus the number of equations in Newton's method is O(n^N - 1).



The remainder of this paper is organized as follows:
We begin with a low-level description of our matrix balancing algorithm in Section <ref> and demonstrate its efficiency in numerical experiments in Section <ref>.
To guarantee the correctness of the algorithm and extend it to tensor balancing, we provide theoretical analysis in Section <ref>.
In Section <ref>, we introduce a generalized log-linear model associated with a partial order structured outcome space, followed by introducing the dually flat Riemannian structure in Section <ref>.
In Section <ref>, we show how to use Newton's method to compute projection of a probability distribution onto a submanifold.
Finally, we formulate the matrix and tensor balancing problem in Section <ref> and summarize our contributions in Section <ref>.



§ THE MATRIX BALANCING ALGORITHM

Given a nonnegative square matrix A = (a_ij) ∈^n × n_≥ 0,
the task of matrix balancing is to find r⃗, s⃗∈^n that satisfy

    (RAS)1⃗ = 1⃗,  
     (RAS)^T1⃗ = 1⃗,

where R = (r⃗) and S = (s⃗).
The balanced matrix A' = RAS is called doubly stochastic, in which each entry a'_ij = a_ij r_i s_j and all the rows and columns sum to one.
The most popular algorithm is the Sinkhorn-Knopp algorithm, which repeats updating r⃗ and s⃗ as r⃗ = 1 / (As⃗) and s⃗ = 1 / (A^T r⃗).
We denote by [n] = {1, 2, …, n} hereafter.

In our algorithm, instead of directly updating r⃗ and s⃗, we update two parameters θ and η defined as

    log p_ij = ∑_i' ≤ i∑_j' ≤ jθ_i'j',  η_ij = ∑_i' ≥ i∑_j' ≥ jp_i'j'

for each i, j ∈ [n], where we normalized entries as p_ij = a_ij / ∑_ij a_ij so that ∑_ij p_ij = 1.
We assume for simplicity that each entry is strictly larger than zero. The assumption will be removed in Section <ref>.



The key to our approach is that we update θ_ij^(t) with i = 1 or j = 1 by Newton's method at each iteration t = 1, 2, … while fixing θ_ij with i,j ≠ 1 so that η_ij^(t) satisfies the following condition (Figure <ref>):

    η_i1^(t) = n - i + 1/n,  η_1j^(t) = n - j + 1/n.

Note that the rows and columns sum not to 1 but to 1 / n due to the normalization.
The update formula is described as

    [ θ_11^(t + 1); θ_12^(t + 1);            ⋮; θ_1n^(t + 1); θ_21^(t + 1);            ⋮; θ_n1^(t + 1) ]
     =
     [ θ_11^(t); θ_12^(t);        ⋮; θ_1n^(t); θ_21^(t);        ⋮; θ_n1^(t) ]
     - J^-1[ η_11^(t) - (n - 1 + 1) / n; η_12^(t) - (n - 2 + 1) / n;                          ⋮; η_1n^(t) - (n - n + 1) / n; η_21^(t) - (n - 2 + 1) / n;                          ⋮; η_n1^(t) - (n - n + 1) / n ],

where J is the Jacobian matrix given as

    J_(ij)(i'j') = ∂η^(t)_ij/∂θ^(t)_i'j' = η_max{i, i'}max{j, j'} - n^2 η_ijη_i'j',

which is derived from our theoretical result in Theorem <ref>.
Since J is a (2n - 1) × (2n - 1) matrix, the time complexity of each update is O(n^3), which is needed to compute the inverse of J.

After updating to θ_ij^(t + 1), we can compute p_ij^(t + 1) and η_ij^(t + 1) by Equation (<ref>).
Since this update does not ensure the condition ∑_ij p_ij^(t + 1) = 1, we again update θ_11^(t + 1) as

    θ_11^(t + 1) = θ_11^(t + 1) - log∑_ijp_ij^(t + 1)

and recompute p_ij^(t + 1) and η_ij^(t + 1) for each i, j ∈ [n].

By iterating the above update process in Equation (<ref>) until convergence, A = (a_ij) with a_ij = n p_ij becomes doubly stochastic.







§ NUMERICAL EXPERIMENTS

We evaluate the efficiency of our algorithm compared to the two prominent balancing methods, the standard Sinkhorn-Knopp algorithm <cit.> and the state-of-the-art algorithm BNEWT <cit.>, which uses Newton's method-like iterations with conjugate gradients.
All experiments were conducted on Amazon Linux AMI release 2016.09 with a single core of 2.3 GHz Intel Xeon CPU E5-2686 v4 and 256 GB of memory.
All methods were implemented in  with the  library and compiled with gcc 4.8.3[An implementation of algorithms for matrices and third order tensors is available at: <https://github.com/mahito-sugiyama/newton-balancing>].
We have carefully implemented BNEWT by directly translating the MATLAB code provided in <cit.> into  with the  library for fair comparison, and
used the default parameters.
We measured the residual of a matrix A' = (a_ij') by the squared norm (A'1⃗ - 1⃗, A'^T1⃗ - 1⃗)_2, where each entry a'_ij is obtained as n p_ij in our algorithm, and ran each of three algorithms until the residual is below the tolerance threshold 10^-6.


Hessenberg Matrix.
The first set of experiments used a Hessenberg matrix, which has been a standard benchmark for matrix balancing <cit.>.
Each entry of an n × n Hessenberg matrix H_n = (h_ij) is given as h_ij = 0 if j < i - 1 and h_ij = 1 otherwise.
We varied the size n from 10 to 5,000, and measured running time (in seconds) and the number of iterations of each method.

Results are plotted in Figure <ref>.
Our balancing algorithm with the Newton's method (plotted in blue in the figures) is clearly the fastest:
It is three to five orders of magnitude faster than the standard Sinkhorn-Knopp algorithm (plotted in red).
Although the BNEWT algorithm (plotted in green) is competitive if n is small, it suddenly fails to converge whenever n ≥ 200, which is consistent with results in the original paper <cit.> where there is no result for the setting n ≥ 200 on the same matrix.
Moreover, our method converges around 10 to 20 steps, which is about three and seven orders of magnitude smaller than BNEWT and Sinkhorn-Knopp, respectively, at n = 100.



To see the behavior of the rate of convergence in detail, we plot the convergence graph in Figure <ref> for n = 20, where we observe the slow convergence rate of the Sinkhorn-Knopp algorithm and unstable convergence of the BNEWT algorithm, which contrasts with our quick convergence.


Trefethen Matrix.
Next, we collected a set of Trefethen matrices from a collection website[<http://www.cise.ufl.edu/research/sparse/matrices/>], which are nonnegative diagonal matrices with primes.
Results are plotted in Figure <ref>, where we observe the same trend as before:
Our algorithm is the fastest and about four orders of magnitude faster than the Sinkhorn-Knopp algorithm.
Note that larger matrices with n > 300 do not have total support, which is the necessary condition for matrix balancing <cit.>, while the BNEWT algorithm fails to converge if n = 200 or n = 300.



§ THEORETICAL ANALYSIS

In the following, we provide theoretical support to our algorithm by formulating the problem as a projection within a statistical manifold, in which a matrix corresponds to an element, that is, a probability distribution, in the manifold.


We show that a balanced matrix forms a submanifold and matrix balancing is projection of a given distribution onto the submanifold, where
the Jacobian matrix in Equation (<ref>) is derived from the gradient of the manifold.



 §.§ Formulation

We introduce our log-linear probabilistic model, where the outcome space is a partially ordered set, or a poset <cit.>.
We prepare basic notations and the key mathematical tool for posets, the Möbius inversion formula, followed by formulating the log-linear model.



  §.§.§ Möbius Inversion

A poset (S, ≤), the set of elements S and a partial order ≤ on S, is a fundamental structured space in computer science.
A partial order “≤” is a relation between elements in S that satisfies the following three properties: For all x, y, z ∈ S,
(1) x ≤ x (reflexivity), (2) x ≤ y, y ≤ x ⇒ x = y (antisymmetry), and (3) x ≤ y, y ≤ z ⇒ x ≤ z (transitivity).
In what follows, S is always finite and includes the least element (bottom) ∈ S; that is, ≤ x for all x ∈ S.
We denote S ∖{} by S^+.

<cit.> introduced the Möbius inversion formula on posets by generalizing the inclusion-exclusion principle.
Let ζ S × S →{0, 1} be the zeta function defined as

    ζ(s, x) = {[          1 if  s ≤ x,;          0 otherwise. ].

The Möbius function μ S × S → satisfies ζμ = I, which is inductively defined for all x, y with x ≤ y as
 
    μ(x, y) = {[                   1          if  x = y,; -∑_x ≤ s < yμ(x, s)          if  x < y,;                   0          otherwise. ].

From the definition, it follows that

    ∑_s ∈ Sζ(s, y)μ(x, s)    = ∑_x ≤ s ≤ yμ(x, s) = δ_xy,
    ∑_s ∈ Sζ(x, s)μ(s, y)    = ∑_x ≤ s ≤ yμ(s, y) = δ_xy

with the Kronecker delta δ such that δ_xy = 1 if x = y and δ_xy = 0 otherwise.
Then for any functions f, g, and h with the domain S such that

    g(x)    = ∑_s ∈ Sζ(s, x) f(s) = ∑_s ≤ x f(s),
    
     h(x)    = ∑_s ∈ Sζ(x, s) f(s) = ∑_s ≥ x f(s),

f is uniquely recovered with the Möbius function:

    f(x) = ∑_s ∈ Sμ(s, x) g(s),  
     f(x) = ∑_s ∈ Sμ(x, s) h(s).

This is called the Möbius inversion formula and is at the heart of enumerative combinatorics <cit.>.



  §.§.§ Log-Linear Model on Posets

We consider a probability vector p on (S, ≤) that gives a discrete probability distribution with the outcome space S.
A probability vector is treated as a mapping p S→ (0, 1) such that ∑_x ∈ S p(x)= 1, where every entry p(x) is assumed to be strictly larger than zero.

Using the zeta and the Möbius functions, let us introduce two mappings θ S→ and η S → as

    θ(x)    = ∑_s ∈ Sμ(s, x) log p(s),
    η(x)    = ∑_s ∈ Sζ(x, s) p(s) = ∑_s ≥ x p(s).

From the Möbius inversion formula, we have

    log p(x)    = ∑_s ∈ Sζ(s, x) θ(s) = ∑_s ≤ xθ(s),
    
     p(x)    = ∑_s ∈ Sμ(x, s) η(s).

They are generalization of the log-linear model <cit.> that gives the probability p(x⃗) of an n-dimensional binary vector x⃗ = (x^1, …, x^n) ∈{0, 1}^n as

    log p(x⃗) = ∑_iθ^i x^i + ∑_i < jθ^ij x^i x^j + ∑_i < j < kθ^ijk x^i x^j x^k + … + θ^1… n x^1 x^2 … x^n - ψ,

where θ⃗ = (θ^1, …, θ^12… n) is a parameter vector, ψ is a normalizer, and η⃗ = (η^1, …, η^12… n) represents the expectation of variable combinations such that

    η^i   = [x^i] = (x^i = 1),
    η^ij   = [x^ix^j] = (x^i = x^j = 1), i < j, …
    η^1… n   = [x^1 … x^n] = (x^1 = … = x^n = 1).

They coincide with Equations (<ref>) and (<ref>) when we let S = 2^V with V = {1, 2, …, n}, each x ∈ S as the set of indices of “1” of x⃗, and the order ≤ as the inclusion relationship, that is, x ≤ y if and only if x ⊆ y.
<cit.> have pointed out that θ⃗ can be computed from p using the inclusion-exclusion principle in the log-linear model.
We exploit this combinatorial property of the log-linear model using the Möbius inversion formula on posets and extend the log-linear model from the power set 2^V to any kind of posets (S, ≤).
<cit.> studied a relevant log-linear model, but the relationship with Möbius inversion formula has not been analyzed yet.




 §.§ Dually Flat Riemannian Manifold

We theoretically analyze our log-linear model introduced in Equations (<ref>), (<ref>) and show that they form dual coordinate systems on a dually flat manifold, which has been mainly studied in the area of information geometry <cit.>.
Moreover, we show that the Riemannian metric and connection of our model can be analytically computed in closed forms.

In the following, we denote by ξ the function θ or η and by ∇ the gradient operator with respect to S^+ = S ∖{}, i.e., (∇ f(ξ))(x) = ∂ f / ∂ξ(x) for x ∈ S^+, and denote by  the set of probability distributions specified by probability vectors, which forms a statistical manifold.
We use uppercase letters P, Q, R, … for points (distributions) in  and their lowercase letters p, q, r, … for the corresponding probability vectors treated as mappings.
We write θ_P and η_P if they are connected with p by Equations (<ref>) and (<ref>), respectively, and abbreviate subscripts if there is no ambiguity.



  §.§.§ Dually Flat Structure

We show that  has the dually flat Riemannian structure induced by two functions θ and η in Equation (<ref>) and (<ref>).
We define ψ(θ) as

    ψ(θ) = -θ() = -log p(),

which corresponds to the normalizer of p.
It is a convex function since we have

    ψ(θ) = log∑_x ∈ Sexp( ∑_ < s ≤ xθ(s) )

from log p(x) = ∑_ < s ≤ xθ(s) - ψ(θ).
We apply the Legendre transformation to ψ(θ) given as

    ϕ(η) = max_θ'(θ'η - ψ(θ')),  θ' η = ∑_x ∈ S^+θ'(x) η(x).

Then ϕ(η) coincides with the negative entropy.

 
 
    ϕ(η) = ∑_x ∈ S p(x) log p(x).



 From Equation (<ref>), we have
 
    θ' η
      = ∑_x ∈ S^+( ∑_ < s ≤ xμ(s, x) log p'(s) ∑_s ≥ x p(s) )
      = ∑_x ∈ S^+ p(x)( log p'(x) - log p'() ).

 Thus it holds that
 
    θ'η - ψ(θ')
         = ∑_x ∈ S p(x) log p'(x).

 Hence it is maximized with p(x) = p'(x).

Since they are connected with each other by the Legendre transformation, they form a dual coordinate system ∇ψ(θ) and ∇ϕ(η) of  <cit.>, which coincides with θ and η as follows.

 
    ∇ψ(θ) = η,  ∇ϕ(η) = θ.



 They can be directly derived from our definitions (Equations (<ref>) and (<ref>)) as

    ∂ψ(θ)/∂θ(x)   = ∑_y ≥ xexp(∑_ < s ≤ yθ(s))/∑_y ∈ Sexp(∑_ < s ≤ yθ(s))
     = ∑_s ≥ x p(s) = η(x),
    ∂ϕ(η)/∂η(x)   = ∂/∂η(x)(θη - ψ(θ)) = θ(x).


Moreover, we can confirm the orthogonality of θ and η as

    [∂/∂θ(x)log p(s) ∂/∂η(y)log p(s)]
        = ∑_s ∈ S[ p(s) ∂/∂θ(x)∑_u ∈ Sζ(u, s)θ(u) ∂/∂η(y)log(∑_u ∈ Sμ(s, u)η(u)) ]
       = ∑_s ∈ S[ p(s) (ζ(x, s) - η(x)) μ(s, y)/p(s) ]
       = ∑_s ∈ Sζ(x, s)μ(s, y) = δ_xy.




The last equation holds from Equation (<ref>), hence the Möbius inversion directly leads to the orthogonality.

The Bregman divergence is known to be the canonical divergence <cit.> to measure the difference between two distributions P and Q on a dually flat manifold, which is defined as

    D[P, Q] = ψ(θ_P) + ϕ(η_Q) - θ_P η_Q.

In our case, since we have
ϕ(η_Q) = ∑_x ∈ S q(x) log q(x) and
θ_Pη_Q - ψ(θ_P) = ∑_x ∈ S q(x) log p(x)
from Theorem <ref> and Equation (<ref>),
it is given as

    D[P, Q]
     = ∑_x ∈ S q(x) logq(x)/p(x),

which coincides with the Kullback–Leibler divergence (KL divergence) from Q to P: D[P, Q] = [Q, P].



  §.§.§ Riemannian Structure

Next we analyze the Riemannian structure on  and show that the Möbius inversion formula enables us to compute the Riemannian metric of .

 
 The manifold (, g(ξ)) is a Riemannian manifold with the Riemannian metric g(ξ) such that for all x, y ∈ S^+
 
    g_xy(ξ) = {[ ∑_s ∈ Sζ(x, s)ζ(y, s)p(s) - η(x)η(y)                            if ξ = θ,;        ∑_s ∈ Sμ(s, x)μ(s, y) p(s)^-1                            if ξ = η. ].



 Since the Riemannian metric is defined as

    g(θ) = ∇∇ψ(θ),   g(η) = ∇∇ϕ(η),

when ξ = θ we have

    g_xy(θ)    = ∂^2/∂θ(x) ∂θ(y)ψ(θ) = ∂/∂θ(x)η(y)
     = ∂/∂θ(x)∑_s ∈ Sζ(y, s)exp(∑_ < u ≤ sθ(u) - ψ(θ))
       = ∑_s ∈ Sζ(x, s)ζ(y, s)p(s) - η(x)η(y).

When ξ = η, it follows that

    g_x, y(η)    = ∂^2/∂η(x) ∂η(y)ϕ(η) = ∂/∂η(x)θ(y)
       = ∂/∂η(x)∑_s ≤ yμ(s, y) log p(s)
     = ∂/∂η(x)∑_s ≤ yμ(s, y) log(∑_u ≥ sμ(s, u) η(u))
       = ∑_s ∈ Sμ(s, x)μ(s, y)/∑_u ≥ sμ(s, u)η(u)
     = ∑_s ∈ Sμ(s, x)μ(s, y) p(s)^-1.


Since g(ξ) coincides with the Fisher information matrix,

    [ ∂/∂θ(x)log p(s) ∂/∂θ(y)log p(s) ]    = g_xy(θ),
    [ ∂/∂η(x)log p(s) ∂/∂η(y)log p(s) ]    = g_xy(η).

Then the Riemannian (Levi–Chivita) connection Γ(ξ) with respect to ξ, which is defined as

    Γ_xyz(ξ) = 1/2( ∂ g_yz(ξ)/∂ξ(x) + ∂ g_xz(ξ)/∂ξ(y) - ∂ g_xy(ξ)/∂ξ(z) )

for all x, y, z ∈ S^+, can be analytically obtained.

 The Riemannian connection Γ(ξ) on the manifold (, g(ξ)) is given in the following for all x, y, z ∈ S^+,
 
    Γ_xyz(ξ) = {[ -1/2∑_s ∈ S(ζ(x, s) - η(x))(ζ(y, s) - η(y))(ζ(z, s) - η(z)) p(s)                                                        if ξ = θ,;                         -1/2∑_s ∈ Sμ(s, x)μ(s, y)μ(s, z) p(s)^-2                                                        if ξ = η. ].



 We have for all x, y, z ∈ S,
 
    ∂ g_y,z(θ)/∂θ(x)
      = ∂/∂θ(x)∑_s ∈ Sζ(y, s)ζ(z, s)p(s) - ∂/∂θ(x)η(y)η(z),

 where
 
    ∂/∂θ(x)∑_s ∈ Sζ(y, s)ζ(z, s)p(s)
         = ∂/∂θ(x)∑_s ∈ Sζ(x, s)ζ(y, s)ζ(z, s) exp(∑_ < u ≤ sθ(u) - ψ(θ))
       = ∑_s ∈ Sζ(x, s)ζ(y, s)ζ(z, s)p(s) - η(x)∑_s ∈ Sζ(y, s)ζ(z, s)p(s)

 and
 
    ∂/∂θ(x)η(y)η(z)
         = ∂η(y)/∂θ(x)η(z) + ∂η(z)/∂θ(x)η(y)
       = η(z)∑_s ∈ Sζ(x, s)ζ(y, s)p(s) + η(y)∑_s ∈ Sζ(x, s)ζ(z, s)p(s) - 2η(x)η(y)η(z).

 It follows that
 
    ∂ g_y,z(θ)/∂θ(x)
      = ∑_s ∈ S(ζ(x, s) - η(x))(ζ(y, s) - η(y))(ζ(z, s) - η(z)) p(s).

 On the other hand,
 
    ∂ g_y,z(η)/∂η(x)   = ∂/∂η(x)∑_s ∈ Sμ(s, y)μ(s, z)p(s)^-1
      = ∂/∂η(x)∑_s ∈ Sμ(s, y)μ(s, z) (∑_u ≥ sμ(s, u) η(s))^-1
       = -∑_s ∈ Sμ(s, x)μ(s, y)μ(s, z) (∑_u ≥ sμ(s, u) η(s))^-2
      = -∑_s ∈ Sμ(s, x)μ(s, y)μ(s, z) p(s)^-2.

Therefore, from the definition of Γ(ξ), it follows that
 
    Γ_x,y,z(θ)
         = 1/2( ∂ g_y,z(θ)/∂θ(x) + ∂ g_x,z(θ)/∂θ(y) - ∂ g_x,y(θ)/∂θ(z) )
       = 1/2∑_s ∈ S(ζ(s, x) - η(x))(ζ(s, y) - η(y))(ζ(s, z) - η(z)) p(s),
    Γ_x,y,z(η)
         = 1/2( ∂ g_y,z(η)/∂η(x) + ∂ g_x,z(η)/∂η(y) - ∂ g_x,y(η)/∂η(z) )
       = -1/2∑_s ∈ Sμ(s, x)μ(s, y)μ(s, z) p(s)^-2.

 







 §.§ The Projection Algorithm

Projection of a distribution onto a submanifold is essential; several machine learning algorithms are known to be formulated as projection of a distribution empirically estimated from data onto a submanifold that is specified by the target model <cit.>.
Here we define projection of distributions on posets and show that Newton's method can be applied to perform projection as the Jacobian matrix can be analytically computed.



  §.§.§ Definition

Let (β) be a submanifold of  such that

    (β) = P ∈ | θ_P(x) = β(x)  for all  x ∈(β)

specified by a function β with (β) ⊆ S^+.
Projection of P ∈ onto (β), called m-projection, which is defined as the distribution P_β∈(β) such that

    {[   θ_P_β(x) = β(x)       if  x ∈(β),; η_P_β(x) = η_P(x) if  x ∈ S^+ ∖(β), ].

is the minimizer of the KL divergence from P to (β):

    P_β   = _Q ∈(β)[P, Q].

The dually flat structure with the coordinate systems θ and η guarantees that the projected distribution P_β always exists and is unique <cit.>.
Moreover, the Pythagorean theorem holds in the dually flat manifold, that is, for any Q ∈(β) we have

    [P, Q] = [P, P_β] + [P_β, Q].

We can switch η and θ in the submanifold (β) by changing [P, Q] to [Q, P], where the projected distribution P_β of P is given as

    {[ θ_P_β(x) = θ_P(x) if  x ∈ S^+ ∖(β),;   η_P_β(x) = β(x)       if  x ∈(β), ].

This projection is called e-projection.

[Boltzmann machine]
 Given a Boltzmann machine represented as an undirected graph G = (V, E) with a vertex set V and an edge set E ⊆{{i, j}| i, j ∈ V}.
 The set of probability distributions that can be modeled by a Boltzmann machine G coincides with the submanifold
 
    _B = P ∈ | θ_P(x) = 0  if  |x| > 2  or  x ∉E,

 with S = 2^V.
 Let P̂ be an empirical distribution estimated from a given dataset.
 The learned model is the m-projection of the empirical distribution P̂ onto _B, where the resulting distribution P_β is given as
 
    {[            θ_P_β(x) = 0  if  |x| > 2  or  x ∉E,;      η_P_β(x) = η_P̂(x) if  |x| = 1  or  x ∈ E. ].





  §.§.§ Computation

Here we show how to compute projection of a given probability distribution.
We show that Newton's method can be used to efficiently compute the projected distribution P_β by iteratively updating P_β^(0) = P as P_β^(0), P_β^(1), P_β^(2), … until converging to P_β.

Let us start with the m-projection with initializing P_β^(0) = P.
In each iteration t, we update θ_P_β^(t)(x) for all x ∈β while fixing η_P_β^(t)(x) = η_P(x) for all x ∈ S^+ ∖(β), which is possible from the orthogonality of θ and η.
Using Newton's method, η_P_β^(t + 1)(x) should satisfy

    (θ_P_β^(t)(x) - β(x)) + ∑_y ∈(β) J_xy(η_P_β^(t + 1)(y) - η_P_β^(t)(y)) = 0,

for every x ∈(β), where J_xy is an entry of the |(β)| × |(β)| Jacobian matrix J and given as

    J_xy = ∂θ_P_β^(t)(x)/∂η_P_β^(t)(y) = ∑_s ∈ Sμ(s, x)μ(s, y) p_β^(t)(s)^-1

from Theorem <ref>.
Therefore, we have the update formula for all x ∈(β) as

    η_P_β^(t + 1)(x)    = η_P_β^(t)(x) - ∑_y ∈(β) J_xy^-1(θ_P_β^(t)(y) - β(y)).

In e-projection, update η_P_β^(t)(x) for x ∈(β) while fixing θ_P_β^(t)(x) = θ_P(x) for all x ∈ S^+ ∖(β).
To ensure η_P_β^(t)() = 1, we add  to (β) and β() = 1.
We update θ_P_β^(t)(x) at each step t as

    θ_P_β^(t + 1)(x)    = θ_P_β^(t)(x) - ∑_y ∈(β)J'_xy^-1(η_P_β^(t)(y) - β(y)),
    J'_xy   = ∂η_P_β^(t)(x)/∂θ_P_β^(t)(y) = ∑_s ∈ Sζ(x, s)ζ(y, s)p_β^(t)(s) - |S|η_P_β^(t)(x)η_P_β^(t)(y).

In this case, we also need to update θ_P_β^(t)() as it is not guaranteed to be fixed.
Let us define

    p_β'^(t + 1)(x)    = p_β^(t)(x)∏_s ∈(β)exp(θ_P_β^(t + 1)(s))/exp(θ_P_β^(t)(s)) ζ(s, x).

Since we have

    p_β^(t + 1)(x)    = exp(θ_P_β^(t + 1)())/exp(θ_P_β^(t)()) p_β'^(t + 1)(x),

it follows that

    θ_P_β^(t + 1)() - θ_P_β^(t)() = - log(exp(θ_P_β^(t)()) +∑_x ∈ S^+p_β'^(t + 1)(x)),

The time complexity of each iteration is O(|(β)|^3), which is required to compute the inverse of the Jacobian matrix.

Global convergence of the projection algorithm is always guaranteed by the convexity of a submanifold (β) defined in Equation (<ref>).
Since (β) is always convex with respect to the θ- and η-coordinates, it is straightforward to see that our e-projection is an instance of the Bregman algorithm onto a convex region, which is well known to always converge to the global solution <cit.>.



§ BALANCING MATRICES AND TENSORS

Now we are ready to solve the problem of matrix and tensor balancing as projection on a dually flat manifold.



 §.§ Matrix Balancing

Recall that the task of matrix balancing is to find r⃗, s⃗∈^n that satisfy (RAS)1⃗ = 1⃗ and (RAS)^T1⃗ = 1⃗ with R = (r⃗) and S = (s⃗)
for a given nonnegative square matrix A = (a_ij) ∈^n × n_≥ 0.

Let us define S as

    S = (i, j) | i, j ∈ [n]  and  a_ij≠ 0,

where we remove zero entries from the outcome space S as our formulation cannot treat zero probability, and give each probability as p((i, j)) = a_ij / ∑_ij a_ij.
The partial order ≤ of S is naturally introduced as

    x = (i, j) ≤ y = (k, l) ⇔ i ≤ j  and  k ≤ l,

resulting in = (1, 1).
In addition, we define ι⃗_k, m for each k ∈ [n] and m ∈{1, 2} such that

    ι⃗_k, m = minx = (i_1, i_2) ∈ S | i_m = k,

where the minimum is with respect to the order ≤.
If ι⃗_k, m does not exist, we just remove the entire kth row if m = 1 or kth column if m = 2 from A.
Then we switch rows and columns of A so that the condition

    ι⃗_1, m≤ι⃗_2, m≤…≤ι⃗_n, m

is satisfied for each m ∈{1, 2}, which is possible for any matrices. 
Since we have

    η(ι⃗_k, m) - η(ι⃗_k + 1, m) = {[ ∑_j = 1^n p((k, j))          if  m = 1,; ∑_i = 1^n p((i, k))           if  m = 2 ].

if the condition (<ref>) is satisfied, the probability distribution is balanced if for all k ∈ [n] and m ∈{1, 2}

    η(ι⃗_k, m) = n - k + 1/n.

Therefore, we obtain the following result.



  
Matrix balancing as e-projection:
Given a matrix A ∈^n × n with its normalized probability distribution P ∈ such that p((i, j)) = a_ij / ∑_ij a_ij.
Define the poset (S, ≤) by Equations (<ref>) and (<ref>) and let (β) be the submanifold of  such that

    (β) = P ∈ | η_P(x) = β(x)  for all  x ∈(β),

where the function β is given as

    (β)    = ι⃗_k, m∈ S | k ∈ [n], m ∈{1, 2},
    β(ι⃗_k, m)    = n - k + 1/n.

Matrix balancing is the e-projection of P onto the submanifold (β), that is, the balanced matrix (RAS) / n is the distribution P_β such that

    {[ θ_P_β(x) = θ_P(x) if  x ∈ S^+ ∖(β),;   η_P_β(x) = β(x)       if  x ∈(β), ].

which is unique and always exists in , thanks to its dually flat structure.
Moreover, two balancing vectors r⃗ and s⃗ are

    exp( ∑_k = 1^i θ_P_β(ι⃗_k, m) - θ_P(ι⃗_k, m)) = {[        r_i if  m = 1,;        a_i if  m = 2, ].

for every i ∈ [n] and r⃗ = r⃗ n / ∑_ij a_ij.▪



 §.§ Tensor Balancing

Next, we generalize our approach from matrices to tensors.
For an Nth order tensor A = (a_i_1 i_2 … i_N) ∈^n_1 × n_2 ×…× n_N and a vector b⃗∈^n_m, the m-mode product of A and b⃗ is defined as

    (A ×_m b⃗)_i_1 … i_m - 1 i_m + 1… i_N = ∑_i_m = 1^n_m a_i_1 i_2 … i_N b_i_m.

We define tensor balancing as follows:
Given a tensor A ∈^n_1 × n_2 ×…× n_N with n_1 = … = n_N = n,
find (N - 1) order tensors R^1, R^2, …, R^N such that

    A' ×_m1⃗ = 1⃗   (∈^n_1 ×…× n_m - 1× n_m + 1×…× n_N)

for all m ∈ [N], i.e., ∑_i_m = 1^n a_i_1 i_2… i_N' = 1, where each entry a'_i_1 i_2… i_N of the balanced tensor A' is given as

    a'_i_1 i_2… i_N = a_i_1 i_2 … i_N ∏_m ∈ [N] R^m_i_1 … i_m - 1 i_m + 1… i_N.

A tensor A' that satisfies Equation (<ref>) is called multistochastic <cit.>.
Note that this is exactly the same as the matrix balancing problem if N = 2.

It is straightforward to extend matrix balancing to tensor balancing as e-projection onto a submanifold.
Given a tensor A ∈^n_1 × n_2 ×…× n_N with its normalized probability distribution P such that

    p(x) = a_i_1 i_2… i_N / ∑_j_1 j_2 … j_N a_j_1 j_2… j_N

for all x = (i_1, i_2, …, i_N).
The objective is to obtain P_β such that ∑_i_m = 1^n p_β( (i_1, …, i_N)) = 1/(n^N - 1) for all m ∈ [N] and i_1, …, i_N ∈ [n].
In the same way as matrix balancing, we define S as

    S = (i_1, i_2, …, i_N) ∈ [n]^N | a_i_1 i_2… i_N≠ 0

with removing zero entries and the partial order ≤ as

    x = (i_1 … i_N) ≤ y = (j_1 … j_N) ⇔∀ m ∈ [N], i_m ≤ j_m.

In addition, we introduce ι⃗_k, m as

    ι⃗_k, m = minx = (i_1, i_2, …, i_N) ∈ S | i_m = k.

and require the condition in Equation (<ref>).



  
Tensor balancing as e-projection:
Given a tensor A ∈^n_1 × n_2 ×…× n_N with its normalized probability distribution P ∈ given in Equation (<ref>).
The submanifold (β) of multistochastic tensors is given as

    (β) = P ∈ | η_P(x) = β(x)  for all  x ∈(β),

where the domain of the function β is given as

    (β) = ι⃗_k, m | k ∈ [n], m ∈ [N]

and each value is described using the zeta function as

    β(ι⃗_k, m) = ∑_l ∈ [n]ζ(ι⃗_k, m, ι⃗_l, m) 1/n^N - 1.

Tensor balancing is the e-projection of P onto the submanifold (β), that is, the multistochastic tensor is the distribution P_β such that

    {[ θ_P_β(x) = θ_P(x) if  x ∈ S^+ ∖(β),;   η_P_β(x) = β(x)       if  x ∈(β), ].

which is unique and always exists in , thanks to its dually flat structure.
Moreover, each balancing tensor R^m is

    R^m_i_1… i_m - 1i_m + 1… i_N = exp(∑_m' ≠ m∑_k = 1^i_m'θ_P_β(ι⃗_k, m') - θ_P(ι⃗_k, m'))

for every m ∈ [N] and R^1 = R^1 n^N - 1 / ∑_j_1… j_N a_j_1… j_N to recover a multistochastic tensor.▪



Our result means that the e-projection algorithm based on Newton's method proposed in Section <ref> converges to the unique balanced tensor whenever (β) ≠∅ holds.



§ CONCLUSION

In this paper, we have solved the open problem of tensor balancing and presented an efficient balancing algorithm using Newton's method.
Our algorithm quadratically converges, while the popular Sinkhorn-Knopp algorithm linearly converges.
We have examined the efficiency of our algorithm in numerical experiments on matrix balancing and showed that the proposed algorithm is several orders of magnitude faster than the existing approaches.

We have analyzed theories behind the algorithm, and proved that balancing is e-projection in a special type of a statistical manifold, in particular, a dually flat Riemannian manifold studied in information geometry.
Our key finding is that the gradient of the manifold, equivalent to Riemannian metric or the Fisher information matrix, can be analytically obtained using the Möbius inversion formula.

Our information geometric formulation can model several machine learning applications such as statistical analysis on a DAG structure.
Thus, we can perform efficient learning as projection using information of the gradient of manifolds by reformulating such models, which we will study in future work.



§ ACKNOWLEDGEMENTS

The authors sincerely thank Marco Cuturi for his valuable comments.
This work was supported by JSPS KAKENHI Grant Numbers JP16K16115, JP16H02870 (MS), JP26120732 and JP16H06570 (HN).
The research of K.T. was supported by JST CREST JPMJCR1502, RIKEN PostK, KAKENHI Nanostructure and KAKENHI JP15H05711.


38


urlstyle

[Agresti(2012)]Agresti12
A. Agresti.
Categorical data analysis.
Wiley, 3 edition, 2012.

[Ahmed et al.(2003)Ahmed, De Loera, and Hemmecke]Ahmed03
M. Ahmed, J. De Loera, and R. Hemmecke.
Polyhedral Cones of Magic Cubes and Squares, volume 25 of
  Algorithms and Combinatorics, pages 25–41.
Springer, 2003.

[Akartunalı and Knight(2016)]Akartunali16
K. Akartunalı and P. A. Knight.
Network models and biproportional rounding for fair seat allocations
  in the UK elections.
Annals of Operations Research, pages 1–19, 2016.

[Amari(2001)]Amari01
S. Amari.
Information geometry on hierarchy of probability distributions.
IEEE Transactions on Information Theory, 470
  (5):0 1701–1711, 2001.

[Amari(2009)]Amari09
S. Amari.
Information geometry and its applications: Convex function and dually
  flat manifold.
In F. Nielsen, editor, Emerging Trends in Visual Computing: LIX
  Fall Colloquium, ETVC 2008, Revised Invited Papers, pages 75–102. Springer,
  2009.

[Amari(2014)]Amari14
S. Amari.
Information geometry of positive measures and positive-definite
  matrices: Decomposable dually flat structure.
Entropy, 160 (4):0 2131–2145, 2014.

[Amari(2016)]Amari16
S. Amari.
Information Geometry and Its Applications.
Springer, 2016.

[Balinski(2008)]Balinski08
M. Balinski.
Fair majority voting (or how to eliminate gerrymandering).
American Mathematical Monthly, 1150 (2):0
  97–113, 2008.

[Censor and Lent(1981)]Censor81
Y. Censor and A. Lent.
An iterative row-action method for interval convex programming.
Journal of Optimization Theory and Applications, 340
  (3):0 321–353, 1981.

[Chang et al.(2016)Chang, Paksoy, and Zhang]Chang16
H. Chang, V. E. Paksoy, and F. Zhang.
Polytopes of stochastic tensors.
Annals of Functional Analysis, 70 (3):0
  386–393, 2016.

[Cui et al.(2014)Cui, Li, and Ng]Cui14
L.-B. Cui, W. Li, and M. K. Ng.
Birkhoff–von Neumann theorem for multistochastic tensors.
SIAM Journal on Matrix Analysis and Applications, 350
  (3):0 956–973, 2014.

[Cuturi(2013)]Cuturi13
M. Cuturi.
Sinkhorn distances: Lightspeed computation of optimal transport.
In Advances in Neural Information Processing Systems 26, pages
  2292–2300, 2013.

[Frogner et al.(2015)Frogner, Zhang, Mobahi, Araya, and
  Poggio]Frogner15
C. Frogner, C. Zhang, H. Mobahi, M. Araya, and T. A. Poggio.
Learning with a Wasserstein loss.
In Advances in Neural Information Processing Systems 28, pages
  2053–2061, 2015.

[Ganmor et al.(2011)Ganmor, Segev, and Schneidman]Ganmor11
E. Ganmor, R. Segev, and E. Schneidman.
Sparse low-order interaction network underlies a highly correlated
  and learnable neural population code.
Proceedings of the National Academy of Sciences, 1080
  (23):0 9679–9684, 2011.

[Gierz et al.(2003)Gierz, Hofmann, Keimel, Lawson, Mislove, and
  Scott]Gierz03
G. Gierz, K. H. Hofmann, K. Keimel, J. D. Lawson, M. Mislove, and D. S. Scott.
Continuous Lattices and Domains.
Cambridge University Press, 2003.

[Idel(2016)]Idel16
M. Idel.
A review of matrix scaling and sinkhorn's normal form for matrices
  and positive maps.
arXiv:1609.06349, 2016.

[Ito(1993)]Ito93
K. Ito, editor.
Encyclopedic Dictionary of Mathematics.
The MIT Press, 2 edition, 1993.

[Knight(2008)]Knight08
P. A. Knight.
The Sinkhorn–Knopp algorithm: Convergence and applications.
SIAM Journal on Matrix Analysis and Applications, 300
  (1):0 261–275, 2008.

[Knight and Ruiz(2013)]Knight13
P. A. Knight and D. Ruiz.
A fast algorithm for matrix balancing.
IMA Journal of Numerical Analysis, 330 (3):0
  1029–1047, 2013.

[Lahr and de Mesnard(2004)]Lahr04
M. Lahr and L. de Mesnard.
Biproportional techniques in input-output analysis: Table updating
  and structural analysis.
Economic Systems Research, 160 (2):0 115–134,
  2004.

[Lamond and Stewart(1981)]Lamond81
B. Lamond and N. F. Stewart.
Bregman's balancing method.
Transportation Research Part B: Methodological, 150
  (4):0 239–248, 1981.

[Livne and Golub(2004)]Livne04
O. E. Livne and G. H. Golub.
Scaling by binormalization.
Numerical Algorithms, 350 (1):0 97–120, 2004.

[Marshall and Olkin(1968)]Marshall68
A. W. Marshall and I. Olkin.
Scaling of matrices to achieve specified row and column sums.
Numerische Mathematik, 120 (1):0 83–90, 1968.

[Miller and Blair(2009)]Miller09
R. E. Miller and P. D. Blair.
Input-Output Analysis: Foundations and Extensions.
Cambridge University Press, 2 edition, 2009.

[Moon et al.(2009)Moon, Gunther, and Kupin]Moon09
T. K. Moon, J. H. Gunther, and J. J. Kupin.
Sinkhorn solves sudoku.
IEEE Transactions on Information Theory, 550
  (4):0 1741–1746, 2009.

[Nakahara and Amari(2002)]Nakahara02
H. Nakahara and S. Amari.
Information-geometric measure for neural spikes.
Neural Computation, 140 (10):0 2269–2316,
  2002.

[Nakahara et al.(2003)Nakahara, Nishimura, Inoue, Hori, and
  Amari]Nakahara03
H. Nakahara, S. Nishimura, M. Inoue, G. Hori, and S. Amari.
Gene interaction in DNA microarray data is decomposed by
  information geometric measure.
Bioinformatics, 190 (9):0 1124–1131, 2003.

[Nakahara et al.(2006)Nakahara, Amari, and Richmond]Nakahara06
H. Nakahara, S. Amari, and B. J. Richmond.
A comparison of descriptive models of a single spike train by
  information-geometric measure.
Neural computation, 180 (3):0 545–568, 2006.

[Parikh(1979)]Parikh79
A. Parikh.
Forecasts of input-output matrices using the R.A.S. method.
The Review of Economics and Statistics, 610
  (3):0 477–481, 1979.

[Parlett and Landis(1982)]Parlett82
B. N. Parlett and T. L. Landis.
Methods for scaling to doubly stochastic form.
Linear Algebra and its Applications, 48:0 53–79,
  1982.

[Rao et al.(2014)Rao, Huntley, Durand, Stamenova, Bochkov, Robinson,
  Sanborn, Machol, Omer, Lander, and Aiden]Rao14
S. S. P. Rao, M. H. Huntley, N. C. Durand, E. K. Stamenova, I. D. Bochkov,
  J. T. Robinson, A. L. Sanborn, I. Machol, A. D. Omer, E. S. Lander, and E. L.
  Aiden.
A 3D map of the human genome at kilobase resolution reveals
  principles of chromatin looping.
Cell, 1590 (7):0 1665–1680, 2014.

[Rota(1964)]Rota64
G.-C. Rota.
On the foundations of combinatorial theory I: Theory of
  Möbius functions.
Z. Wahrseheinlichkeitstheorie, 2:0 340–368, 1964.

[Sinkhorn(1964)]Sinkhorn64
R. Sinkhorn.
A relationship between arbitrary positive matrices and doubly
  stochastic matrices.
The Annals of Mathematical Statistics, 350
  (2):0 876–879, 06 1964.

[Sinkhorn and Knopp(1967)]Sinkhorn67
R. Sinkhorn and P. Knopp.
Concerning nonnegative matrices and doubly stochastic matrices.
Pacific Journal of Mathematics, 210 (2):0
  343–348, 1967.

[Solomon et al.(2015)Solomon, de Goes, Peyré, Cuturi, Butscher,
  Nguyen, Du, and Guibas]Solomon15
J. Solomon, F. de Goes, G. Peyré, M. Cuturi, A. Butscher, A. Nguyen, T. Du,
  and L. Guibas.
Convolutional Wasserstein distances: Efficient optimal
  transportation on geometric domains.
ACM Transactions on Graphics, 340 (4):0
  66:1–66:11, 2015.

[Soules(1991)]Soules91
G. W. Soules.
The rate of convergence of sinkhorn balancing.
Linear Algebra and its Applications, 150:0 3–40,
  1991.

[Sugiyama et al.(2016)Sugiyama, Nakahara, and Tsuda]Sugiyama2016ISIT
M. Sugiyama, H. Nakahara, and K. Tsuda.
Information decomposition on structured space.
In 2016 IEEE International Symposium on Information Theory,
  pages 575–579, July 2016.

[Wu and Michor(2016)]Wu16
H.-J. Wu and F. Michor.
A computational strategy to adjust for copy number in tumor Hi-C
  data.
Bioinformatics, 320 (24):0 3695–3701, 2016.





