Tensor Balancing on Statistical Manifold
Mahito Sugiyama (NII, JST PRESTO), Hiroyuki Nakahara (RIKEN BSI), Koji Tsuda (UTokyo, RIKEN AIP, NIMS)
The 34th International Conference on Machine Learning (ICML 2017), August 6–11, 2017

Resultse-Projection and m-Projection

Log-Linear Model on Poset• Balancing of higher order (more than two) tensors is firstly(theoretically) achieved• Let (S, ≤) be a poset (= DAG)

Each x ∈ S has a triple:(p(x), θ(x), η(x))

Q e-projection• A probability vector p∶S → (0, 1)
s.t. ∑ x∈S p(x) = 1• A fast balancing algorithm with quadratic convergence using
Newton’s method (an existing algorithm is linear convergence)– (Normalized) weight for each node

• [Theory] We provide dually flat Riemannian manifold of
probability distributions with the structured outcome space

p21 p22 p23 p24

r2

p31 p32 p33 p34
p41 p42 p43 p44
a1

a2

a3

s≤x

Realize balancing as “projection”~10,000x faster!!

PBalancing = projectionr4a4

diag(r)Pdiag(s)n×n• Matrix balancing: Given a nonnegative matrix P = (p i j ) ∈ R+ ,nT(RPS) 1 = 1102
10–10 2
10–4
1050

500

5000

nNewton (proposed)

10108
106

Residual

Number of iterations

Running time (sec.)–– Applications: input-output analysis, Hi-C data analysis, the Sudoku
puzzle, and Wasserstein metric approximation– Standard balancing algorithm: Sinkhorn–Knopp algorithm

104

10–1–

10 3102

10–5

10

10–7
1050

Sinkhorn

n500

5000

BNEWT

NewtonSinkhorn
BNEWT0500 1000 1500 2000 2500 3000
Number of iterationsilog p(x) = ∑ θ x + ∑011i+θ001ij ii< j

1...n 1 2jθ x x +...
nx x . . . x − ψ,iiijii

η = E[x ] = Pr(x = 1),
jijη = E[x x ] = Pr(x = x = 1), . . .

• θ and η form a dual coordinate system: ∇ψ(θ) = η, ∇ϕ(η) = θ– ψ(θ) = −θ(⊥) = − log p(⊥), ϕ(η) = ∑ x∈S p(x) log p(x)– ψ(θ) and ϕ(η) are connected via the Legendre transformation:′(θη−ψ(θ)),ϕ(η) = max′θη(x) = β(x)+∀x ∈ B = S / AB

• e-projection and m-projection can be computed
by Newton’s Method• Each step of Newton’s method in e-projection:⎡⎢⎡⎢⎤⎥⎤⎥⋮⋮⎢⎢⎢⎢ (t)⎥⎥⎥⎥⎢⎢⎢⎢η (x) − β(x)⎥⎥⎥⎥⋮⎢⎢ Pβ⎥⎥ + J ⎢⎢ (t+1)⎥⎥ = 0,(t)⎢⎢θ⎢⎢⎥⎥⎥⎥(y)−θ(y)⋮⎢⎢ Pβ⎢⎢⎥⎥⎥⎥Pβ⎢⎢⎢⎢⎥⎥⎥⎥⋮⋮⎣⎣⎦⎦– J is the ∣B∣ × ∣B∣ Jacobian matrix given as

Dually Flat Structure′′
R = diag(r), S = diag(s), each entry is given as p i j = p i j r i s j

104010i000

r1a3p13 + r2a3p23 + r3a3p33 + r4a3p43 = 1
r2a1p21 + r2a2p22 + r2a3p23 + r2a4p24 = 1 Use geometric structure
in Information Geometry

find r, s ∈ R s.t.
(RPS)1 = 1 and100

101Ps≥x• Our model is generalization of
the log-linear model on binary vectors
nwith x ∈ {0, 1} = S :

111
110

r3A

• We introduce θ∶S → R and η∶S → R asMatrix (Tensor) Balancing
r1

Rm-projection

log p(x) = ∑ θ(s), η(x) = ∑ p(s)p11 p12 p13 p14

θ(x) = α(x)∀x ∈ A′θη=∑′x∈S\{⊥}θ (x)η(x)• The gradients: g(θ) = ∇∇ψ(θ) = ∇η, g(η) = ∇∇ϕ(η) = ∇θ⎧∂η(x)
⎪⎪(θ)
g==
∑
ζ(x,s)ζ(y,s)p(s)−η(x)η(y)⎪xy⎪∂θ(y) s∈S⎪⎪⎪⎨⎪⎪∂θ(x)−1⎪⎪g(η)==∑µ(s,x)µ(s,y)p(s)
⎪xy⎪
∂η(y)
⎪s∈S
⎩

– Zeta function: ζ(s, x) = 1 if s ≤ x , ζ(s, x) = 0 otherwise– Möbius function: µ∶S × S → Z satisfying ζ µ = I– The manifold (S , g(ξ)) is a Riemannian manifold withthe set S of probability vectors and the Riemannian metric g(ξ)

Jx y =(t)∂η Pβ (x)(t)∂θ Pβ (y)= ∑ ζ(x, s)ζ(y, s)p β(t)(t)(t)(s) − η Pβ (x)η Pβ (y)s∈Sfor each x, y ∈ BBalancing = e-Projection
p11 p12 p13 p14
p21 p22 p23 p24

η11
η21

p31 p32 p33 p34
p41 p42 p43 p44
Matrix balancing is achieved if:
η11 = 4, η21 = 3, η31 = 2, η41 = 1
η11 = 4, η12 = 3, η13 = 2, η41 = 1

η31
η41

432

η123

η132

η14

θ22

θ23θ24

θ32

θ33

θ34

θ42

θ43

θ44

1
Change η

1Fix θ