日本オペレーションズリサーチ学会

March 17, 2022

思う

私の数理工学

来嶋秀治
九州大学 大学院システム情報科学研究院

研究略歴
略歴
2002. 3. 東大 工 計数工学科 数理工学コース 卒
2004. 3. 同院 情報理工 数理情報学専攻 修士了2007. 3. 同 博士了
博士論文「MCMC法における近似精度保証と完璧サンプリング法」
2007. 4.– 2010.3. 京大 数理解析研 助教
2010. 4.- 現在 九大 システム情報科学研究院 准教授
2016.10 – 2020.3 JSTさきがけ研究員
領域「社会的課題の解決に向けた数学と諸分野の協働 (國府総括)」
研究分野
数理工学, 理論計算機科学
キーワード
 確率過程: ランダムウォーク, MCMC（マルコフ連鎖モンテカルロ）
 離散数学: 劣モジュラ関数，グラフ
 分散計算: 自律分散ロボット，自己安定2機関誌「オペレーションズ・リサーチ」, 334巻, 11号, 616頁, 1989. (伏見先生)

数理工学者・技術者は，広い工学的技能・常識，古典的近代的数学の基本的な考え方に通じているにと
どまらず，これを自由に駆使して工学の基礎的な問
題にも現実の実践的な問題にも深く切り込むことができなければならない．
---伏見正則

「数理工学 ------ 「工学」諸問題を数理する！」(文部科学省科学技術政策研究4
所 編著：「数学イノベーション」 第７章(2007年)の 元原稿) (藤重先生)「数理工学」とは，端的に述べると，「工学」諸問
題を数理する学問である．
---藤重悟今日の話題I. 私の数理工学: 理論編 (20分: 研究紹介)
II. 私の数理工学: 実践編 (15分)
III. 私の数理工学: 開拓編 (15分)

5I. 私の数理工学: 理論編
良い研究をして、良い論文を書く。

私の数理工学: 理論編 (研究宣伝) 完璧サンプリング法
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012]
 確率と計算 決定性RW [RSA 2015, TCS 2017, SIDMA 2018]
 高次元体積の決定性近似 [Algorithmica 2016, TCS 2020]

 対数優モジュラ分布
 劣モジュラ関数の変数変換 [SIDMA 2020] 任意の有限分配束は𝑀♮凹関数の最小化元族に同型 [OR lett. 2021]

 動的環境下のアルゴリズム論
 増えるクーポン収集 [SODA 2021]7

私の数理工学: 理論編 (研究宣伝)8 完璧サンプリング法 (学生～数理解析研)
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012][1] S. Kijima and Tomomi Matsui, Polynomial time perfect sampling algorithm for two 確率と計算
rowed
contingency tables, Random Structures & Algorithms, 29:2 (2006), 243--256.

[2] S. 
Kijimaand Tomomi
Matsui,Approximation
algorithm
決定性RW
[RSA
2015, TCS
2017, SIDMA
2018] and perfect sampler for closed
Jackson networks with single servers, SIAM Journal on Computing, 38:4 (2008), 1484--1503.
 高次元体積の決定性近似 [Algorithmica 2016, TCS 2020][3] S. Kijima and Tomomi Matsui, Randomized approximation scheme and perfect sampler for
closedJackson networks with multiple servers, Annals of Operations Research, 162 (2008),
 対数優モジュラ分布35--55.
 劣モジュラ関数の変数変換 [SIDMA 2020][4] Tomomi Matsui, Mitsuo Motoki, Naoyuki Kamatani and S. Kijima, Polynomial time
 任意の有限分配束は𝑀
凹関数の最小化元族に同型[OR lett.
2021]
approximate
or perfect samplers for ♮discretized
Dirichlet distribution, Japan
JournalofIndustrial
and Applied Mathematics, 27:1 (2010), 91--123.
 動的環境下のアルゴリズム論[5] S. Kijima, Sampling from log-super/submodular distributions, in Proc. 7th Hungarian 増えるクーポン収集[SODA 2021] and Its Applications (2011), 227--236.
Japanese
Symposium on Discrete Mathematics

9

例: 2行分割表のランダム生成
2元分割表
 各セルには非負整数が入る （与えられた）周辺和を満たす

1218
5 4 3 7 5 6 30

5 4 3 0 0 0 12
0 0 0 7 5 6 18
5 4 3 7 5 6 30
状態A

0 0 0 1 5 6 12
5 4 3 6 0 0 18
5 4 3 7 5 6 30
状態C

4 3 1 3 1 0 12
1 1 2 4 4 6 18
5 4 3 7 5 6 30
状態B

問題
Given: 周辺和
出力: 分割表の一様ランダム生成

10

例: 2行分割表のランダム生成
2元分割表
 各セルには非負整数が入る （与えられた）周辺和を満たす

1218
5 4 3 7 5 6 30与えられた周辺和を満たす2行分割表の個数を求める問題
#P完全 (NP困難) [‘97 Dyer, Kannan, & Mount]問題
Given: 周辺和
出力: 分割表の一様ランダム生成

マルコフ連鎖を用いた
サンプリング法

11

例: 2行分割表に対するマルコフ連鎖 [K & Matsui ‘06]
𝑗 𝑗+1

K-M chain
1. 𝑗 ∈ {1, … , 𝑛 − 1} を 一様ランダムに選ぶ．
𝑥1,𝑗 𝑥1,𝑗+1
2.
𝑥2,𝑗 𝑥2,𝑗+1 を実行可能な中から一様ランダムに選ぶ．

2 3 5
1 4 5
3 7 104 3 2 3 0 0 12
1 1 1 4 5 6 18
5 4 3 7 5 6 30
1/4
1/4

+

+𝑘 −𝑘
−𝑘 +𝑘

3 2
0 5

2 3
1 4

1/4

1/4

1 4
2 3

定理
このマルコフ連鎖 の定常分布は一様分布．

0 5
3 2

12

MCMC法のアイデア
1. 所望の分布を定常分布にもつマルコフ連鎖を設計する．
2. 十分な回数推移させて，定常分布からサンプリングする．
状態空間

何回?
x
x
x

・・・

x

x

 (漸近的に)所望の分布に従うサンプリングを実現

問題点は何か？
「何回推移させれば十分か？」近似サンプリング法
 収束スピード(mixing time)の算定
完璧サンプリング法
 マルコフ連鎖の推移シミュレーションを工夫 無限回の推移の結果を出力 ( 定常分布に厳密に従う)• 過去からのカップリング [Propp & Wilson 1996]

13

14

過去からのカップリング[Propp & Wilson 1996]
初期状態を知る必要がない !!

4



1

6

3

2



5

s1

s1

s1

s1

s1

s1

s1
uniquestate at 0

s2

s2

s2

s2

s2

s2

ss22

s3

s3

s3

s3

s3

s3

s3 coalesce

–6

–5

–4

–3

–2

–1

0

時刻

15

2行分割表の完璧サンプリング法
2元分割表
 各セルには非負整数が入る （与えられた）周辺和を満たす

1218
5 4 3 7 5 6 30

[K & Matusi 2006]
2行分割表の一様分布に対して，過去からのカップリングを用いたO 𝑛3 log 𝑁 時間の完璧サンプリング法を設計．
[1] S. Kijima and Tomomi Matsui, Polynomial time perfect samplingalgorithm for two-rowed contingency tables, Random Structures &Algorithms, 29:2 (2006), 243--256.

問題
Given: 周辺和
出力: 分割表の一様ランダム生成

私の数理工学: 理論編 (研究紹介)16 完璧サンプリング法 (学生～数理解析研)
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012][1] S. Kijima and Tomomi Matsui, Polynomial time perfect sampling algorithm for two-rowed 確率と計算
contingency
tables, Random Structures & Algorithms, 29:2 (2006), 243--256.

[2] S. 
Kijimaand Tomomi
Approximationand perfect sampler for closed
決定性RW[RSAMatsui,
2015, TCS
2017, SIDMAalgorithm
2018]Jackson networks with single servers, SIAM Journal on Computing, 38:4 (2008), 1484--1503.
 高次元体積の決定性近似 [Algorithmica 2016, TCS 2020][3] S. Kijima and Tomomi Matsui, Randomized approximation scheme and perfect sampler
forJackson networks with multiple servers, Annals of Operations Research, 162closed
対数優モジュラ分布(2008), 35--55.
 劣モジュラ関数の変数変換 [SIDMA 2020][4] Tomomi Matsui, Mitsuo Motoki, Naoyuki Kamatani and S. Kijima, Polynomial time
 任意の有限分配束は𝑀
凹関数の最小化元族に同型[OR lett.
2021]
approximate
or perfect samplers for ♮discretized
Dirichlet distribution, Japan
JournalofIndustrial
and Applied Mathematics, 27:1 (2010), 91--123.
 動的環境下のアルゴリズム論[5] S. Kijima, Sampling from log-super/submodular distributions, in Proc. 7th Hungarian 増えるクーポン収集[SODA 2021] and Its Applications (2011), 227--236.
Japanese
Symposium on Discrete Mathematics

17

待ち行列ネットワーク• 顧客はノードでサービスを受ける• サービス終了後，顧客は次のノードへ確率的に移動する．

どれくらい待つのだろう?
 待ち行列の長さ: 顧客 (ジョブ) ノードでの滞在時間: ノード (サーバ)

Everybody waits for Space Mountain, the Haunted Mansion, Jungle cruise, and so on …

閉ジャクソンネットワーク (単一サーバ)
𝐾: ネットワーク中のジョブの総数
1, 2, … , 𝑛 : ノードの集合
𝜇 = 𝜇1 , 𝜇2 , … , 𝜇𝑛 , 𝜽 = 𝜃1 , 𝜃2 , … , 𝜃𝑛 : パラメータ

1. 閉ネットワーク (= ジョブの総数 K は一定)
2. 各ノードでは 先着順 (First Come First Served).
3. 各ノードには1個のサーバ (単一サーバモデル).
4. ノード i でのサービス時間は期待値i の指数分布.
5. ノード i でサービス終了したジョブは, 確率 Pij でノード j へ.
セントラルサーバモデルなどに応用

18

閉ジャクソンネットワーク (単一サーバ) の均衡分布
𝐾: ネットワーク中のジョブの総数
1, 2, … , 𝑛 : ノードの集合
𝜇 = 𝜇1 , 𝜇2 , … , 𝜇𝑛 , 𝜽 = 𝜃1 , 𝜃2 , … , 𝜃𝑛 : パラメータ
ネットワーク中のジョブの均衡分布

単体中の
整数格子点

状態空間:
確率関数:

ただし

正規化定数

積形式解
(product form solution)19

閉ジャクソンネットワーク (単一サーバ) の均衡分布
𝐾: ネットワーク中のジョブの総数
1, 2, … , 𝑛 : ノードの集合
𝜇 = 𝜇1 , 𝜇2 , … , 𝜇𝑛 , 𝜽 = 𝜃1 , 𝜃2 , … , 𝜃𝑛 : パラメータ
ネットワーク中のジョブの均衡分布

単体中の
整数格子点

[K & Matusi 2008]
状態空間:

閉ジャクソンネットワークの積形式解に対して，過去からの確率関数:
カップリングを用いたO 𝑛3 log 𝑁 時間の完璧サンプリング法を設計．
[2] S. Kijima and Tomomi Matsui, Approximation
algorithm and perfect
積形式解
sampler for closed Jackson networks with single servers, SIAM Journal
ただし
on Computing, 38:4 (2008), 1484--1503.(product form solution)[3] S. Kijima and Tomomi Matsui, Randomized approximation schemeand perfect sampler for closed Jackson networks with multiple servers,
Annals of Operations Research, 162 (2008), 35--55.

正規化定数

20

私の数理工学: 理論編 (研究紹介)21 完璧サンプリング法 (学生～数理解析研)
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012][1] S. Kijima and Tomomi Matsui, Polynomial time perfect sampling algorithm for two-rowed 確率と計算
contingency
tables, Random Structures & Algorithms, 29:2 (2006), 243--256.

[2] S. 
Kijimaand Tomomi
Matsui,Approximation
algorithm
決定性RW
[RSA
2015, TCS
2017, SIDMA
2018] and perfect sampler for closed
Jackson networks with single servers, SIAM Journal on Computing, 38:4 (2008), 1484--1503.
 高次元体積の決定性近似 [Algorithmica 2016, TCS 2020][3] S. Kijima and Tomomi Matsui, Randomized approximation scheme and perfect sampler for
closedJackson networks with multiple servers, Annals of Operations Research, 162 (2008),
 対数優モジュラ分布35--55.
 劣モジュラ関数の変数変換 [SIDMA 2020][4] Tomomi Matsui, Mitsuo Motoki, Naoyuki Kamatani and S. Kijima, Polynomial time
 任意の有限分配束は𝑀
凹関数の最小化元族に同型[OR lett.
2021]
approximate
or perfect samplers for ♮discretized
Dirichlet distribution, Japan
JournalofIndustrial
and Applied Mathematics, 27:1 (2010), 91--123.
 動的環境下のアルゴリズム論[5] S. Kijima, Sampling from log-super/submodular distributions, in Proc. 7th Hungarian 増えるクーポン収集[SODA 2021] and Its Applications (2011), 227--236.
Japanese
Symposium on Discrete Mathematics

完璧サンプリングと単調優モジュラ分布

定理 [5]
分配束上の reversible Hasse walk が 単調な更新関数 を持つ
 定常分布が対数優モジュラ[1--4]よりも素朴なマルコフ連鎖．
特に[4](離散化Dirichlet分布)は(一般に)対数優モジュラではない．
[5] S. Kijima, Sampling from log-super/submodular distributions, in Proc. 7th HungarianJapanese Symposium on Discrete Mathematics and Its Applications (2011), 227--236.

22私の数理工学: 理論編 (研究宣伝) 完璧サンプリング法
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012]
 確率と計算 決定性RW [RSA 2015, TCS 2017, SIDMA 2018][6]S. Kijima
and Toshio Nemoto, [Algorithmica
On randomized
approximation
for finding
高次元体積の決定性近似
2016,
TCS 2020]
a level ideal of a poset and the generalized median stable matchings,
Mathematics
対数優モジュラ分布of Operations Research, 37:2 (2012), 356--371.
 劣モジュラ関数の変数変換 [SIDMA 2020] 任意の有限分配束は𝑀♮凹関数の最小化元族に同型 [OR lett. 2021]

 動的環境下のアルゴリズム論
 増えるクーポン収集 [SODA 2021]23

安定結婚問題[Gale & Shapley ‘62]

2012年ノーベル経済学賞
Shapley and Rotha. Preference list

favor

Men’s lists

m1
m2
m3
m4
m5

w1
w2
w3
w4
w5

w2
w1
w4
w3
w3

w4
w3
w2
w5
w1

w3
w4
w5
w1
w2

w5
w5
w1
w2
w4

Women’s lists

favor

w1
w2
w3
w4
w5

m3
m2
m3
m4
m1

men
1

women
1

2

2

3

3

4

4

5

5

m2
m3
m5
m1
m4

m1
m1
m4
m3
m5

m4
m4
m2
m5
m2

m5
m5
m1
m2
m3

24

25

安定マッチング“駆け落ち”を誘引しないマッチング
a. Preference list
Men’s lists

m1
m2
m3
m4
m5

w1
w2
w3
w4
w5

Women’s lists

w2
w1
w4
w3
w3

w4
w3
w2
w5
w1

w3
w4
w5
w1
w2

w5
w5
w1
w2
w4

w1
w2
w3
w4
w5

men
1

women
1

2

2

3

3

4

4

5

5

m2
m3
m5
m1
m4

m1
m1
m4
m3
m5

m3
m2
m3
m4
m1

m4
m4
m2
m5
m2

m5
m5
m1
m2
m3

26

安定マッチングは一つとは限らない
a. Preference list
Men’s lists

m1
m2
m3
m4
m5

w1
w2
w3
w4
w5

Women’s lists

w2
w1
w4
w3
w3

w4
w3
w2
w5
w1

w3
w4
w5
w1
w2

w5
w5
w1
w2
w4

w1
w2
w3
w4
w5

b. Stable matchings (from men’s view)1 2 3 4 5 6 7 8
m1 1 2 1 2 1 4 2 4
m2 2 1 2 1 2 1 1 1
m3 3 3 4 4 4 2 4 2
m4 4 4 3 3 5 3 5 5
m5 5 5 5 5 3 5 3 3

m2
m3
m5
m1
m4

m1
m1
m4
m3
m5

m3
m2
m3
m4
m1

m4
m4
m2
m5
m2

m5
m5
m1
m2
m3

men
1

women
1

2

2

3

3

4

4

5

5

[Conway (Knuth ’77,‘91)]

安定マッチングは分配束をなす
8
41253

6

b. Stable matchings (from men’s view)27
1 2 3 4 5 6 7 8
m1 1 2 1 2 1 4 2 4
m2 2 1 2 1 2 1 1 1
m3 3 3 4 4 4 2 4 2
m4 4 4 3 3 5 3 5 5
m5 5 5 5 5 3 5 3 3

7

41235

favor for men

favor for women

21453

Every finite distributive lattice is a set of stable matchings
4
5
任意の有限分配束は安定マッチングに同型
12453

21435

[Blair 1984]
2

3

21345

12435
1

(Hasse 図)

12345

男女の間で
の“公平性”

[Teo & Sethuraman 1998]

abcd
41253

abd
21453

メディアン安定マッチングbd
男女間で“公平”な安定マッチング
ab

1
2
3
4
5
1

1
2
3
4
5
2
1

1
2
4
3
5
3

2
1
4
3
5
4

3

2
1
4
3
5
5
4

2
1
4
5
3
6

4
1
2
5
3
7

7

favor for men

m1
m2
m3
m4
m5

favor for women

メディアン安定マッチング

abc
41235

28

c. GMSMs
4
1
2
5
3
8

8

12453
21435
Teo, Chung-Piaw; Sethuraman, Jay (1998), "The geometry of fractional stablematchings and its applications", Mathematics of Operations Research, 23 (4): 874–891.

a

b

21345

12435

12345

男女の間で
の“公平性”

[6] S. Kijima and Toshio Nemoto, On randomized
29
approximation for finding a level ideal of a poset
メディアン安定マッチングの計算量
and the generalized median stable matchings,
Mathematics of Operations Research, 37:2
定理 [Cheng ‘08]
(2012), 356--371.

i = (N) の時, 「i 番目の LI は要素 e を含むか？」は#P困難.
定理 [Cheng ‘08]
i = O(log |P|) の時, i 番目の LI は多項式時間で見つかる.
NP問題？

定理 [K, Nemoto ‘12]
「イデアル S は i 番目の LI か？」は,i = O(N1/c)の場合でも #P困難．但し c  1 は任意の定数.
定理 [K, Nemoto ‘12]
i = O(|P|c) の時, i 番目の LI は多項式時間で見つかる.
O(log log N)O((log N)c)

poly timepoly time((log N)c)&o(N1/c)

未解決

O(N1/c)

O(N)

#P困難

#P-hardi

私の数理工学: 理論編 (研究宣伝) 完璧サンプリング法
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012]
 確率と計算 (2010-)
 決定性RW [RSA 2015, TCS 2017, SIDMA 2018]
 高次元体積の決定性近似 [Algorithmica 2016, TCS 2020]

S.対数優モジュラ分布[7]
Kijima, Kentaro Koga and Kazuhisa Makino, Deterministic random walks on finite
graphs,
Structures & Algorithms, 46:4
(July2020]2015), 739--761.
 Random
劣モジュラ関数の変数変換[SIDMA♮凹関数の最小化元族に同型[8] Takeharu
Shiraga, Yukiko Yamauchi,
S. Kijima and Masafumi Yamashita, 任意の有限分配束は𝑀[ORDeterministic
lett. 2021]
random walks for rapidly mixing chains, SIAM Journal on Discrete Mathematics, 32:3
 動的環境下のアルゴリズム論(2018),
2180--2193.
 増えるクーポン収集
[SODA
2021]
[9] Takeharu
Shiraga, Yukiko Yamauchi,
Shuji
Kijima and Masafumi Yamashita, Total
variation discrepancy of deterministic random walks for ergodic Markov chains,

Theoretical Computer Science, 699 (7 November 2017), 63--74.

30

31

確率と計算

「計算にランダムネスは必要か？」32

ランダムウォーク
トークンがグラフ上をランダムウォークする.
 0: 初期分布 (トークンは確率0vで頂点vに居る)
 P: 推移確率行列 (確率Puvで頂点u頂点vに移動)
 時刻tの確率分布t:=0Pt (頂点vに居る確率 (t)v )
.3

.2
d

a
.4

e

.1

b

c

33

ランダムウォーク (複数トークンによる分布の近似)
N個のトークンがグラフ上を独立にランダムウォークする.
 0: 初期配置 (0  0/N)
 P: 推移確率行列 (確率Puvで頂点u頂点vに移動)
 時刻tの期待配置t := 0Pt (t  t/N)(Nが非常に大きければ, )トークンを約2:4:1:3の割合で隣接点にばらまいている.

.3

.2
d

a

1/3
.4

1/3

.1
決定性過程でも大差ない?

e
1/3

b

c

34

決定性RW (Propp機械; rotor-router)
N個のトークンがグラフ上を移動する.
 0: 初期配置 (0 = 0) : “rotor router” (比率Puvで頂点u頂点vに“順次”移動) 時刻tの配置t (t  t ???)
a=<d,b,c,a,a,a,d,c,c,c>
.2
.3
d
a
.4
.1

省略

t=0
e=<c,d,e>
1/3
1/3
e
1/3

b

c

(𝑇)

主結果

𝜒𝑤 : #tokens @ 𝑤 ∈ 𝑉, @ time 𝑇 > 0 in detRW 35(𝑇)

𝜇𝑤 : E[#tokens] @ 𝑤 ∈ 𝑉, @ time 𝑇 > 0 in RW

定理 [K, Koga, Makino 10+]省略
対応する推移確率行列 𝑃 の固有値がすべて非負ならば,任意の多重有向グラフ, 任意の初期状態, 任意のrotor-router,
任意の頂点 𝑤, 任意の時刻 𝑡 について,𝑇𝑇
𝜒𝑤 − 𝜇𝑤 ≤ 4𝑚𝑛 − O(𝑚)但し, 𝑛, 𝑚 は多重グラフの頂点数,枝数.

Remark「推移確率行列Pの固有値がすべて非負」という条件.
 reversible lazy Markov chainはこの条件を満たす.
 MCMC法で使われるマルコフ連鎖
 +Pが対称  Pは半正定値行列

(𝑇)𝜒𝑤 : #tokens @ 𝑤 ∈ 𝑉, @ time 𝑇 > 0 in detRW 36

主結果

(𝑇)

𝜇𝑤 : E[#tokens] @ 𝑤 ∈ 𝑉, @ time 𝑇 > 0 in RW

定理 [K, Koga, Makino 10+]省略
対応する推移確率行列 𝑃 の固有値がすべて非負ならば,任意の多重有向グラフ, 任意の初期状態, 任意のrotor-router,
任意の頂点 𝑤, 任意の時刻 𝑡 について,𝑇𝑇
𝜒𝑤 − 𝜇𝑤 ≤ 4𝑚𝑛 − O(𝑚)但し, 𝑛, 𝑚 は多重グラフの頂点数,枝数.
トークン数を増やせば，
deterministic random walk (決定性過程)は(𝑇)
𝑇
Remark𝜒𝑤
𝜇𝑤
1 𝑁→∞
ランダムウォークをよく模倣する．
−=O
0「推移確率行列Pの固有値がすべて非負」という条件.
𝑁𝑁𝑁
 reversible lazy Markov chainはこの条件を満たす.
 MCMC法で使われるマルコフ連鎖Rem.
 +Pが対称  Pは半正定値行列
𝜇 (𝑇)
(𝑇)
≃ 𝜋 : 時刻 𝑇 の分布

𝑁

私の数理工学: 理論編 (研究宣伝) 完璧サンプリング法
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012]
 確率と計算 (2010-)
 決定性RW [RSA 2015, TCS 2017, SIDMA 2018]
 高次元体積の決定性近似 [Algorithmica 2016, TCS 2020]

 対数優モジュラ分布
劣モジュラ関数の変数変換[10] EiAndoand S. Kijima, An FPTAS for the [SIDMA
volume 2020]
computation of 0-1 knapsack polytopes♮凹関数の最小化元族に同型
based on approximate
convolution, Algorithmica,
76:4 (December 2016),[OR1245--1263.
任意の有限分配束は𝑀
lett. 2021][11]
Ando and S. Kijima, An FPTAS for the volume of some V-polytopes --- It is hard to Ei動的環境下のアルゴリズム論
compute the volume of the intersection of two cross-polytopes, Theoretical Computer 833:12増えるクーポン収集
2021]
Science,(September 2020), [SODA87--106.

37高次元多面体の体積計算の決定性近似
Lovász [1986]“membership oracleで与えられた 𝑛次元凸体の体積に対して，
近似比 1.999𝑛 の決定性多項式時間アルゴリズムは存在しない．”
乱択近似 (1980～)

 Jerrum, Valiant, Vazirani (1986): 近似数え上げ≃ランダムサンプリング Toda (1991): PH ⊆ P #P (Gödel賞)
 Dyer, Frieze, Kanna (1991): 高次元凸体体積の乱択近似可能 (Fulkerson賞)
 Jerrum, Sinclair, Vigoda (2004): パーマネントの乱択近似可能 (Fulkerson賞)

決定性近似 (話者の研究 2014～)
 ナップサック多面体の体積は決定性近似可能(FPTAS) [Algorithmica16]• H多面体: 線形不等式系で記述される多面体
 ‘双対’ナップサック多面体の体積は決定性近似可能(FPTAS) [TCS 20]• V多面体: 端点集合で記述される多面体•

𝑛次元L1球2個の交差領域(resp. union)の体積計算は#P困難

38

私の数理工学: 理論編 (研究宣伝) 完璧サンプリング法
 分割表 [RSA 2006]
待ち行列ネットワーク[SICOMP
2008,andAORMasafumi
2008] Yamashita, Finding[12]
JunpeiNakashima, Yukiko Yamauchi,
S. Kijima
submodularity hidden in symmetric difference, SIAM Journal on Discrete Mathematics,

メディアン安定マッチング [MathOR 2012]
34:1 (2020), 571--585.

確率と計算
[13]
Tomohito Fujii and Shuji Kijima, Every finite distributive lattice is isomorphic to the♮-concave set function, Operations Research Letters, 49:1 (January
minimizer
set of an M
 決定性RW[RSA
2015, TCS 2017, SIDMA 2018]2021), 1--4.
 高次元体積の決定性近似 [Algorithmica 2016, TCS 2020]

 対数優モジュラ分布 (2016-2020)
 劣モジュラ関数の変数変換 [SIDMA 2020] 任意の有限分配束は𝑀♮凹関数の最小化元族に同型 [OR lett. 2021]

 動的環境下のアルゴリズム論
 増えるクーポン収集 [SODA 2021]

39

40

凸性はアフィン変換で保存される
定義：凸関数関数𝑓: ℝ𝑛 → ℝに対し，任意の𝒙, 𝒚 ∈ ℝ𝑛 と任意の𝜆 ≥ 0に対して
𝜆𝑓 𝒙 + 1 − 𝜆 𝑓 𝒚 ≥ 𝑓(𝜆𝒙 + 1 − 𝜆 𝒚)が成り立つとき，𝑓は凸関数という．

𝑓 𝑥 = 𝑥2
𝑓 𝑥 = 1 − exp

凸関数

凸ではない関数

𝑥2
−
2

41

問題

𝑓 𝑥, 𝑦 = 2𝑥 2 + 2𝑥𝑦 + 5𝑦 2は凸関数か？

42問題

𝑓 𝑥, 𝑦 = 2𝑥 2 + 2𝑥𝑦 + 5𝑦 2は凸関数か？
解答
変数変換
𝑥 ≔ 𝑠 + 2𝑡𝑦 ≔𝑠−𝑡を施す．
𝑔 𝑠, 𝑡 = 𝑓(𝑠 + 2𝑡, 𝑠 − 𝑡)
= 2 𝑠 + 2𝑡 2 + 2 𝑠 + 2𝑡 𝑠 − 𝑡 + 5 𝑠 − 𝑡 2= 2𝑠 2 + 8𝑠𝑡 + 8𝑡 2 + 2𝑠 2 + 2𝑠𝑡 − 4𝑡 2 + 5𝑠 2 − 10𝑠𝑡 + 5𝑡 2
= 9𝑠 2 + 9𝑡 2

𝒈(𝒔, 𝒕)は凸．

43

問題

𝑓 𝑥, 𝑦 = 𝑥 2 + 4𝑥𝑦 + 3𝑦 2は凸関数か？

44問題

𝑓 𝑥, 𝑦 = 𝑥 2 + 4𝑥𝑦 + 3𝑦 2は凸関数か？
解答
変数変換
𝑥 ≔ 𝑠 + 3𝑡
𝑦 ≔𝑠−𝑡を施す．
𝑔 𝑠, 𝑡 = 𝑓(𝑠 + 3𝑡, 𝑠 − 𝑡)
= 𝑠 + 3𝑡 2 + 4 𝑠 + 3𝑡 𝑠 − 𝑡 + 3 𝑠 − 𝑡 2= 𝑠 2 + 6𝑠𝑡 + 9𝑡 2 + 4𝑠 2 + 8𝑠𝑡 − 12𝑡 2 + 3𝑠 2 − 6𝑠𝑡 + 3𝑡 2
= 8𝑠 2 + 8𝑠𝑡𝒈(𝒔, 𝒕)は凸ではない．
e.g., 𝑔 1, −2 = 8 − 16 = −8, 𝑔 −1,2 = 8 − 16 = −8,𝑔1
1(1, −2) + (−1,2)
2
2

1
2

1
2

= 𝑔 0,0 = 0 > 𝑔 1, −2 + 𝑔 −1,2 .

凸性はアフィン変換で保存される
定義：凸関数関数𝑓: ℝ𝑛 → ℝに対し，任意の𝒙, 𝒚 ∈ ℝ𝑛 と任意の𝜆 ≥ 0に対して
𝜆𝑓 𝒙 + 1 − 𝜆 𝑓 𝒚 ≥ 𝑓(𝜆𝒙 + 1 − 𝜆 𝒚)が成り立つとき，𝑓は凸関数という．

定理 cf. [Rockafellar]ℎ: ℝ𝑛 → ℝ𝑛 はアフィン写像とする．
関数𝑓: ℝ𝑛 → ℝが凸のとき𝑔 ≔ 𝑓 ∘ ℎは凸関数．
(i.e., ℎ = 𝐴𝒙 + 𝒃: 𝐴 ∈ ℝ𝑛×𝑛 , 𝒃 ∈ ℝ𝑛 のとき，

𝑔 𝒙 = 𝑓(𝐴𝒙 + 𝒃)は凸関数)

45

凸性はアフィン変換で保存される
定義：凸関数関数𝑓: ℝ𝑛 → ℝに対し，任意の𝒙, 𝒚 ∈ ℝ𝑛 と任意の𝜆 ≥ 0に対して
𝜆𝑓 𝒙 + 1 − 𝜆 𝑓 𝒚 ≥ 𝑓(𝜆𝒙 + 1 − 𝜆 𝒚)離散関数のアナロジー？

が成り立つとき，𝑓は凸関数という．

劣モジュラ関数

定理 cf. [Rockafellar]ℎ: ℝ𝑛 → ℝ𝑛 はアフィン写像とする．
関数𝑓: ℝ𝑛 → ℝが凸のとき𝑔 ≔ 𝑓 ∘ ℎは凸関数．
(i.e., ℎ = 𝐴𝒙 + 𝒃: 𝐴 ∈ ℝ𝑛×𝑛 , 𝒃 ∈ ℝ𝑛 のとき，

𝑔 𝒙 = 𝑓(𝐴𝒙 + 𝒃)は凸関数)

46定義: 劣モジュラ
集合関数 𝑓: 0,1 𝑉 → ℝ に対して，Φ𝑓 : 0,1 𝑉 × 0,1 𝑉 → ℝを
Φ𝑓 𝑋, 𝑌 ≔ 𝑓 𝑋 + 𝑓 𝑌 − 𝑓 𝑋 ∪ 𝑌 − 𝑓 𝑋 ∩ 𝑌とする．
• 集合関数𝑓が劣モジュラとは，任意の𝑋, 𝑌 ∈ 0,1 𝑉 に対して
Φ𝑓 𝑋, 𝑌 ≥ 0が成り立つこととする．
定義: 対称差変換

47

𝑓 1,1,1 = 1対称差写像𝜎𝑠 : 0,1 𝑉 → 0,1 𝑉 を
𝜎𝑠 𝑋 ≔ 𝑋 ⊕ 𝑆 (𝑋 ∈ 0,1 𝑉 )とする．

𝑓 1,1,0 = 2

𝑓 1,0,1 = 0

𝑓 0,1,1 = 2

𝑓 1,0,0 = 1

𝑓 0,1,0 = 1

𝑓 0,0,1 = 1集合関数𝑔: 0,1 𝑉 → ℝ を𝑔 ≔ 𝑓 ∘ 𝜎𝑆を 𝑆 による 𝑓 の対称差変換と呼ぶ．
𝑓 0,0,0 = 0定義: 劣モジュラ
集合関数 𝑓: 0,1 𝑉 → ℝ に対して，Φ𝑓 : 0,1 𝑉 × 0,1 𝑉 → ℝを
Φ𝑓 𝑋, 𝑌 ≔ 𝑓 𝑋 + 𝑓 𝑌 − 𝑓 𝑋 ∪ 𝑌 − 𝑓 𝑋 ∩ 𝑌とする．
• 集合関数𝑓が劣モジュラとは，任意の𝑋, 𝑌 ∈ 0,1 𝑉 に対して
Φ𝑓 𝑋, 𝑌 ≥ 0が成り立つこととする．

48

𝑔 ≔ 𝑓 ∘ 𝜎(1,0,0) の場合

定義: 対称差変換

𝑔 0,1,1 = 1対称差写像𝜎𝑠 : 0,1 𝑉 → 0,1 𝑉 を
𝜎𝑠 𝑋 ≔ 𝑋 ⊕ 𝑆 (𝑋 ∈ 0,1 𝑉 )とする．

𝑔 0,1,0 = 2

𝑔 0,0,1 = 0

𝑔 1,1,1 = 2

𝑔 0,0,0 = 1

𝑔 1,1,0 = 1

𝑔 1,0,1 = 1集合関数𝑔: 0,1 𝑉 → ℝ を𝑔 ≔ 𝑓 ∘ 𝜎𝑆を 𝑆 による 𝑓 の対称差変換と呼ぶ．
𝑔 1,0,0 = 0

一般に，対称差変換は劣モジュラ性を保存しない．

定義: 劣モジュラ
集合関数 𝑓: 0,1 𝑉 → ℝ に対して，Φ𝑓 : 0,1 𝑉 × 0,1 𝑉 → ℝを
Φ𝑓 𝑋, 𝑌 ≔ 𝑓 𝑋 + 𝑓 𝑌 − 𝑓 𝑋 ∪ 𝑌 − 𝑓 𝑋 ∩ 𝑌とする．
• 集合関数𝑓が劣モジュラとは，任意の𝑋, 𝑌 ∈ 0,1 𝑉 に対して
Φ𝑓 𝑋, 𝑌 ≥ 0が成り立つこととする．

49

𝑔 ≔ 𝑓 ∘ 𝜎(1,0,1) の場合

定義: 対称差変換

𝑔 0,1,0 = 1対称差写像𝜎𝑠 : 0,1 𝑉 → 0,1 𝑉 を
𝜎𝑠 𝑋 ≔ 𝑋 ⊕ 𝑆 (𝑋 ∈ 0,1 𝑉 )とする．

𝑔 0,1,1 = 2

𝑔 0,0,0 = 0

𝑔 1,1,0 = 2

𝑔 0,0,1 = 1

𝑔 1,1,1 = 1

𝑔 1,0,0 = 1集合関数𝑔: 0,1 𝑉 → ℝ を𝑔 ≔ 𝑓 ∘ 𝜎𝑆を 𝑆 による 𝑓 の対称差変換と呼ぶ．
𝑔 1,0,1 = 0

劣モジュラ性を保存する非自明な対称差変換も存在．

50

定義: 劣モジュラ
集合関数 𝑓: 0,1 𝑉 → ℝ に対して，Φ𝑓 : 0,1 𝑉 × 0,1 𝑉 → ℝを
Φ𝑓 𝑋, 𝑌 ≔ 𝑓 𝑋 + 𝑓 𝑌 − 𝑓 𝑋 ∪ 𝑌 − 𝑓 𝑋 ∩ 𝑌とする．
• 集合関数𝑓が劣モジュラとは，任意の𝑋, 𝑌 ∈ 0,1 𝑉 に対して
Φ𝑓 𝑋, 𝑌 ≥ 0が成り立つこととする．

𝑔 ≔ 𝑓 ∘ 𝜎(1,0,1) の場合[Nakashima et al.]𝑔 0,1,0 = 1
定義: 対称差変換
劣モジュラ性を保存する対称差変換の特徴づけ．
対称差写像𝜎𝑠 : 0,1 𝑉 → 0,1 𝑉 を
狭義劣モジュラ関数ならO 𝑑 時間で発見．
𝑉
𝜎𝑠 𝑋 ≔ 𝑋 ⊕ 𝑆 (𝑋 [12]∈ 0,1
) Nakashima, Yukiko Yamauchi,
Junpei
Kijima
𝑔 S.
0,0,0
=and
0 Masafumi
𝑔 1,1,0 = 2とする．

𝑔 0,1,1 = 2
Yamashita, Finding submodularity hidden in symmetric difference,
SIAM Journal on Discrete Mathematics, 34:1 (2020), 571--585.

集合関数𝑔: 0,1 𝑉 → ℝ を𝑔 ≔ 𝑓 ∘ 𝜎𝑆

𝑔 0,0,1 = 1

𝑔 1,1,1 = 1

を 𝑆 による 𝑓 の対称差変換と呼ぶ．
𝑔 1,0,1 = 0

劣モジュラ性を保存する非自明な対称差変換も存在．

𝑔 1,0,0 = 1

51

任意の有限分配束は安定マッチングに同型

答え: 対数M#凸分布でもサンプリングは依然難しい.

II. 私の数理工学: 実践編
使ってこその数理工学

応用(お手伝い)研究でも、貢献してきました。
 学習理論 [ISAAC2011, ALT2012, 他][14] Daiki Suehiro, Kohei Hatano, S. Kijima, Eiji Takimoto and Kiyohito Nagano,
Online prediction under submodular constraints, Proc. ALT 2012, 260--274.
[15] Shota Yasutake, Kohei Hatano, S. Kijima, Eiji Takimoto and Masayuki Takeda,
Online linear optimization over permutations, Proc. ISAAC 2011, 534--543.

 分散計算論 [JACM 2017, TOCT 2021, 他][16] Akihiro Monde, Yukiko Yamauchi, Shuji Kijima and Masafumi Yamashita, Can askywalker localize the midpoint of a rope?, ACM Transactions on Computation
Theory, 13:3 (2021), Article 17, 1--23.
[17] Yukiko Yamauchi, Taichi Uehara, Shuji Kijima and Masafumi Yamashita, Plane
formation by synchronous mobile robots in the three dimensional Euclidean space,
Journal of the ACM, 64:3 (2017), Article 16.

53いつも論文化がうまくいくわけではない。
 単調CFTP ⇔ 対数優モジュラ分布
 [5] S. Kijima, Sampling from log-super/submodular distributions, in Proc. 7th HungarianJapanese Symposium on Discrete Mathematics and Its Applications (2011), 227--236.

 最大ポピュラーマッチングの構成員は一意 Mizuki Hirakawa, Yukiko Yamauchi, *S. Kijima, Masafumi Yamashita, On the structure of
popular matchings in the stable marriage problem – who can join a popular matching?, MATCHUP 2015, Glasgow, UK, April 16-18, 2015 (April 17).
 平川瑞樹, The Family of Vertex Sets of Popular Matchings in The Stable Marriage Problem,
九大システム情報 修士論文, 2015.

 大規模系統樹の作成 (マルウェア分類)
 *米巧，山内由紀子, 来嶋秀治, 山下雅史, 大規模マルウェアデータ群に対する系統樹推定, 日本オペレーションズ・リサーチ学会 2016年春季研究発表会, 慶應大学矢上キャン
パス（東京都新宿区）, 2016年3月17日-18日(3月18日).
 米巧, 大規模マルウェアデータ群に対する系統樹推定, 九大システム情報 修士論文,
2016.

 ストリーミングデータのO log log 𝑛 領域サンプリング [18] Masatora Ogata, Yukiko Yamauchi, Shuji Kijima and Masafumi Yamashita,A randomized
algorithm for finding frequent elements in streams using O(loglog N) space,in Proc. ISAAC
2011, 514--523.

54産官学連携の一プロジェクトの一チームの一研究課題

大規模マルウェアデータ群に対する系統樹推定
1九州大学 システム情報科学府 / 2九州先端科学技術研究所○米 巧1,2, 山内 由紀子1, 来嶋 秀治1,2, 山下 雅史1,2

謝辞
本研究は, 総務省「国際連携によるサイバー攻撃の予知技術の研究開発」の支援を受けた.
55

数理工学: 実践編
Cilibrasi & Vitanyiでマルウェア分類とかどうでしょう？

56

正規圧縮距離を用いたデータ分類
R. Cilibrasi, P. Vitanyi, Clusterring by compression, IEEE
Trans Information Theory, 51(2005), 1523--1545.
CilibrasiとVitanyiによるデータ分類概要
データ群
正規圧縮距離データ間の距離を定義
4点法(quartet method)

距離の木距離近似(系統樹作成)57

正規圧縮距離を用いたデータ分類
R. Cilibrasi, P. Vitanyi, Clusterring by compression, IEEE
Trans Information Theory, 51(2005), 1523--1545.
CilibrasiとVitanyiによるデータ分類概要
データ群
正規圧縮距離データ間の距離を定義
4点法(quartet method)

距離の木距離近似(系統樹作成)58

正規圧縮距離R. Cilibrasi, P. Vitanyi, Clusterring by compression, IEEE
Trans Information Theory, 51(2005), 1523--1545.

データ間の距離(類似度)を圧縮度に基づいて定義
Cf. コルモゴロフ複雑度
M.Li, X. Chen, X.Li, B. Ma, and P. Vitanyi, The similarity metric,
IEEE Transactions on Information Therory, 50, 3250--3264, 2004

max {𝑐 𝑥𝑦 − 𝑐 𝑥 , 𝑐 𝑦𝑥 − 𝑐(𝑦)}NCD(𝑥, 𝑦)＝max {𝑐 𝑥 , 𝑐(𝑦)}𝑐 𝑥 : 文字列xを圧縮プログラムCを用いて圧縮した文字列長
𝑥𝑦 : 文字列xと文字列yを連結した文字列文字列で表現されるデータであれば適用可能

59

正規圧縮距離を用いたデータ分類
R. Cilibrasi, P. Vitanyi, Clusterring by compression, IEEE
Trans Information Theory, 51(2005), 1523--1545.
CilibrasiとVitanyiによるデータ分類概要
データ群
正規圧縮距離データ間の距離を定義
4点法(quartet method)

距離の木距離近似(系統樹作成)60

系統樹推定R. Cilibrasi, P. Vitanyi, Clusterring by compression, IEEE
Trans Information Theory, 51(2005), 1523--1545.

距離行列を木距離で近似する．

0.2∼0.12345
0.166112

0.1∼
0.2∼

0.3∼

0.1∼

61

系統樹推定R. Cilibrasi, P. Vitanyi, Clusterring by compression, IEEE
Trans Information Theory, 51(2005), 1523--1545.

入力 : データ間の距離行列D
出力 : 各データに対応する葉をもつ系統樹T0.2∼
0.12345
0.166112

0.1∼
0.2∼

0.3∼

0.1∼

62

系統樹推定R. Cilibrasi, P. Vitanyi, Clusterring by compression, IEEE
Trans Information Theory, 51(2005), 1523--1545.

入力 : データ間の距離行列D
出力 : 各データに対応する葉をもつ系統樹T0.2∼
0.123450.1∼

0.166112

0.2∼

0.3∼

0.1∼

63

マルウェア分類に使ってみた．

数理工学: 実践編
Cilibrasi & Vitanyiでマルウェア分類とかどうでしょう？
こんな感じです。とりあえず50件でやってみました。

データ数𝑛 = 5 × 104 ．

あのね、データは5万件あるんだよ。

 距離行列作成の時間/領域計算量 O 𝑛2 = 2.5 × 109※ 便宜上，ここでは文字列圧縮にかかる時間は1とする．

 系統樹作成の時間計算量
 4点法 [Cilibrasi &Vitányi 05]: ヒューリスティクス (NP困難な最小化問題)
 近隣結合法[Saitou & Nei 87]: O 𝑛3 = 1.25 × 1014
[Yone et al. ]距離行列計算+系統樹作成: O 𝑠 3 + 𝑠𝑛 = 109 + 5 × 107
𝑠 = 103 はサンプルサイズ(パラメータ)

64高速化のための提案手法

概要
1. 一部のデータに対し, 系統樹を作成.
n’ : ベースサイズ (n’≪ n)𝑛𝑛`

3. 距離行列再計算．

…, 𝑛1, 2, 3, …

2. 残りのデータをベース系統樹上で分割．
1
2
3

4. 頂点の縮約．
5. 系統樹作成.

𝑛65

提案手法のイメージの部分で系統樹を作成，残りのデータをその系統樹上で分割し，ブロックごとに系統樹を作成する．

𝑎2

𝑎3
𝑎1

𝑎1

𝑎2
𝑎4
𝑎3
𝑎4
従来手法では計算していたの部分を計算せずに系統樹を作成．

66

67

私の思う数理工学: 実践編
一論文一主張
実践編を論文にするのは難しい．
• 「論文」は「ポイント(=原理)」が明確な方が評価されやすい．
 マルウェア分類はまだ論文に近づいて行った方．

• 一言で解決した話も多い (暗黙知)．
 本当は，こちらが理想なのでは？
• 実践は実際に使われることで貢献する方が健全なのでは?
 成果発表: 特許，ハード/ソフト開発，起業 etc.
 評価指標: DL数，売上高

実践編で私が大事に思うこと．
• サイズ感覚．
 データサイズ，要望応答時間，精度(粒度)．
 使うべき技術が変わる．

• 易しいアルゴリズム．

数理工学の出番 根拠は深く．アルゴリズム解析自体は難しくてよい．

III. 私の思う数理工学: 開拓編
拓く数理工学私の思う数理工学: 開拓編〇〇をどうにかしたい。
だけど、何をすればよいのかわからない。

69

私の数理工学: 理論編 (研究宣伝) 完璧サンプリング法
 分割表 [RSA 2006]
 待ち行列ネットワーク [SICOMP 2008, AOR 2008] メディアン安定マッチング [MathOR 2012]
 確率と計算[20]
S. Kijima,
Nobutaka
Shiraga,
How many vertices does a
決定性RW[RSA Shimizu,
2015, TCSTakeharu
2017, SIDMA
2018]
random walk miss in a network with moderately increasing the number of vertices?,
SODA
高次元体積の決定性近似[Algorithmica 2016, TCS 2020]
Proc.
2021, 106--122,

 対数優モジュラ分布

 劣モジュラ関数の変数変換 [SIDMA 2020] 任意の有限分配束は𝑀♮凹関数の最小化元族に同型 [OR lett. 2021]

 動的環境下のアルゴリズム論
 増えるクーポン収集 [SODA 2021]70

1. 学術的背景 (未発達のシーズ)動的ネットワーク上のアルゴリズム解析
起点: 「現実のネットワークは刻々と変化する」
 ネットワーク工学で特に注目 [Kuhn & Oshman SIGACT News 2011]
 [Augustine et al. SIGACT News 2016]
 [Michail & Spirakis CACM 2018]「しかし、理論的なことはほとんどわかっていない。」
 海外の大型プロジェクト
 “NETWORKS” (蘭 2014-2023)

 “DYNAMIC MARCH” (ERC Starting Grant 2016-2021)
今日は動的グラフ上のRWの話を中心に72

ランダムウォークとアルゴリズム
 RWとアルゴリズム
 ネットワーク探索
cf. 深さ優先探索，幅優先探索，
universal traversal sequence Cover time
 全探索にかかる(期待)時間
 [Aleliunas et al. 1979]
 単純RW 𝑡cov = O 𝑚𝑛 Ω 𝑛3 : lollipop
 [Ikeda et al. 2003, Nonaka et al. 2010]
 “𝜷-RW”, Metropolis walk 𝑡cov = O 𝑛2 log 𝑛 Ω(𝑛2 log 𝑛): spider [David & Feige 2017] local-min-degree RW 𝑡cov = O 𝑛2

⇒ グラフが時間変化するとき、全探索にかかる時間は？

73動的グラフ上のランダムウォーク（2010—2015頃の来嶋）
木場孝輔，山内由紀子, 来嶋秀治, 山下雅史, 動的
グラフ上のランダムウォークの到達時間と全訪問時間, 数理解析研講究録, 1849, 100—104, 2013.

 解析が難しい。
想像はつくけど、面倒そう。
 [L. Saloff-Coste and J. Zuniga 2009, 2011]

 頂点集合が変化する場合はどうなんだろう？
 [C. Cooper and A. Frieze 2003]

𝑐 stepごとに1頂点増える =>𝑐𝑛 をcover𝑐+1

 動的グラフ上のRWの黎明期

 [Avin, et al. 2008], [Clementi et al. 2008], [Baumann et al. 2009] etc.
 Workshop "Random Walks on Random Graphs and Applications,"
Eindhoven, The Netherlands, April 14-16, 2015.
 （当時のマイブームは決定性ランダムウォーク）74

動的ネットワーク上のアルゴリズム解析 --- なぜ難しい？

 解析には高い技術が必要．
 静的ネットワーク上のアルゴリズム解析より“難しい”．

 「何を」解析すれば良い？: 数理モデリング

モデルができる前には「何もない」



モデルができてみれば「当たり前」

数理工学

解析可能な数理モデル75

増えるクーポン収集問題*来嶋秀治 (九大), 清水伸高(東工大)，白髪丈晴 (東工大)*本研究はJSPS科研費 JP17K19982の助成を受けたものです。

クーポン収集• ビックリマンシール• プロ野球チップス
• ガチャ

77

• P●KEM●Nカード100種類．コンプまでに何枚引けば良い？

• 𝑛種類のカードが裏返しにおいてある．1枚選んでめくり，種類を確認したら元に戻して，よく混ぜる．全ての種類を確認するまでに何回めくればよいか？
• 年賀状のお年玉くじ下2桁．00から99までそろえるには年賀状を何枚もらえばよいか？
• 一万円札を千円札10枚に両替する．千円札の記番号下2桁を00から99までそろえるには，
いくら両替すればよいか？

1

2

3

4

5

6

7

8

78

クーポン収集 𝑋𝑖 ∈ {1, … , 𝑛} (𝑖 = 1,2, …)は独立とし，Pr[𝑋𝑖 = 𝑘] =1𝑛

𝑘 = 1, … , 𝑛 ．

 𝑇 = min 𝑡 𝑋1 , … , 𝑋𝑡 = 1, … , 𝑛 をコンプ回数(completion time)と呼ぶ．

Q. E 𝑇 を求めよ．

79

クーポン収集 𝑋𝑖 ∈ {1, … , 𝑛} (𝑖 = 1,2, …)は独立とし，Pr[𝑋𝑖 = 𝑘] =1𝑛

𝑘 = 1, … , 𝑛 ．

 𝑇 = min 𝑡 𝑋1 , … , 𝑋𝑡 = 1, … , 𝑛 をコンプ回数(completion time)と呼ぶ．

Q. E 𝑇 を求めよ．
𝑇𝑘 ≔ min 𝑡

𝑋1 , … , 𝑋𝑡

= 𝑘 :𝑘種コンプ回数

𝑆𝑘 ≔ 𝑇𝑘 − 𝑇𝑘−1 :𝑘 − 1種コンプしてから新種が出るまでの回数
Claim. E 𝑆𝑘 =

𝑛
．
𝑛−𝑘+1

 𝑘 − 1種コンプした状態で，新種の出る確率 𝑝𝑘 =
 従ってE 𝑆𝑘 =𝑛−(𝑘−1)．
𝑛1𝑛=(幾何分布の期待値)．
𝑝𝑘
𝑛−𝑘+1

𝑇 = 𝑇𝑛 = σ𝑛𝑘=1 𝑆𝑘 に注意して，≃519 (𝑛 = 100の時)

𝑛

𝑛

𝑛

𝑘=1

𝑘=1

𝑘 =1

𝑛
1
E 𝑇 = E 𝑇𝑛 = ෍ E 𝑆𝑘 = ෍= 𝑛 ෍ ≃ 𝑛 ln 𝑛 + 0.577
𝑛−𝑘+1
𝑘′
′

80

クーポン収集 𝑋𝑖 ∈ {1, … , 𝑛} (𝑖 = 1,2, …)は独立とし，Pr[𝑋𝑖 = 𝑘] =1𝑛

𝑘 = 1, … , 𝑛 ．

 𝑇 = min 𝑡 𝑋1 , … , 𝑋𝑡 = 1, … , 𝑛 をコンプ回数(completion time)と呼ぶ．

Q. E 𝑇 を求めよ．
Q. Pr 𝑇 ≥ 2𝑛 ln 𝑛 はどれくらいか？

81

クーポン収集 𝑋𝑖 ∈ {1, … , 𝑛} (𝑖 = 1,2, …)は独立とし，Pr[𝑋𝑖 = 𝑘] =1𝑛

𝑘 = 1, … , 𝑛 ．

 𝑇 = min 𝑡 𝑋1 , … , 𝑋𝑡 = 1, … , 𝑛 をコンプ回数(completion time)と呼ぶ．

Q. E 𝑇 を求めよ．
Q. Pr 𝑇 ≥ 2𝑛 ln 𝑛 はどれくらいか？
Chevyshevの不等式を使うとPr 𝑋 ≥ 2E 𝑋

𝑛2 𝜋 2
2
Var 𝑋
𝜋
6
≤
≤
=
2
2
E𝑋
𝑛 ln 𝑛6 ln 𝑛 2≃0.078 (𝑛 = 100の時)

増えるクーポン収集
毎日1回ガチャを引く．
• 毎日1種ずつ増えるとき，コンプ回数は？

82

83

増えるクーポン収集
Day (𝒊)

1日目

2日目

3日目

4日目

5日目

6日目

7日目

8日目

9日目

種類

1

2

4

𝟏
𝟏

𝟏
𝟐

5
𝟏
𝟓

6
𝟏
𝟔

7

𝐏𝐫 𝑿𝒊 = 𝒌

3
𝟏
𝟑

8
𝟏
𝟖

9
𝟏
𝟗

1

1

2

𝟏
𝟒

1

2

3

𝟏
𝟕

1

2

3

4

5

6

7

8

増えるクーポン収集
毎日1回ガチャを引く．
• 毎日1種ずつ増えるとき，コンプ回数は？
 コンプ不可能84

増えるクーポン収集
毎日1回ガチャを引く．
• 毎日1種ずつ増えるとき，コンプ回数は？
 コンプ不可能• 2日に1回新種リリースされたら，コンプ回数は？

85

86

増えるクーポン収集
Day

1日目

2日目

3日目

4日目

5日目

6日目

7日目

8日目

9日目

種類

1

1

2

2

4

𝟏
𝟏

𝟏
𝟏

𝟏
𝟐

𝟏
𝟐

3
𝟏
𝟑

4

𝐏𝐫 𝐗 𝐢 = 𝐤

3
𝟏
𝟑

𝟏
𝟒

𝟏
𝟒

5
𝟏
𝟓

2

1

2

3

1

1

1

2

1

増えるクーポン収集
毎日1回ガチャを引く．
• 毎日1種ずつ増えるとき，コンプ回数は？
 コンプ不可能• 2日に1回新種リリースされたら，コンプ回数は？
 コンプ不可能87

増えるクーポン収集
毎日1回ガチャを引く．
• 毎日1種ずつ増えるとき，コンプ回数は？
 コンプ不可能• 2日に1回新種リリースされたら，コンプ回数は？
 コンプ不可能• 10日に1回リリースされたら？
 コンプ不可能(種類数1000の時を考えてみよ)88ゆっくりと増えるクーポン収集
毎日1回ガチャを引く．
• 𝑛種類目がリリースされて𝑛日後に新種がリリース．
 Q. コンプ回数は？

89

90

ゆっくりと増えるクーポン収集
Day

1日目

2日目

3日目

4日目

5日目

6日目

7日目

8日目

9日目

種類

1

2

2

4

𝟏
𝟐

𝟏
𝟐

3
𝟏
𝟑

4

𝟏
𝟏

3
𝟏
𝟑

4

𝐏𝐫 𝑿𝒊 = 𝒌

3
𝟏
𝟑

𝟏
𝟒

𝟏
𝟒

𝟏
𝟒

1期

2期

1

1

3期

2

4期

1

2
3

1

2

3

4

ゆっくりと増えるクーポン収集
毎日1回ガチャを引く．
• 𝑛種類目がリリースされて𝑛日後に新種がリリース．
 Q. コンプ回数は？
 A. コンプ不可能(新種は出続けるから)

91

92

ゆっくりと増えるクーポン収集
毎日1回ガチャを引く．
• 𝑛種類目がリリースされて𝑛日後に新種がリリース．
 Q. コンプ回数は？
 A. コンプ不可能(新種は出続けるから)
 Q. 𝒏期末にはどのくらいコンプされる？
1. o 𝑛
2. 𝑐𝑛
3. 𝑛 − 𝑐 𝑛
4. 𝑛 − 𝑐 log 𝑛
5. 𝑛 − 𝑐𝑐は適当な定数93

ゆっくりと増えるクーポン収集
毎日1回ガチャを引く．
• 𝑛種類目がリリースされて𝑛日後に新種がリリース．
 Q. コンプ回数は？
 A. コンプ不可能(新種は出続けるから)
 Q. 𝒏期末にはどのくらいコンプされる？
1. o 𝑛
2. 𝑐𝑛
3. 𝑛 − 𝑐 𝑛
4. 𝑛 − 𝑐 log 𝑛

5. 𝒏 − 𝒄 (正解)

𝑐は適当な定数毎日1回ガチャを引く(一様ランダム)94
ゆっくりと増えるクーポン収集 𝔡 𝑛 : 𝑛期の日数
𝑈𝑛 : 𝑛期末の未収集アイテムの種類数

命題
𝔡 𝑛 = 𝑛のとき，E 𝑈𝑛 < 1．

毎日1回ガチャを引く(一様ランダム)95
ゆっくりと増えるクーポン収集 𝔡 𝑛 : 𝑛期の日数
𝑈𝑛 : 𝑛期末の未収集アイテムの種類数

命題
𝔡 𝑛 = 𝑛のとき，E 𝑈𝑛 < 1．
証明

1 (𝑛期末にアイテム 𝑖 が未収集)
 ℰ𝑖,𝑛 ≔ ൝
0 (𝑛期末にアイテム 𝑖 が収集済み)(𝑖 = 1,2, … , 𝑛)

 𝑈𝑛 = σ𝑛𝑖=1 ℰ𝑖,𝑛

 アイテム𝑛が𝑛期末に未収集の確率𝑛1
Pr[ℰ𝑛,𝑛 = 1] = 1 −< e−1𝑛 アイテム 𝑖 (𝑖 ≤ 𝑛)が 𝑛期末に未収集の確率
𝑖
𝑖+1
𝑛
𝑛+1−𝑖
1
1
1
1
Pr ℰ𝑖,𝑛 = 1 = 1 −
1−
… 1−
<
𝑖
𝑖+1
𝑛
e
1



1 𝑛+1−𝑖1111𝑛𝑛
E 𝑈𝑛 = σ𝑖=1 Pr ℰ𝑖,𝑛 < σ𝑖=1
= + 2 + ⋯ + 𝑛 < e1 =< 0.582
e
e
e
e
e−1
1−
e

クーポン収集からランダムウォークへ

・・・
クーポンのcompletion time = 完全グラフ上のRWのcover time96

頂点の増えるグラフ上のランダムウォーク [SODA 2021]
RWoGG: 𝔡,

𝐺𝑖

∞
𝑖=1

,

𝑃𝑖

∞
𝑖=1

定理(一般上界)𝔡 𝑖 ≥ 𝑐𝑡hit (𝑖) (𝑐 ≥ 1)のとき，E 𝑈 = 𝑂 1 ．
𝑛→∞𝔡 𝑖 𝑖→∞
さらに∞なら，E 𝑈𝑛0．
𝑡ℎ𝑖𝑡 (𝑖)定理(lazy reversible)
𝑃 𝑖 はlazy reversibleとする．
𝑡hit 𝑖𝑖𝛾≥ かつ𝔡 𝑖𝑡mix (𝑖)𝑐3𝑐𝑡hit 𝑖≥𝑖𝛾

8𝑛𝛾(𝑐 > 0)のとき，E 𝑈𝑛 ≤+ 32．
𝑐S. Kijima, N. Shimizu, T. Shiraga, How many vertices does a random walk miss in a
network with moderately increasing the number of vertices?, in Proc. SODA 2021,
106—122.

97増えるクーポン収集問題まとめ

参入お待ち
してます．

98

 増えるクーポン収集問題
 頂点がゆっくり増加するグラフ上のRWの非被覆具合
 期間𝔡と未訪問数E[𝑈𝑛 ]の関係
今後の課題

 𝐄 𝑼𝒏 = 𝐎 𝟏 の必要条件: 𝔡 𝑖 =? Ω 𝑡hit 𝑖 “Cover time”以外の指標
 ランダムウォークからアルゴリズムへ「動的環境下のアルゴリズム解析」
 とは？ “入力”が静的でない???
計算機科学では「問題」は「入力」と「質問」からなる:
例: 素数判定問題
入力: 𝑧 ∈ ℕ
質問: 𝑧は素数か?

まとめ

100

私の数理工学
私の数理工学: 理論編
 良い研究をして、良い論文を書く
II. 私の数理工学: 実践編
 使ってこその数理工学
III. 私の数理工学: 開拓編
 数理で工学を拓く
I.

これからも精進いたします．ご鞭撻のほど，お願いいたします．
理論，実践，開拓，いずれでも，ご興味のある方，お声がけください．ぜひ，一緒に研究しましょう．

Fin.