モジュラリティの上界値算出 - 上智大学 情報理工学部 情報理工学科

2011 年度 卒業論文
モジュラリティの上界値算出
指導教員 宮本 裕一郎 准教授
上智大学 理工学部 情報理工学科
A0878365 宮内 敦史
概要
コミュニティ検出結果の評価関数として，Newman らにより提案されたモジュラリティが知られている．
モジュラリティ最大化問題は NP-困難であるため，これまでに多くのヒューリスティクスが提案されてきた．
しかし一方で，厳密解法や精度保証付き近似解法はほとんど提案されていない．さらに，モジュラリティの上
界値算出に特化した手法は全く提案されていない．そのため，提案されたヒューリスティクスの性能評価は十
分に行えていないと言える．本論文では，モジュラリティの上界値算出に特化した手法を提案する．そして，
従来の手法では扱いきれなかった規模のネットワークに対して，その上界値を求める．
目次
1
序論
1
2
従来研究
2
2.1
自然な定式化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.2
疎な定式化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2.3
列生成法によるアプローチ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
制約式制限法
6
3.1
アイデアと単純な手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7
3.2
制約式制限法の工夫
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8
3.3
計算実験 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
3
変数制限法
12
4.1
アイデアと単純な手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
13
4.2
変数制限法の工夫 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
4.3
計算実験 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
結論と今後の課題
17
4
5
謝辞
17
参考文献
18
付録 A 本文中の定理の証明
19
付録 B 列生成法における補助問題の導出と上界値の計算手順
24
1 序論
近年，ソーシャル・ネットワーキング・サービスの普及 [4] などにより，ネットワーク解析が注
目されている．その中でも特に，ネットワークを密な部分ネットワークに分割する手法は，盛んに
研究されている [8]．このような密な部分ネットワークをコミュニティと呼び，ネットワークをコ
ミュニティに分割することをコミュニティ検出という．図 1 として，その例を挙げる．
図 1.
コミュニティ検出結果の例．左のネットワークのコミュニティ検出結果が，右のネット
ワークである．同じ色の頂点が同一のコミュニティに属している．
本論文では，無向グラフ G = (V, E) におけるコミュニティ検出を考える．グラフ G の頂点数
を n := |V |，枝数を m := |E| とする．2004 年に，コミュニティ検出結果の“良さ”を表す，モ
ジュラリティという指標が提案された [11]．正確には，モジュラリティは，頂点集合 V の分割
C = {V1 , V2 , . . . , Vk } を引数とする関数である．分割 C の要素であるコミュニティは非空かつ互い
∪k
に素であり， i=1 Vi = V を満たす．ここで，隣接行列の (i, j) 成分を Aij ，頂点 i の次数を di と
おく．また，頂点 i が属すコミュニティを Ci とし，δ をクロネッカーの記号とする．すると，モ
ジュラリティは
(
)
1 ∑∑
di dj
Q(C) =
Aij −
δ(Ci , Cj )
2m
2m
(1)
i∈V j∈V
と表される．モジュラリティは必ず 1 以下であり，値が高いほどコミュニティ検出の結果が良い
とされる．そして多くの場合，コミュニティ検出問題はモジュラリティ最大化問題として定義され
る．すなわち，入力として無向グラフ G = (V, E) が与えられ，最大のモジュラリティを実現する
1
頂点集合 V の分割 C を返す問題である．ここで，入力としてコミュニティの個数は与えられてい
ないことに注意する．本論文でも，コミュニティ検出問題として，モジュラリティ最大化問題を採
用する．
2008 年に，モジュラリティ最大化問題は NP-困難であることが示された [5]．それ以前から，最
適解を得るのは困難な問題であると予想されていたため，多くのヒューリスティクスが提案されて
きた [8]．最も有名なものとして，枝の媒介性を利用した GN アルゴリズム [11] が挙げられる．こ
のヒューリスティクスが，コミュニティ検出分野の発展を促したと言える．また，性能面で優れた
ものとして，ネットワークの階層性を利用したアルゴリズム [3] が挙げられる．この方法により，1
億頂点を超える規模のグラフに対しても，高速かつ良好なコミュニティ検出が可能となった．
ヒューリスティクスの性能評価を行う際，次の二つの尺度が考えられる．すなわち，検出に要す
る時間と，実現されたモジュラリティと最適値との差である．しかし，後者の尺度では，現在は詳
しい評価を行えていないと言える．なぜならば，厳密解法や精度保証付き近似解法がほとんど提案
されておらず [2, 6]，計算実験に使用されている大規模なネットワークでは，どの程度のモジュラ
リティが実現できるのかわかっていないからである．
本論文の目的は，提案されたヒューリスティクスの厳密な評価を可能にすることである．これを
達成するため，モジュラリティの上界値算出に特化した手法を提案する．提案手法の特徴は，高速
かつ省メモリであり，大規模なネットワークにも適用可能なことである．
ここで，本論文の構成を説明する．第 2 節では，モジュラリティの上界値算出に関連する従来研
究について触れ，それらの問題点を指摘する．第 3 節では，制約式制限法を提案し，計算実験によ
りその性能を評価する．さらに，第 4 節では，制約式制限法を利用した変数制限法を提案し，計算
実験によりその性能を評価する．最後に，第 5 節で本論文の結論と今後の課題を述べる．
2 従来研究
前節で述べた通り，モジュラリティの上界値算出に特化した手法は，これまでに提案されていな
い．しかし，厳密解法やヒューリスティクスの設計を目標として，その過程でモジュラリティの上
界値を求めている研究は見られる [1, 2]．よって，本節では，モジュラリティの上界値算出に関連
する従来研究について議論する．第 3 節および第 4 節の提案手法に直結するのは，2.2 節の疎な定
式化のみである．しかしながら，これは 2.1 節の自然な定式化から得られたものであるため，そ
れについても説明する．また，列生成法によるアプローチは代表的な手法であるため，2.3 節で触
れる．
2
2.1 自然な定式化
本節では，コミュニティ検出問題を整数計画問題として定式化する．そして，それを線形緩和
し，モジュラリティの上界値を得る定式化を与える．
頂点集合を V := {1, 2, . . . , n} とし，qij := Aij −
di dj
2m
とする．そして，頂点 i と頂点 j が同
一のコミュニティに属すとき 1，そうでないとき 0 をとる変数を xij とおく．任意の頂点 i に対し
て xii = 1 であるから，xii は変数として採用せず，定数 C :=
1
2m
∑n
i=1 qii
を目的関数値に加え
る．さらに，任意の頂点対 i, j に対して，xij = xji と qij = qji が成り立つ．従って，変数として
xij (i < j) のみを採用すると，コミュニティ検出問題は
1 ∑∑
qij xij + C
m i=1 j>i
n
(IP) : max.
n
(2)
s. t. xij + xjk − xik ≤ 1
for all i < j < k,
(3)
xij − xjk + xik ≤ 1
for all i < j < k,
(4)
− xij + xjk + xik ≤ 1
for all i < j < k,
(5)
xij ∈ {0, 1}
for all i < j,
(6)
(n)
3
3 = Θ(n ) 本ある．定式化
として定式化できる．変数は Θ(n2 ) 個，0-1 制約以外の制約式は 3
(IP) において，式 (6) を xij ∈ [0, 1] に置き換えることで，線形緩和問題
1 ∑∑
qij xij + C
m i=1 j>i
n
(LP) : max.
n
(7)
s. t. xij + xjk − xik ≤ 1
for all i < j < k,
(8)
xij − xjk + xik ≤ 1
for all i < j < k,
(9)
− xij + xjk + xik ≤ 1
for all i < j < k,
(10)
xij ∈ [0, 1]
for all i < j,
(11)
を得る．この (LP) を解くことで，モジュラリティの上界値が得られる．特に，顕著なコミュニ
ティ構造をもつグラフに対しては，最適値に近い上界値を得られることがわかっている．ただし，
頂点数が 200 を超える規模のグラフでは，実行時間とメモリ使用量の両面で，解くのが困難である．
2.2 疎な定式化
最近，定式化 (IP) に比べて制約式の本数を大幅に減らしながらも，同一の最適解集合をもつ定式
化が提案された [6]．ここで，頂点 i の隣接点の集合を N (i) とし，N (i, j) := N (i) ∪ N (j) \ {i, j}
3
とする．すると，提案された定式化は
1 ∑∑
qij xij + C
m i=1 j>i
n
(IPsparse ) : max.
s. t.
n
(12)
xki + xkj − xij ≤ 1
for all k (∈ N (i, j)) < i < j,
(13)
xik + xkj − xij ≤ 1
for all i < k (∈ N (i, j)) < j,
(14)
xik + xjk − xij ≤ 1
for all i < j < k (∈ N (i, j)) ,
(15)
xij ∈ {0, 1}
for all i < j,
(16)
と書ける．*1
目的関数は定式化 (IP) のそれと全く変わっていない．変わったのは，制約式の本数である．定
式化 (IP) では，制約式は Θ(n3 ) 本あった．これに対して，定式化 (IPsparse ) では，制約式の本数
の上界値が
n ∑
n
∑
(di + dj ) = (n − 1)
i=1 j>i
n
∑
di = O(mn)
(17)
i=1
で与えられる．コミュニティ検出の対象となる多くのネットワークは疎であり，枝数 m が頂点数
n に比例する．従って，応用上のほとんどの場面では，制約式は O(n2 ) 本であると評価できる．そ
して，次の定理が成り立つ [6]．
定理 1. (IP) と (IPsparse ) は同一の最適解集合をもつ．
証明については付録 A を参照されたい．この定理により，(IPsparse ) を解くことで，コミュニティ
検出問題の厳密解（およびモジュラリティの最適値）を得られることがわかる．定式化 (IPsparse )
において，式 (16) を xij ∈ [0, 1] に置き換えることで，線形緩和問題
1 ∑∑
qij xij + C
m i=1 j>i
n
(LPsparse ) : max.
s. t.
n
(18)
xki + xkj − xij ≤ 1
for all k (∈ N (i, j)) < i < j,
(19)
xik + xkj − xij ≤ 1
for all i < k (∈ N (i, j)) < j,
(20)
xik + xjk − xij ≤ 1
for all i < j < k (∈ N (i, j)) ,
(21)
xij ∈ [0, 1]
for all i < j,
(22)
を得る．そして，次の定理が成り立つ [6]．
定理 2. (LP) と (LPsparse ) は同一の最適解集合をもつ．
先ほどと同じく，証明については付録 A を参照されたい．この定理により，(LPsparse ) で得ら
れるモジュラリティの上界値は，(LP) で得られる値と等しいことがわかる．そして，実験的には
*1
文献 [6] の定式化に同値な変形を施し，細部に若干の修正を加えている．
4
(LP) よりも (LPsparse ) の方が高速に解けることがわかっている．しかし，1000 頂点を超える規
模のグラフに対しては，コミュニティ構造が自明である場合を除いて，非常に時間がかかる．さら
に，変数が Θ(n2 ) 個あるため，数千頂点規模のグラフでは，メモリ使用量の面で問題がある．
2.3 列生成法によるアプローチ
Aloise らにより，コミュニティ検出問題に対する列生成法を利用した厳密解法が提案された [2]．
本節では，その手法の中で，モジュラリティの上界値算出に必要な部分のみを説明する．
本節を通して，頂点集合 V のべき集合から空集合を除いた集合族 P := 2V \ ∅ を分割集合とい
い，その要素を分割要素と呼ぶ．従って，分割要素は P := 2n − 1 個あり，P = {1, 2, . . . , P } と
ラベル付けされているとする．ここで，頂点 i が分割要素 p に入っているとき 1，そうでないとき
0 である定数を aip とおく．また，分割要素 p をコミュニティとして採用するとき 1，そうでない
とき 0 をとる変数を xp とおく．すると，コミュニティ検出問題は
1 ∑∑∑
aip ajp qij xp + C
m p=1 i=1 j>i
P
(IPCG ) : max.
s. t.
P
∑
n
n
aip xp = 1
(23)
for all 1 ≤ i ≤ n,
(24)
for all 1 ≤ p ≤ P,
(25)
p=1
xp ∈ {0, 1}
として定式化できる．この定式化は，先に見た定式化 (IP) と比べて，かなり不器用に見える．な
ぜならば，変数が 2n − 1 個と n の指数に比例しているからである．しかしながら，列生成法の適
用を考えると，自然な定式化であると言える．定式化 (IPCG ) において，式 (25) を線形緩和して，
線形計画問題
1 ∑∑∑
aip ajp qij xp + C
m p=1 i=1 j>i
P
(LPCG ) : max.
s. t.
P
∑
n
n
aip xp = 1
(26)
for all 1 ≤ i ≤ n,
(27)
for all 1 ≤ p ≤ P,
(28)
p=1
xp ≥ 0
を得る．ここで，式 (28) が xp ∈ [0, 1] ではないことに注意する．式 (27) により，任意の分割要素
p に対して xp ≤ 1 であることは保証されているから，どちらの条件も同値である．この (LPCG )
を解くことで，モジュラリティの上界値が得られる．そこで，列生成法により解くことを考える．
定式化 (LPCG ) の双対問題の暫定最適解を y ∗ = (y1∗ , y2∗ , . . . , yn∗ )T とおくと，生成する列を決定す
5
る補助問題は
(QP) : min.
n
∑
i=1
s. t.
1 ∑∑
qij ai aj
m i=1 j>i
n
yi∗ ai −
ai ∈ {0, 1}
n
for all 1 ≤ i ≤ n,
(29)
(30)
として定式化できる．定式化 (QP) では，y1∗ , y2∗ , . . . , yn∗ は定数であり，a1 , a2 , . . . , an が変数であ
ることに注意する．補助問題 (QP) の導出の詳細については，付録 B を参照されたい．さらに，そ
の後の上界値算出の手順も付録 B に記した．
ところで，(QP) は変数が n 個でかつ係数が完全に密の二次整数計画問題であり，頂点数が 100
を超える規模のネットワークでは，計算時間の面で解くのが困難である．そこで，文献 [2] では，
新たに変数 xe (e ∈ E) と変数 D を導入し，定式化 (QP) は
(QP∗ ) : min.
n
∑
yi∗ ai −
i=1
s. t.
D=
∑ xe ( D )2
+
m
2m
(31)
e∈E
n
∑
di ai ,
(32)
i=1
x e ≤ ai
for all e = {i, j} ∈ E,
(33)
x e ≤ aj
for all e = {i, j} ∈ E,
(34)
xe ∈ {0, 1}
for all e = {i, j} ∈ E,
(35)
と書き換えられている．定式化 (QP∗ ) では，変数が m + n + 1 個，0-1 制約以外の制約式が 2m + 1
本ある．定式化 (QP∗ ) の最適解は (QP) の最適解を含んでおり，かつ等しい最適値が得られる．
よって，補助問題として (QP∗ ) を採用しても良い．そして，実験的には (QP) よりも (QP∗ ) の方
が高速に解けることがわかっている．しかしながら，頂点数が 500 を超える規模のネットワークで
は，計算時間の面で解くのが困難である．なお，この書き換えには，文献 [12] で議論されているコ
ミュニティ検出問題の定式化法が利用されている．
列生成法の問題点として，列生成を繰り返しても，最適値が同じ値であり続けることが度々起こ
ることが挙げられる．特に，コミュニティ検出のように最適解が 0 を多く含む問題においては，最
適値への収束が遅いことが知られている．このような問題を解決するため，Aloise ら [2] は文献 [7]
の手法を適用しているが，ここでは触れない．
3 制約式制限法
本節では，モジュラリティの上界値を求める制約式制限法を提案する．ここで提案する制約式制
限法には，I から III までの三つの手法がある．いずれの手法も，2.2 節の定式化 (LPsparse ) と，文
献 [9] で議論されている切除平面法のテクニックを基礎としている．はじめに，3.1 節で制約式制
6
限法のアイデアを示し，制約式制限法 I と制約式制限法 II を提案する．3.2 節では，それらを計算
速度と省メモリ性の両面で改善した手法である制約式制限法 III を提案し，それが制約式制限法 II
の拡張であることを見る．そして，3.3 節で計算実験によりそれらの性能を評価する．
3.1 アイデアと単純な手法
2.2 節の疎な定式化による制約式の減少は劇的であった．定式化 (LP) では Θ(n3 ) 本あった制約
式が，定式化 (LPsparse ) では疎なグラフにおいて O(n2 ) 本にまで減少した．しかしながら，頂点
数が 1000 程度のグラフにおいては，相当大規模な線形計画問題を解くことになってしまう．そこ
で，加える制約式の本数を制限し，(LPsparse ) の部分問題を解くことを考える．このとき，全ての
制約式を採用せずとも，目的関数も含めて同等の実行可能領域が得られるならば，(LPsparse ) を単
純に解いたときと等しい上界値が得られる．これこそが，制約式制限法のアイデアである．ここ
で，アイデアから自明な制約式制限法 I を提案する．
制約式制限法 I：
Step1. 定式化 (LPsparse ) を無制約で解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を得る．
Step2. 定式化 (LPsparse ) の全ての制約式を走査する．x̄∗ がその全てを満たしていれば，
モジュラリティの上界値として z̄ ∗ を返して終了する．そうでなければ，x̄∗ が満たさ
ない全ての制約式を加えて解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を更新する．そして，
Step2 を繰り返す．
定式化 (LPsparse ) では，1000 頂点程度のグラフであれば全ての制約式の走査が可能となったた
め，それを利用している．この手法で最終的に加える制約式の本数は，(LPsparse ) の制約式の本数
よりも少なくなると考えられる．そして，扱う線形計画問題が小規模に抑えられることで，単純
に (LPsparse ) を解くよりも高速に上界値が求められると期待できる．しかも，その最適性は保たれ
ており，得られる上界値は (LPsparse ) を単純に解いたときと一致する．次に，制約式制限法 I の
Step1 に変更を施した制約式制限法 II を提案する．
制約式制限法 II：
Step1. 定式化 (LPsparse ) について，k ∈ N (i, j) を k ∈ N (i) ∩ N (j) \ {i, j} とした制約
式のみを全て採用して解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を得る．
Step2. 定式化 (LPsparse ) の全ての制約式を走査する．x̄∗ がその全てを満たしていれば，
モジュラリティの上界値として z̄ ∗ を返して終了する．そうでなければ，x̄∗ が満たさ
ない全ての制約式を加えて解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を更新する．そして，
Step2 を繰り返す．
7
制約式制限法 II では，定式化 (LPsparse ) の制約式の中でも厳しいであろうものを，最初に全て
加えている．これによって，より速く最適解に近づき，実行時間が短くなると期待できる．
3.2 制約式制限法の工夫
制約式制限法 I および制約式制限法 II では，暫定最適解が満たしていない制約式を全て加えてい
た．しかし，そのような制約式の中には，暫定最適解が大きく満たしていないものから僅かに満た
していないものまでが混在している．ここでは，制約式制限法のアイデアを基本としながら，暫定
最適解が“最も満たしていない”制約式を中心に加えていくことを考える．
定式化 (LPsparse ) の制約式 xik + xkj − xij ≤ 1 (i < k ∈ N (i, j) < j) に注目する．以下の議論
は，制約式 xki + xkj − xij ≤ 1 (k ∈ N (i, j) < i < j) と xik + xjk − xij ≤ 1 (i < j < k ∈ N (i, j))
でも同様に行えるため，それらについての説明は省く．ここで，暫定最適解 x̄∗ が制約式 constr を
満たしていない度合いを δ(x̄∗ , constr) := x̄∗ik + x̄∗kj − x̄∗ij − 1 とし，不満足度と呼ぶことにする．
すると，任意の暫定最適解 x̄∗ と任意の制約式 constr に対して δ(x̄∗ , constr) ∈ [−2, 1] である．そ
して，δ(x̄∗ , constr) > 0 のとき，暫定最適解 x̄∗ は制約式 constr を満たしていないことがわかる．
さらに，δ(x̄∗ , constr) が大きい制約式 constr ほど，それを加えたときに暫定最適解 x̄∗ への影響
が大きいと考えられる．そこで，次の手法について考える．
制約式制限法 III (仮)：
Step1. 定式化 (LPsparse ) について，k ∈ N (i, j) を k ∈ N (i) ∩ N (j) \ {i, j} とした制約式
のみを全て採用して解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を得る．
Step2. 定式化 (LPsparse ) の全ての制約式について，不満足度 δ(x̄∗ , constr) を計算する．
全ての制約式について δ(x̄∗ , constr) ≤ 0 であれば，モジュラリティの上界値として z̄ ∗
を返して終了する．そうでなければ，最大の δ(x̄∗ , constr) を実現した制約式を 1 本加
えて解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を更新する．そして，Step2 を繰り返す．
上の Step2 において，最大の δ(x̄∗ , constr) を実現した制約式が複数本ある場合には，その中の
1 本をランダムに加えることとする．この手法は，Step1 を初期操作とみなすと，暫定最適解が最
も満たしていない制約式を 1 本ずつ加える単純な手法である．
ただ，この手法には大きな問題点が二つある．一つ目は，制約式を 1 本だけ加えるために，毎回
O(mn) 本の制約式を走査していることである．これにより，計算時間が膨大になることが予想さ
れる．二つ目は，最大の不満足度をもつ制約式が多数存在する場合が多く，1 本ずつ採用する意味
があまりないことである．
これらの問題点を解消するには，なるべく不満足度の大きい制約式から，複数本ずつ加えていく
必要がある．これを実現する手法を，制約式制限法 III として提案する．
8
制約式制限法 III：
Step1. 定式化 (LPsparse ) について，k ∈ N (i, j) を k ∈ N (i) ∩ N (j) \ {i, j} とした制約式
のみを全て採用して解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を得る．
Step2. 定式化 (LPsparse ) の全ての制約式について，不満足度 δ(x̄∗ , constr) を計算する．
全ての制約式について δ(x̄∗ , constr) ≤ 0 であれば，モジュラリティの上界値として z̄ ∗
を返して終了する．そうでなければ，閾値 t ∈ (0, 1] 以上の不満足度 δ(x̄∗ , constr) を実
現した全ての制約式を加えて解き，暫定最適解 x̄∗ と暫定最適値 z̄ ∗ を更新する．さら
に，閾値 t も更新する．そして，Step2 を繰り返す．
制約式制限法 III において，常に閾値 t = 0 とした手法が制約式制限法 II である．つまり，制約
式制限法 III は制約式制限法 II の自然な拡張である．上手く閾値を設定することで，高速性と省メ
モリ性の両面で優れた手法になると期待できる．
3.3 計算実験
制約式制限法の性能を評価するため，計算実験を行った．はじめに，計算実験の環境を表 1 に示
す．次に，計算実験に使用したネットワーク（問題例）の情報を表 2 に挙げる．
表 1.
計算実験の環境
OS
Mac OS 10.6.8
プロセッサ
2.4 GHz Intel Core 2 Duo
メモリ
4 GB 1067 MHz DDR3
ソルバー
Gurobi Optimizer 4.6.1
プログラミング言語
Python 2.6.1
従来手法と制約式制限法 I の計算結果を，表 3 に示す．定式化 (LP) を単純に解く手法では，
300 頂点のネットワークで，3600 秒以内に問題の入力ができないことを確認した．一方，定式化
(LPsparse ) を単純に解く手法では，1500 頂点のネットワークでも，3600 秒以内に問題の入力がで
きることを確認した．しかし，問題番号 8 以降の大規模なネットワークでは，問題の入力に多くの
時間がかかり，精度の良い上界値は得られなかった．
制約式制限法 I では，最終的に加えた制約式が (LPsparse ) の制約式の 10 % 程度にとどまり，扱
う線形計画問題が小規模に抑えられた．これにより，多くのネットワークにおいて実行時間が短縮
された．さらに，問題番号 8 以降の大規模なネットワークに対しては，計算を止めた 3600 秒の時
点で (LPsparse ) を単純に解く手法よりも精度の良い上界値が得られた．これは，問題の入力にかか
9
表 2.
計算実験に用いたネットワーク．問題番号 1 から 4 のネットワークは，計算実験に
頻繁に使用される実ネットワークである．問題番号 5 以降は，次数分布がべき乗則に従うス
ケールフリー・グラフである．これらは，Python のパッケージである Networkx 2.7.2 の関数
scale free graph(n) により作成した．また，下界値とは，モジュラリティの下界値である．問
題番号 1 から 4 のネットワークの下界値は，文献 [2] に記されている最適値である．問題番号 5
以降のネットワークの下界値は，文献 [3] のアルゴリズムを基本とした手法で求めた．
問題番号
ネットワーク名
n
m
下界値
1
Zachary’s karate club
34
78
0.420
2
Dolphin’s social network
62
159
0.529
3
Books about U.S. politics
105
441
0.527
4
American College Football
115
613
0.605
5
Scale-free 300
300
466
0.534
6
Scale-free 500
500
841
0.515
7
Scale-free 700
700
1194
0.528
8
Scale-free 1000
1000
1670
0.534
9
Scale-free 1200
1200
2137
0.518
10
Scale-free 1500
1500
2586
0.532
る時間が大幅に短くなり，より早く計算を始められたことによると考えられる．
次に，制約式制限法 II の計算結果を，表 4 に示す．時間内に計算が終了した全てのネットワー
クに対して，制約式制限法 II では制約式制限法 I よりも高速に上界値が求められた．これは，定式
化 (LPsparse ) の中でも厳しい制約式を最初に全て加えたことで，より速く最適解に近づけたこと
によると考えられる．また，多くのネットワークに対して，最終的に加えた制約式の本数が減少し
た．これは，厳しい制約式を最初に全て加えたことで，それと同等の実行可能領域を得るためによ
り多くの制約式を加えることを防いだためであると考えられる．問題番号 8 以降のネットワークに
対しては，制約式制限法 I よりも最終的に加えた制約式の本数が増えている．これは，計算を止め
た 3600 秒の時点で，より真の計算終了に近づいていたためであると考えられる．なぜならば，制
約式制限法 II で得られた上界値の精度が，制約式制限法 I で得られた上界値の精度を上回っている
からである．従って，高速性と省メモリ性の両面で，制約式制限法 II の方が優れた手法である．
最後に，制約式制限法 III の計算結果を，表 5 に示す．手法の内部で用いる閾値 t の設定方法と
して，t = 1 からスタートし，t ← t − 0.1 と更新を続けた．そして，t = 0 となるときに更新を止
め，それ以降は δ(x̄∗ , constr) > 0 を満たす制約式を全て加えるようにした．制約式制限法 III で
は，問題番号 6 までの小規模なネットワークにおいて，制約式制限法 II よりも実行時間が長くなっ
10
表 3.
計算結果（その 1）．2 段になっているセルの上段は，得られたモジュラリティの上界値
と計算時間である．つまり，その手法の精度と高速性を見ることができる．OM とは，3600 秒
以内に問題の入力ができなかったことを示している．さらに，問題の入力ができても 3600 秒以
内に計算が終了しなかった場合は OT と表示し，3600 秒の時点で得られた上界値を記した．ま
た，下段の括弧内は，左が最終的に加えた変数の個数（後に提案する手法を評価する際に必要）
，
右が最終的に加えた制約式の本数である．つまり，その手法の省メモリ性を示している．
問題番号
n
(LP)
(LPsparse )
制約式制限法 I
1
34
0.420 / 1s
(0.6k, 18k)
0.420 / 0s
(0.6k, 4k)
0.420 / 0s
(0.6k, 1k)
2
62
0.531 / 6s
(1.9k, 113k)
0.531 / 1s
(1.9k, 18k)
0.531 / 1s
(1.9k, 2k)
3
105
0.528 / 83s
(5.5k, 562k)
0.528 / 20s
(5.5k, 86k)
0.528 / 5s
(5.5k, 9k)
4
115
0.605 / 123s
(6.6k, 741k)
0.605 / 23s
(6.6k ,133k)
0.605 / 6s
(6.6k, 11k)
5
300
OM
(45k, 13M)
0.564 / 53s
(45k, 262k)
0.564 / 29s
(45k, 25k)
6
500
OM
(125k, 62M)
0.546 / 580s
(125k, 791k)
0.546 / 367s
(125k, 87k)
7
700
OM
(245k, 171M)
0.568 / 2476s
(245k, 1602k)
0.568 / 2969s
(245k, 176k)
8
1000
OM
(500k, 499M)
0.592 / OT
(500k, 3185k)
0.576 / OT
(500k, 284k)
9
1200
OM
(719k, 862M)
0.630 / OT
(719k, 4948k)
0.592 / OT
(719k, 295k)
10
1500
OM
(1124k, 1684M)
0.778 / OT
(1124k, 7446k)
0.599 / OT
(1124k, 422k)
た．しかし，問題番号 7 のネットワークでは，制約式制限法 II の実行時間を大きく下回った．ま
た，最終的に加えた制約式の本数は最も少なく，省メモリ性においても高い性能を示している．特
に，問題番号 10 のネットワークに対して，制約式制限法 II よりも加えた制約式の本数は 10 万本
以上少ないにも関わらず，等しい上界値が得られた．従って，これまでに見た手法の中では，高速
性と省メモリ性の両面で，制約式制限法 III が最も優れている．
11
表 4.
計算結果（その 2）．表の見方は，表 3 に同じである．
問題番号
n
(LPsparse )
制約式制限法 I
制約式制限法 II
1
34
0.420 / 0s
(0.6k, 4k)
0.420 / 0s
(0.6k, 1k)
0.420 / 0s
(0.6k, 1k)
2
62
0.531 / 1s
(1.9k, 18k)
0.531 / 1s
(1.9k, 2k)
0.531 / 1s
(1.9k, 2k)
3
105
0.528 / 20s
(5.5k, 86k)
0.528 / 5s
(5.5k, 9k)
0.528 / 3s
(5.5k, 9k)
4
115
0.605 / 23s
(6.6k ,133k)
0.605 / 6s
(6.6k, 11k)
0.605 / 4s
(6.6k, 11k)
5
300
0.564 / 53s
(45k, 262k)
0.564 / 29s
(45k, 25k)
0.564 / 27s
(45k, 26k)
6
500
0.546 / 580s
(125k, 791k)
0.546 / 367s
(125k, 87k)
0.546 / 271s
(125k, 86k)
7
700
0.568 / 2476s
(245k, 1602k)
0.568 / 2969s
(245k, 176k)
0.568 / 2488s
(245k, 169k)
8
1000
0.592 / OT
(500k, 3185k)
0.576 / OT
(500k, 284k)
0.575 / OT
(500k, 289k)
9
1200
0.630 / OT
(719k, 4948k)
0.592 / OT
(719k, 295k)
0.590 / OT
(719k, 296k)
10
1500
0.778 / OT
(1124k, 7446k)
0.599 / OT
(1124k, 422k)
0.593 / OT
(1124k, 483k)
4 変数制限法
前節の制約式制限法では，加える制約式を制限し，扱う線形計画問題を小規模に抑えた．表 3 から
表 5 の通り，制約式制限法の各種では，大規模なネットワークを中心に，最終的に加えた制約式の本
数よりも変数の個数の方が多くなった．つまり，省メモリ性の観点からは，変数の個数が問題になっ
ている．本節では，採用する変数をも制限し，高速性および省メモリ性をさらに向上させた変数制限
法を提案する．ここで提案する変数制限法には，I と II の二つの手法がある．はじめに，4.1 節で変
数制限法のアイデアを示し，単純な手法である変数制限法 I を提案する．4.2 節では，それを精度面
で改善した変数制限法 II を提案する．そして，4.3 節で計算実験によりそれらの性能を評価する．こ
こで，定式化 (LPsparse ) の制約式 xik +xkj ≤ xij +1 (i < k ∈ N (i, j) < j) に対する本節の全ての議
12
表 5.
計算結果（その 3）．表の見方は，表 3 に同じである．
問題番号
n
(LPsparse )
制約式制限法 II
制約式制限法 III
1
34
0.420 / 0s
(0.6k, 4k)
0.420 / 0s
(0.6k, 1k)
0.420 / 0s
(0.6k, 1k)
2
62
0.531 / 1s
(1.9k, 18k)
0.531 / 1s
(1.9k, 2k)
0.531 / 2s
(1.9k, 2k)
3
105
0.528 / 20s
(5.5k, 86k)
0.528 / 3s
(5.5k, 9k)
0.528 / 8s
(5.5k, 9k)
4
115
0.605 / 23s
(6.6k ,133k)
0.605 / 4s
(6.6k, 11k)
0.605 / 11s
(6.6k, 11k)
5
300
0.564 / 53s
(45k, 262k)
0.564 / 27s
(45k, 26k)
0.564 / 38s
(45k, 24k)
6
500
0.546 / 580s
(125k, 791k)
0.546 / 271s
(125k, 86k)
0.546 / 285s
(125k, 77k)
7
700
0.568 / 2476s
(245k, 1602k)
0.568 / 2488s
(245k, 169k)
0.568 / 1717s
(245k, 145k)
8
1000
0.592 / OT
(500k, 3185k)
0.575 / OT
(500k, 289k)
0.575 / OT
(500k, 215k)
9
1200
0.630 / OT
(719k, 4948k)
0.590 / OT
(719k, 296k)
0.582 / OT
(719k, 278k)
10
1500
0.778 / OT
(1124k, 7446k)
0.593 / OT
(1124k, 483k)
0.593 / OT
(1124k, 376k)
論は，制約式 xki +xkj ≤ xij +1 (k ∈ N (i, j) < i < j) と xik +xjk ≤ xij +1 (i < j < k ∈ N (i, j))
でも同様に成り立つことを断っておく．
4.1 アイデアと単純な手法
はじめに，孤立点はモジュラリティに影響を与えないため，存在しないと仮定する．任意の頂点
対 {i, j} ∈
/ E について Aij = 0 であるため，qij = Aij −
di dj
2m
< 0 が成り立つ．従って，定式化
(LPsparse ) において変数 xij ({i, j} ∈
/ E) は採用せず，これらの変数が関わる制約式を一切加えな
ければ，それを解くことでモジュラリティの上界値が得られる．これこそが，変数制限法のアイデ
アである．
しかし，xij ({i, j} ∈
/ E) を全く採用しないと，ほとんど自明な上界値しか得られない．そこで，
13
目的関数値を下げる方向に強く働く変数 xij ({i, j} ∈
/ E) を採用したい．つまり，最終的に 0 にな
らないか，あるいは制約式上で他の変数に影響を与える変数 xij ({i, j} ∈
/ E) を見つけたい．
こ こ で ，頂 点 i と 頂 点 j の グ ラ フ 上 で の 距 離 を d(i, j) と 表 す こ と に す る ．い ま ，変 数
xij (d(i, j) = 2) と，それが関わる制約式 xik + xkj ≤ xij + 1 (i < k ∈ N (i, j) < j) について考え
る．d(i, j) = 2 であるから，頂点 i と j は N (i, j) の要素に必ず共通の隣接点をもつ．つまり，
{i, k} ∈ E かつ {k, j} ∈ E をみたす頂点 k が存在する．いま，頂点対 {i, k} ∈ E と {k, j} ∈ E
について Aik = Akj = 1 であるから，多くの場合 qik > 0 かつ qkj > 0 が成り立つ．従って，制
約式が xij ∈ [0, 1] (i < j) のみであれば，たいていは xik = xkj = 1 となる．ところが，Aij = 0
により qij < 0 であるから，制約式の右辺にある変数 xij は 0 にしておきたい．従って，変数
xij (d(i, j) = 2) と制約式 xik + xkj ≤ xij + 1 (i < k ∈ N (i, j) < j) を加えることで，xij が 0 に
ならないか，あるいは xik = xkj = 1 とならないかのどちらかが起こると期待される．こうなれ
ば，目的関数値は減少する．よって，変数 xij (d(i, j) = 2) のほとんど全てが，上で見つけたいと
考えている変数である．ここまでの議論から，以下に変数制限法 I を提案する．
変数制限法 I：
Step1. 定式化 (LPsparse ) において，d(i, j) ≤ 2 を満たす変数 xij と，それらのみが関わ
る制約式を全て採用した部分問題 (sub LPsparse ) をつくる．
Step2. 部分問題 (sub LPsparse ) に対して制約式制限法 III を適用する．
変数制限法 I では，採用する変数の個数が大幅に減少し，それによって最終的に加える制約式の
本数も減少すると考えられる．さらに，扱う線形計画問題が小規模に抑えられることで，計算時間
の短縮も期待できる．しかし，それで得られる上界値は，(LPsparse ) を単純に解いたときに得られ
る上界値と等しい保証はなく，精度面での問題がある．
4.2 変数制限法の工夫
変数制限法では，変数 xij (d(i, j) ≤ 2) のみを採用することで，扱う線形計画問題を小規模に
抑えた．しかしながら，その代償として，得られる上界値の精度が低下すると予想された．これ
を改善するには，変数 xij (d(i, j) ≥ 3) を採用しなければならない．しかし，実ネットワークの
スモールワールド性 [10] により，変数 xij (d(i, j) = 3) を最初から全て採用すると，変数制限が
不十分になると考えられる．そこで，変数制限法 I を終えてから，その最適解を利用して，変数
xij (d(i, j) = 3) を部分的に採用することを考える．
いま，変数 xij (d(i, j) = 3) と，それが関わる制約式 xik + xkj ≤ xij + 1 (i < k ∈ N (i, j) < j)
について考える．すると，d(i, j) = 3 かつ k ∈ N (i, j) により，{i, k} ∈ E と {k, j} ∈ E のどち
らか一方のみが成り立つ．ここで，{i, k} ∈ E が成り立つと仮定する．以下の議論は，{k, j} ∈ E
14
であると仮定しても，同様に成り立つ．すると，Aik = 1 であるため，多くの場合 qik > 0 であ
る．従って，変数制限法 I で最終的に得られた解 x̄∗ において，たいていは x̄∗ik = 1 となっている．
よって，もし x̄∗kj > 0 であれば，制約式の左辺は 1 より大きいと期待される．ところが，Aij = 0
により qij < 0 であるから，制約式の右辺にある変数 xij は 0 にしておきたい．ここで，変数制限
法 I で最終的に解 x̄∗ が得られたとき，
(a) d(i, j) = 3 かつ {i, k} ∈ E かつ x̄∗kj > t (∈ [0, 1))，
(b) d(i, j) = 3 かつ {k, j} ∈ E かつ x̄∗ik > t (∈ [0, 1))
のどちらか一方を満たす変数 xij と制約式 xik + xkj ≤ xij + 1 (i < k ∈ N (i, j) < j) を加えると
する．すると，これまでの議論により，xij が 0 にならないか，x̄∗ik または x̄∗kj の値が変化するか
のどちらかが起こると期待できる．こうなれば，変数制限法 I で得られた上界値は減少し，すなわ
ち精度は向上する．ここまでの議論から，以下に変数制限法 II を提案する．
変数制限法 II：
Step1. 定式化 (LPsparse ) において，d(i, j) ≤ 2 を満たす変数 xij と，それらのみが関わ
る制約式を全て採用した部分問題 (sub LPsparse ) をつくる．
Step2. 部分問題 (sub LPsparse ) に対して制約式制限法 III を適用し，最適解 x̄∗ と最適値
z̄ ∗ を得る．以前に Step3 を実行したならば，上界値として z̄ ∗ を返して終了する．
Step3. 閾値 t ∈ [0, 1) を定める．上の (a) と (b) のどちらか一方を満たす全ての変数 xij
を (sub LPsparse ) に加える．そして，Step2 に戻る．
変数制限法 II は，変数制限法 I を終えた後に操作を行うため，実行時間は単純に長くなる．しか
し，上述の通り，変数制限法 I で得られた上界値の精度の改善が期待できる．
4.3 計算実験
変数制限法の性能を評価するため，計算実験を行った．計算実験の環境は表 1 の通りであり，計
算実験に用いたネットワークも表 2 の通りである．
計算実験の結果を，表 6 に示す．変数制限法の各種では，その内部で制約式制限法 III を実行し
ているが，そこでの閾値の設定方法は 3.3 節のものと同じである．また，変数制限法 II の Step3
における閾値は，t = 0.5 とした．
時間内に計算が終了した全てのネットワークに対して，変数制限法の各種では制約式制限法 III
よりも高速に上界値が求められた．また，採用した変数の個数が減少し，それによって最終的に加
えた制約式の本数も減少した．よって，変数制限法の各種は，高速性と省メモリ性の両面で優れた
手法である．
15
表 6.
計算結果（その 4）．表の見方は，表 3 に同じである．表中の (∗) は，変数制限法 I が
3600 秒以内に計算を終えず，変数制限法 II に固有の手順である Step3 が始められなかったこ
とを示している．つまり，(∗) の部分の計算結果は，変数制限法 I に同じである．
問題番号
n
制約式制限法 III
変数制限法 I
変数制限法 II
1
34
0.420 / 0s
(0.6k, 1k)
0.420 / 0s
(0.3k, 1k)
0.420 / 0s
(0.4k, 1k)
2
62
0.531 / 2s
(1.9k, 2k)
0.568 / 1s
(0.6k, 1k)
0.534 / 1s
(1.3k, 2k)
3
105
0.528 / 8s
(5.5k, 9k)
0.562 / 3s
(2.0k, 5k)
0.529 / 5s
(4.1k, 8k)
4
115
0.605 / 11s
(6.6k, 11k)
0.617 / 4s
(2.9k, 7k)
0.605 / 6s
(5.4k, 9k)
5
300
0.564 / 38s
(45k, 24k)
0.579 / 13s
(11k, 13k)
0.569 / 19s
(15k, 17k)
6
500
0.546 / 285s
(125k, 77k)
0.558 / 88s
(38k, 45k)
0.549 / 180s
(48k, 57k)
7
700
0.568 / 1717s
(245k, 145k)
0.583 / 311s
(55k, 66k)
0.573 / 886s
(76k, 91k)
8
1000
0.575 / OT
(500k, 215k)
0.584 / 1743s
(128k, 153k)
0.577 / OT
(157k, 185k)
9
1200
0.582 / OT
(719k, 278k)
0.593 / OT
(141k, 170k)
(∗)
10
1500
0.593 / OT
(1124k, 376k)
0.600 / OT
(264k, 301k)
(∗)
しかしながら，変数制限法 I で得られた上界値の精度は悪く，必ずしも有用な手法であるとは言
えない．この問題を解決するための変数制限法 II では，多少の高速性を犠牲にして，より精度の良
い上界値が得られた．その上界値は，制約式制限法 III で得られた上界値と，最大で 1 % 程度しか
変わらなかった．よって，高速性と省メモリ性を考えると，変数制限法 II は優れた手法である．
ところで，変数制限法 I では扱う線形計画問題がかなり小規模に抑えられたものの，時間内に解
ける問題のサイズは制約式制限法 III と比べて大きく変わらなかった．これは，定式化 (LPsparse )
と定式化 (sub LPsparse ) とで，問題の本質的な難しさが変わっていないことによると考えられる．
コミュニティ検出問題はグラフを分割する問題の一種であり，どの枝を切るかを決定する問題であ
ると捉えられる．つまり，変数 xij ({i, j} ∈ E) の値を決定することが問題の本質的な部分である
16
と言える．従って，変数 xij (d(i, j) ≤ 2) のみを採用した問題 (sub LPsparse ) においてそれを決定
できれば，全ての変数を扱う問題 (LPsparse ) でも決定できると考えられる．この計算結果および考
察は，コミュニティ検出問題だけでなくグラフ分割問題などに対しても，重要な示唆を与えている．
5 結論と今後の課題
本論文では，モジュラリティの上界値算出について議論した．そして，定式化 (LPsparse ) と切除
平面法のテクニックを基礎とした手法を提案した．
第一の提案手法は，制約式制限法である．制約式制限法の各種では，定式化 (LPsparse ) を単純に
解くのではなく，暫定最適解が満たしていない制約式を順次加えていった．その結果，最終的に加
えた制約式は (LPsparse ) の制約式の 10 % 程度にとどまり，扱う線形計画問題が小規模に抑えられ
た．これにより，実行時間も短縮された．よって，制約式制限法は，高速性と省メモリ性の両面で
優れた手法である．さらに，最終的に得られる上界値は，定式化 (LPsparse ) を単純に解いたときに
得られる上界値と一致することが保証されている．
第二の提案手法は，変数制限法である．制約式だけでなく採用する変数をも制限することで，
制約式制限法よりも扱う線形計画問題が小規模に抑えられた．得られる上界値の精度が最大で 1
% 程度低下したが，高速性と省メモリ性の向上を考えると，優れた手法である．また，変数制限法
の計算結果に対する考察では，コミュニティ検出問題の本質的な難しさがどこにあるかに迫った．
上記の提案手法により，コミュニティ構造が自明でない 1500 頂点のネットワークにおいて，精度
の良い上界値が得られた．ヒューリスティクスの性能評価を考えると，これは大きな成果である．
最後に，今後の課題を述べる．本論文の提案手法では，定式化 (LP) を解いて得られる上界値よ
りも精度の良い上界値は得られない．よって，コミュニティ構造が顕著でないグラフでは，最適値
からいくらか離れた値しか得られない．ヒューリスティクスのより精密な性能評価のためには，コ
ミュニティ構造が顕著でないグラフでも精度の良い上界値を得られる手法が求められる．
謝辞
本論文を作成するにあたり，宮本裕一郎准教授には，終始温かく指導していただきました．心か
ら感謝いたします．
参考文献
[1] Gaurav Agarwal and David Kempe. Modularity-maximizing graph communities via mathematical programming. The European Physical Journal B - Condensed Matter and Com-
17
plex Systems, 66(3):409–418, December 2008.
[2] Daniel Aloise, Sonia Cafieri, Gilles Caporossi, Pierre Hansen, Sylvain Perron, and Leo
Liberti. Column generation algorithms for exact modularity maximization in networks.
Physical Review E, 82:046112, Oct 2010.
[3] Vincent D. Blondel, Jean-Loup Guillaume, Renaud Lambiotte, and Etienne Lefebvre.
Fast unfolding of communities in large networks. Journal of Statistical Mechanics, page
P10008, 2008.
[4] Danah M. Boyd and Nicole B. Ellison. Social network sites: Definition, history, and
scholarship. Journal of Computer-Mediated Communication, 13(1):article 11, October
2007.
[5] Ulrik Brandes, Daniel Delling, Marco Gaertler, Robert Görke, Martin Hoefer, Zoran
Nikoloski, and Dorothea Wagner. On modularity clustering. IEEE Transactions on
Knowledge and Data Engineering, 20:172–188, 2008.
[6] Thang N. Dinh and My T. Thai. Finding community structure with performance guarantees in complex networks. CoRR, abs/1108.4034, 2011.
[7] Olivier du Merle, Daniel Villeneuve, Jacques Desrosiers, and Pierre Hansen. Stabilized
column generation. Discrete Mathematics, 194(1-3):229–237, 1999.
[8] Santo Fortunato. Community detection in graphs. Physics Reports, 486(3-5):75 – 174,
2010.
[9] Martin Grötschel and Yoshiko Wakabayashi. A cutting plane algorithm for a clustering
problem. Mathematical Programming, 45:59–96, August 1989.
[10] Stanley Milgram. The small-world problem. Psychology Today, 1(1):61–67, 1967.
[11] M. E. J. Newman and Michelle Girvan. Finding and evaluating community structure in
networks. Physical Review E, 69:026113, Feb 2004.
[12] Gang Xu, Sophia Tsoka, and Lazaros G. Papageorgiou. Finding community structures in
complex networks using mixed integer optimisation. The European Physical Journal B Condensed Matter and Complex Systems, 60(2):231–239, 2007.
18
付録 A
本文中の定理の証明
ここでは，2.3 節の定理 1 と定理 2 の証明を記す．以下の証明は，文献 [6] で書かれているものに
若干の補足を加えているだけである．従って，そこで考えられているコミュニティ検出問題の自然
な定式化を挙げておく．また，本文中で述べた通り，定式化 (IPsparse ) と定式化 (LPsparse ) には，
提案された形から若干の変更を施してある．これについても，証明で考えられている形の定式化を
挙げておく．しかしながら，どちらの定式化も本文中のものと同値であるため，本文中の定理を証
明したことになる．
はじめに，頂点 i と頂点 j が同一のコミュニティに属すとき 0，そうでないとき 1 をとる変
数を dij とおく．任意の頂点 i に対して dii = 0 であるから，dii は変数として採用せず，定数
C :=
1
2m
∑n
i=1 qii
を目的関数値に加える．さらに，任意の頂点対 i, j に対して，dij = dji と
qij = qji が成り立つ．従って，変数として dij (i < j) のみを採用すると，コミュニティ検出問題は
1 ∑∑
(IP) : max.
qij (1 − dij ) + C
m i=1 j>i
(36)
s. t. dij + djk − dik ≥ 0
for all i < j < k,
(37)
dij − djk + dik ≥ 0
for all i < j < k,
(38)
− dij + djk + dik ≥ 0 for all i < j < k,
(39)
dij ∈ {0, 1}
(40)
n
n
for all i < j,
として定式化できる．簡単な計算により，式 (36) の目的関数は
1 ∑∑
max. −
qij dij
m i=1 j>i
n
n
(41)
と同値であることがわかる．定式化 (IP) において，式 (40) を dij ∈ [0, 1] に置き換えることで，線
形緩和問題
1 ∑∑
qij dij
m i=1 j>i
n
(LP) : max. −
n
(42)
s. t. dij + djk − dik ≥ 0
for all i < j < k,
(43)
dij − djk + dik ≥ 0
for all i < j < k,
(44)
− dij + djk + dik ≥ 0 for all i < j < k,
(45)
dij ∈ [0, 1]
(46)
for all i < j,
を得る．付録 A を通して，定式化 (IP) と定式化 (LP) とは，いま定義したものを指す．
19
定式化 (IP) に対応する疎な定式化は
1 ∑∑
qij dij
m i=1 j>i
n
(IPsparse ) : max. −
n
(47)
s. t. dki + dkj ≥ dij
for all k (∈ N (i, j)) < i < j,
(48)
dik + dkj ≥ dij
for all i < k (∈ N (i, j)) < j,
(49)
dik + djk ≥ dij
for all i < j < k (∈ N (i, j)) ,
(50)
dij ∈ {0, 1}
for all i < j,
(51)
と書ける．定式化 (IPsparse ) において，式 (51) を dij ∈ [0, 1] に置き換えることで，線形緩和問題
1 ∑∑
(LPsparse ) : max. −
qij dij
m i=1 j>i
n
n
(52)
s. t. dki + dkj ≥ dij
for all k (∈ N (i, j)) < i < j,
(53)
dik + dkj ≥ dij
for all i < k (∈ N (i, j)) < j,
(54)
dik + djk ≥ dij
for all i < j < k (∈ N (i, j)) ,
(55)
dij ∈ [0, 1]
for all i < j,
(56)
を得る．先ほどと同様に，付録 A を通して，定式化 (IPsparse ) と定式化 (LPsparse ) とは，いま定義
したものを指す．
以下の証明の中で，擬似距離なる概念を使うため，その定義をする．関数 d(i, j) = dij が
(1) d(i, j) ≥ 0,
(2) d(i, i) = 0 (i ̸= j である i, j に対しても，d(i, j) = 0 となり得る),
(3) d(i, j) = d(j, i),
(4) d(i, j) ≤ d(i, k) + d(k, j),
を満たすとき，d = {dij } は擬似距離であるという．全ての変数 dij が 0 以上 1 以下を満たす解 d
が擬似距離であることは，d が定式化 (LP) の実行可能解であるための必要十分条件である．
最後に，以下の証明における注意を記す．dij のように変数には添字が付いているが，添字の大
小関係は考慮していない．しかし一方で，定式化 (IP)・(LP)・(IPsparse )・(LPsparse ) では，i < j
を満たす変数 dij のみを採用している．よって，変数 dij について考えるときには，i < j のとき
には dij のまま，j < i のときには dji に置き換えて考える必要がある．
定理 1. (IP) と (IPsparse ) は同一の最適解集合をもつ．
証明. (IP) の最適解が (IPsparse ) の実行可能解であることと，(IPsparse ) の最適解が (IP) の実行可
能解であることを示せばよい．
(IPsparse ) の任意の制約式は (IP) の制約式であるため，(IP) の最適解は (IPsparse ) の実行可能解
である．よって，前者は示された．
20
以下では，後者を示す．はじめに，(IPsparse ) の最適解を d とおく．d が (IP) の実行可能解であ
ることを示すには，d が擬似距離であることを言えばよい．ここで，G = (V, E) は孤立点をもた
ないと仮定して，一般性を失わない．dij = 0 である全ての頂点対 i, j に対してのみ枝 {i, j} をお
いたグラフ Gd = (V, Ed ) をつくる．Cd = {Cd1 , Cd2 , . . . , Cdl } を Gd における連結成分の頂点集合
を要素にもつ集合族とする．つまり，Cdt は t 番目の連結成分の頂点集合である．
補題 1. 全ての連結成分 Cdt は，G = (V, E) において連結な部分グラフを誘導する．
証明. 背理法による．はじめに，ある連結成分 Cdt が G において非連結な部分グラフを誘導すると
仮定する．従って，G において，Cdt をその間に枝を一切もたない頂点部分集合 S と T に分割で
きる．
ここで，P (A, B) := {{i, j} | i ∈ A, j ∈ B} とする．(IPsparse ) の最適解 d から，全ての
{i, j} ∈ P (S, T ) について d′ij = 1 とした新たな解 d′ を構成する．任意の {i, j} ∈ P (S, T )
について Aij = 0 かつ di , dj > 0 より，qij = Aij −
di dj
2m
< 0 が成り立つ．従って，任意の
{i, j} ∈ P (S, T ) について d′ij = 1 とすることで，目的関数値は増加する．ところで，dij = 0 であ
る {i, j} ∈ P (S, T ) が一つもなければ，Gd において Cdt が連結成分とならない．よって，dij = 0
である {i, j} ∈ P (S, T ) は少なくとも一つは存在する．
補題 2. 上で構成した解 d′ は，(IPsparse ) の全ての制約式を満たす．
証明. 任意の制約式 dik + dkj ≥ dij = 0 (k ∈ N (i, j), {i, j} ∈ P (S, T )) において，d′ij = 1 とした
ときに，d′ik + d′kj ≥ d′ij = 1 が成り立つことを示せばよい．すなわち，d′ik = d′kj = 0 でないこと
を示す．
ここで，d′ik = d′kj = 0 であると仮定すると，k ∈ N (i, j) について k ∈ S∩T が成り立つ．すると，
S ∩ T = ∅ であるから，そうした k は存在しないことになる．しかし，dik + dkj ≥ dij (k ∈ N (i, j))
が (IPsparse ) の制約式であるため，k は存在する．ゆえに，矛盾である．
従って，最適解 d から，より高い目的関数値を実現する解 d′ を構成できた．これは矛盾であ
る．
残りは，Gd の各連結成分 Cdt に関して，i, j ∈ Cdt ならば dij = 0 であることを示すだけである．
そこで，クリーク拡大を繰り返し適用することで示していく．つまり，拡大の各段階において，ク
リーク内の全ての頂点対 i, j について dij = 0 であることを仮定する．そして，隣接点を一つ加え
ることでクリークを拡大し，新しくできたクリークにおいても，全ての頂点対 i, j について dij = 0
が成り立つことを示す．以下では，連結成分 Cdt について議論する．他の連結成分においても同様
に議論すればよい．
21
初期操作：はじめに，dij = 0 を満たすグラフ G の枝 {i, j} ∈ E が Cdt 内に存在することを示
す．もし存在しないならば，dij = 0 をみたす Cdt 内の全ての頂点対 i, j に関して Aij = 0 であり，
qij = Aij −
di dj
2m
< 0 となる．いかなる制約式も破ることなく dij = 0 をみたす全ての頂点対 i, j
について dij = 1 とでき，最適解よりも高い目的関数値を実現する実行可能解を構成できる．これ
は矛盾である．よって，いま存在が保証された dij = 0 をみたす枝 {i, j} ∈ E の端点を，頂点数 2
の初期クリークとして採用する．
拡大操作：現在のクリークを Kt とする．Kt = Cdt であれば証明は終わる．そうでなければ，枝
{u, w} が G にも Gd にも属すような頂点 u ∈ Kt と頂点 w ∈ Cdt \ Kt が存在する．（全ての頂点対
/ Ed であるならば，Gd において Kt と Cdt \ Kt は異な
{u, w} ∈ P (Kt , Cdt \ Kt ) について {u, w} ∈
る連結成分となる．これは Cdt が連結であることに矛盾する．そこで，Gd において枝をもつ頂点
対 {u, w} ∈ P (Kt , Cdt \ Kt ) に注目する．その全ての頂点対が G において枝をもたないとすると，
G において Cdt は連結な部分グラフを誘導しない．これは矛盾である．）各頂点 v ∈ Kt \ {u} に関
して dvu + duw ≥ dvw (u ∈ N (v, w)) は (IPsparse ) の制約式であり，Kt の性質より dvu = 0 であ
る．さらに，{u, w} ∈ Ed より duw = 0 であるから，全ての v ∈ Kt について dvw = 0 が成り立
つ．よって，頂点 w ∈ Cdt \ Kt を Kt に加えると，全ての頂点対 i, j ∈ Kt について dij = 0 であ
るという性質を保ちながら，クリークを拡大できる．
頂点集合 Cdt の要素数は高々 n であるため，拡大の操作は Kt = Cdt となったところで終結す
る．
定理 2. (LP) と (LPsparse ) は同一の最適解集合をもつ．
証明. (LP) の最適解が (LPsparse ) の実行可能解であることと，(LPsparse ) の最適解が (LP) の実行
可能解であることを示せばよい．
(LPsparse ) の任意の制約式は (LP) の制約式であるため，(LP) の最適解は (LPsparse ) の実行可
能解である．よって，前者は示された．
以下では，後者を示す．はじめに，(LPsparse ) の最適解を d とおく．この d が (LP) の実行可
能解であることを示す．ここで，全ての枝 {i, j} ∈ E に対して重み wij = dij を与える．頂点対
{i, j} ∈
/ E については wij = ∞ と定める．ここで，枝重み w 上での頂点 i と j の距離を d′ij とす
ると，
(1) 任意の頂点対 i, j について d′ij ≥ dij であり，任意の {i, j} ∈ E について d′ij = dij である，
(2) 任意の頂点対 i, j について d′ij = minnk=1 {d′ik + d′kj } である
という二つの主張が成り立つ．主張 (2) は d′ が擬似距離であることと同値である．主張 (1) の
前半を定式化 (LPsparse ) の三角不等式（制約式）を適用することで示す．頂点 i と j が G の異
なる連結成分に属すとき，d′ij = ∞ であるから d′ij > dij が成り立つ．従って，頂点 i と j は
22
同一の連結成分に属すと仮定する．d′ij は G 上の重み w に関する頂点 i と j の距離であるから，
距離 d′ij = du0 u1 + du1 u2 + · · · + dul−1 ul である道 i = u0 , u1 , . . . , ul = j が存在する．任意の
k = 1, 2, . . . , l について {uk−1 , uk } ∈ E であるから，三角不等式を繰り返し適用すると
dij ≤ du0 u1 + du1 ul ≤ du0 u1 + du1 u2 + du2 ul
≤ · · · ≤ du0 u1 + du1 u2 + · · · + dul−1 ul = d′ij
(57)
となり，主張 (1) の前半は導かれる．また，{i, j} ∈ E であれば d′ij ≤ dij が成り立つため，主張
(1) の後半も導かれる．また，主張 (2) は d′ij の定義より明らかである．
ここで，いま構成した解 d′ の各変数は 1 以下とは限らないことに注意する．そこで，d∗ij :=
min{d′ij , 1} とする．すると，d′ij ≥ 1 のとき d∗ij = 1 ≥ dij であり，d′ij < 1 のとき d∗ij = d′ij ≥ dij
であるから，上の主張 (1) の前半は d∗ij においても成り立つ．また，{i, j} ∈ E について d∗ij = d′ij
であるから，主張 (1) の後半も d∗ij においても成り立つ．さらに，解 d∗ もまた擬似距離である．
なぜならば，d∗ik + d∗kj = min{d′ik , 1} + min{d′kj , 1} ≥ min{d′ik + d′kj , 1} ≥ min{d′ij , 1} = d∗ij が
成り立つからである．
いま，全ての頂点対 i, j について dij = d∗ij であれば，(LPsparse ) の最適解 d は (LP) の全ての
制約式を満たすから，証明は終わる．そうでなければ，上の議論により，ある頂点対 i, j について
d∗ij > dij が成り立つと仮定できる．この (LPsparse ) の実行可能解 d∗ が最適解 d よりも高い目的
関数値を実現することを示し，矛盾を導く．
全ての {i, j} ∈ E について d∗ij = dij であり，ある頂点対 {i, j} ∈
/ E について qij = Aij −
0 かつ d∗ij > dij が成り立つ．従って，解 d∗ は d よりも高い目的関数値を実現する．
23
di dj
2m
<
付録 B
列生成法における補助問題の導出と上界値の計算手順
コミュニティ検出問題の線形緩和問題
1 ∑∑∑
aip ajp qij xp + C
m p=1 i=1 j>i
P
(LPCG ) : max.
s. t.
P
∑
n
n
aip xp = 1
(58)
for all 1 ≤ i ≤ n,
(59)
for all 1 ≤ p ≤ P,
(60)
p=1
xp ≥ 0
を列生成法により解く際の補助問題を導出する．問題 (LPCG ) の双対問題は
(DPCG ) : min.
s. t.
n
∑
i=1
n
∑
yi + C
(61)
1 ∑∑
aip ajp qij
m i=1 j>i
n
aip yi ≥
i=1
n
for all 1 ≤ p ≤ P,
(62)
として定式化できる．定式化 (DPCG ) では，変数が n 個，制約式が 2n − 1 本ある．ところで，
(LPCG ) は自明な解をもつため，実行可能である．ゆえに，強双対定理により，(LPCG ) と (DPCG )
の最適値は一致するため，(DPCG ) の最適値はモジュラリティの上界値である．よって，(DPCG )
を解いてモジュラリティの上界値を得ることを考える．
定式化 (DPCG ) について，制約式を k ≤ 2n − 1 本だけ採用して解く．このとき，暫定最適解
y ∗k = (y1∗k , y2∗k , . . . , yn∗k )T と暫定最適値 z ∗k が得られたとする．この暫定最適解 y ∗k が“最も満
たしていない”制約式を見つける問題は
(QP) : min.
n
∑
i=1
s. t.
1 ∑∑
−
qij ai,k+1 aj,k+1
m i=1 j>i
n
yi∗k ai,k+1
ai,k+1 ∈ {0, 1}
n
for all 1 ≤ i ≤ n,
(63)
(64)
と書ける．この問題 (QP) こそが，所望の補助問題である．本文でも述べたが，(QP) では
y1∗k , y2∗k , . . . , yn∗k は定数であり，a1,k+1 , a2,k+1 , . . . , an,k+1 が変数であることに注意する．（QP）
の最適値が 0 以上であれば，暫定最適解 y ∗k は (DPCG ) の真の最適解であるから，その最適値は
z ∗k であることがわかる．ゆえに，モジュラリティの上界値として z ∗k を返して終了する．それに
対して，(QP) の最適値が 0 未満であれば y ∗k は (DPCG ) の真の最適解ではないから，z ∗k も真の
最適値とは断定できない．従って，k + 1 本目の制約式として
n
∑
1 ∑∑
ai,k+1 yi ≥
ai,k+1 aj,k+1 qij
m i=1 j>i
i=1
n
24
n
(65)
を (DPCG ) に加えて，k ← k + 1 として反復を繰り返す．生成した列の本数 k は 2n − 1 まで増加
し得るが，実験的にはそれよりも圧倒的に少なく終わることが知られている．以上が，列生成法に
よりモジュラリティの上界値を求めるための計算手順である．
25