Nov. 20–23, 2019
IBIS 2019

隣接代数と双対平坦構造を用いた学習杉山 麿人（国立情報学研究所，JST さきがけ研究者）

Matrix Balancing
p₁₁ p₁₂p₂₁ p₂₂

1/41

Matrix Balancing
Find r and s:(Make doublystochasticmatrix)=r₁ 0

p₁₁ p₁₂s₁ 00 r₂

p₂₁ p₂₂

0 s₂r₁s₁ p₁₁ r₁s₂ p₁₂
r₂s₁ p₂₁ r₂s₂ p₂₂∑j r₁sj p₁j = 1∑j r₂sj p₂j = 1∑i ris₁ pi₁ = 1 ∑i ris₂ pi₂ = 1

1/41

Sinkhorn-Knopp Algorithm• Alternately rescale all rows and columnsof a matrix 𝑃 to sum to 1
• Commonly used to compute entropy-regularized
Optimal transport (Wasserstein distance)

2/41

Revisit Matrix Balancing [ICML2017]
p₁₁p₁₂p₁₃

p₂₁p₂₂p₂₃

p₃₁p₃₂p₃₃
3/41

Introduce 𝜼η₁₂ η₁₃
η₁₁
p₁₁ p₁₂ p₁₃
η₂₁
p₂₁

p₂₂

p₂₃

η₃₁
p₃₁

p₃₂

p₃₃
4/41

Introduce 𝜼η₁₂ η₁₃
η₁₁
p₁₁ p₁₂ p₁₃
η₂₁
p₂₁

p₂₂

p₂₃

η₃₁
p₃₁

p₃₂

p₃₃
4/41

Introduce 𝜼η₁₂ η₁₃
η₁₁
p₁₁ p₁₂ p₁₃
η₂₁
p₂₁

p₂₂

p₂₃

η₃₁
p₃₁

p₃₂

p₃₃
4/41

Introduce 𝜼η₁₂ η₁₃
η₁₁
p₁₁ p₁₂ p₁₃
η₂₁
p₂₁

p₂₂

p₂₃

η₃₁
p₃₁

p₃₂

p₃₃
4/41

Introduce 𝜼
32121

η₂₁
p₂₁

p₂₂

p₂₃

η₃₁
p₃₁

p₃₂

p₃₃

η₁₂ η₁₃
η₁₁
p₁₁ p₁₂ p₁₃

4/41

Introduce 𝜽
p₁₁

p₁₂

p₁₃

p₂₁

p₂₂

p₂₃
θ₂₂ θ₂₃

p₃₁

p₃₂

p₃₃
θ₃₂ θ₃₃

θij = log pij– log pi–₁j – log pij–₁+ log pi–₁j–₁

5/41

Introduce 𝜽
p₁₁

p₁₂

p₁₃

p₂₁

p₂₂

p₂₃
θ₂₂ θ₂₃

p₃₁

p₃₂

p₃₃
θ₃₂ θ₃₃

θij = log pij– log pi–₁j – log pij–₁+ log pi–₁j–₁

5/41

Introduce 𝜽
p₁₁

p₁₂

p₁₃

p₂₁

p₂₂

p₂₃
θ₂₂ θ₂₃

p₃₁

p₃₂

p₃₃
θ₃₂ θ₃₃

θij = log pij– log pi–₁j – log pij–₁+ log pi–₁j–₁

5/41

Introduce 𝜽p₁₁

p₁₂

p₁₃

p₂₁

p₂₂

p₂₃
θ₂₂ θ₂₃

p₃₁

p₃₂

p₃₃
θ₃₂ θ₃₃

θij = log pij– log pi–₁j – log pij–₁+ log pi–₁j–₁

5/41

Balancing as Constraints on 𝜼 and 𝜽
32121

η₂₁
p₂₁

p₂₂

η₃₁
p₃₁

p₂₃
θ₂₂ θ₂₃p₃₂

p₃₃
θ₃₂ θ₃₃

η₁₂ η₁₃
η₁₁
p₁₁ p₁₂ p₁₃

θij = log pij– log pi–₁j – log pij–₁+ log pi–₁j–₁

Matrix balancing ⇔Satisfy ηi₁ = η₁i = 3 – i + 1
with keeping all θij
6/41

Natural Gradient• Given 𝑃 ∈ ℝ𝑛×𝑛 , introduce (𝜃, 𝜂) as∑∑∑∑log 𝑝𝑖𝑗 =𝜃 𝑘𝑙 , 𝜂 𝑖𝑗 =𝑝𝑘𝑙
𝑘≤𝑖 𝑙≤𝑗𝑘≥𝑖 𝑙≥𝑗• Let 𝐼 = {11, 12, … , 1𝑛, 21, … , 𝑛1}, 𝜽 = (𝜃 𝑖 )𝑇𝑖∈𝐼 , 𝜼 = (𝜂 𝑖 )𝑇𝑖∈𝐼• Using Fisher information matrix 𝐺 ∈ ℝ|𝐼|×|𝐼| s.t.
𝑔(𝑖𝑗)(𝑘𝑙) = 𝜂 max{𝑖,𝑘} max{𝑗,𝑙} − 𝜂 𝑖𝑗 𝜂 𝑘𝑙 , update formula is𝜽 next = 𝜽 − 𝐺 −1 (𝜼 − 𝜼correct )
7/41

108
106
104
102
100
100 500n106

Running time (sec.)

Number of iteration

Results on Hessenberg Matrix
1041021005000

10–2100 500
n> x1000
fasterNatural
BNEWT
Sinkhorn

5000
8/41

Introduce Partial Order Structure(p(s), θₛ, ηₛ)
p₁₁

p₁₂

p₁₃

p₂₁

p₂₂

p₂₃

p₃₁

p₃₂

p₃₃

s

9/41

Partially Ordered Sets (Posets)2{a,b,c}{a,b,c}ℕ²

Network{a,b} {a,c} {b,c}{a}{b} {c}
∅

10/41

Incidence Algebra• Incidence algebra is defined over a poset (𝑆, ≤)– (Closed) Interval [𝑎, 𝑏] = {𝑠 ∈ 𝑆 ∣ 𝑎 ≤ 𝑠 ≤ 𝑏}

• Members of the incidence algebra arefunctions 𝑓(𝑎, 𝑏) from intervals [𝑎, 𝑏] to a scalar with(𝑓 + 𝑔)(𝑎, 𝑏) = 𝑓(𝑎, 𝑏) + 𝑔(𝑎, 𝑏)∑(𝑓𝑔)(𝑎, 𝑏) =𝑓(𝑎, 𝑥)𝑔(𝑥, 𝑏) (convolution)
𝑎≤𝑥≤𝑏

11/41

Analogy to Matrix Multiplication• For [𝑎, 𝑏], define 𝒇, 𝒈 ∈ ℝ|[𝑎,𝑏]| as⎡𝑔(𝑎, 𝑎)⎤⎡𝑓(𝑎, 𝑎)⎤𝒇 = ⎢ ⋮ ⎥, 𝒈 = ⎢ ⋮ ⎥⎣𝑔(𝑏, 𝑎)⎦⎣𝑓(𝑎, 𝑏)⎦• For 𝑓 and 𝑔,(𝑓𝑔)(𝑎, 𝑏) = 𝒇𝑇 𝒈

12/41

Special Elements• Delta function 𝛿:
1 if 𝑎 = 𝑏
𝛿(𝑎, 𝑏) = {0 otherwise• Zeta function 𝜁:𝜁(𝑎, 𝑏) = {(integral)

1 if 𝑎 ≤ 𝑏0 otherwise• Möbius function 𝜇 = 𝜁 −1 : 𝜁𝜇 = 𝛿(differential)
13/41

Möbius Inversion Formula• Given a poset 𝑆, for any functions 𝑓, 𝑔 ∶ 𝑆 → ℝ,the Möbius inversion formula is given as∑∑∑⎧ 𝑔(𝑥) =𝜁(𝑠, 𝑥)𝑓(𝑠) =𝜁(𝑠, 𝑥)𝑓(𝑠) =𝑓(𝑠)⎪𝑠∈𝑆𝑠≤𝑆𝑠≤𝑥∑∑⎨ 𝑓(𝑥) =𝜇(𝑠, 𝑥)𝑔(𝑠) =𝜇(𝑠, 𝑥)𝑔(𝑠)⎪𝑠∈𝑆𝑠≤𝑆⎩
14/41

E.g.1: Inclusion-Exclusion Principle
A∪B∪CABCA∩B A∩C B∩CA∩B∩C|A∪B∪C|= |A|+|B|+|C|–|A∪B|–|A∪C|–|B∪C|+|A∩B∩C|f (X) = |X| g(X) = |X \ ⋃ Y|f (X) = ∑Y≤X g(X)g(X) = ∑Y≤X μ(Y,X) f (Y)

Y⊂X

15/41

E.g.2: Divisibility• Divisibility poset: 𝑎 ≤ 𝑏 ⇐⇒ 𝑏|𝑎 (𝑎 divides 𝑏)• The Möbius function: 𝑛 = 𝑏∕𝑎 and(−1)𝑘𝜇(𝑛) = {0if 𝑛 = 𝑝1 𝑝2 … 𝑝𝑘 for 𝑘 distinct primes
otherwise– 𝜇(𝑎, 𝑏) = 𝜇(𝑏|𝑎), the Möbius function in number theory• The Riemann zeta function 𝜁 is given by∑∞1∕𝜁(𝑠) =𝜇(𝑛)∕𝑛𝑠𝑛=1

16/41

Log-Linear Model on Poset [ICML2017]∑• For probability 𝑝∶𝑆 → (0, 1) with 𝑥∈𝑆 𝑝(𝑥) = 1,introduce 𝜽 and 𝜼 as∑𝜃𝑥 =𝜇(𝑠, 𝑥) log 𝑝(𝑠),𝑠∈𝑆∑∑𝜂𝑥 =𝜁(𝑥, 𝑠)𝑝(𝑠) =𝑝(𝑠)𝑠∈𝑆𝑠≥𝑥• From the Möbius inversion formula, log-linear model is:∑∑
log 𝑝(𝑥) =𝜁(𝑠, 𝑥)𝜃 𝑠 =𝜃𝑠𝑠∈𝑆𝑠≤𝑥

17/41

Log-Linear Model on Poset(p(x), θx, ηx)
x

18/41

Log-Linear Model on Poset(p(x), θx, ηx)log p(x) = ∑s≤x θsx

18/41

Log-Linear Model on Poset(p(x), θx, ηx)log p(x) = ∑s≤x θsxηx = ∑s≥x p(s)18/41

Exponential Family• The log-linear model on posets belongs tothe exponential family• 𝜽 : Natural parameter• 𝜼 : Expectation parameter

19/41

Binary Log-Linear Model [AAAI2019]2{a,b,c}{a,b,c}(= Boltzmann machine){a,b} {a,c} {b,c}{a}{b} {c}
∅

20/41

Binary Log-Linear Model [AAAI2019]{a,b,c}2{a,b,c}{a,b} {a,c} {b,c}{a}{b} {c}∅(= Boltzmann machine)log p(x) = ∑s≤x θs
ηx = ∑s≥x p(s)20/41Binary Log-Linear Model [AAAI2019]{0, 1}3111(= Boltzmann machine)

110 101 011
100 010 001
000
21/41

Binary Log-Linear Model [AAAI2019]{0, 1}3111

110 101 011(= Boltzmann machine)log p(x) = ∑s≤x θs= –ψ + ∑i θixi + ∑i,jθijxixj + ···100 010 001 For x with xi₁ = ··· = xik = 1,ηx = ∑s≥x p(s)= 𝔼[xi₁···xik] = Pr(xi₁ = ··· = xik = 1)
000

21/41

Dually Flat Structure• Let 𝜓(𝜃) = −𝜃(⊥)(convex, partition function)∑Legendre transformation
𝜓(𝜃) ,,,,,,,,,,,,,,,,,,,,,,→ 𝜙(𝜂) =𝑝(𝑥) log 𝑝(𝑥)𝑥∈𝑆

• (𝜓(𝜃), 𝜙(𝜂)) leads to dually flat coordinate system (𝜃, 𝜂):𝜕𝜓(𝜃) = 𝜂 𝑥∇𝜓(𝜃) = 𝜂,𝜕𝜃𝑥𝜕
∇𝜙(𝜂) = 𝜃,𝜙(𝜂) = 𝜃 𝑥
𝜕𝜂 𝑥

22/41

Riemannian Metric (Fisher Information)∑𝜕 𝜕𝜕
𝜓(𝜃) =𝜂𝑦 =𝜁(𝑥, 𝑠)𝜁(𝑦, 𝑠)𝑝(𝑠) − 𝜂 𝑥 𝜂 𝑦𝜕𝜃 𝑥 𝜕𝜃 𝑦𝜕𝜃𝑥𝑠∈𝑆∑𝜕 𝜕𝜕
𝜙(𝜂) =
𝜃𝑦 =𝜇(𝑠, 𝑥)𝜇(𝑠, 𝑦)𝑝(𝑠)−1𝜕𝜂 𝑥 𝜕𝜂𝑦𝜕𝜂𝑥𝑠∈𝑆
𝜕𝜕log 𝑝(𝑠)] = 𝛿(𝑥, 𝑦)
𝔼𝑠 [log 𝑝(𝑠)𝜕𝜂 𝑦𝜕𝜃 𝑥23/41

Riemannian Metric (Fisher Information)∑𝜕𝜕
𝔼𝑠 [log 𝑝(𝑠)log 𝑝(𝑠)] =𝜁(𝑥, 𝑠)𝜁(𝑦, 𝑠)𝑝(𝑠) − 𝜂 𝑥 𝜂 𝑦𝜕𝜃 𝑥𝜕𝜃 𝑦𝑠∈𝑆∑𝜕𝜕
𝔼𝑠 [log 𝑝(𝑠)log 𝑝(𝑠)] =𝜇(𝑠, 𝑥)𝜇(𝑠, 𝑦)𝑝(𝑠)−1𝜕𝜂 𝑥𝜕𝜂 𝑦𝑠∈𝑆𝜕𝜕
log 𝑝(𝑠)] = 𝛿(𝑥, 𝑦)
𝔼𝑠 [log 𝑝(𝑠)𝜕𝜂 𝑦𝜕𝜃 𝑥24/41

Mixed Coordinate System• Many problems are formulated as coordinate mixing

P = ( θ1, θ2, ..., θi–1, θi, θi+1, ..., θₙ )Q = ( η1, η2, ..., ηi–1, θi, θi+1, ..., θₙ )R = ( η1, η2, ..., ηi–1, ηi, ηi+1, ..., ηₙ )e-projection(MLE)
m-projection

Pythagorean theorem:(Q is always unique)
KL(P, R) = KL(P, Q) + KL(Q, R)

25/41

Mixed Coordinate System (Example)• Many problems are formulated as coordinate mixing

P = ( 0 , 0 , ..., 0 , 0 , 0, ..., 0 )

Q = ( η1, η2, ..., ηi–1, 0 , 0, ..., 0 )

R = ( η1, η2, ..., ηi–1, ηi, ηi+1, ..., ηₙ )Uniform dist.
Empirical dist.

Pythagorean theorem:(Q is always unique)
KL(P, R) = KL(P, Q) + KL(Q, R)

26/41

Two Submanifolds
Q = (η1,η2,...,ηi–1,θi,θi+1,...,θₙ)Fix
P = (θ1,θ2,...,θi–1,θi,θi+1,...,θₙ)(Model manifold)e-projection
= MLE
Fix
R = (η1,η2,...,ηi–1,ηi,ηi+1,...,ηₙ) (Data manifold)
27/41

Gradient methods for e-projection• e-projection is convex optimization• Gradient descent (first-order):𝜽 next ← 𝜽 − 𝜀(𝜼 − 𝜼̂ target )• Natural gradient (second-order)
𝜽 next ← 𝜽 − 𝐺 −1 (𝜼 − 𝜼̂ target )– 𝐺 is Fisher information matrix w.r.t. 𝜃• Coordinate descent [IBIS2019]28/41

Boltzmann Machine Training
BoltzmannmachineSamplespace 1111
1110

1100 1010
1000 0100
0000

29/41

Boltzmann Machine Training
BoltzmannmachineSamplespace 1111

θ=011101100 1010
1000 0100
0000

η=η29/41

Matrix Balancing
3x3 matrix
as poset:

η=3–i+1

θ=θ30/41

Legendre Decomposition [NeurIPS 2018]
3x3x3 tensor
as poset:

31/41

Legendre Decomposition [NeurIPS 2018]
3x3x3 tensor
as poset:
θ=0
η=η

31/41

Introducing Hidden Variables in BM
Samplespace 1111
BM withhidden variables
1110
1100 1010
1000 0100
0000

32/41

Introducing Hidden Variables in BM
Samplespace 1111
BM withhidden variables1110

η0101 = ?

1100 1010
1000 0100
0000

32/41

Introducing Hidden Variables in BM
Samplespace 1111
BM withhidden variables
1110
1100 1010

θ=0
η=η1000 0100
Need to estimate η0000 via EM
nonconvex
32/41

Introducing Hidden States
BoltzmannmachineSamplespace 1111
1110

1100 1010
1000 0100
0000

33/41

Introducing Hidden StatesBoltzmann
machineSamplespace 1111
1110

1100 1010
1000 0100
0000

Hiddenstates
34/41

Introducing Hidden States
BoltzmannmachineSamplespace 1111

Optimization isstill convex!

1110

1100 1010
1000 0100
0000

Hiddenx states
ηx35/41

Introducing Hidden States
BoltzmannmachineSamplespace 1111

θ=011101100 1010
1000 0100
0000

η=η36/41

Blind Source Separation [arXiv]
Mixing Source Received
e.g. image

37/41

Blind Source Separation [arXiv]
Mixing Source Received
e.g. image

η=η

θ=037/41

Results of BSS
SourceReceived
Reconst.

Method RMSE
IGBSS0.27032
FastICA 0.43630
NMF0.62195
DicLearn 0.37167
38/41

Relationship to Homology• Möbius function is Euler characteristics• Consider order complex ∆(𝑆) of a poset 𝑆 with ⊥, ⊤ ∈ 𝑆(i) Vertices of ∆(𝑆) are elements of 𝑆(ii) Faces of ∆(𝑆) are chains of 𝑆

• For the Euler characteristic 𝜒(∆(𝑆)),𝜇(⊥, ⊤) = 𝜒(∆(𝑆)) + 1– Two spaces are homotopy equivalent⇒ Euler characteristics are the same

39/41

Summary: Recipe for Poset-LogLinear
1. Treat the target as a poset
2. Introduce the log-linear model on poset
3. Formulate the objective as coordinate mixing
4. Solve it by a gradient method(This slide is at https://mahito.nii.ac.jp/)40/41

Acknowledgment• Tsuda, K. (UTokyo), Nakahara, H. (RIKEN CBS)• Yamada, R. (KyotoU), Mimura, K. (HiroshimaCU)• Luo, S., Azizi, L. (USydney)• Hayashi, S., Matsushima, S. (UTokyo)• Borgwardt, K. and his lab members (ETH Zürich)• Lab members at NII
41/41