欲張り法の続きです。今回はグラフ関連のアルゴリズムを紹介します。それから、ナップザック問題を例題に「動的計画法」を説明します。
今回は閉路が存在しない無向グラフを考えてみましょう。次の図を見てください。
Ａ●─────●Ｂ Ａ● ●Ｂ
│ ／│ ＼ ／
│ ／ │ ＼ ／
│ ／ │ Ｅ●──●Ｆ
│ ／ │ ／ ＼
│／ │ ／ ＼
Ｃ●─────●Ｄ Ｃ● ●Ｄ
(1) 閉路が存在 (2) 閉路の無いグラフ
図 1 : 無向グラフ
無向グラフの場合、頂点を自分自身と連結するような長さ 3 以上の経路を「閉路」といいます。上図 (1) では、A - B - C - A と B - C - D - B と A - B - C - D - A という 3 つの閉路が存在します。A - C - A のような経路は閉路として考えません。これに対して、(2) では閉路は存在しません。
いま、n 個の頂点を全てつなぐようなグラフを考えると、最低 n - 1 本の辺が必要になります。最小の辺の本数で全ての頂点を接続したグラフを「展張木 (spanning tree)」と呼ぶます。(2) も展張木となっていますが、6 個の頂点を全てつなげばいいので、このほかにも展張木となる接続パターンが多数存在します。このとき、辺に重みを与えて、その合計が最小となる展張木を MST (Minimum Spanning tree) といいます。
簡単な例を示しましょう。次の図を見てください。
Ａ───８───Ｂ
／ ＼ ／ ＼
５ １ ６ ２
／ ＼ ／ ＼
Ｃ───６───Ｄ───４───Ｅ
＼ ／ ＼ ／
４ ３ ５ ３
＼ ／ ＼ ／
Ｆ───７───Ｇ
図 2 : 経路図
上の経路図は、各地点をケーブルでつなぐ場合にかかるコストを表しています。各地点をケーブルでつなぐ場合、7 個の頂点を全てつなげばいいので、展張木として考えることができます。
たとえば、図 3 に示すような展張木を考えることができます。この場合、最も安いコストでケーブルを接続する経路は MST で表すことができます。上の 2 つの例は MST ではありません。
MST を求める方法では「プリム (R. C. Prim) のアルゴリズム」と「クラスカル (J. B. Kruskal) のアルゴリズム」の 2 つが有名です。この 2 つは、「欲張り法」に基づいたアルゴリズムです。最初に、プリムのアルゴリズムから説明しましょう。
Ａ──────Ｂ Ａ───────Ｂ
＼ ＼
＼ ＼
＼ ＼
Ｃ──────Ｄ──────Ｅ Ｃ──────Ｄ──────Ｅ
／ ／
／ ／
／ ／
Ｆ──────Ｇ Ｆ───────Ｇ
総コスト＝３０ 総コスト＝２９
図 3 : 展張木(spanning tree)の例
基本的な考え方は、MST に属する頂点とそうでないものとに分け、まだ MST に属していない頂点から、最もコストの安く済む頂点を付け加えていく、という欲張り法になります。最初に選ぶ頂点はどれでもかまいません。どの頂点でも最後は MST に属するで、どの頂点からスタートしても大丈夫です。
問題は最小のコストになる頂点の選び方です。MST に属する頂点が複数あり、まだ選んでいない頂点も複数ある場合、その中から最小コストの辺を選ばなくてはいけません。このあたりの工夫が、このアルゴリズムを実装する上でのポイントになります。
それでは、具体的にプログラミングしていきましょう。グラフは隣接行列を使って表すことにします。
リスト : 隣接行列
weight = [
[INF, 8, 5, 1, INF, INF, INF],
[ 8, INF, INF, 6, 2, INF, INF],
[ 5, INF, INF, 6, INF, 4, INF],
[ 1, 6, 6, INF, 4, 3, 5],
[INF, 2, INF, 4, INF, INF, 3],
[INF, INF, 4, 3, INF, INF, 7],
[INF, INF, INF, 5, 3, 7, INF]
]
INF はつながっていないことを示すグローバル変数で、値は MAX_VALUE - 1 (0xffffff) としました。あとは各辺の重みを設定するだけです。
次に、グローバル変数を定義します。
リスト : グローバル変数の定義 # 定数 MAX_VALUE = 0x1000000 INF = MAX_VALUE - 1 MAX_SIZE = 7 # グローバル変数 closest = [0] * MAX_SIZE low_cost = [0] * MAX_SIZE
プリムのアルゴリズムの場合、まだ MST に属していない頂点の中から、最小コストの辺とつながっている頂点を選択しなくてはなりません。それを管理するために、その頂点につながっている辺の中で最小コストを変数 low_cost に格納し、対応する頂点を変数 closest に格納します。頂点を選択するたびに、この 2 つのグローバル変数を更新していくことで、プリムのアルゴリズムは動作します。
それでは、本題である探索アルゴリズムを見ていきましょう。
リスト : 探索（プリムのアルゴリズム）
def search():
# 初期化
for i in range(MAX_SIZE):
low_cost[i] = weight[0][i]
# 頂点の選択
for i in range(1, MAX_SIZE):
min = low_cost[1]
k = 1
for j in range(2, MAX_SIZE):
if min > low_cost[j]:
min = low_cost[j]
k = j
print('{:c} - {:c}'.format(k + ord('A'), closest[k] + ord('A')))
low_cost[k] = MAX_VALUE
for j in range(1, MAX_SIZE):
if low_cost[j] < MAX_VALUE and low_cost[j] > weight[k][j]:
low_cost[j] = weight[k][j]
closest[j] = k
最初に、グローバル変数 low_cost を初期化します。頂点 A からスタートするので、low_cost には A と各頂点を結ぶ辺の重みをセットします。closest には接続されている頂点をセットしますが、最初は頂点 A を対象とするので 0 に初期化しておきます。
次の for ループで頂点を一つずつ選択していきます。その次の for ループで、low_cost の中から最小コストの辺とつながっている頂点を選びます。変数 k に選択した頂点をセットし、変数 min にそのコストをセットします。頂点 A からスタートするので、最初に B の値を k と min にセットして、C から G までの値と比較します。
最小コストの頂点を見つけたら、それと対応する頂点 colsest[k] をいっしょに表示します。そして、low_cost[k] に MAX_VALUE をセットして、頂点 k が選択されて MST に属したことを表します。
最後の for ループで low_cost と closest を更新します。今、頂点 k が MST に加えられたので、その頂点とほかの頂点の接続を調べて low_cost の値を更新します。まず、調べる頂点が MST に属していないことを確認します。それから、重み weight[k][j] と最小コスト low_cost[j] を比較して、重みの方が小さければ low_cost[j] と closest[j] の値を更新します。これで low_cost には常に最小コストが格納されていることになります。
実際に実行してみると、結果は次のようになりました。
実行結果
D - A Ａ Ｂ
F - D ＼ ＼
C - F １ ５
E - D ＼ ＼
B - E Ｃ Ｄ───４───Ｅ
G - E ＼ ／ ／
３ ２ ６
＼ ／ ／
Ｆ Ｇ
【注意】数字は順番を表す
図 4 : 「プリムのアルゴリズム」実行結果
この場合、総コストは 13 になります。プリムのアルゴリズムは、頂点の個数を N とすると N - 1 回のループの中に N - 1 回のループがあるので、実行時間は N の 2 乗に比例します。したがって N が大きくなると、このアルゴリズムでは効率が不十分になります。そのような場合、次に説明するクラスカルのアルゴリズムを使った方がよいでしょう。クラスカルのアルゴリズムは、辺の本数を M とすると実行時間は M * log2 M に比例します。
次はクラスカルのアルゴリズムを説明します。基本的な考え方は、辺が一本もない状態から始めて、一本ずつ辺を追加していくことで MST を構成します。このとき、重みの小さい辺から加えていく「欲張り法」を使うのですが、単純に辺を追加していくだけでは、途中で閉路ができてしまって展張木の条件を満たさなくなる可能性があります。そこで、「MST にならないような辺は追加しない」という条件を付け加えます。
具体的には、次のように考えます。
Ａ●──●Ｂ Ａ●──●Ｂ Ａ●──●Ｂ
│ ＼ │
│ ＼│
Ｃ● ●Ｄ Ｃ● ●Ｄ Ｃ● ●Ｄ
（１） （２） （３）
Ａ●──●Ｂ Ａ●──●Ｂ
│
│
Ｃ●──●Ｄ Ｃ●──●Ｄ
（４） （５）
図 5 : 辺の選択
たとえば、上図 (1) のように辺 A - B を選んだとします。これを一つのグループとして管理します。このとき、各頂点に自分の属するグループ番号をつけておきます。最初は頂点に別々の番号をつけておき、辺を選んだときに両端のグループ番号をチェックします。違う番号であれば、その辺を選択して頂点のグループ番号を統一します。同じグループ番号ならば、閉路を形成するので、その辺は選択しません。
たとえば、(2) の場合では AB と D は違うグループなので、辺 B - D を追加できます。そして、D を AB と同じグループ番号に統一します。次に、(3) のように辺 A - D を追加しようとした場合、A と D は同じグループ番号だから、この辺を追加することはできないことがわかります。
(4) のように、まず辺 A - B を選択し、次に辺 C - D を選択した場合は、AB と CD は違うグループ番号になります。したがって、(5) のように B と D では番号が違うので、辺 B - D を選択することができます。そして、頂点の個数 n に対して n - 1 本の辺を選択した時点で、MST の作成は終了です。
それでは、クラスカルのアルゴリズムを実装してみましょう。このアルゴリズムでは、重みが最小である辺を探す処理で、MST に属していない辺を線形探索するようでは、その真価を発揮することはできません。ここは、最小の重みを持つ辺を、簡単に取り出せるようなデータ構造を採用すべきでしょう。この条件にぴったりのデータ構造が「ヒープ」です。ヒープについては拙作のページ「二分木とヒープ」をお読みください。
それから、もう一つ考えておかなくてはならない問題があります。それは、二つのグループを一つにまとめる併合（マージ）処理です。いちばん簡単な方法は、次のようにプログラムできます。
リスト : グループの併合 (力任せバージョン)
def merge(g1, g2):
for i in range(NODE_SIZE):
if group[i] == g2: group[i] = g1
配列 group に各頂点のグループ番号を格納します。グループ g1 と g2 を併合する場合、g2 に属する頂点を全て g1 に書き換えるだけです。この方法はとても簡単ですが、頂点の個数が多くなると効率が悪くなるのが欠点です。
そこで、今回はグループを Python の set で表すことにしましょう。set のメソッド update() を使うと、グループを簡単に併合することができます。ただし、このままでは頂点からグループ番号を簡単に求めることができません。たとえば、頂点 3 のグループ番号を調べるには、グループ番号 0 から順番に set の要素をチェックする必要があります。
この欠点を補うため配列 group を用意し、そこに頂点のグループ番号を格納することにします。グループを一つに併合するときに group を書き換えておけば、頂点が属するグループ番号を group から直ぐに求めることができます。
それでは、プログラムを作りましょう。辺の重みを基準にヒープを作成するので、隣接行列を使わずに、ずばり「辺」そのものをデータ構造として定義します。
リスト : 辺の定義
class Edge:
def __init__(self, p1, p2, weight):
self.p1 = p1
self.p2 = p2
self.weight = weight
def __gt__(x, y): return x.weight > y.weight
def __le__(x, y): return x.weight <= y.weight
クラス名は Edge としました。インスタンス変数 p1 と p2 は両端の頂点を格納します。変数 weight は、その辺の重みを格納します。それから、Edge に特殊メソッド __gt__() と __le__() を定義します。これで「二分木とヒープ」で作成したプライオリティキュー (PQueue) をそのまま利用することができます。
次は辺のデータを定義します。
リスト : 辺のデータ (p1, p2, weight)
edge_data = [
(0, 1, 8), (0, 2, 5), (0, 3, 1),
(1, 3, 6), (1, 4, 2), (2, 3, 6),
(2, 5, 4), (3, 4, 4), (3, 5, 3),
(3, 6, 5), (4, 6, 3), (5, 6, 7)
]
タプルの要素 (p1, p2, weight) は Edge のインスタンス変数 p1, p2, weight と同じで、両端の頂点と重みを表しています。
次はグループ用のグローバル変数を定義します。
リスト : グループ用グローバル変数
group = list(range(NODE_SIZE))
group_set = [{x} for x in range(NODE_SIZE)]
group は各頂点のグループ番号を格納します。最初は、全て異なる番号を割り当てておきます。ここでは、自分自身の番号で初期化しています。処理が進むとグループ番号は書き換えられていき、最後には一つのグループ番号に統合されます。group_set はグループを表す set を格納します。最初は頂点番号とグループ番号が一致するので、set は頂点番号で初期化します。
次は Edge のオブジェクトを生成してプライオリティキューを初期化する関数 make_edge() を作ります。
リスト : Edge のオブジェクトを生成
def make_edge():
global edges
edges = PQueue()
for p1, p2, w in edge_data:
edges.push(Edge(p1, p2, w))
グローバル変数 edge_data から Edge のオブジェクトを生成して、プライオリティキュー (ヒープ) に追加します。プライオリティキューはグローバル変数 edges にセットします。あとは、Edge(p1, p2, w) で Edge のオブジェクトを生成して、プライオリティキューのメソッド push() で追加します。これで edges には重みを基準にしたヒープが構成されます。
それでは、探索本体を行う関数 search() を作りましょう。
リスト : クラスカルのアルゴリズムによる探索
def search():
i = 0
while i < NODE_SIZE - 1:
e = edges.pop()
if group[e.p1] != group[e.p2]:
merge(e.p1, e.p2)
print('{:c} - {:c}'.format(e.p1 + ord('A'), e.p2 + ord('A')))
i += 1
while ループで辺を一つずつ選択していきます。変数 i は選択した辺の個数を管理し、i が NODE_SIZE - 1 になったら while ループを終了します。まず、edges から最小の重みを持つ辺をメソッド pop() で取り出して変数 e にセットします。それから、両端のグループ番号をチェックします。違う番号であれば、その辺を選択することができます。関数 merge を呼び出してグループ番号を一つにまとめます。そのあと、選択した辺を print() で表示し、変数 i の値を一つ増やします。PQueue を使っているので、プログラムはとても簡単になりました。
次は、グループの併合を行う関数 merge() を作ります。
リスト : グループを併合する
def merge(node1, node2):
gs1 = group_set[group[node1]]
gs2 = group_set[group[node2]]
if len(gs1) <= len(gs2):
gs1.update(gs2)
gs3 = gs2
g = group[node1]
else:
gs2.update(gs1)
gs3 = gs1
g = group[node2]
for x in gs3:
group[x] = g
gs3.clear()
頂点 node1 が属するグループを gs1 に、頂点 node2 が属するグループを gs2 にセットします。そして、個数の少ないグループを多いほうに併合します。これは set のメソッド update を使うと簡単です。個数の少ないグループを gs3 に、多いグループの番号を g にセットします。最後に、gs3 のグループ番号を g に書き換えて、gs3 をメソッド clear で空にします。
最後に、make_edge() でデータを作成してから search() を呼び出します。
実行結果は次のようになりました。
実行結果
A - D Ａ Ｂ
B - E ＼ ＼
D - F １ ２
E - G ＼ ＼
C - F Ｃ Ｄ───６───Ｅ
D - E ＼ ／ ／
５ ３ ４
＼ ／ ／
Ｆ Ｇ
【注意】数字は順番を表す
図 6 : 「クラスカルのアルゴリズム」実行結果
プリムのアルゴリズムと違い、辺を追加して一つのグループにまとめていく様子がよくわかると思います。
それでは、プリムとクラスカルのアルゴリズムを証明してみましょう。これらのアルゴリズムでは、展張木が作られるのは自明なので、これが MST になることを示します。
頂点を 2 つの集合 U と V に分けます。この 2 つの集合を連結する辺 e は、展張木の特性から 1 本しかありません。もし、辺 e が MST に属するのならば、その重みは U と V を接続する辺の中で、最小でなければなりません。
この証明は簡単です。辺 f を含む MST が存在すると仮定します。辺 f の重みが辺 e よりも大きいのであれば、その重みの総和は辺 e を含む MST よりも大きくなるので、仮定に反することになります。つまり、2 つのグループを連結する辺の中で、最小の重みを持つ辺を使って展張木を構成すれば、それは MST になるのです。
┌─────── ＭＳＴ ───────┐
│ Ｕ Ｖ │
│ ┌───┐ ┌───┐ │
│ │ │ ｅ │ │ │
│ │ Ａ─┼───────┼─Ｘ │ │
│ │ ・ │ │ ・ │ │
│ │ ・ │ │ ・ │ │
│ │ ・ │ ｆ │ ・ │ │
│ │ Ｂ─┼───────┼─Ｙ │ │
│ │ │ │ │ │
│ └───┘ └───┘ │
└───────────────────┘
図 7 : アルゴリズムの証明
プリムのアルゴリズムの場合、MST に属する頂点とそうでない頂点に分けて、その中で最小の重みを持つ辺を選択しているので、構成された展張木は MST になります。クラスカルのアルゴリズムでは 2 つの違うグループ同士を統合するときに、最小の重みを持つ辺を使って連結しているので、これも MST となるのです。
クラスカルのアルゴリズムは Union-Find というデータ構造を使うともっと簡単に実装することができます。Union-Find の詳しい説明は拙作のページ「Union-Find」をお読みください。
プログラムは次のようになります。
リスト : クラスカルのアルゴリズムによる探索
def search():
i = 0
u = UnionFind(NODE_SIZE)
while i < NODE_SIZE - 1:
e = edges.pop()
if u.find(e.p1) != u.find(e.p2):
u.union(e.p1, e.p2)
print('{:c} - {:c}'.format(e.p1 + ord('A'), e.p2 + ord('A')))
i += 1
最初に UnionFind で NODE_SIZE 個の集合を生成します。メソッド find() は集合を表す値を返します。返り値が同じ値であれば、e.p1 と e.p2 は同じ集合に属していることがわかります。そうでなければ、メソッド union() で e.p1 と e.p2 を併合します。Union-Find は find() と union() の処理を高速に行うことができるデータ構造なので、今回作成した merge() よりも効率的に処理することができます。
#
# prim.py : プリムのアルゴリズム
#
# Copyright (C) 2007-2022 Makoto Hiroi
#
# 定数
MAX_VALUE = 0x1000000
INF = MAX_VALUE - 1
MAX_SIZE = 7
# グローバル変数
closest = [0] * MAX_SIZE
low_cost = [0] * MAX_SIZE
# 隣接行列
weight = [
[INF, 8, 5, 1, INF, INF, INF],
[ 8, INF, INF, 6, 2, INF, INF],
[ 5, INF, INF, 6, INF, 4, INF],
[ 1, 6, 6, INF, 4, 3, 5],
[INF, 2, INF, 4, INF, INF, 3],
[INF, INF, 4, 3, INF, INF, 7],
[INF, INF, INF, 5, 3, 7, INF]
]
# 探索
def search():
# 初期化
for i in range(MAX_SIZE):
low_cost[i] = weight[0][i]
# 頂点の選択
for i in range(1, MAX_SIZE):
min = low_cost[1]
k = 1
for j in range(2, MAX_SIZE):
if min > low_cost[j]:
min = low_cost[j]
k = j
print('{:c} - {:c}'.format(k + ord('A'), closest[k] + ord('A')))
low_cost[k] = MAX_VALUE
for j in range(1, MAX_SIZE):
if low_cost[j] < MAX_VALUE and low_cost[j] > weight[k][j]:
low_cost[j] = weight[k][j]
closest[j] = k
# 実行
search()
#
# krus.py : クラスカルのアルゴリズム
#
# Copyright (C) 2007-2022 Makoto Hiroi
#
from pqueue import *
# 定数
NODE_SIZE = 7
EDGE_SIZE = 12
# 辺のデータ (p1, p2, weight)
edge_data = [
(0, 1, 8), (0, 2, 5), (0, 3, 1),
(1, 3, 6), (1, 4, 2), (2, 3, 6),
(2, 5, 4), (3, 4, 4), (3, 5, 3),
(3, 6, 5), (4, 6, 3), (5, 6, 7)
]
# グループ用変数
group = list(range(NODE_SIZE))
group_set = [{x} for x in range(NODE_SIZE)]
# 辺の定義
class Edge:
def __init__(self, p1, p2, weight):
self.p1 = p1
self.p2 = p2
self.weight = weight
def __gt__(x, y): return x.weight > y.weight
def __le__(x, y): return x.weight <= y.weight
# 辺のデータを作成
def make_edge():
global edges
edges = PQueue()
for p1, p2, w in edge_data:
e = Edge(p1, p2, w)
edges.push(e)
# 併合
def merge(node1, node2):
gs1 = group_set[group[node1]]
gs2 = group_set[group[node2]]
if len(gs1) <= len(gs2):
gs1.update(gs2)
gs3 = gs2
g = group[node1]
else:
gs2.update(gs1)
gs3 = gs1
g = group[node2]
for x in gs3:
group[x] = g
gs3.clear()
# 探索
def search():
i = 0
while i < NODE_SIZE - 1:
e = edges.pop()
if group[e.p1] != group[e.p2]:
merge(e.p1, e.p2)
print('{:c} - {:c}'.format(e.p1 + ord('A'), e.p2 + ord('A')))
i += 1
# 実行
make_edge()
search()
#
# krus.py : クラスカルのアルゴリズム (Union-Find 版)
#
# Copyright (C) 2012-2022 Makoto Hiroi
#
from pqueue import *
from unionfind import *
# 定数
NODE_SIZE = 7
EDGE_SIZE = 12
# 辺のデータ (p1, p2, weight)
edge_data = [
(0, 1, 8), (0, 2, 5), (0, 3, 1),
(1, 3, 6), (1, 4, 2), (2, 3, 6),
(2, 5, 4), (3, 4, 4), (3, 5, 3),
(3, 6, 5), (4, 6, 3), (5, 6, 7)
]
# 辺の定義
class Edge:
def __init__(self, p1, p2, weight):
self.p1 = p1
self.p2 = p2
self.weight = weight
def __gt__(x, y): return x.weight > y.weight
def __le__(x, y): return x.weight <= y.weight
# 辺のデータを作成
def make_edge():
global edges
edges = PQueue()
for p1, p2, w in edge_data:
e = Edge(p1, p2, w)
edges.push(e)
# 探索
def search():
i = 0
u = UnionFind(NODE_SIZE)
while i < NODE_SIZE - 1:
e = edges.pop()
if u.find(e.p1) != u.find(e.p2):
u.union(e.p1, e.p2)
print('{:c} - {:c}'.format(e.p1 + ord('A'), e.p2 + ord('A')))
i += 1
# 実行
make_edge()
search()
動的計画法は難しそうな名前がついていますが、これに惑わされてはいけません。動的計画法は、大きな問題を小さな問題に分けて、それを一つ一つ解いていくことで大きな問題を解く方法です。
問題によっては、小さな問題に分割していくと同じ小問題が何回も現れる場合があります。この場合、同じ問題を何回も解くよりも、その解を表などに保存しておいて、必要なときにその表から答を求めた方が、効率良く問題を解くことができるはずです。
どうせ小問題を解かなければならないのであれば、はじめから必要になりそうな小問題を解いて表を埋めておいたほうが、プログラムを作りやすい場合もあります。このように、与えられた問題を解くために小問題の表を埋めてしまう、というのが「動的計画法」の基本的な考え方です。
簡単な例題として、組み合わせの数 \({}_n \mathrm{C}_r\) を求めるプログラムを作ってみましょう。\({}_n \mathrm{C}_r\) を求めるには、次の公式を使えば簡単です。
皆さんお馴染みの公式ですね。(1) と (2) の公式を使うと簡単に (高速に) 答えを求めることができます。ただし、(3) の公式をそのままプログラムすると二重再帰になるので、大きな値を求めると時間がかかってしまいます。実際にプログラムを作って確かめてみましょう。
リスト : 組み合わせの数
import time
def comb(n, r):
if r == 0 or n == r:
return 1
else:
return comb(n - 1, r) + comb(n - 1, r - 1)
if __name__ == '__main__':
s = time.time()
print(comb(22, 11))
print(time.time() - s)
s = time.time()
print(comb(24, 12))
print(time.time() - s)
s = time.time()
print(comb(26, 13))
print(time.time() - s)
705432 0.2252349853515625 2704156 0.7691750526428223 10400600 2.870633840560913
このように、公式 (3) を使うと時間がかかります。公式 (1), (2) を使えば高速に答えを求めることができますが、今回は動的計画法の例題として、あえてこのプログラムの高速化に挑戦してみましょう。
公式からわかるように、\({}_n \mathrm{C}_r\) の値は \({}_{n-1} \mathrm{C}_r\) と \({}_{n-1} \mathrm{C}_{r-1}\) を足したものです。n = 0 から順番に組み合わせの数を求めて表に格納しておけば、n が大きな値でも簡単に求めることができるはずです。プログラムは次のようになります。
リスト : 組み合わせの数 (動的計画法[1])
import time
def comb_dp_sub(n, r):
global comb_table
if r == 0 or n == r:
return 1
else:
return comb_table[n - 1][r] + comb_table[n - 1][r - 1]
def comb_dp(n, r):
global comb_table
comb_table = [[0] * (n + 1) for _ in range(n + 1)]
for i in range(n + 1):
for j in range(i + 1):
comb_table[i][j] = comb_dp_sub(i, j)
return comb_table[n][r]
if __name__ == '__main__':
s = time.time()
print(comb_dp(30, 15))
print(time.time() - s)
s = time.time()
print(comb_dp(100, 50))
print(time.time() - s)
大域変数 comb_table に値を格納する 2 次元配列をセットします。関数 comb_dp_sub() は comb_table から組み合わせの数を求めます。組み合わせの数を (n, r) で表すことにすると、関数 comb_dp() は (0, 0) から順番に、(1, 0), (1, 1), (2, 0), (2, 1), (2, 2) ... と組み合わせの数を求めて comb_table にセットします。ようするに「パスカルの三角形」を作っていくわけです。最後に comb_table から (n, r) の値を求めて返します。
それでは実際に試してみましょう。
155117520 0.0004627704620361328 100891344545564193334812497256 0.001954317092895508
動的計画法の効果はとても高いですね。なお、表は二次元配列ではなく一次元配列 (ベクタ) で済ますこともできます。たとえば、(6, 3) を求めてみましょう。次の図を見てください。
0 1 2 3 4 5 6
-------------------------
0 [ 1 1 1 1 1 1 1 ]
1 [ 1 1 1 1 1 1 1 ]
＼|
2 [ 1 2 1 1 1 1 1 ]
＼|＼|
3 [ 1 3 3 1 1 1 1 ]
＼|＼|＼|
4 [ 1 4 6 4 1 1 1 ]
＼|＼|＼|＼
5 [ 1 5 10 10 5 1 1 ]
＼|＼|＼|＼|＼|
6 [ 1 6 15 20 15 6 1 ]
図 8 : パスカルの三角形
最初にベクタの内容を 1 に初期化します。n = 0, 1 の場合はこのままで大丈夫です。あとは図のように、隣の要素を足し算するだけです。3 番目の要素の値 20 が (6, 3) の値になります。プログラムは次のようになります。
リスト : 組み合わせの数 (動的計画法[2])
def comb_dp1(n, r):
table = [1] * (n + 1)
for i in range(1, n + 1):
for j in range(i - 1, 0, -1):
table[j] += table[j - 1]
return table[r]
ベクタの値を書き換えていくので、ベクタの後方から計算していくことに注意してください。前方から計算すると、値がおかしくなります。
なお、関数 comb() はメモ化関数を使って高速化することもできます。次のリストを見てください。
リスト : 組み合わせの数 (メモ化)
# メモ化関数
def memoize(f):
table = {}
def func(*args):
if not args in table:
table[args] = f(*args)
return table[args]
return func
# メモ化
comb = memoize(comb)
プログラムの説明は拙作のページ「再帰定義」をお読みください。
動的計画法とメモ化関数は、どちらも表を使って計算を高速化する方法です。動的計画法は小さな問題の解を積み上げていく、つまりボトムアップな方法なのに対し、メモ化関数は大きな問題を小さな問題に分割していく、トップダウンな方法 (分割統治法) といえます。
もうひとつ、数値計算の例を取り上げましょう。整数 n を 1 以上の自然数の和で表すことを考えます。これを「整数の分割」といいます。整数を分割するとき、同じ自然数を何回使ってもかまいませんが、並べる順序が違うだけのものは同じ分割とします。簡単な例を示しましょう。次の図を見てください。
─┬─ 6 : 6
│
├─ 5 ─ 1 : 5 + 1
│
├─ 4 ┬ 2 : 4 + 2
│ │
│ └ 1 ─ 1 : 4 + 1 + 1
│
├─ 3 ┬ 3 : 3 + 3
│ │
│ ├ 2 ─ 1 : 3 + 2 + 1
│ │
│ └ 1 ─ 1 ─ 1 : 3 + 1 + 1 + 1
│
├─ 2 ┬ 2 ┬ 2 : 2 + 2 + 2
│ │ │
│ │ └ 1 ─ 1 : 2 + 2 + 1 + 1
│ │
│ └ 1 ─ 1 ─ 1 ─ 1 : 2 + 1 + 1 + 1 + 1
│
└─ 1 ─ 1 ─ 1 ─ 1 ─ 1 ─ 1 : 1 + 1 + 1 + 1 + 1 + 1
図 9 : 整数 6 の分割
6 の場合、分割の仕方は上図のように 11 通りあります。この数を「分割数」といいます。分割の仕方を列挙する場合、整数 n から k 以下の整数を選んでいくと考えてください。まず、6 から 6 を選びます。すると、残りは 0 になるので、これ以上整数を分割することはできません。次に、6 から 5 を選びます。残りは 1 になるので、1 を選ぶしか方法はありません。
次に、4 を選びます。残りは 2 になるので、2 から 2 以下の整数を分割する方法になります。2 から 2 を選ぶと残りは 0 になるので 2 が得られます。1 を選ぶと残りは 1 になるので、1 + 1 が得られます。したがって、4 + 2, 4 + 1 + 1 となります。同様に、6 から 3 を選ぶと、残りは 3 から 3 以下の整数を選ぶ方法になります。
6 から 2 以下の整数を選ぶ方法は、残り 4 から 2 以下の整数を選ぶ方法になり、そこで 2 を選ぶと 2 から 2 以下の整数を選ぶ方法になります。1 を選ぶと 4 から 1 以下の整数を選ぶ方法になりますが、これは 1 通りしかありません。最後に 6 から 1 を選びますが、これも 1 通りしかありません。これらをすべて足し合わせると 11 通りになります。
整数 n を k 以下の整数で分割する総数を求める関数を p(n, k) とすると、p(n, k) は次のように定義することができます。
たとえば、p(6, 6) は次のように計算することができます。
p(6, 6) => p(0, 6) + p(6, 5)
=> 1 + p(1, 5) + p(6, 4)
=> 1 + 1 + p(2, 4) + p(6, 3)
=> 1 + 1 + 2 + 7
=> 11
p(2, 4) => p(-2, 4) + p(2, 3)
=> 0 + p(-1, 3) + p(2, 2)
=> 0 + 0 + p(0, 2) + p(2, 1)
=> 0 + 0 + 1 + 1
=> 2
p(6, 3) => p(3, 3) + p(6, 2)
=> p(0, 3) + p(3, 2) + p(4, 2) + p(6, 1)
=> 1 + p(1, 2) + p(3, 1) + p(2, 2) + p(4, 1) + 1
=> 1 + 1 + 1 + p(0, 2) + p(2, 1) + 1 + 1
=> 1 + 1 + 1 + 1 + 1 + 1 + 1
=> 7
分割数を求める関数 partition_number() は、関数 p(n, k) を使うと次のようにプログラムすることができます。
リスト : 分割数
import time
def part_num(n, k):
if n < 0 or k < 1:
return 0
if n <= 1 or k == 1:
return 1
else:
return part_num(n - k, k) + part_num(n, k - 1)
def partition_number(n):
return part_num(n, n)
if __name__ == '__main__':
s = time.time()
print(partition_number(50))
print(time.time() - s)
s = time.time()
print(partition_number(60))
print(time.time() - s)
s = time.time()
print(partition_number(70))
print(time.time() - s)
204226 0.08218860626220703 966467 0.3387720584869385 4087968 1.6849095821380615
関数 part_num() は p(n, k) の定義をそのままプログラムしただけです。ただし、このプログラムは二重再帰で何度も同じ値を求めているため実行速度はとても遅くなります。part_num() をメモ化すると高速化することができますが、大きな値を計算すると Python のスタックがオーバーフローしてしまいます。
動的計画法を使うと、大きな値でも高速に計算することができます。次の図を見てください。
k 1 : [1, 1, 1, 1, 1, 1, 1] 2 : [1, 1, 1+1=2, 1+1=2, 2+1=3, 2+1=3, 3+1=4] => [1, 1, 2, 2, 3, 3, 4] 3: [1, 1, 2, 1+2=3, 1+3=4, 2+3=5, 3+4=7] => [1, 1, 2, 3, 4, 5, 7] 4: [1, 1, 2, 3, 1+4=4, 1+5=6, 2+7=9] => [1, 1, 2, 3, 5, 6, 9 5: [1, 1, 2, 3, 5, 1+6=7, 1+9=10] => [1, 1, 2, 3, 5, 7, 10] 6: [1, 1, 2, 3, 5, 7, 10+1=11] => [1, 1, 2, 3, 5, 7, 11]
大きさ n + 1 のベクタを用意します。ベクタの添字が n を表していて、p(n, 1) から順番に値を求めていきます。p(n, 1) の値は 1 ですから、ベクタの要素は 1 に初期化します。次に、p(n, 2) の値を求めます。定義により p(n, 2) = p(n - 2, 2) + p(n, 1) なので、2 番目以降の要素に n - 2 番目の要素を加算すれば求めることができます。あとは、k の値をひとつずつ増やして同様の計算を行えば p(n, n) の値を求めることができます。
プログラムは次のようになります。
リスト : 分割数 (動的計画法)
def partition_number2(n):
table = [1] * (n + 1)
for k in range(2, n + 1):
for m in range(k, n + 1):
table[m] += table[m - k]
return table[n]
説明をそのままプログラムしただけなので、とくに難しいところはないと思います。それでは実際に試してみましょう。1000 の分割数を求めてみました。
24061467864032622473692149727991 0.06751608848571777
動的計画法の効果はとても高いですね。
最後に、簡単な例題として「ナップザック問題」を取り上げます。ナップザック (knapsack) とは辞書を引いてみると、ランドセルのような背中にせおう四角形の袋や箱のことを意味します。ここでは物を入れる袋と簡単に考えてください。
ここで、ナップザックの中に品物を詰め込むことを考えてみます。一つのナップザックと複数の品物が与えられたとき、袋に詰めた品物の合計金額が最大になるような選び方を求めることが「ナップザック問題」です。ここでは、同じ品物を何個も選んでもいいのですが、ナップザックの大きさをオーバーしてはいけません。
実はこの「ナップザック問題」が「NP 問題」なのです。世の中にはさまざまな問題が山積していますが、スーパーコンピュータを使っても解くのに数億年かかる、というような難問が「NP 問題」です。
これは、厳密に解を求めようとすると、全ての場合について総当たりで調べるしか方法がなく、データ数が多くなると時間がかかるため、現実的な時間では解答を出すことができないというものです。品物の詰め方が難問の一つ、といわれてもピンとこないと思いますが、ナップザック問題は品物の種類が増えるにしたがって、その組み合わせ方が爆発的に増えるのです。
ところが、幸いなことに「ナップザック問題」は実用的には解決済みの問題と考えられています。とくに有名なのが「動的計画法」を用いた解法です。ナップザックと品物の大きさを整数値に限定すれば、動的計画法を用いることで厳密解を求めることができるのです。
それでは具体的に、ナップザック問題に動的計画法を適用してみましょう。ナップザックの大きさは 10 で、次の 3 種類の品物を詰め込むことにします。
(Ａ) 大きさ ４ ：金額 ６ (Ｂ) 大きさ ３ ：金額 ４ (Ｃ) 大きさ １ ：金額 １
まず、大きさが 0 から 10 までのナップザックを用意します。これらのナップザックに品物を順番に詰め込んで、その合計金額を配列に格納しておきます。この配列は品物を詰め込んでいない状態 (金額は全て 0) に初期化します。
最初に品物 A を詰め込みます。このとき、小さなナップザックから順番に詰め込んでいきます。
Ａを詰め込む
↓
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│／│／│／│／│／│／│／│／│／│／│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│０│０│０│０│０│０│０│０│０│０│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
↓ ↓
金額[0] + 6 = 6 > 0 金額[4] に 6 をセット
図 10 : ナップザック (1)
品物 A が入る大きさのナップザックから詰め込みます。品物を入れるときは、それより小さいナップザックには、その時点での最適な値が決定されていると考えます。ナップザック (4) に品物 A を詰め込む場合、A の大きさだけ空いているナップザック (0) の状態に詰め込めば、ナップザック (4) の最適な値を求めることができるはずです。このように、前に計算された値を使うところが動的計画法の特徴なのです。
具体的には、金額[0] に A の金額 6 を足した値を計算し、金額[4] より大きくなれば金額[4] をその値に更新します。もし、金額[4] より小さいのであれば金額[4] は更新しません。つまり、品物Ａはナップザック (4) には詰め込まないのです。ほかの組み合わせの方が正解だというわけです。
詰め込んだ品物を記憶しておくため、もう一つ配列を用意して、そこに追加した品物の種類を格納しておきます。ナップザックの中身全てを記憶しておく必要はありません。この配列を使って、あとからナップザックの中身を求めることができます。
次に、品物 B を詰め込んでいきます。
Ｂを追加する
↓
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│／│／│Ｂ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│０│０│４│６│６│６│６│12│12│12│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
↓ ↓
金額[1] + 4 = 4 < 6 金額[4] は更新しない
図 11 : ナップザック (2)
まず、ナップザック (3) に B が詰め込まれます。これは品物 A の場合と同じですね。次に、ナップザック (4) に B を詰めようとします。その値を計算すると 4 となり、金額[4] の値 6 より小さいので、B は詰め込みません。ナップザック (5) の場合も同様です。
Ｂを追加する
↓
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│／│／│Ｂ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│Ａ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│０│０│４│６│６│６│６│12│12│12│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
↓ ↓
金額[3] + 4 = 8 < 6
金額[6] を 8 に更新し、選択[6] にＢをセットする
図 12 : ナップザック (2)
次はナップザック (6) に B を詰めます。値を計算すると 8 になり、今度は金額[6] の値 6 より大きくなります。つまり、A を詰め込むよりも B を詰め込む方が金額が高くなるのです。金額[6] と選択[6] の値を更新します。ナップザック (7) の場合も同様ですね。
あとは、順番に同じことを繰り返して、配列の値を更新していきます。そして、品物 C を最後まで詰め込むと、次のようになります。
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│Ｃ│Ｃ│Ｂ│Ａ│Ａ│Ｂ│Ｂ│Ａ│Ａ│Ｂ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
Ｃ Ｃ
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
金額│０│１│２│４│６│６│８│10│12│12│14│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴─┘
７ 13
図 13 : ナップザック (3)
このときの金額[10] の 14 が答となります。この状態からナップザックに詰め込まれた品物を求めます。
０ １ ２ ３ ４ ５ ６ ７ ８ ９ 10
┌─┬─┬─┬─┬─┬─┬─┬─┬─┬─┬─┐
選択│／│Ｃ│Ｃ│Ｂ│Ａ│Ｃ│Ｂ│Ｂ│Ａ│Ｃ│Ｂ│
└─┴─┴─┴─┴─┴─┴─┴─┴─┴─┴┼┘
↑ │ │ │
└───────┴←────┴←────┘
図 14 : ナップザック (4)
まず、選択[10] にセットされた品物を取り出します。この場合は B ですね。次は、10 からＢの大きさ 3 を引いた 7 のときに追加された品物を取り出します。この場合も B ですね。同様に、7 から 3 を引いた 4 のときに追加された品物を求めます。これはＡですね。4 から A の大きさを引くと 0 になるので、これ以上品物は入っていません。したがって、ナップザックの中には A が 1 個、B が 2 個入っていることがわかります。
それでは、実際にプログラミングしましょう。まず、品物を表すクラスを定義します。
リスト : 品物の定義
class Item:
def __init__(self, price, size, count):
self.price = price
self.size = size
self.count = count
インスタンス変数 price が金額で、size が大きさを表します。count は、ナップザックの中の品物を数えるときに使用します。
次は品物を表すデータを設定します。
リスト : データの設定 knap_size = 10 item_data = [(4, 6), (3, 4), (1, 1)] item_list = [Item(price, size, 0) for size, price in item_data]
kanp_size はナップザックの大きさを格納します。item_data は品物のデータで、タプルは (size, price) を表しています。item_data からクラス Item のオブジェクトを生成し、それを配列 (Python のリスト) に格納して item_list にセットします。
次は、ナップザック問題の解を求める関数 knapsack() を作ります。
リスト : ナップザック問題の解法
def knapsack():
gain = [0] * (knap_size + 1)
choice = [None] * (knap_size + 1)
for item in item_list:
size = item.size
price = item.price
j = size
while j <= knap_size:
space = j - size
new_value = price + gain[space]
if gain[j] < new_value:
gain[j] = new_value
choice[j] = item
j += 1
# 結果を表示する
i = knap_size
while choice[i]: # 修正 (2012/01/22)
item = choice[i]
item.count += 1
i -= item.size
for item in item_list:
print('price = {:d}, size = {:d}, count = {:d}'.format(item.price, item.size, item.count))
最初に、金額を表す配列 gain と選択した品物を格納する配列 choice を用意します。gain は 0 で初期化します。次の for ループで、品物を一つずつ item_list から取り出します。品物の金額と大きさを price と size にセットし、次の while ループで、size から knap_size まで gain と choice を更新していきます。
まず、size 分だけ空けたナップザックの大きさを space にセットします。次に、品物を追加した場合の金額を new_value にセットします。この new_value が gain[j] より大きくなれば、gain と choice を更新します。
最後にナップザックの中身を表示します。while ループで、ナップザックに詰め込まれた品物の個数を種類毎に数えます。まず choice に格納された品物を item にセットします。品物の個数 item.count を一つ増やして、i から size を引いて、次に詰め込まれた品物へ移ります。あとは次の for ループで、item_list にセットされたデータを表示するだけです。
それでは実行結果を示します。
price = 6, size = 4, count = 1 price = 4, size = 3, count = 2 price = 1, size = 1, count = 0
正常に動作していますね。興味のある方はいろいろなデータで試してみてください。