（その１）（その２）では，時間をかけたあげく最終的には素朴な解決法に落ち着いたのですが，その後，意外な事実がわかりました．高次元図形の場合はマッピングをもってしても求めたい領域の１部しか描けないのです．
マッピング解は求める領域を完全ではないにしてもほぼカバーできるものと思い込んでいたので，これは衝撃的な発見でした．（その１）のマッピング解は４重ループにより計算したものですが，７重ループにしたところでこの事情は変わらないでしょう．
すぐ頭の中に「モーザーのパラドックス」が思い浮かびました．すなわち，人間の直観や勘は１次元から高々３次元までの低次元の世界ではよく働くものの高次元の世界では直観は働きにくく，しかも次元が大きくなるに従い解の配位は非常に複雑となり，われわれが３次元空間でイメージするものとは大きく異なってくるという「教訓」です．→［補］
４次元以上の高次元についてのイメージはあまり働かないのがふつうであって，結局，ある高次元図形を網羅的に描きたいならば，多面体の頂点の座標をすべて求めそれらを正しく結ぶしかありません．今回のコラムではシンプレックス法を応用して多面体のすべての頂点を求めてみることにします．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】シンプレックス法
線形計画法は最小の費用の下で得られるべき利益を最大にすることを目的とした最適化の１つの方法です．商業上重要というだけでなく，数学のいろいろな分野でも使われます．
線形計画法におけるシンプレックス法は，有限個の線形不等式
Ａｘ≧ｂ， maximize ｃｘ
の解の集合から最適解を求めるもので，ダンツィックにより１９４７年から開発が進められました．その手順はまず実行可能な頂点を見つけ，そこから多面体の辺に沿って，線形の目的関数ｃｘが増加するように進んでいくというものです．
ｄ次元の場合，多面体の各頂点からはｄ本の辺がでているのですが，しかし，目的関数が増加するような次の辺をどうやって選べばよいかは現在でも完全には解決されていない問題となっています（多項式時間アルゴリズムがあるか？）．
また，一口にシンプレックス法といっても，データが得られた状況によってそれぞれに最適な手法が用意されていて，それらを適宜選択し使い分ける必要が出てきます．シンプレックス法の教科書はこれまで多数出版されていますが，たいていはｘが非負の場合のみを解説しています．しかし，ここではおきまりのｘが常に非負の値をとるようにした制約条件ではなく，下限Ｌと上限Ｕの間の値をとるようにした制約条件の場合（有界変数法）を考えることになります．
変数ｘが非負の値をとるようにした制約条件の場合，非基底変数の値は０ですが，下限Ｌと上限Ｕの間の値をとるようにした制約条件の場合，非基底変数の値はＬかＵのどちらかとなります．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
具体的には原点を通る３つの超平面
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
で定めた７次元矩形領域を定義域とするポリトープのすべての頂点を求める問題をシンプレックス法を応用して解こうというものです．
しかし，このままの形では扱いにくく，スラック変数を導入して，
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]
の形にします．係数行列Ａの中には単位行列Ｅが含まれるのですが，この形を正準形を呼びます．
シンプレックス法の計算過程では基底変数の入れ替えがあり，その都度，係数行列
Ａ＝［Ｂ＋Ｃ］
に基底逆行列Ｂ^(-1)をかけて
Ｂ^(-1)Ａ＝［Ｂ^(-1)Ｂ＋Ｂ^(-1)Ｃ］＝［Ｅ＋Ｂ^(-1)Ｃ］
として正準形に直します．
*INVERSE.BASE.MATRIX:
基底行列から出される列，基底行列へ加入する列はそれぞれ，
*XS.DETECT:
*XR.DETECT:
で検出されます．しかし，基底変数の変更があるごとに基底逆行列を計算していたのでは手間がかかりますから，実際のプログラムでは改訂シンプレックス法を利用して逆行列計算を行わなくてもいいようにしています．
*PIVOT.MATRIX:
最初の頂点が求まると，あとは有界変数法のサブルーチンを利用して芋づる式にすべての頂点が求まります．したがって，最初の頂点を確実に求めることが重要な問題になるのですが，そのために２段階シンプレックス法の第１段階を利用して，
c=[a11+a21+a31,a12+a22+a32,･･･,a17+a27+a37]
とおいています．
*DATAIN1:
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】幅優先探索
以下，
［参］フバータル「線形計画法」啓学出版
の幅優先探索に則って記しますが，係数行列Ａがｍ×（ｎ＋ｍ）行列のとき，
ｘkが基底のとき ｅk＝１
ｘkが非基底で下限Ｌのとき ｅk＝０
ｘkが非基底で上限Ｕのとき ｅk＝２
とするベクトル
ｅ^0＝［ｅ1，ｅ2，・・・，ｅn+m］
は可能基底分割を表していて，これが決まれば頂点のひとつが一意に決定されます．
そして，ｅ^0には最大でｎ個の近傍が存在するのですが，ｅ^0のすべての近傍
ｅ^1，ｅ^2，・・・，ｅ^n
をつくる，次にｅ^1のすべての近傍をつくる，ｅ^2のすべての近傍をつくる，・・・という作業を続けていって，もうこれ以上新たな近傍がつくれなくなったとき，すべての頂点が求まったことになります．
その際，つくられた近傍がすでにリストにあるのがどうかを調べる必要が生じます．同じものがなければリストに登録し，あれば登録しない・・・．しかし，データ数が多くなるとリストの検索には非常に時間がかかってしまいます．そこで，可能基底分割ｅ^kが得られるごとに２分挿入法（binary insertion）で挿入位置を求めます．これで結果的にリストをソートしたことになるので，検索時間を短縮できるというわけです．
17430 GOSUB *INDEX.LIST:IF SW=1 THEN 17480
18190 *INDEX.LIST:
18200 GOSUB *INDEX.LIST.SUB
18210 L=1:R=MM:SW=0
18220 WHILE L<=R AND SW=0
18230 O=(L+R)\2
18240 IF E0$<EE$(O) THEN R=O-1
18250 IF E0$>EE$(O) THEN L=O+1
18260 IF E0$=EE$(O) THEN SW=1
18270 WEND
18280 IF SW=1 THEN RETURN
18290 '
18300 FOR J=MM TO L STEP -1
18310 EE$(J+1)=EE$(J)
18320 NEXT J
18330 EE$(L)=E0$
18340 RETURN
なお，ここでは平行２ｎ胞体のｍ次元断面を求めているので，原点に対して対称な点も求める頂点となります．検索時間をさらに短縮させるためには，インデックスリストのたとえば"20120112"とその０←→２を反転させた"021021110"は同じ点と見なすことにして，リストには"021021110"の方だけを登録することにします．
細かいことですが，そのためには
ｘkが基底のとき ｅk＝０
ｘkが非基底で下限Ｌのとき ｅk＝−１
ｘkが非基底で上限Ｕのとき ｅk＝＋１
とするよりも
ｘkが基底のとき ｅk＝１
ｘkが非基底で下限Ｌのとき ｅk＝０
ｘkが非基底で上限Ｕのとき ｅk＝２
としたほうが便利です．
18090 *INDEX.LIST.SUB:
18100 E1$="":E2$=""
18110 FOR I=1 TO NV
18120 E1$=E1$+MID$(STR$(E1(I)),2)
18130 E2$=E2$+MID$(STR$(2-E1(I)),2)
18140 NEXT I
18150 IF E1$<E2$ THEN E0$=E1$ ELSE E0$=E2$
18160 PRINT E0$
18170 RETURN
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【３】実行例
頂点数が約５００個のポリトープの輪郭が描かれています．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【４】ムーア・ペンローズの一般逆行列・再考
（その１）で述べたように
Ｆ＝（Ｊ’）~τ
式で与えられるＦをもう一度Ｊ’を用いてτ空間に戻すと，もとの値に戻らないという現象が起こります．
τ→Ｆでは直線が点に退化して射影され，Ｆ→τでは点が直線に射影されます．しかし，線形空間ですから１対１の対応はとれていて，異なるものに射影されるということはないと思われたのですが，実際に起こるのです．
このことより，コラム「ｎ次元楕円をｍ次元空間に投影する」で作成したムーア・ペンローズの一般逆行列を用いたプログラムを実際に動かしてみると，思ったほどには機能せず，アイディア倒れに終わった原因もムーア・ペンローズの一般逆行列に問題があるためだとわかりました．
本シリーズ「平行２ｎ胞体のｍ次元断面」で作成したプログラムはムーア・ペンローズの一般逆行列を用いているのですべて失敗作といわざるを得ません．とはいえ，冗長でない場合は正方行列となるので，真の値を出力してくれるプログラムになっています．それでご勘弁願います．（冗長な場合でも１対１の対応であるため，元の値と何らかの関係（幾何学的な関係）はあるはずで，それをもとにプログラムを修正することが考えられたのですが，この時点ではノーアイディアでした．）
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
また，このプログラムでは
（Ｊ’）~＝（ＪＪ’）^(-1)Ｊ
を用いているのですが，一般逆行列の計算ではしばしばより複雑な
a) X~=X'*X*(X'*X*X'*X)^(-1)*X'
b) X~=X'*(X*X')^(-1)*X*(X'*X)^(-1)*X'
が用いられます．
ａ，ｂいずれの方法を用いるにせよ，計算誤差が累積することが考えられました．このように一般逆行列の計算では桁落ち誤差は避けられませんし，
a) X~=X'*X*(X'*X*X'*X)^(-1)*X'
b) X~=X'*(X*X')^(-1)*X*(X'*X)^(-1)*X'
のどちらかを適切に選ばないと解が発散したりするわけですから，かなり危ない綱渡りをしていることになります．
とくに，ａは行列要素が大きい場合には発散するのですが，佐賀大学の佐々木誠先生のデータのように要素が小さい場合はｂよりもａのほうが計算安定性に優れているようでした．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【補】モーザーのパラドックス
ｎ次元ユークリッド空間において，１辺の長さが１の立方体[-1/2,1/2]^nをｎ次元単位立方体といいます．その体積は１ですが，もっとも離れた２頂点を結ぶ対角線の長さはｎ次元ユークリッド空間の距離の定義から
√（１^2＋１^2＋・・・＋１^2）＝√ｎ
となります．したがって，次元ｎが大きくなると対角線の長さはどんどん大きくなり，身長１７０ｃｍの人間はおろか，ついには地球でさえ含むことができるようになります．
辺の長さが４の正方形に４つの単位円板を詰めると，４つの円板で囲まれた部分に，第５の小さな円を入れることができます．また，辺の長さが４の立方体の８つのカドに単位球を８個詰めると，中にできる隙間に第９の小さな球を入れることができます．ピタゴラスの定理によって第５の円，第９の球の半径はそれぞれ√２−１，√３−１だとわかります．
これと同じことを４次元以上の空間で行うことができます．もはやイメージすることは不可能ですが，１辺の長さが４の４次元超立方体の１６個のカドに１６個の単位球を詰めると，中の隙間には半径√４−１＝１の４次元超球（すなわち単位球）が入ります．同様に，１辺の長さが４のｎ次元超立方体の２n個のカドに単位球を詰めると，中の隙間に半径√ｎ−１のｎ次元超球が詰められるのです．
しかし，ここの驚きが潜んでいます．たとえば，ｎ＝９の場合，中に詰められるｎ次元超球の半径は√９−１＝２であり，この球は外側の立方体の表面に接してしまい，ｎ＞９だとはみ出してしまうのです．この驚くべき結論は，日常生活ではありえないだけに面食らってしまいます．
一方，球に相当するｎ次元の図形を超球と呼びます．ｎ次元単位超球{x1^2+x2^2+･･･+xn^2≦1}の体積をＶnとすると，Ｖ1=2（直径），Ｖ2=π（面積），Ｖ3=4π/3（体積）はご存知でしょう．ｎ次元単位球はどんなに次元が高くても，長さが２より大きな線分を含むことはできません．
したがって，ｎ＝２，３，４では単位立方体（対角線の長さ√ｎ）は単位球体の中に含まれますが，ｎ≧５でははみ出る部分があり，次元とともにはみ出る部分が増えていきます．単位球体の直径は次元によらず２なのです．
次元とともにはみ出る部分が増えているのですが，球の詰め込みに関するこのはみ出し現象は，モーザーのパラドックスとして知られているものです．この逆説は，人間の直観や勘は３次元までの世界では働きますが，４次元以上の高次元についてはあまり働かないという例として，しばしば引き合いに出されます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
