大規模な制約なし最小化問題に対する コーダル部分グラフを用いた

特別研究報告書

大規模な制約なし最小化問題に対する コーダル部分グラフを用いた

スパース準ニュートン法

指導教員 福嶋雅夫 教授      山下信雄 助教授     

京都大学工学部情報学科 数理工学コース

平成１５年４月入学 平成１９年３月卒業

黒川 典俊

平成１９年１月３１日提出

大規模な制約なし最小化問題に対する

コーダル部分グラフを用いたスパース準ニュートン法

黒川 典俊

摘要

本報告書では，大規模な制約なし最小化問題に対する準ニュートン法を取り扱う．準ニュートン法が中小 規模の問題に対して有効な解法であることはよく知られているが，各反復で更新される近似ヘッセ行列は 密になるため，大規模な問題に対して準ニュートン法を適用するには何らかの工夫が必要である．大規模な 問題では，目的関数のヘッセ行列は一般に疎な行列となる．その疎性を利用した準ニュートン更新の手法と して，Matrix Completion Quasi-Newton method (行列補完を用いた準ニュートン法，MCQN法)が提案され ている．

ヘッセ行列の疎構造に対応するグラフ（疎構造グラフ）がコーダルグラフとなるとき，MCQN法で更新 される近似行列は，疎な三角行列の積として表すことができる．疎構造グラフがコーダルグラフとならな いとき，従来の

MCQN

問題によっては，ヘッセ行列の非ゼロ要素数と比べて近似行列の非ゼロ要素数が大幅に増え，必要とする記 憶量と計算時間が増大してしまう．

上記の問題を解決するため，本報告書では疎構造グラフのコーダルな部分グラフを用いて準ニュートン 更新を行う手法を提案する．この手法を用いると，近似行列の非ゼロ要素数はヘッセ行列の非ゼロ要素数以 下に抑えることができるので，実装に必要な記憶量と各反復の計算時間を減らすことができる．提案手法 が速く収束するためには，元の疎構造グラフの枝をより多く含むようなコーダル部分グラフを用いること が望ましい．本報告書では，そのような部分グラフを求める工夫についても考察する．また，いくつかの数 値実験結果から提案手法の有効性を議論する．

目 次

1

1

2

2

2.1

. . . . 2

2.2

. . . . 4

2.3

. . . . 5

2.4

. . . . 6

3

8 3.1

. . . . 8

3.2

MCQN

. . . . 9

3.2.1 Ext-MCQN

. . . . 9

3.2.2 Ext-MCQN

. . . . 11

3.3

MCQN

. . . . 12

3.3.1 Del-MCQN

. . . . 12

3.3.2 Del-MCQN

. . . . 14

4

16 5

19 A

22 A.1 MCS

. . . . 22

A.2 RIP

. . . . 22

B

22

C 4

2

24

1 序論

次の制約なし最小化問題を考える．

minimize f (x)

subject to x ∈ R

(1.1)

ここで，

f : R

→ R

2

n

f

0

制約なし最小化問題の解法としては，最急降下法，ニュートン法，準ニュートン法，共役勾配法などがあ る．その中でも準ニュートン法は，実装が容易であり収束の性質も良いので，広く用いられている

[5, 9, 14]．

準ニュートン法は，目的関数のヘッセ行列（もしくはその逆行列）の近似行列を用いて，点列を生成する反 復法である．近似行列の更新公式としては

BFGS

DFP

比較的少ない反復回数で問題

(1.1)

DFP

2

大規模な最小化問題では，一般にヘッセ行列が疎になる．そのため，ヘッセ行列の疎性を利用して近似行 列を更新すれば，記憶容量を大幅に減らすことができる．ヘッセ行列の疎性を利用した近似行列の更新規則 として，行列補完を用いた準ニュートン法

(Matrix Completion Quasi-Newton method,

MCQN

[15]．MCQN

公式

(BFGS

DFP

の

(行列式を最大にする)

関連した性質を持つとき，部分行列の

(行列式を最大にする)

一般にヘッセ行列の疎構造に対応するグラフ（以下，疎構造グラフ）は必ずしもコーダルグラフとは限 らない．そこで，

[15]

MCQN

(Ext-MCQN

)

Ext-MCQN

Ext-MCQN

法は大規模な問題には非常に有効である． 理論的な面からも

Ext-MCQN

DFP

そこで，本報告書では，疎構造グラフのコーダル部分グラフを用いる手法

(Del-MCQN

法は

Ext-MCQN

一方，Del-MCQN法ではヘッセ行列の情報を省いてしまうため，MCQN法がもつ高速な収束性が失われ てしまう可能性がある．そのため，コーダル部分グラフを求める際，疎構造グラフから枝をできるだけ削ら ないようにすることが望ましい．そのような問題は，疎構造グラフの枝数最大のコーダル部分グラフを求 める問題（以後，最大コーダル部分グラフ問題と呼ぶ）として定式化できる．最大コーダル部分グラフ問題 は

NP

[8]，それを近似的に解く実用的なヒューリスティックアルゴリズムがいくつか提案

[2, 13]．本報告書では，Xue

[13]

て

Del-MCQN

本報告書の構成は，以下のとおりである．まず，2節において行列の疎性と関連の深いコーダルグラフの 性質について簡単にまとめる．Del-MCQN法を実装する際に必要となる，最大コーダル部分グラフ問題に

1コーダル拡張については，2章で説明する．

対するアルゴリズムについてもこの節で紹介する．次に，3.1,3.2節で準ニュートン法，ならびに

Ext-MCQN

(Del-MCQN

Del-MCQN

2

2 行列の疎性とコーダルグラフ

2.1

本節では，コーダルグラフに関する基本的な事柄をまとめる．本節で述べるコーダルグラフの性質は，主

に

[1], [4]

[1]

ダルグラフとクリーク木について丁寧に解説している論文である．一方，

[4]

本報告書を通して，Vを頂点の集合，E

⊆ V × V

= (V, E)

v ∈ V

v) < E

本報告書で用いるグラフに関する用語を，以下で定義する．

定義

2.1 (基本的な用語)

• 2

u, v ∈ V

(u, v) ∈ E

• v ∈ V

Adj

(v) = {u ∈ V | (u, v) ∈ E}

Adj(v)

• v ∈ V

v

•

2

• 2

G = (V, E)

G

= (V

, E

)

⊆ V

E

⊆ E

G

• V

G = (V, E)

V

V

E

:=

E ∩ (V

× V

)

G

= (V

, E

)

• C

G = (V, E)

V

G = (V, E)

G

2

C

•

v ∈ V

Adj(v)

G = (V, E)

v

•

•

2

2すなわち，deg(v)=|Adj_G(v)|である．

3G⁰はGの部分グラフになっていることに注意

4すなわち，異なるi,j∈Cのすべてのペアに対して(i,j)∈E

定義

2.2 (コーダルグラフ)

G = (V, E)

4

G

図

1:

図

1

(a)

(b)

4

コーダルグラフの最も基本的な性質は，次の

2

性質

1

性質

2 G = (V, E)

∈ V

\ {v

}

性質

1,2

\ {v

}

v

) ∩ {v

, v

, · · · , v

}

i = 1, 2, · · · , n − 1

(v

, v

, · · · v

)（ただし，n = |V|）を行うことができる．この順序を完全消去順序(perfect elimination ordering,

PEO

性質

3

G = (V, E)

PEO

注意

1 PEO

2

2, 3, 4, 5, 6, 7)

PEO

(3, 2, 4, 6, 7, 1, 5)

PEO

コーダルグラフの極大クリークは，PEOを用いて次のように簡単に列挙することができる．G

= (V, E)

, v

, · · · , v

)

G

PEO

} ∪ Adj(v

)

, v

, · · · , v

}

G = (V, E)

{v

} ∪ Adj(v

) ∩ {v

, v

, · · · , v

}

, i = 1, 2, · · · , n

G

性質

4 G = (V, E)

, · · · , v

)

G

PEO

{C

| r = 1, 2, · · · , l}

C

= {v

} ∪ Adj(v

) ∩ {v

, v

, · · · , v

}

, i = min

v_j∈C_r

j

性質

4

l

n

2

2, 3, 4, 5, 6, 7)

PEO

= {1, 5, 7}, C

= {2, 3, 4, 6}, C

= {4, 6, 7}, C

= {5, 6, 7}

さらに，コーダルグラフの極大クリーク族は，次の条件を満たすように添え字をつけることができること が知られている．

図

2:

5 G = (V, E)

= 1, 2, · · · , l − 1

∃s ≥ r + 1 : C

∩ (C

∪ C

∪ · · · ∪ C

) ( C

(2.1)

{C

| r = 1, 2, · · · , l}

性質

5

Running Intersection Property（以降 RIP

2

{C

, C

, C

, C

}

2.2

本節では，行列の疎性とコーダルグラフの関係について述べる．

まず，n

× n

A

E

E := {(i, j) | A

, 0}

A

V, E

{1, 2, · · · , n}, E := E \ {(i, i) | i = 1, · · · , n}

V

E

G = (V, E)

G = (V , E)

A

例

1

A

3(a)

∗

A

, 0，他は A

= 0

A

3(b)

図

3:

A

図

3(b)

G = (V, E)

定義

2.3

G = (V, E)

•

G

= (V, F)

G

(chordal extension)

かつ

F ⊇ E

5したがって，G=(V,E)はループを持たない．

6行列Aの非対角成分において非ゼロな箇所は，グラフG=(V,E)の隣接行列の非対角成分において非ゼロな箇所と一致すること に注意しておく．

•

G

= (V , F)

G

(chordal deletion)

かつ

F ⊆ E

すなわち，一般のグラフ

G = (V, E)

G

なお，グラフ

G

また，コーダル拡張，コーダル縮小は一意ではない．

例

2

3(b)

G

4(a),(b)

G

(3, 4), (5, 7), (6, 4), (6, 7)

G

(1, 7)

図

4:

3(b)

G = (V, E)

(a)，コーダル縮小 (b)

コーダル拡張は，疎行列のコレスキー分解と密接に関連することが知られている

[11]．行列 A

A = LL

L

A

= 0

L

, 0

0

0

(fill-in)

L

G

= (V, F)（ただし，V := {1, 2, · · · , n}, F := {(i, j) | L

, 0 or L

, 0}, F := F \ {(i, i) | i = 1, 2, · · · , n}）

で表現すると，（コレスキー分解の計算過程で数値的な打ち消しが起こらないと仮定すれば）グラフ

G

行列

A

G = (V , E)

[4]．

充填が起こると，行列

L

A

A

なるべく充填の数を少なくすることが求められる⁷．上で述べたことから，行列

A

A

一方，コーダル縮小は，一般のグラフに含まれる最大のクリークを見つける問題（最大クリーク問題）や 疎な線形方程式を解くアルゴリズムへの応用に関連して研究されている

[2, 13]．一般のグラフに対する最

NP

しかし，与えられたグラフの枝数最小のコーダル拡張を求める問題は

NP

る

[8]．同様に，与えられたグラフの枝数最大のコーダル縮小を求める問題も NP

[8]．そのため，

これらの問題に対する発見的手法がいくつか提案されている．その方法を次節以降で説明する．

2.3

枝数最小のコーダル拡張を求めるヒューリスティクスとしては,最小次数法（最小次数順序法，Minimum

degree ordering algorithm）を用いた手法が有名である．最小次数法は，コレスキー分解における充填の発生

7充填の数は，コレスキー分解の計算過程における行（または列）の順序に依存することが知られている．したがって，充填の数が 少なくなるようにあらかじめ行列Aの行（列）を並べ替えておく．

を観察することにより考案されたグラフの頂点のオーダリング手法である．その基本的な考え方は，「隣接 頂点の少ない順に頂点に番号を振る」ということである．具体的には，次数の小さい頂点をグラフから除去 し，その除去された部分グラフを作る．そして，頂点がすべて除去されるまで同様の操作を繰り返すという アルゴリズムである⁸．その手続きを正確に述べると次のようになる．

¶ ³

最小次数法

Step 1 : V

:= V, E

:= E, S := ∅, k := 1

Step 2 : G

:= (V

, E

)

G

v

(S , v

)

Step 3 : G

v

G

:= (V

, E

)

Step 4 : k = k + 1

> |V|

2

µ ´

このアルゴリズムを適用すると，頂点の順番は

S = (v

, · · · , v

)

i

i

最小次数順序

S

PEO

Parter

[10]

S = (v

, v

, · · · , v

)

¶ ³

順番

S

G = (V, E)

Step 1 : k := 1

= G

= E

Step 2 : k > n

v

G

Step 3 : Adj

(v

)

E

Step 4 : k := k + 1

2

µ ´

このアルゴリズムにより得られるグラフ

G

= (V, E

)

PEO

G

2.4

一般のグラフに含まれる枝数最大のコーダル部分グラフを見つける問題

(最大コーダル部分グラフ問題)

NP

いま，グラフ

G = (V, E)

{v

, v

, · · · , v

}

:=

{v

| j ≥ k}

V

E

⊆ V

× V

G

= (V

, E

)

G

まず，1頂点のみからなるグラフ

G

:= ({v

}, ∅)

= (V

, E

)

G

, v

, v

, · · · , v

}

G

= (V

, E

)

Step 0 : G

:= G

Step 1 : G

V

v

Step 2 : v

G

C ⊆ V

v

E

8グラフにおいて次数の少ない頂点は，対応する行列において非ゼロ要素の少ない列に対応する．

9与えられたグラフの枝数最大のコーダル縮小を求める問題は，与えられたグラフに含まれる枝数最大のコーダル部分グラフを求 める問題と等価であることに注意．

¶ ³

【条件】

V

(v

) := Adj

(v

) ∩ V

(v

)

G

µ ´

こうして得られたグラフ

G

k = n, n − 1, n − 2, · · · , 1

G

= (V

, E

)

= V

G

G

(v

, v

, · · · , v

)

PEO

3

G

G

である．

枝数のできるだけ多いコーダル縮小を求めるため，頂点の選び方にさまざまな工夫がなされている．ここ では，Xueのアルゴリズム

[13]

Suc

(v

) := {v

| j > i and v

∈ Adj

(v

)}

とする．

¶ ³

Edge-maximal chordal subgraph (Xue [13])

:

G = (V, E)

出力

:

G

= (V, E

)

Step 0 : k := n

:= ∅， E

:= ∅， peo := ∅， U := V

∈ V

= ∅, s(v) = 0

Step 1 : k = 1

:= (V

, E

)

peo = (v

, · · · , v

)

PEO

Step 2 : s(v) = max{s(u) | u ∈ U}

v ∈ U

1

V

:= V

∪{v}， peo := (v, peo)， U := U \{v}

とする. E^k−1

:= E

∪ {(v, u) | u = t(v) or u ∈ Adj

(v) ∩ Suc

(t(v))}

v

v

Step 3 : ∀u ∈ Adj

(v) ∩ U

:= 1 + |Suc

(v) ∩ Adj

(u)|

もし，ru

≥ s(u)

v， s(u) := r

k − 1

µ ´

注意

2 Step2

max{s(u) | u ∈ U}

G

[13]．

Xue

G

定理

2.1 ([13],Theorem 3.1) G = (V, E)

Xue

G

= (V, E

)

peo = (v

, · · · , v

)

PEO

G

なお，Xueのアルゴリズムの計算量は，O(|V|

+ P

deg(v) + P

deg

(v)) ∼ O(∆ × |E|)

[13]．

∆ = max

deg(v) | v ∈ V

10すなわち，グラフG^k=(V^k,E^k)の頂点のうち，与えられたグラフG=(V,E)においてvk−1と隣接した頂点の集合をV^k(vk−1)と する．

3 コーダル部分グラフを用いた準ニュートン法

3.1

この節では，準ニュートン法とその問題点について述べる．

準ニュートン法は,目的関数のヘッセ行列（もしくはその逆行列）の近似行列を用いて点列を生成する反 復法である．xkを現在の反復点とし，Bkを目的関数のヘッセ行列

∇

f (x

)

は正定値対称行列であるとしよう．準ニュートン法では，まず探索方向

p

p

= −B

∇ f (x

)

で定める¹¹．次に，pk方向に直線探索を行って，次回の反復点

x

x

= x

+ α

p

により定める．ただし，αkは非負の実数（ステップ幅）である．

行列

B

H

= B

H

H

y

= s

を満たすように更新される．ただし，

s

= x

− x

, y

= ∇ f (x

) − ∇ f (x

) (3.1)

BFGS

DFP

DFP

H

[3]．

min

ψ(H

HH

1 2

k

) subject to Hy

= s

H = H

H 0

(3.2)

ただし，ψ

: R

→ R

ψ(A) = trace(A) − ln det(A)

で定義される狭義凸関数である．また，H

0 (H 0)

特に，Hが正定値行列で，条件：

s

y

> 0 (3.3)

を満たすとき¹³，問題

(3.2)

H

[3]．

H

= H

− H

y

(H

y

)

y

H

y

+ s

s

s

y

(3.4)

また，BFGS法は，次式で与えられる更新公式（BFGS公式）を用いて近似行列の更新を行う方法である．

H

= H

− H

y

s

+ s

(H

y

)

s

y

+

 

 1 + y

H

y

s

y

 

 s

s

s

y

(3.5) DFP

BFGS

H

(s

)

f (x)

H

O(n

)

n

11探索方向pkは，必ず目的関数の降下方向になる．なぜなら，Bkの正定値性より∇f (xk)^Tpk=−∇f (xk)B⁻¹_k ∇f (xk)<0が成り立つ からである．

12このとき，探索方向pkはpk=−Hk∇f (xk)と行列とベクトルの積で計算できる．

13fが強凸関数のとき，条件(3.3)は任意の2点xk,xk+1に対して成り立つ．また．f が非凸であるとき，ステップ幅αkがWolfe の条件を満たすようにすれば，条件(3.3)は満たされることが知られている[9]．

3.2

MCQN

3.2.1 Ext-MCQN

次元数

n

(3.2)

min

ψ(H

HH

1 2

k

) subject to Hy

= s

H = H

H 0

(H

)

= 0 ∀(i, j) < F

(3.6)

ここで，集合

F

⊆ V ×V

F ≈ E := {(i, j) |

x ∈ R

[∇

f (x)]

, 0 }

ただし，V

= {1, 2, · · · , n}

i) ∈ F, i = 1, 2, · · · , n (2) (i, j) ∈ F ⇒ ( j, i) ∈ F

F

F = E

残念なことに，問題

(3.2)

(3.6)

(3.6)

H

[15]．

¶ ³

(MCQN

Step 1:

(3.4)

(3.5)

H

(H

)

, ∀(i, j) ∈ F

Step 2: (H

)

, (i, j) ∈ F

min

ψ(H

HH

1 2

k

) subject to H = H

H

= (H

)

∀(i, j) ∈ F H 0

(H

)

= 0 ∀(i, j) < F

(3.7)

の最適解を

H

µ ´

問題

(3.7)

(3.6)

Hy

= s

= (H

)

, ∀(i, j) ∈ F

(3.7)

(3.6)

= V × V

(3.7)

(3.2)

H

H

問題

(3.7)

= H

, ∀(i, j) ∈ F

j) < F

H

B

を参照．）問題

(3.7)

せることが知られている．

¶ ³

【コーダル条件】 

頂点集合を

V

F := F \ {(i, i) | i = 1, · · · , n}

G

= (V, F)

µ ´

問題の疎性を利用するためには

F = E

G = (V, E)

[15]

F

14グラフG⁰にはループはないことに注意．

提案されている．

条件

1 F ⊇ E．

条件

2

F

(3.8) (3.8)

2

F

G = (V, E)

G

= (V, F)

G

G

MCQN

F \ E

F

G

G

2.3

G

さて，コーダル条件を満たす集合

F

(3.7)

1

補題

3.1 ( [15],Theorem 2) s

y

> 0

)

= 0, ∀(i, j) < F

(3.7)

maxmize det(H)

subject to H

= H

∀(i, j) ∈ F H = H

H 0

(3.9)

さらに，この解は唯一つであり，解として得られた行列

H

)

= 0, ∀(i, j) < F

注意

3

(3.9)

, ∀(i, j) ∈ F

F

(3.9)

G

= (V, F)

2

{C

| r = 1, 2, · · · , l}

• F = S

r=1

C

× C

• RIP（式 (2.1)）

= 1, 2, · · · , l − 1

∃s ≥ r + 1 : C

∩ (C

∪ C

∪ · · · C

) ( C

を満たす．

この極大クリーク族を用いて，{Cr

}

{S

| r = 1, 2, · · · , l}, {U

| r = 1, 2, · · · , l}

S

:= C

\ (C

∪ C

∪ · · · ∪ C

), r = 1, · · · , l (3.10) U

:= C

∩ (C

∪ C

∪ · · · ∪ C

), r = 1, · · · , l (3.11)

= ∪

S

S

∩ S

= ∅

, S

, · · · , S

1, 2, · · · , n

G

P

15定義により，S1の各要素は非空であり，S1に属する各頂点はG⁰の単体的頂点である．さらに，Siの各要素は，Si∪Si+1∪ · · ·Sl

により導かれるG⁰の部分グラフの単体的頂点である．したがって，S1,S2,· · ·,Slの順にそれぞれの要素を取り出して並べた順序は G⁰のPEOになる．こうして求めたPEOが(1,2,· · ·,n)となるように頂点の番号を付け替えることになる．