1. No category

アルゴリズムの進化
線形計画法の場合
前田英次郎
わかみず会
2014.6.18
アルゴリズムの進化（１）
• 出発
– 問題の提示
– モデルの作成
– アルゴリズムの考案
– プログラムの作成
• 計算機の制約
– 遅い，もっとメモリーを
– アルゴリズム，プログラムの工夫
アルゴリズムの進化（２）
• 計算機の進化
– 速くなった，メモリー増えた，２次記憶ついた
– アルゴリズムの工夫
• 問題の大型化
– 問題の本当の顔が少し見えてくる
– 新しい顔に対応したアルゴリズム
• 以下、これの繰り返し
線形計画法
線形計画問題とは、
１次不等式の制約 の下で
１次式の目的関数 を最大にせよ
というものである。
制約条件
å
n
j =1
a ij x j £ bi (i = 1～ m )
x j ³ 0 ( j = 1～ n )
目的関数
z=
å
n
j =1
cjxj
例１(生産計画)
Ｐ社では、ｍ 種類の原料から、ｎ 種類の製品を製造して いる。
１単位の製品ｊには 原料ｉを a ij 単位使用し、 c j の利益を生む。
原料ｉの使用量は b i 単位を越えることはで きない。
利益を最大にする製造 計画をつくれ。
製品ｊの製造量を x j 単位とする。制約条件 は
å
n
j =1
a ij x j £ bi
(i = 1～ m )
x j ³ 0 ( j = 1～ n )
目的関数は
z=
å
となる。
n
j =1
cjxj
例２（輸送問題）
Ｑ社では、ｍヶ所の工 場で製造した商品を、
全国ｎヶ所の営業所に 送って販売している。
工場ｉから営業所ｊへ の輸送単価は c ij 、
工場ｉの搬出量を a i 、営業所ｊの受取量を b jとして、
総輸送費を最小にする 輸送計画をつくれ。
工場ｉから営業所ｊへ の輸送量を x ijとすると、制約条件は
(i = 1～ m )
å j =1 x ij = a i n
( j = 1～ n )
å i =1 x ij = b j m
x ij ³ 0 (i = 1～ m , j = 1～ n )) 総輸送費
z =
å å
m
n
i =1
j =1
を最小にせよ。
c ij x ij
例３（在庫管理）
Ｒ工場では、ある原料の今後ｎ期の在庫計画を作成する所である。
第ｉ期の原料の価格をci、使用量をu i、在庫単価をv、
初期在庫をs0、倉庫の容量をＷとする。
ｎ期末までの原料購入費と在庫費の和を最小にする在庫プランをつくれ。
第ｉ期の購入量をpi、期末の在庫量をsiとする。制約条件は
si -1 + pi - u i = si
0 £ si £ W , pi ³ 0　(i = 1～n)
目的関数は
z = åi =1 (c i pi + vsi )
n
の最小化である。
モデルの標準形
不等式の制約条件
å
n
j =1
a ij x j £ bi に slack 変数 s iを入れて、等式にする 。
å
n
j =1
a ij x j + s i = bi si ³ 0
slack 変数も普通の変数とし て扱うと、モデルは
制約条件
å j =1 aij x j = bi n
xj ³ 0
目的関数
z = å j =1 c j x j
n
となる。
行列／ベクトル表記
ベクトル表
max
示で
z = cx
subject to
Ax = b
x³0
と書くことにする。
a 12
æ a 11
ç
a 22
ç a 21
A=ç


ç
ç
è a m1 a m 2
c = (c 1 c 2 




cn )
a 1n ö
÷
a 2n ÷
 ÷
÷
a mn ø÷
æ x1
ç
ç x2
x =ç

ç
ç
è xn
ö
÷
÷
÷
÷
÷
ø
æ b1
ç
ç b2
b =ç

ç
ç
è bm
ö
÷
÷
÷
÷
÷
ø
基底（basis）
単体法の第１ステップ は基底をつくることで ある。
ｍ個の変数を選んで基 底とする。
基底に入った変数は基 底変数 (basic variable )、
それ以外の変数は非基 底変数 (nonbasic variable )という。
基底変数のベクトルを x B ，非基底変数のベクト ルを x Nと書く。
係数行列 Aの x B 部分だけを集めて B、x N 部分は Nとする。
z の係数 c の x B 部分を c B , x N 部分を c Nとする。
モデルは
max z = c B x B + c N x N
sub. to
Bx B + Nx N = b
x B ³ 0、x N ³ 0
となる。
従属変数、独立変数
max z = c B x B + c N x N
sub . to
Bx B + Nx N = b
x B ³ 0、x N ³ 0
x B = B -1 (b - Nx N ) = B -1b - B -1 Nx N
z = c B ( B -1b - B -1 Nx N ) + c N x N = c B B -1b + (c N - c B B -1 N ) x N
x B と z を x N で表した形になった。 x N = 0 とおくと、
x B = B -1b、z = c B B -1b となる。
B -1 が存在しないと始めら れないので、基底を修 正する。
slack 変数を入れれば B -1 は作れる。
単体法(simplex method)
もし、B -1b ³ 0 であれば、x N = 0、x B = B -1b は制約を満している、
実行可能解、許容解 (feasible solution)である。目的関数
z = c B B -1b + (c N - c B B -1 N ) x N
の変数 x N は値が 0 、係数 (c N - c B B -1 ) の正の要素に対応する 変数を
増加させると、 z も増える。その中の一 つ x p を任意に選ぶ。
x B = B -1b - ( B -1 a p ) x p
x p を増やすと x B の値が変わる。B -1 a p に正の要素があれば、 x B の要素
で負になるものが出る 。非負制約を最初に破 ろうとする変数を基底 から
外し、x p を基底に入れる。
BとNを修正し、反復を続け る。(c N - c B B -1 N )に正の要素がなくなれ ば
最適解 (optimal solution)が得られている。
線形計画法
最初の基底で B -1b ³ 0 でなければ、目的関数を
負の変数の和を最大化するものに変えて単体法を行う。
負の変数がなくなれば本来の目的関数に戻す。
負の変数を根絶できなければ、制約を満たす解はない。
基底に入れる変数の選び方は反復数に影響する。
最適解(optimal solution)
単体法は、許容解から出発、z が増加する方向に移動し、増加する
方向がなくなれば、そこが最適解、という山登り法の１つである。
集合に属する任意の２点を結ぶ線分がその集合に含まれるとき、
その集合を凸集合という。ＬＰの許容領域は凸集合である。
２点 x P , xQを結ぶ線分上の点は、x = px P + qxQ (p + q = 1, p, q ³ 0)
と表される。
Ax P = AxQ = b　＝＞　Ax = pAx P + qAxQ = ( p + q )b = b
単体法の終点より良い点があれば、それに向かう方向は
目的関数を改善する方向である。動ける方向はx N の凸結合なので、
z に正の係数があるはずである。
Simplex Tableau
初期の単体法に、従属変数／独立変数の見方はなかった。
éc B B -1b  c - c B B -1 Aù
ê
ú
 ú
ê  +
-1
ê B -1b
ú

B
A
ë
û
のような表を作って、要素を全部計算した。この表を
Simplex Tableauという。
基 底の入れ替えは、掃き出し計算で、表を修正する。
私が付き合った計算機（１）
• 1963～１９６６ 東京大学計算センター
– OKITAC 5090
– 4000 words, 1 word = 10進12桁 = 48 bits
– 入力：紙テープ、出力：ラインプリンター
• 何と言ってもメモリー不足
– m=30, n= 100 (3000 words) 程度が限界
Ｃｏｍｐａｃｔ Ｔａｂｌｅａｕ
• Simplex Tableau の基底変数の部分は単位行
列である。これをメモリーに置く必要はない。
• 基底に入った変数の場所に基底から出た変
数の情報を入れればいい。
– 表の各列に入っている変数を記録する。
– 各行にその行に１を持つ基底変数を記録する。
• これをCompact Tableau という。
私が付き合った計算機（２）
• 1966～1967 三菱原子力電子計算部
– IBM 7090
– 32000 words 1 word = 36 bits
– 入力：パンチカード，出力：ラインプリンター
– 外部記憶：磁気テープ
• 日本に３台しかない、大型計算機だと聞いた。
• ＬＰ９０というＬＰパッケージがあった。
疎行列(sparse matrix) （１）
• モデルをつくって解こうとすれば、Ａ行列の大
部分の要素が０であることに気付く。
• 入力データは、非零要素だけを (i , j, aij) の
組で表すようになる。
• 番号では、モデルの修正が難しい。
• 変数と制約に名前を付ける。
– (変数名 制約名 要素の値）
疎行列（２）
• モデル記述言語が開発されるまでには、相当
な年月が必要だった。
• 非零要素が少なければ、解きやすいだろう。
疎の度合いを示しておこう。
– density = 非零要素数／（ｍ×ｎ）
– ＬＰソフトにこの density が表示されるようになっ
た。
タブローは必要か
単体法の手順は
１．目的関数の係数を計算し、　（c - c B B -1 A）
正のものから基底に入れる変数 x p を選ぶ。
２．x p の列を計算し、（B -1 A p ）　基底を出る変数をみつける。
３．基底を入れ替える。
必要な物は c - c B B -1 A と B -1 A p、タブローの１行と１ 列だけ。
c B B -1 Aではまず c B B -1を計算する。密なｍベクトルである。
A は疎行列だから、(c B B -1 ) Aの計算量はオーダーｎと考えてよい。
この部分の計算量は c B B -1とB -1 A p が主要部分と考えられる。
タブローの全要素を書き換えるのとは雲泥の差である。
Product Form（１）
基底の入れ替えでは、 B の１列が A p に変わる。
[
]
B -1 A p = t a1' a 2'  a m' としよう。
'
新しい基底行列を B とすると、
é1
ê
ê
ê0
ê
B -1 B ' = ê0
ê0
ê
ê
ê
ë0
 0
 
a1'

 1 a r' -1
 0 a r'
 0 a r' +1
0 
 0

a m'
0  0ù
ú
 0 ú
0  0ú
ú
0  0ú
0  0ú
ú
  ú
0  0ûú
( B -1 B ' ) -1 B -1 = (( B ' ) -1 B ) B -1 = ( B ' ) -1
Product Form（２）
-1
é1 
é1  0 a1' 0  0ù
ê
ú
ê







ê 
ú
ê
ê0 
ê0  1 a r' -1 0  0ú
ê
ú
ê
'
( B -1 B ' ) -1 = =
0
0
0
0


a
ê0 
ú
ê
r
ê0 
ê0  0 a '
0  0ú
r +1
ê
ú
ê

  ú
ê 
ê  
ê0 
ú
ê0  0 a '
0
0

m
ë
û
ë
( B ' ) -1 になる。
この行列をB -1 の前にかければ、
0

- a1' / ar'

1 - a r' -1 / a r'
0
1 / a r'
0 - a r' +1 / a r'


0
- a m' / ar'
0  0ù
ú
  ú
0  0ú
ú
0  0ú
0  0ú
ú
  ú
0  0ûú
これは第r 列であることと、その要素を記憶すれば済む。
反復ごとにこの掃き出し行列を追加するのが、Product Formである。
Product Form（３）
最初にB -1 があるとしれをP0 と書こう。
反復ごとに掃き出し行列 P1 , P2、, P3 が追加される。
B -1 = P?  P3 P2 P1 P0
掃き出し行列は２次記憶に置くことになる。
当時の２次記憶は磁気テープだった。
記録を読みだすには巻き戻さなければならない。
c B B -1 = c B P?  P3 P2 P1 P0 , B -1 A p = P?  P3 P2 P1 P0 A p
左の計算では Pを逆順に使い、右では作成順に使う。
悩ましいが、何かうまい工夫があったかどうか。
再逆転(reinversion)
反復が進み、掃き出し 行列が増えると、計算 時間が増え、
計算誤差もたまってく る。 B -1 を直接求め直すことに なる。
これを再逆転（ reinversio n ）という。
B は疎行列である。工夫 次第で計算時間も、 B -1 つまり P0の
ためのメモリーも格段 に違ってくる。 P0 の出来が良ければ
以後の計算時間が少な くてすみ、作業量が少 なければ、
発生する計算誤差も少 なくて済む。　B -1 も掃き出し行列の積で 表すので、ピボット要 素の選択が
鍵である。
三角化
B の行と列を並べ替えて 、下三角行列にできれ ば、掃き出しで
非零要素が増えること はない。
第１行の非零要素は b11だけだから、第１列の 掃き出しでは、
非零要素は増えない。 第２行の非零要素は b22 だけだから、
第２列の掃き出しでも 非零要素は増えない。 以下同様で、
最後まで増えない。
完全三角化は無理だが 、生きている非要素が １個の行を
見つけることで、左上 隅に三角部分をつくれ る。
同様に非零要素が１個 の列を見つけて右下隅 を三角にする。
ピボット順
左上と右下の三角をとると、中央に正方形が残る。
ピボットをうまく選び、非零要素の増え方を抑えたい。
いろいろなヒューリスチックが使用されていただろうが、
マルコヴィッツの方法というのがあった。
生きている非零要素の少ない行を選ぶ、といったもの
ではなかったか。記憶が曖昧である。
私が付き合った計算機（３）
• 1967～1968 三菱原子力電子計算部
– IBM 360
– 512000 bytes 1 word = 4 bytes
– 入力：パンチカード，出力：ラインプリンター
– ２次記憶：磁気ディスク
ＬＵ分解
ガウスの掃き出しが、非零要素を増やさないことは知られていた。
左上隅から対角要素をピボットにして、まず下三角部分を掃き出し、
上三角行列を残す。第２段階は下から掃き上げる。
掃き上げるときには非零要素は増えない。
第１段階での増え方を抑えるためにピボットをうまく選びたい。
画期的な方法が現れた。P0 はＬＵ分解で求める。反復ごとの
B -1 の修正はＬＵ分解の部分に手を入れる、というものである。
反復ごとのデータの増え具合は少ないらしい。
残念ながら、ＬＰアルゴリズムを追いかける状況ではなく、
詳細については知らない。
完全準三角化
上下の三角を取った後 の部分の処理に決定版 が出た。
ｍ個の非零要素を各列 各行に１つすつ入るよ うに選び、
それらが対角要素にな るように、行と列を並 べ替える。
対角要素を点とするグ ラフを考える。 b ji ¹ 0 なら、
点 i から点 j への 線を引く。３つの場合 がある。
１．点 i から点 j に行き、帰ってこられ る。
２．点 i から点 j に行けるが、帰ってこ られない。
３．点 i から点 j には行けず、点 j から点 i にも行けない。
１ならば、点 i と点 j は同じブロックに属す る。
２ならば、点 i の属するブロックを 点 j のブロックの上に置く 。
３ならば、点 i と点 j は無関係。位置もどう でもいい。
それぞれのブロックは 対角線上の正方形のブ ロックになる。
これがもっとも細かい ブロック分割である。
Density(非零要素密度)
• Density = 非零要素数／（ｍ×ｎ）
• あるとき、この density についてちょっと考え
てみた。
• なんだ、そうだったのか。何十年か何も考え
なかったのだな。
• Density は無意味。Column density を使おう。