■平行体の体積とグラミアン(その71)
4次元以上の高次元についてのイメージはあまり働かないのがふつうである.結局,ある高次元図形を網羅的に描きたいならば,多面体の頂点の座標をすべて求めそれらを正しく結ぶしかない.今回のコラムではシンプレックス法を応用して多面体のすべての頂点を求めることを考えてみたい.
===================================
【1】シンプレックス法
線形計画法は最小の費用の下で得られるべき利益を最大にすることを目的とした最適化の1つの方法です.商業上重要というだけでなく,数学のいろいろな分野でも使われます.
線形計画法におけるシンプレックス法は,有限個の線形不等式
Ax≧b, maximize cx
の解の集合から最適解を求めるもので,ダンツィックにより1947年から開発が進められました.その手順はまず実行可能な頂点を見つけ,そこから多面体の辺に沿って,線形の目的関数cxが増加するように進んでいくというものです.
d次元の場合,多面体の各頂点からはd本の辺がでているのですが,しかし,目的関数が増加するような次の辺をどうやって選べばよいかは現在でも完全には解決されていない問題となっています(多項式時間アルゴリズムがあるか?).
また,一口にシンプレックス法といっても,データが得られた状況によってそれぞれに最適な手法が用意されていて,それらを適宜選択し使い分ける必要が出てきます.シンプレックス法の教科書はこれまで多数出版されていますが,たいていはxが非負の場合のみを解説しています.しかし,ここではおきまりのxが常に非負の値をとるようにした制約条件ではなく,下限Lと上限Uの間の値をとるようにした制約条件の場合(有界変数法)を考えることになります.
変数xが非負の値をとるようにした制約条件の場合,非基底変数の値は0ですが,下限Lと上限Uの間の値をとるようにした制約条件の場合,非基底変数の値はLかUのどちらかとなります.
===================================
具体例をあげると,原点を通る3つの超平面
a11x1+a12x2+a13x3+a14x4+a15x6+a16x6+a17x7=0
a21x1+a22x2+a23x3+a24x4+a25x6+a26x6+a27x7=0
a31x1+a32x2+a33x3+a34x4+a35x6+a36x6+a37x7=0
のカーネルで,上限・下限を
li≦xi≦si i=1~7
で定めた7次元矩形領域を定義域とするポリトープのすべての頂点を求める問題をシンプレックス法を応用して解こうというものです.
しかし,このままの形では扱いにくく,スラック変数を導入して,
a11x1+a12x2+a13x3+a14x4+a15x6+a16x6+a17x7 +x8 =0
a21x1+a22x2+a23x3+a24x4+a25x6+a26x6+a27x7 +x9 =0
a31x1+a32x2+a33x3+a34x4+a35x6+a36x6+a37x7 +x10 =0
[a11,a12,a13,a14,a15,a16,a17,1,0,0] [0]
A=[a21,a22,a23,a24,a25,a26,a27,0,1,0] b=[0]
[a31,a32,a33,a34,a35,a36,a37,0,0,1] [0]
の形にします.係数行列Aの中には単位行列Eが含まれるのですが,この形を正準形を呼びます.
シンプレックス法の計算過程では基底変数の入れ替えがあり,その都度,係数行列
A=[B+C]
に基底逆行列B^(-1)をかけて
B^(-1)A=[B^(-1)B+B^(-1)C]=[E+B^(-1)C]
として正準形に直します.
基底行列から出される列,基底行列へ加入する列を検出しますが,基底変数の変更があるごとに基底逆行列を計算していたのでは手間がかかりますから,実際のプログラムでは改訂シンプレックス法を利用して逆行列計算を行わなくてもいいように工夫します.
最初の頂点が求まると,あとは有界変数法のサブルーチンを利用して芋づる式にすべての頂点が求まります.したがって,最初の頂点を確実に求めることが重要な問題になるのですが,そのために2段階シンプレックス法の第1段階を利用して,
c=[a11+a21+a31,a12+a22+a32,・・・,a17+a27+a37]
とおきます.
===================================
【2】幅優先探索
以下,
[参]フバータル「線形計画法」啓学出版
の幅優先探索に則って記しますが,係数行列Aがm×(n+m)行列のとき,
xkが基底のとき ek=1
xkが非基底で下限Lのとき ek=0
xkが非基底で上限Uのとき ek=2
とするベクトル
e^0=[e1,e2,・・・,en+m]
は可能基底分割を表していて,これが決まれば頂点のひとつが一意に決定されます.
そして,e^0には最大でn個の近傍が存在するのですが,e^0のすべての近傍
e^1,e^2,・・・,e^n
をつくる,次にe^1のすべての近傍をつくる,e^2のすべての近傍をつくる,・・・という作業を続けていって,もうこれ以上新たな近傍がつくれなくなったとき,すべての頂点が求まったことになります.
その際,つくられた近傍がすでにリストにあるのがどうかを調べる必要が生じます.同じものがなければリストに登録し,あれば登録しない・・・.しかし,データ数が多くなるとリストの検索には非常に時間がかかってしまいます.そこで,可能基底分割e^kが得られるごとに2分挿入法(binary insertion)で挿入位置を求めます.これで結果的にリストをソートしたことになるので,検索時間を短縮できるというわけです.
===================================