Tensor Balancing on Statistical Manifold

Mahito Sugiyama 1 2 Hiroyuki Nakahara 3 Koji Tsuda 4 5 6Abstract
We solve tensor balancing, rescaling an N th order nonnegative tensor by multiplying N tensors of order N − 1 so that every fiber sums toone. This generalizes a fundamental process ofmatrix balancing used to compare matrices in awide range of applications from biology to economics. We present an efficient balancing algorithm with quadratic convergence using Newton’s method and show in numerical experimentsthat the proposed algorithm is several orders ofmagnitude faster than existing ones. To theoretically prove the correctness of the algorithm,
we model tensors as probability distributions ina statistical manifold and realize tensor balancing as projection onto a submanifold. The key toour algorithm is that the gradient of the manifold,used as a Jacobian matrix in Newton’s method,
can be analytically obtained using the Mö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 learningmodels such as weighted DAGs and Boltzmann
machines.

1. Introduction
Matrix balancing is the problem of rescaling a given square
nonnegative matrix A ∈ Rn×n
≥0 to a doubly stochastic matrix RAS, where every row and column sums to one, bymultiplying two diagonal matrices R and S. This is afundamental process for analyzing and comparing matrices in a wide range of applications, including input-output
analysis in economics, called the RAS approach (Parikh,
1979; Miller & Blair, 2009; Lahr & de Mesnard, 2004),
seat assignments in elections (Balinski, 2008; Akartunalı &1
National Institute of Informatics 2 JST PRESTO 3 RIKEN
Brain Science Institute 4 Graduate School of Frontier Sciences,
The University of Tokyo 5 RIKEN AIP 6 NIMS. Correspondence
to: Mahito Sugiyama <mahito@nii.ac.jp>.

Proceedings of the 34 th International Conference on Machine
Learning, Sydney, Australia, PMLR 70, 2017. Copyright 2017by the author(s).

Statistical manifold 𝓢(dually flat Riemannian manifold)

Given tensor A

Probability
distribution P
Tensor balancing

Projection

Multistochastic tensor A’

Submanifold 𝓢(β)

Every fibersums to 1

Projected
distribution Pβ

Figure 1. Overview of our approach.

Knight, 2016), Hi-C data analysis (Rao et al., 2014; Wu &Michor, 2016), the Sudoku puzzle (Moon et al., 2009), andthe optimal transportation problem (Cuturi, 2013; Frogner
et al., 2015; Solomon et al., 2015). An excellent review ofthis theory and its applications is given by Idel (2016).
The standard matrix balancing algorithm is the SinkhornKnopp algorithm (Sinkhorn, 1964; Sinkhorn & Knopp,
1967; Marshall & Olkin, 1968; Knight, 2008), a special
case of Bregman’s balancing method (Lamond & Stewart,1981) that iterates rescaling of each row and column untilconvergence. The algorithm is widely used in the above
applications due to its simple implementation and theoretically guaranteed convergence. However, the algorithm
converges linearly (Soules, 1991), which is prohibitively
slow for recently emerging large and sparse matrices. Although Livne & Golub (2004) and Knight & Ruiz (2013)
tried to achieve faster convergence by approximating eachstep of Newton’s method, the exact Newton’s method withquadratic 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
N th order nonnegative tensor to a multistochastic tensor,
in which every fiber sums to one, by multiplying (N − 1)thorder N tensors. There are some results about mathematical properties of multistochastic tensors (Cui et al., 2014;
Chang et al., 2016; Ahmed et al., 2003). However, there
is no result for tensor balancing algorithms with guaranteed convergence that transforms a given tensor to a multistochastic tensor until now.

Tensor Balancing on Statistical Manifold

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 1),
which is a statistical manifold and known to be the essential structure for probability distributions in informationgeometry (Amari, 2016). Using a partially ordered outcome space, we generalize the log-linear model (Agresti,2012) used to model the higher-order combinations of binary variables (Amari, 2001; Ganmor et al., 2011; Nakahara & Amari, 2002; Nakahara et al., 2003), which allows
us to model tensors as probability distributions in the statistical manifold. The remarkable property of our model isthat the gradient of the manifold can be analytically computed using the Möbius inversion formula (Rota, 1964), theheart of combinatorial mathematics (Ito, 1993), which enables us to directly obtain the Jacobian matrix in Newton’s
method. Moreover, we show that (n − 1)N entries for thesize nN of a tensor are invariant with respect to one of thetwo coordinate systems of the statistical manifold. Thus
the number of equations in Newton’s method is O(nN −1 ).
The remainder of this paper is organized as follows: We
begin with a low-level description of our matrix balancing
algorithm in Section 2 and demonstrate its efficiency in numerical experiments in Section 3. To guarantee the correctness of the algorithm and extend it to tensor balancing, we
provide theoretical analysis in Section 4. In Section 4.1, we
introduce a generalized log-linear model associated with apartial order structured outcome space, followed by introducing the dually flat Riemannian structure in Section 4.2.
In Section 4.3, we show how to use Newton’s method tocompute projection of a probability distribution onto a submanifold. Finally, we formulate the matrix and tensor balancing problem in Section 5 and summarize our contributions in Section 6.

2. The Matrix Balancing Algorithm
Given a nonnegative square matrix A = (aij ) ∈ Rn×n
≥0 , thetask of matrix balancing is to find r, s ∈ Rn that satisfy(RAS)1 = 1,

(RAS)T 1 = 1,

(1)where R = diag(r) and S = diag(s). The balanced matrix
A′ = RAS is called doubly stochastic, in which each entry
a′ij = aij ri sj 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) ands = 1/(AT r). We denote by [n] = {1, 2, . . . , n} hereafter.
In our algorithm, instead of directly updating r and s, weupdate two parameters θ and η defined as∑∑∑∑pi′ j ′
(2)
θi′ j ′ , ηij =log pij =i′ ≤i j ′ ≤ji′ ≥i j ′ ≥jMatrixa11 a12 a13 a14

Constraints for balancing
η11 η12 η13 η14a21 a22 a23 a24

η21 θ22 θ23 θ24

a31 a32 a33 a34

η31 θ32 θ33 θ34

a41 a42 a43 a44

η41 θ42 θ43 θ44 Invariant

Figure 2. Matrix balancing with two parameters θ and η.

for each∑ i, j ∈ [n], where∑ we normalized entries as pij =aij / ij aij so that ij pij = 1. We assume for simplicity that each entry is strictly larger than zero. The assumption will be removed in Section 5.
(t)

The key to our approach is that we update θij with i = 1
or j = 1 by Newton’s method at each iteration t = 1, 2, . . .
(t)while fixing θij with i, j ̸= 1 so that ηij satisfies the following condition (Figure 2):(t)ηi1 = (n − i + 1)/n,(t)

η1j = (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
 (t)
 (t+1)   (t) θ11θ11η11 − (n − 1 + 1)/n .   . ..

 ..   .. 
.


  

 (t)

 (t+1)   (t) θ1n  θ1n −1 η1n − (n − n + 1)/n , (3) (t+1)  =  (t)  − J  (t) η21 − (n − 2 + 1)/n θ21  θ21 



  


 ..   .. ..

 .   . 

.
(t)
(t)
(t+1)
ηn1 − (n − n + 1)/nθn1θn1where J is the Jacobian matrix given as(t)J(ij)(i′ j ′ ) =

∂ηij

(t)

∂θi′ j ′= ηmax{i,i′ } max{j,j ′ } −n2 ηij ηi′ j ′ , (4)which is derived from our theoretical result in Theorem 3.
Since J is a (2n−1)×(2n−1) matrix, the time complexityof each update is O(n3 ), which is needed to compute theinverse of J.
(t+1)

(t+1)

(t+1)

After updating to θij , we can compute pij
and ηij
by Equation (2). Since this update does not ensure the∑ (t+1)(t+1)
condition= 1, we again update θ11
asij pij
∑
(t+1)
(t+1)
(t+1)
(t+1)
θ11
= θ11
− log ij pij
and recompute pij(t+1)

and ηij

for each i, j ∈ [n].

By iterating the above update process in Equation (3) untilconvergence, A = (aij ) with aij = npij becomes doubly
stochastic.

3. Numerical Experiments
We evaluate the efficiency of our algorithm compared to thetwo prominent balancing methods, the standard SinkhornKnopp algorithm (Sinkhorn, 1964) and the state-of-the-art

104
102
10

104
102
10
10–2
10–4

10

50

500
5000
n
Newton (proposed)

1050

Sinkhorn

n

500

5000

BNEWT

Figure 3. Results on Hessenberg matrices. The BNEWT algorithm (green) failed to converge for n ≥ 200.

105

109

Running time (sec.)

106

Number of iterations

108

Running time (sec.)

Number of iterations

Tensor Balancing on Statistical Manifold

106
103
10
20

100

200
300
n
Newton (proposed)

103
101
10–1
10–3
20

100

Sinkhorn

n

200

300

BNEWT

Figure 5. Results on Trefethen matrices. The BNEWT algorithm(green) failed to converge for n ≥ 200.

Residual

10
10–1
–

10 3
NewtonSinkhorn
BNEWT

10–5
10–7
0

500

1000
1500
2000
Number of iterations

2500

3000

Figure 4. Convergence graph on H20 .

algorithm BNEWT (Knight & Ruiz, 2013), 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 C++ with the Eigen library andcompiled with gcc 4.8.31 . We have carefully implemented
BNEWT by directly translating the MATLAB code provided in (Knight & Ruiz, 2013) into C++ with the Eigen
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′T 1 − 1)∥2 , where each entry a′ij is obtained as npij in our algorithm, and ran eachof 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 (Parlett & Landis, 1982; Knight &Ruiz, 2013). Each entry of an n × n Hessenberg matrix
Hn = (hij ) is given as hij = 0 if j < i − 1 and hij = 1
otherwise. We varied the size n from 10 to 5, 000, andmeasured running time (in seconds) and the number of iterations of each method.
Results are plotted in Figure 3. Our balancing algorithmwith the Newton’s method (plotted in blue in the figures)
1

An implementation of algorithms for matrices and thirdorder tensors is available at: https://github.com/
mahito-sugiyama/newton-balancing

is clearly the fastest: It is three to five orders of magnitude
faster than the standard Sinkhorn-Knopp algorithm (plottedin 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 theoriginal paper (Knight & Ruiz, 2013) 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 BNEWTand Sinkhorn-Knopp, respectively, at n = 100.
To see the behavior of the rate of convergence in detail, weplot the convergence graph in Figure 4 for n = 20, where
we observe the slow convergence rate of the SinkhornKnopp 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 website2 , which are nonnegative diagonal matrices with primes. Results are plotted in
Figure 5, where we observe the same trend as before: Ouralgorithm 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, whichis the necessary condition for matrix balancing (Knight &Ruiz, 2013), while the BNEWT algorithm fails to converge
if n = 200 or n = 300.

4. Theoretical Analysis
In the following, we provide theoretical support to our algorithm by formulating the problem as a projection withina statistical manifold, in which a matrix corresponds to anelement, that is, a probability distribution, in the manifold.
We show that a balanced matrix forms a submanifold andmatrix balancing is projection of a given distribution onto
the submanifold, where the Jacobian matrix in Equation (4)is derived from the gradient of the manifold.
2
http://www.cise.ufl.edu/research/sparse/
matrices/Tensor Balancing on Statistical Manifold

4.1. Formulation
We introduce our log-linear probabilistic model, where theoutcome space is a partially ordered set, or a poset (Gierz
et al., 2003). We prepare basic notations and the key mathematical tool for posets, the Möbius inversion formula, followed by formulating the log-linear model.

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 Möbius functions, let us introduce
two mappings θ : S → R and η : S → R as∑
θ(x) =µ(s, x) log p(s),(6)
s∈S

4.1.1. M ÖBIUS I NVERSION
A poset (S, ≤), the set of elements S and a partial order≤ on S, is a fundamental structured space in computerscience. 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 + .
Rota (1964) introduced the Möbius inversion formula on
posets by generalizing the inclusion-exclusion principle.
Let ζ : S × S → {0, 1} be the zeta function defined as{
1 if s ≤ x,
ζ(s, x) =0 otherwise.
The Möbius function µ : S ×S → Z satisfies ζµ = I, whichis inductively defined for all x, y with x ≤ y as
if x = y,
 1 ∑
−
µ(x,
s)
if x < y,µ(x, y) =x≤s<y0otherwise.

∑

x≤s≤y

ζ(x, s)µ(s, y) =

s∈S

∑

(5)

x≤s≤y

with the Kronecker delta δ such that δxy = 1 if x = y andδxy = 0 otherwise. Then for any functions f , g, and h withthe domain S such that∑∑g(x) =
ζ(s, x)f (s) =f (s),
h(x) =

∑

s∈S

s≤x

ζ(x, s)f (s) =∑f (s),

s≥x

f is uniquely recovered with the Möbius function:∑∑f (x) =µ(s, x)g(s), f (x) =µ(x, s)h(s).
s∈Sζ(x, s)p(s) =

s∈S

∑

p(s).

From the Möbius inversion formula, we have∑∑log p(x) =ζ(s, x)θ(s) =
θ(s),
s∈S

p(x) =

(7)

s≥x

∑

(8)

s≤x

µ(x, s)η(s).

(9)

s∈S

They are generalization of the log-linear model (Agresti,2012) that gives the probability p(x) of an n-dimensional
binary vector x = (x1 , . . . , xn ) ∈ {0, 1}n as∑∑∑θijk xi xj xk
log p(x) =θ i xi +θij xi xj +
i

i<j

i<j<k

+ · · · + θ1...n x1 x2 . . . xn − ψ,
where θ = (θ1 , . . . , θ12...n ) is a parameter vector, ψ is anormalizer, and η = (η 1 , . . . , η 12...n ) represents the expectation of variable combinations such thatη i = E[xi ] = Pr(xi = 1),
η 1...n = E[x1 . . . xn ] = Pr(x1 = · · · = xn = 1).

µ(s, y) = δxys∈S∑η ij = E[xi xj ] = Pr(xi = xj = 1), i < j, . . .

From the definition, it follows that∑∑
ζ(s, y)µ(x, s) =µ(x, s) = δxy ,s∈S

η(x) =

They coincide with Equations (8) and (7) when we let
S = 2V with V = {1, 2, . . . , n}, each x ∈ S as the setof indices of “1” of x, and the order ≤ as the inclusion relationship, that is, x ≤ y if and only if x ⊆ y. Nakahara
et al. (2006) have pointed out that θ can be computed fromp using the inclusion-exclusion principle in the log-linear
model. We exploit this combinatorial property of the loglinear model using the Möbius inversion formula on posets
and extend the log-linear model from the power set 2V to
any kind of posets (S, ≤). Sugiyama et al. (2016) studied arelevant log-linear model, but the relationship with Möbius
inversion formula has not been analyzed yet.
4.2. Dually Flat Riemannian Manifolds∈SThis is called the Möbius inversion formula and is at theheart of enumerative combinatorics (Ito, 1993).
4.1.2. L OG -L INEAR M ODEL ON P OSETS
We consider a probability vector p on (S, ≤) that gives a
discrete probability distribution with the outcome space S.

We theoretically analyze our log-linear model introduced in
Equations (6), (7) and show that they form dual coordinate
systems on a dually flat manifold, which has been mainly
studied in the area of information geometry (Amari, 2001;
Nakahara & Amari, 2002; Amari, 2014; 2016). Moreover,
we show that the Riemannian metric and connection of ourmodel can be analytically computed in closed forms.

Tensor Balancing on Statistical Manifold

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 S theset of probability distributions specified by probability vectors, which forms a statistical manifold. We use uppercase
letters P, Q, R, . . . for points (distributions) in S and theirlowercase 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 (6) and (7), respectively,
and abbreviate subscripts if there is no ambiguity.
4.2.1. D UALLY F LAT S TRUCTURE
We show that S has the dually flat Riemannian structure
induced by two functions θ and η in Equation (6) and (7).
We define ψ(θ) as
ψ(θ) = −θ(⊥) = − log p(⊥),(10)

which corresponds to the normalizer of p. It is a convex
function since we have∑∑ψ(θ) = log
exp θ(s) 
∑

⊥<s≤x

x∈S

from log p(x) = ⊥<s≤x θ(s) − ψ(θ). We apply the Legendre transformation to ψ(θ) given as
(
)
∑
′
′
′
φ(η) = max
θ
η
−
ψ(θ
)
,
θ
η
=
θ′ (x)η(x). (11)
′
θ

x∈S +

Proof. They can be directly derived from our definitions(Equations (6) and (11)) as
(∑
)
∑
exp
θ(s)
∑
y≥x
⊥<s≤y
∂ψ(θ)
(∑
) =
=∑
p(s) = η(x),
∂θ(x)
exp
θ(s)
y∈S

(

⊥<s≤y

s≥x

)

∂φ(η)
∂
=
θη − ψ(θ) = θ(x).
∂η(x)
∂η(x)
Moreover, we can confirm the orthogonality of θ and η as[
] ∑∂ log p(s) ∂ log p(s)E=
ζ(x, s)µ(s, y) = δxy .
∂θ(x)
∂η(y)
s∈S

The last equation holds from Equation (5), hence theMöbius inversion directly leads to the orthogonality.
The Bregman divergence is known to be the canonical divergence (Amari, 2016, Section 6.6) 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 1
and Equation (12), it is given as∑q(x)D [P, Q] =q(x) log
,
p(x)
x∈S

Then φ(η) coincides with the negative entropy.

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

Theorem 1 (Legendre dual).
∑φ(η) =
p(x) log p(x).

4.2.2. R IEMANNIAN S TRUCTURE
Next we analyze the Riemannian structure on S and showthat the Möbius inversion formula enables us to compute
the Riemannian metric of S.

x∈S

Proof. From Equation (5), we have


∑
∑
∑

θ′ η =µ(s, x) log p′ (s)p(s) 
⊥<s≤x

x∈S +

=

∑

Theorem 3 (Riemannian metric). The manifold (S, g(ξ))
is a Riemannian manifold with the Riemannian metric g(ξ)
such that for all x, y ∈ S +
]
∑[


ζ(x, s)ζ(y, s)p(s) − η(x)η(y) if ξ = θ,s∈Sgxy (ξ) = ∑µ(s, x)µ(s, y)p(s)−1if ξ = η.


s≥x
′

′

p(x) ( log p (x) − log p (⊥) ) .

x∈S +Thus it holds thatθ′ η − ψ(θ′ ) =∑p(x) log p′ (x).

x∈S

Proof. Since the Riemannian metric is defined asHence it is maximized with p(x) = p′ (x).

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

Since they are connected with each other by the Legendre
transformation, they form a dual coordinate system ∇ψ(θ)
and ∇φ(η) of S (Amari, 2016, Section 1.5), which coincides with θ and η as follows.
Theorem 2 (dual coordinate system).
∇ψ(θ) = η,∇φ(η) = θ.

s∈S

(12)

(13)

g(η) = ∇∇φ(η),

when ξ = θ we have∂2∂
ψ(θ) =
η(y)∂θ(x)∂θ(y)
∂θ(x)


∑
∂ ∑
ζ(y, s) exp θ(u) − ψ(θ)
=∂θ(x)gxy (θ) =

s∈S

⊥<u≤s

Tensor Balancing on Statistical Manifold=∑

ζ(x, s)ζ(y, s)p(s) − |S|η(x)η(y).

s∈SWhen ξ = η, it follows that∂2∂φ(η) =
θ(y)∂η(x)∂η(y)
∂η(x)


∑
∂ ∑
=
µ(s, y) log µ(s, u)η(u)
∂η(x)
s≤y
u≥s
∑
=
µ(s, x)µ(s, y)p(s)−1 .

gxy (η) =s∈SSince g(ξ) coincides with the Fisher information matrix,
[
]
∂
∂
E
log p(s)log p(s) = gxy (θ),
∂θ(x)
∂θ(y)
[
]
∂
∂
E
log p(s)log p(s) = gxy (η).
∂η(x)
∂η(y)
Then the Riemannian (Levi–Chivita) connection Γ(ξ) withrespect to ξ, which is defined as()
1 ∂gyz (ξ) ∂gxz (ξ) ∂gxy (ξ)Γxyz (ξ) =
+
−
2
∂ξ(x)
∂ξ(y)
∂ξ(z)
for all x, y, z ∈ S , can be analytically obtained.
+

Theorem 4 (Riemannian connection). The Riemannian
connection Γ(ξ) on the manifold (S, g(ξ)) is given in thefollowing for all x, y, z ∈ S + ,

)(
)
1 ∑(


ζ(x,
s)
−
η(x)
ζ(y,
s)
−
η(y)

 2

)

s∈S (

ζ(z, s) − η(z) p(s)if ξ = θ,Γxyz (ξ) =



1∑


µ(s, x)µ(s, y)µ(s, z)p(s)−2 if ξ = η.
−

 2
s∈S

Proof. Connections Γxyz (θ) and Γxyz (η) can be obtainedby directly computing ∂gyz (θ)/∂θ(x) and ∂gyz (η)/∂η(x),respectively.
4.3. 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 thetarget model (Amari, 2016). Here we define projection of
distributions on posets and show that Newton’s method canbe applied to perform projection as the Jacobian matrix canbe analytically computed.
4.3.1. D EFINITION
Let S(β) be a submanifold of S such that
S(β) = {P ∈ S | θP (x) = β(x), ∀x ∈ dom(β)} (14)specified by a function β with dom(β) ⊆ S + . Projectionof P ∈ S onto S(β), called m-projection, which is definedas the distribution Pβ ∈ S(β) such that{θPβ (x) = β(x)if x ∈ dom(β),ηPβ (x) = ηP (x) if x ∈ S + \ dom(β),
is the minimizer of the KL divergence from P to S(β):
Pβ = argmin DKL [P, Q].
Q∈S(β)The dually flat structure with the coordinate systems θ andη guarantees that the projected distribution Pβ always exists and is unique (Amari, 2009, Theorem 3). Moreover,the Pythagorean theorem holds in the dually flat manifold,that is, for any Q ∈ S(β) we have
DKL [P, Q] = DKL [P, Pβ ] + DKL [Pβ , Q].
We can switch η and θ in the submanifold S(β) by changing DKL [P, Q] to DKL [Q, P ], where the projected distribution Pβ of P is given as{
θPβ (x) = θP (x) if x ∈ S + \ dom(β),ηPβ (x) = β(x) if x ∈ dom(β),
This projection is called e-projection.
Example 1 (Boltzmann machine). Given a Boltzmann machine represented as an undirected graph G = (V, E) witha vertex set V and an edge set E ⊆ {{i, j} | i, j ∈ V }.
The set of probability distributions that can be modeled bya Boltzmann machine G coincides with the submanifold
S B = {P ∈ S | θP (x) = 0 if |x| > 2 or x ̸∈ E},with S = 2V . Let P̂ be an empirical distribution estimated from a given dataset. The learned model is the mprojection of the empirical distribution P̂ onto S 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.
4.3.2. C OMPUTATION
Here we show how to compute projection of a given probability distribution. We show that Newton’s method can beused to efficiently compute the projected distribution Pβ by
(0)
(0)
(1)
(2)
iteratively updating Pβ = P as Pβ , Pβ , Pβ , . . . untilconverging to Pβ .
(0)

Let us start with the m-projection with initializing Pβ=(t)

P . In each iteration t, we update θPβ (x) for all x ∈ domβ(t)while fixing ηPβ (x) = ηP (x) for all x ∈ S + \ dom(β),
which is possible from the orthogonality of θ and η. Using(t+1)
Newton’s method, ηPβ (x) should satisfy
(
) ∑
(
)
(t)
(t+1)
(t)
θPβ (x) − β(x) +Jxy ηPβ (y) − ηPβ (y) = 0,y∈dom(β)

Tensor Balancing on Statistical Manifold

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

Jxy =

∂θPβ (x)(t)∂ηPβ (y)=∑µ(s, x)µ(s, y)pβ (s)−1
(t)

s∈S

from Theorem 3. Therefore, we have the update formulafor all x ∈ dom(β) as
(
)
∑
(t+1)
(t)
(t)
−1
ηPβ (x) = ηPβ (x) −JxyθPβ (y) − β(y) .
y∈dom(β)(t)

In e-projection, update ηPβ (x) for x ∈ dom(β) while fix(t)ing θPβ (x) = θP (x) for all x ∈ S + \ dom(β). To ensure(t)ηPβ (⊥) = 1, we add ⊥ to dom(β) and β(⊥) = 1. We(t)update θPβ (x) at each step t as)
∑ −1 ( (t)
(t+1)
(t)
J ′ xy ηPβ (y) − β(y) ,θPβ (x) = θPβ (x) −y∈dom(β)(t)∂ηPβ (x) ∑(t)
J ′ xy ==ζ(x, s)ζ(y, s)pβ (s)(t)∂θPβ (y) s∈S(t)(t)− |S|ηPβ (x)ηPβ (y).
(t)

In this case, we also need to update θPβ (⊥) as it is notguaranteed to be fixed. Let us define
(
)
(t+1)
exp θPβ (s)∏′(t+1)(t)() ζ(s, x).
pβ(x) = pβ (x)(t)
s∈dom(β) exp θPβ (s)Since we have

)
(
(t+1)
exp θPβ (⊥)(t+1)() p′(t+1)
pβ
(x) =
(x),
β
(t)
exp θPβ (⊥)

it follows that
(t+1)

(t)

(⊥) − θPβ (⊥)()()∑ ′(t+1)(t)= − log exp θPβ (⊥) +
pβ
(x) ,
θPβ

x∈S +

The time complexity of each iteration is O(|dom(β)|3 ),which is required to compute the inverse of the Jacobianmatrix.
Global convergence of the projection algorithm is always
guaranteed by the convexity of a submanifold S(β) definedin Equation (14). Since S(β) 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 (Censor & Lent, 1981).

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

5.1. Matrix Balancing
Recall that the task of matrix balancing is to find r, s ∈ Rn
that satisfy (RAS)1 = 1 and (RAS)T 1 = 1 with R =diag(r) and S = diag(s) for a given nonnegative squaren×nmatrix A = (aij ) ∈ R≥0.
Let us define S as
S = {(i, j) | i, j ∈ [n] and aij ̸= 0},(15)where we remove zero entries from the outcome space S as
our formulation cannot treat zero∑probability, and give eachprobability as p((i, j)) = aij / ij aij . The partial order≤ of S is naturally introduced as
x = (i, j) ≤ y = (k, l) ⇔ i ≤ j and k ≤ l,(16)
resulting in ⊥ = (1, 1). In addition, we define ιk,m foreach k ∈ [n] and m ∈ {1, 2} such that
ιk,m = min{ x = (i1 , i2 ) ∈ S | im = k },where the minimum is with respect to the order ≤. If ιk,mdoes not exist, we just remove the entire kth row if m = 1
or kth column if m = 2 from A. Then we switch rows andcolumns of A so that the condition
ι1,m ≤ ι2,m ≤ · · · ≤ ιn,m(17)
is satisfied for each m ∈ {1, 2}, which is possible for any
matrices. Since we have {∑np((k, j)) if m = 1,
∑j=1
η(ιk,m ) − η(ιk+1,m ) =n
if m = 2i=1 p((i, k))if the condition (17) is satisfied, the probability distribution
is balanced if for all k ∈ [n] and m ∈ {1, 2}
n−k+1.
η(ιk,m ) =nTherefore, we obtain the following result.
Matrix balancing as e-projection:
Given a matrix A ∈ Rn×n with its normalized probabil∑
ity distribution P ∈ S such that p((i, j)) = aij / ij aij .
Define the poset (S, ≤) by Equations (15) and (16) and let
S(β) be the submanifold of S such that
S(β) = {P ∈ S | ηP (x) = β(x) for all x ∈ dom(β)},
where the function β is given asdom(β) = {ιk,m ∈ S | k ∈ [n], m ∈ {1, 2}},n−k+1.
β(ιk,m ) =nMatrix balancing is the e-projection of P onto the submanifold S(β), that is, the balanced matrix (RAS)/n is the
distribution Pβ such that{θPβ (x) = θP (x) if x ∈ S + \ dom(β),ηPβ (x) = β(x) if x ∈ dom(β),
which is unique and always exists in S, thanks to its duallyflat structure. Moreover, two balancing vectors r and s are
( i
) {
∑
ri if m = 1,exp
θPβ (ιk,m ) − θP (ιk,m ) =ai if m = 2,
k=1∑
for every i ∈ [n] and r = rn/ ij aij .
■Tensor Balancing on Statistical Manifold

5.2. Tensor Balancing
Next, we generalize our approach from matrices to tensors.
For an N th order tensor A = (ai1 i2 ...iN ) ∈ Rn1 ×n2 ×···×nN
and a vector b ∈ Rnm , the m-mode product of A and b isdefined as
nm∑(A ×m b)i1 ...im−1 im+1 ...iN =ai1 i2 ...iN bim .
im =1

We define tensor balancing as follows: Given a tensor A ∈Rn1 ×n2 ×···×nN with n1 = · · · = nN = n, find (N − 1)order tensors R1 , R2 , . . . , RN such thatA′ ×m 1 = 1 (∈ Rn1 ×···×nm−1 ×nm+1 ×···×nN ) (18)∑n
for all m ∈ [N ], i.e., im =1 a′i1 i2 ...iN = 1, where eachentry a′i1 i2 ...iN of the balanced tensor A′ is given as∏a′i1 i2 ...iN = ai1 i2 ...iN
Rim1 ...im−1 im+1 ...iN .
m∈[N ]A tensor A′ that satisfies Equation (18) is called multistochastic (Cui et al., 2014). Note that this is exactly thesame 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 ∈ Rn1 ×n2 ×···×nN with its normalized probability
distribution P such that/∑p(x) = ai1 i2 ...iN(19)
aj1 j2 ...jN
j1 j2 ...jN

for all x = ∑(i1 , i2 , . . . , iN ). The objective is to obtainnPβ such that im =1 pβ ((i1 , . . . , iN )) = 1/(nN −1 ) for allm ∈ [N ] and i1 , . . . , iN ∈ [n]. In the same way as matrix
balancing, we define S as{}
S = (i1 , i2 , . . . , iN ) ∈ [n]N ai1 i2 ...iN ̸= 0
with removing zero entries and the partial order ≤ as
x = (i1 . . . iN ) ≤ y = (j1 . . . jN ) ⇔ ∀m ∈ [N ], im ≤ jm .
In addition, we introduce ιk,m asιk,m = min{ x = (i1 , i2 , . . . , iN ) ∈ S | im = k }.
and require the condition in Equation (17).
Tensor balancing as e-projection:
Given a tensor A ∈ Rn1 ×n2 ×···×nN with its normalized
probability distribution P ∈ S given in Equation (19). The
submanifold S(β) of multistochastic tensors is given as
S(β) = {P ∈ S | ηP (x) = β(x) for all x ∈ dom(β)},
where the domain of the function β is given asdom(β) = { ιk,m | k ∈ [n], m ∈ [N ] }and each value is described using the zeta function as∑
1β(ιk,m ) =ζ(ιk,m , ιl,m ) N −1 .
nl∈[n]Tensor balancing is the e-projection of P onto the submanifold S(β), that is, the multistochastic tensor is the distribution Pβ such that{θPβ (x) = θP (x) if x ∈ S + \ dom(β),ηPβ (x) = β(x) if x ∈ dom(β),
which is unique and always exists in S, thanks to its duallyflat structure. Moreover, each balancing tensor Rm is
Rim1 ...im−1 im+1 ...iN


im′
∑ ∑
= exp
θPβ (ιk,m′ ) − θP (ιk,m′ )
m′ ̸=m k=1

for every m ∈ [N ] and R1 = R1 nN −1 /to recover a multistochastic tensor.

∑

j1 ...jN aj1 ...jN■

Our result means that the e-projection algorithm based on
Newton’s method proposed in Section 4.3 converges to theunique balanced tensor whenever S(β) ̸= ∅ holds.

6. 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 ouralgorithm in numerical experiments on matrix balancing
and showed that the proposed algorithm is several ordersof magnitude faster than the existing approaches.
We have analyzed theories behind the algorithm, andproved that balancing is e-projection in a special type ofa statistical manifold, in particular, a dually flat Riemannian manifold studied in information geometry. Our key
finding is that the gradient of the manifold, equivalent toRiemannian metric or the Fisher information matrix, can beanalytically obtained using the Mö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 researchof K.T. was supported by JST CREST JPMJCR1502,
RIKEN PostK, KAKENHI Nanostructure and KAKENHI
JP15H05711.

Tensor Balancing on Statistical Manifold

References
Agresti, A. Categorical data analysis. Wiley, 3 edition,
2012.
Ahmed, M., De Loera, J., and Hemmecke, R. Polyhedral
Cones of Magic Cubes and Squares, volume 25 of Algorithms and Combinatorics, pp. 25–41. Springer, 2003.
Akartunalı, K. and Knight, P. A. Network models andbiproportional rounding for fair seat allocations in the
UK elections. Annals of Operations Research, pp. 1–19,
2016.

the National Academy of Sciences, 108(23):9679–9684,
2011.
Gierz, G., Hofmann, K. H., Keimel, K., Lawson, J. D., Mislove, M., and Scott, D. S. Continuous Lattices and Domains. Cambridge University Press, 2003.
Idel, M. A review of matrix scaling and sinkhorn’s normal
form for matrices and positive maps. arXiv:1609.06349,
2016.
Ito, K. (ed.). Encyclopedic Dictionary of Mathematics. The
MIT Press, 2 edition, 1993.

Amari, S. Information geometry on hierarchy of probability distributions. IEEE Transactions on Information
Theory, 47(5):1701–1711, 2001.

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

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

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

Amari, S. Information geometry of positive measuresand positive-definite matrices: Decomposable dually flat
structure. Entropy, 16(4):2131–2145, 2014.
Amari, S. Information Geometry and Its Applications.
Springer, 2016.
Balinski, M. Fair majority voting (or how to eliminate gerrymandering). American Mathematical Monthly, 115(2):97–113, 2008.
Censor, Y. and Lent, A. An iterative row-action method forinterval convex programming. Journal of Optimization
Theory and Applications, 34(3):321–353, 1981.
Chang, H., Paksoy, V. E., and Zhang, F. Polytopes ofstochastic tensors. Annals of Functional Analysis, 7(3):386–393, 2016.
Cui, L.-B., Li, W., and Ng, M. K. Birkhoff–von Neumanntheorem for multistochastic tensors. SIAM Journal on
Matrix Analysis and Applications, 35(3):956–973, 2014.
Cuturi, M. Sinkhorn distances: Lightspeed computation of
optimal transport. In Advances in Neural Information
Processing Systems 26, pp. 2292–2300, 2013.
Frogner, C., Zhang, C., Mobahi, H., Araya, M., and Poggio, T. A. Learning with a Wasserstein loss. In Advancesin Neural Information Processing Systems 28, pp. 2053–
2061, 2015.
Ganmor, E., Segev, R., and Schneidman, E. Sparse loworder interaction network underlies a highly correlated
and learnable neural population code. Proceedings ofLahr, M. and de Mesnard, L. Biproportional techniques
in input-output analysis: Table updating and structural
analysis. Economic Systems Research, 16(2):115–134,
2004.
Lamond, B. and Stewart, N. F. Bregman’s balancing
method. Transportation Research Part B: Methodological, 15(4):239–248, 1981.
Livne, O. E. and Golub, G. H. Scaling by binormalization.
Numerical Algorithms, 35(1):97–120, 2004.
Marshall, A. W. and Olkin, I. Scaling of matrices to achieve
specified row and column sums. Numerische Mathematik, 12(1):83–90, 1968.
Miller, R. E. and Blair, P. D. Input-Output Analysis: Foundations and Extensions. Cambridge University Press, 2
edition, 2009.
Moon, T. K., Gunther, J. H., and Kupin, J. J. Sinkhorn
solves sudoku. IEEE Transactions on Information Theory, 55(4):1741–1746, 2009.
Nakahara, H. and Amari, S. Information-geometric measure for neural spikes. Neural Computation, 14(10):2269–2316, 2002.
Nakahara, H., Nishimura, S., Inoue, M., Hori, G., andAmari, S. Gene interaction in DNA microarray data isdecomposed by information geometric measure. Bioinformatics, 19(9):1124–1131, 2003.
Nakahara, H., Amari, S., and Richmond, B. J. A comparison of descriptive models of a single spike train by
information-geometric measure. Neural computation, 18(3):545–568, 2006.

Tensor Balancing on Statistical Manifold

Parikh, A. Forecasts of input-output matrices using the
R.A.S. method. The Review of Economics and Statistics,
61(3):477–481, 1979.
Parlett, B. N. and Landis, T. L. Methods for scaling to doubly stochastic form. Linear Algebra and its Applications,
48:53–79, 1982.
Rao, S. S. P., Huntley, M. H., Durand, N. C., Stamenova,
E. K., Bochkov, I. D., Robinson, J. T., Sanborn, A. L.,
Machol, I., Omer, A. D., Lander, E. S., and Aiden, E. L.
A 3D map of the human genome at kilobase resolution
reveals principles of chromatin looping. Cell, 159(7):1665–1680, 2014.
Rota, G.-C. On the foundations of combinatorial theory I:
Theory of Möbius functions. Z. Wahrseheinlichkeitstheorie, 2:340–368, 1964.
Sinkhorn, R. A relationship between arbitrary positive matrices and doubly stochastic matrices. The Annals of
Mathematical Statistics, 35(2):876–879, 06 1964.
Sinkhorn, R. and Knopp, P. Concerning nonnegative matrices and doubly stochastic matrices. Pacific Journal of
Mathematics, 21(2):343–348, 1967.
Solomon, J., de Goes, F., Peyré, G., Cuturi, M., Butscher,A., Nguyen, A., Du, T., and Guibas, L. Convolutional
Wasserstein distances: Efficient optimal transportation
on geometric domains. ACM Transactions on Graphics,
34(4):66:1–66:11, 2015.
Soules, G. W. The rate of convergence of sinkhorn balancing. Linear Algebra and its Applications, 150:3–40,
1991.
Sugiyama, M., Nakahara, H., and Tsuda, K. Information
decomposition on structured space. In 2016 IEEE International Symposium on Information Theory, pp. 575–
579, July 2016.
Wu, H.-J. and Michor, F. A computational strategy to adjust
for copy number in tumor Hi-C data. Bioinformatics, 32(24):3695–3701, 2016.