Colloquium

NO.1894　　　 ダイクストラ法と領域認識　　　　2010.12.26.　　ＤＤＴ

皆さん、ダイクストラ法ってご存知ですか？。知りませんよね。知ってたらマニア過ぎます^^。ダイクストラ法とは、道路ネットワークが与えられた時に、ある交差点から別の交差点にいたる最短経路を探索する標準的（古典的）方法です。主に交通計画分野にいる学生か講師陣しか知らないと思います。ダイクストラさんが1959年に提出しました。
道路ネットワークとは、図-1にようなものを指します。図-1は、札幌や京都のような碁盤の目の市街地に対応した、道路ネットワークです。

道路ネットワークでは、考える範囲の全ての交差点と、一つの交差点に隣接する全ての交差点のデータ、および隣り合う交差点間の道路距離が与えられます。
図-1のように規則性のある場合は、例えば点（交差点）Ａから点Ｂにいたる最短経路（赤）は、(北に3ブロック，東に2ブロック)進む任意の経路とすぐわかりますが、ネットワークに規則性がないと、そうは行きません。いきおいトライandエラーの探索という事になります。
すぐ思いつくのは、道路ネットワークの全ての可能な通り方を多分木で全探索して経路長順に並べ、ＡからＢにいたる経路だけを抽出し最短を選ぶというやり方です。一つの交差点に隣接する全ての交差点のデータと、隣り合う交差点間の道路距離が与えられるので、これは原理的には可能です。ところが少し概算すると、このやり方は現実的に実行不可能なのがわかります。
多分木全探索を図-1に試すとします。交差点は概ね3分木になっています。すなわち、ある道路から交差点に進入した場合、3通りの分かれ道というわけです。3通りの内の一つを選ぶと、また3通りの分かれ道が現れます。もし交差点が100個あるとすれば（縦横に道路10本）、超概算で3×3×・・・×3＝3¹⁰⁰の探索数になります。 3¹⁰⁰は、ほぼ10⁵⁰のオーダーです。前もやりましたが、一回の探索を1GHzのCPU一回の動作で完了できるとしても、10³⁰年ほどかかる計算になります。もちろん10⁵⁰個の中には、同じブロックをぐるぐる巡回したり、Ａから出発しなかったり、Ｂに行かなかったり、通過経路に戻ったりするのもあるので、制約条件をつけて大幅に数を減らす事は可能です。例えば、Ａから出発するものだけに限れば、1/100になるのは明らかです。
それでもその結果、例え探索数が100億分の1になったとしても焼け石に水です。10²⁰年です。図-1に示したように、わざわざ遠回りする経路（青）も全て含まれるからです。交差点100個（縦横に道路10本程度）なんて、ちょっと大きな町なら持ってるはずです。
という訳でダイクストラさんは、最短経路探索を効率化するダイクストラ法を考えました。ダイクストラ法は帰納法で証明できますが、探索効率を上げるための手順が追加されていて、わかりにくいものになっています。そこで、ダイクストラ法ほど効率は良くない、次の基本ダイクストラ（仮称です）を考えます。

１．基本ダイクストラの証明
最初に出発点Ａと目的点Ｂを決め、Ａから行ける隣接点C₁⁽¹⁾，C₂⁽¹⁾，･･･，C_a⁽¹⁾への経路長を計算します。それらをL₁⁽¹⁾，L₂⁽¹⁾，･･･，L_a⁽¹⁾とすれば、これはＡで分岐する道路の道路長そのものです。また任意のC_i⁽¹⁾が、Ａからの最短経路上にあるのは明らかです。
　経路長L_i⁽¹⁾を、経路長リストＬ⁽¹⁾＝{Ｌ_i⁽¹⁾}として記録します。
　経路P_i⁽¹⁾＝Ａ→C_i⁽¹⁾を、 経路リストＰ⁽¹⁾＝{P_i⁽¹⁾}として記録します。
次に、各C_i⁽¹⁾から行けるC_j⁽²⁾を全て集めて C₁⁽²⁾，C₂⁽²⁾，･･･， C_b⁽²⁾とし、それらへの経路長をL_j⁽²⁾とします。 L_j⁽²⁾は、L_i⁽¹⁾にそれぞれの道路長を足すだけです。 経路P_j⁽²⁾＝Ａ→C_i⁽¹⁾→C_j⁽²⁾も作り、 それぞれリストに加えます。

　　Ｌ⁽²⁾＝{L_i⁽¹⁾，L_i⁽²⁾}
　　Ｐ⁽²⁾＝{P_i⁽¹⁾，P_i⁽²⁾}

上記のＬ，Ｐは省略して書いてますが、 本当はＬ⁽²⁾＝{L_i⁽¹⁾}∪{L_i⁽²⁾}などの意味です。
{C_i⁽²⁾}には、C_i⁽¹⁾へ逆行するものも含まれますし、 逆行しなくても他のC_j⁽²⁾に一致するものも出る可能性もあります。 ありますが、とりあえず気にしません。 重要なのは{C_i⁽²⁾}が、2回のステップでＡから行ける全ての点を網羅している事です。
Ｐ⁽²⁾に含まれる一つの経路P_i^(k)に注目します（k＝1，2）。 P_i^(k)はC_i^(k)への、最短経路とは限りません。何故なら、k回以下のステップmでC_i^(k)達するP_j^(m)があって、 L_i^(k)よりL_j^(m)の方が短いかもしれないからです。
また(k＋1)回以上のステップnでC_i^(k)達する、より短いP_j⁽ⁿ⁾が、この先出るかもしれないからです。しかしＬ^(k)の最小値L_p^(k)に対応する、経路P_p^(k)とその端点C_p^(k)については違います。
L_p^(k)はk回以下のステップで生成される経路長の最小なので、k回以下のステップに限れば、経路P_p^(k)は、端点C_p^(k)への最短経路です。 それより小さな経路長が、k回以下のステップにはないからです。 さらに、(k＋1)回以上のステップでC_p^(k)に達する経路長は、 全てL_p^(k)より大きい事が言えます。 (k＋1)回以上のステップで到達できる点は全て、k回以下で行ける点を通過しなければならないからです。 ここで、{C_i^(k)}が、k回のステップでＡから行ける全ての点を網羅している事と、L_p^(k)がＬ^(k)の最小値である事が効きます。 以上のステップを、Ｂに到達するまで繰り返します。 いつか必ず達するのは明らかです。
ステップの繰り返しの終了条件は、 経路Pp(k)の端点C_p^(k)が、C_p^(k)＝Ｂとなる事です。

もう帰納法に持ち込む準備は整いましたよね？（整いました！・・・何か某謎掛け芸人みたいですね^^）。なので証明抜きで、次の基本定理（だと思います）をあげます。

［ダイクストラ法の基本定理］
　Ｌ^(k)＝{L_i⁽¹⁾，L_i⁽²⁾，･･･，L_i^(k)}の最小値L_p^(k)に対応する経路P_p^(k)は、 その端点C_p^(k)への最短経路。ただし(k-1)以下のステップで最小値となったものは、最小値評価から除外する。

［証明］
　帰納法による。k＝1の時は、Ｐ⁽¹⁾＝{P_i⁽¹⁾}が全て最短経路なので、帰納法は成立する。
［証明終］

この基本定理さえやってしまえば、後は探索をもっと効率化するだけです。基本ダイクストラの不経済な点は以下です。
　　(B-1)逆行する。
　　(B-2)ループする（以前のステップで探索した点に戻る）。
　　(B-3)既に最短経路がわかった点も探索・評価する。
　　(B-4)より短い経路があるのに、そちらに探索を絞らない。
　　(B-5)一回のステップで、不要に遠い点まで探索する。

(B-1)〜(B-5)を解消すると、次のわかりにくいダイクストラ法になります。

[０]初期状態
　　(1)出発点Ａと目的点Ｂを決める。
　　(2)Ａから到達できる隣接点Ｃ₁，Ｃ₂，･･･，Ｃ_sへの経路長と経路をリストする。
　　(3)リストの最小値となる点Ｃ_mを見出す。経路Ａ→Ｃ_mが、Ｃ_mへの最短経路なのは明らか。
　　(4)Ｃ_mを探索点から除外する。
　　　（最短経路のわかった点は探索しない）
　　(5)Ｃ_mを経路長と経路リストから除外し、別途用意した最短経路リストに、Ａ→Ｃ_mを記録する。
　　　（最短経路のわかった点は評価しない）

[１]以後の手順
　　(1)出発点をＣ_mとする。
　　　（一回のステップで、不要に遠い点まで探索しない）
　　(2)Ｃ_mから到達できる隣接点Ｄ₁，Ｄ₂，･･･，Ｄ_tへの、Ａからの経路長と、Ａからの経路をリストする。
　　　（差分を足すだけ）
　　　　ただしＤ₁，Ｄ₂，･･･，Ｄ_tの中に、経路Ａ→･･･→Ｃ_mの点は含まない。
　　　（逆行とループ禁止）
　　(3)Ｄ₁，Ｄ₂，･･･，Ｄ_tの中に、Ｃ_iと一致するものがあれば、既存の経路長と現在の経路長を比較し、短かいほうで、Ａからの経路長リストとＡからの経路リストを置き換える。そうでないものは、たんにリストに追加する。
　　　（より短い経路を優先）
　　(4)Ｃ₁，Ｃ₂，･･･，Ｃ_s，Ｄ₁，Ｄ₂，･･･，Ｄ_uの経路長の中から、最小のＤ_mを選ぶ。
　　　（基本定理）
　　(5)Ｄ_mを探索点から除外する。
　　　（最短経路のわかった点は探索しない）
　　(6)Ｄ_mを経路長と経路リストから除外し、最短経路リストに、Ａ→･･･→Ｃ_m→Ｄ_mを記録する。
　　　（最短経路のわかった点は評価しない）

[2]終了条件
　　(1)Ｄ_mがＢなら終了する（目的点に達した）。
　　(2)探索点がないなら終了する（Ｂが最後に探索された）。
　　(3)(1)でも(2)でもないなら、Ｃ_m＝Ｄ_mとして[1]-(1)に戻る。
　基本ダイクストラは、一回の探索ステップごとにＢへの最短経路か否かを判定できる点を除いて、じつは多分木全探索と同じです。しかし判定アルゴリズムがあるので、多分木探査の最小ステップでＢへの最短経路が得られ、そこで計算を打ち切れます。これが、ダイクストラ法が効率的である基本的な理由です。 　(B-1)，(B-2)の解消は、最短経路探索ならば必ず満たすべき制約なので、これらで本質的に速くなる訳ではありませんが、多少速くなるのは確実です。 　(B-3)の解消は、ダイクストラ法を本質的に速くします。(B-3)を解消すれば、一回の探索ステップごとに探索点が一個づつ減っていくので、必要探索量の上限が1/nづつ複利で減る事態を招きます。ここでnは、各点の分岐数の平均程度の整数です。多分木探索で探索数が、n^m（m：探索点数）の複利で増えるのと逆の状況です。ダイクストラ法は、後半に行くほどスピードアップします。 　(B-4)，(B-5)の解消も。セットで本質的な加速です。アルゴリズムとしては、[0]-(3)または[2]-(3)から[1]-(1)への移行と、[1]-(3)のリスト更新として実装されます。 　[1]で選択されたＤm点への経路は常に、現ステップでの「最最短経路」です。そこに次の探索ステップを絞ります。そうするとダイクストラ法では、k回目のステップで、k回のステップで到達できる点を「全て網羅はしない」事になります。ただし、k回以下のステップでは、必ず網羅されているステップが存在します。例えばk＝1は、明らかに網羅です。 　現ステップでの「最最短経路」を選ぶので、「最最短経路」は、網羅されたステップの点の側から決まっていきます。結果として、Ａ点からの近さの順に最短経路が決まります。これによって不要に遠い点まで探索する事が防がれます。ＢがＡから10番目に近い点なら、10ステップで到達できる訳です。 　このようなダイクストラ法の計算量はm²に比例します（証明できません）。ですが、個人的に比較できるアルゴリズムとしては、m元連立一次方程式を解く、ガウスの消去法があります。これの計算量は、m³に比例です。現行の標準市販PCのCPU、1GHz×quad Coreなら、m＝10000であっても1000本の最短経路決定に、長くて数秒と予想できます。 　ダイクストラ法は、「言われてみれば当たり前」の事をやっています。しかしそれを、ここまで整理して定式化したダイクストラさんは、やっぱりすごいと思います。

２．有限要素法
自分の本当の専門は、構造屋で橋梁屋です。それがどうしてダイクストラ法なんかを持ち出したかと言うと、久しぶりに有限要素法解析を頼まれたからです。指定された解析ソフトは馴染み深いもので、最初はそんなに気にしてませんでした。ところがその解析ソフトを久しぶりに使用したところ、恐ろしく不便なのに気づきました（切れかけました）。昔は不便さに、気づいていなかったのです。という訳で、有限要素法を概説させて頂きます（←全く理由になってない^^;）。

有限要素法が適用される典型的な状況は、例えば正方形の鉄の塊があって、それが自重でどのように変形するのかを計算したい場合です。図-2の碁盤の目で区切られた正方形が、今回の鉄の塊です。碁盤の目の事をメッシュと言い、目の一つ一つを要素と言います。
有限要素法の前提は、じつは一般的な物理的仮定と数学的原理にあります。いま正方形の変形状況を知りたいので、図-2にx-y座標を入れ、正方形の各点のx方向，y方向の変位dx，dyが、(x，y)の関数としてわかればOKです。

　　　

つまりd＝(d_x，d_y)は、変形前と変形後のある点の座標の差d＝(x＋dx，y＋dy)−(x，y)です。
有限要素法の最初の前提は、次の経験事実と数学的原理です。

　　(1)物理量は、みな微分可能である。
　　(2)微分可能なら、局所的に線形近似（比例配分）可能。これは微分の意味そのもの。

(1)より(d_x，d_y)は、微分可能です。そして「局所的に」なので、図-2の要素が十分小さいとして、(dx，dy)を要素上で考えると結論として、要素の変形に対して次の歪み分解が得られます。

局所的には要素の変形は、x，y方向の伸びに対応する変位、せん断変形に対応する変位、回転に対応する変位の、たんなる足し算になります。足し算なるのは、(d_x，d_y)が微分可能で局所的に線形近似可能だからです。
材料は変形に対して抵抗します。例えばx方向に伸びたら、ある力でx方向に縮もうとします（バネを想像して下さい）。通常の有限要素法では、ここで2つ目の物理的仮定の登場です。

　　(3)変形は、微小である。

これは大抵成り立ちます。例えば何10万tというタンカーでも、自重で微小にしか変形せず自立できるので、建造中に乾ドックに置いておけます。
変形に応答する抵抗も物理量です。(1)より微分できます。しかも変形は微小です。(2)を言い換えると「微分可能なら、微小入力に対する応答は、入力に比例する」です。つまり歪み分解の各変形モードに対して、本当にバネを想定できるという事です。歪み分解の変形モードのうち、要素を本当に変形させているのは、最初の3つです。それで要素には、次のようなバネが内蔵されていると考えます。

歪み分解の図からわかるように、これらのバネの伸び縮みは、要素の四つ角の点の変位が支配します。バネ定数は、鉄なら鉄の変形にふさわしいように材質に合わせて調整し、実際の未知数は、要素の四つ角の変位という離散化量です。
以上の手順で、図-2の碁盤の目の目を全てバネモデルでおきかえ、物体への外力を与えると、バネを縮ませる（伸ばす）四つ角の変位に関する連立一次方程式が得られます。それを解いて全ての要素の四つ角変位を得たら、再び要素が小さい事を利用して（(1)，(2)を利用して）、四つ角変位の値から要素内部の点の変位を、線形近似で求めます。これで終了です。
有限要素法は(1)，(2)，(3)という一般的原理に基づいているので、ほとんどの場合で現象に相当するバネモデルが存在し、ほとんどあらゆる分野に適用可能な近似計算法です（宇宙論にも）。とは言え「メッシュを切らなきゃ始まらない」・・・事は、わかって頂けると思います。
物体などの形状（解析領域）が単純なら、図-2のような碁盤の目メッシュで済みますが、形状が複雑だとそうはいきません。じっさい有限要素法で一番手がかかるのは、実はこの「メッシュ切り作業」なのです。

３．領域認識

図-3は、とある山の断面です。山の形はラジコン機による航空測量の結果と思われます。図中のＡ，Ｂ，Ｃの領域では、土や岩の固さが違います。Ａ，Ｂ，Ｃの境のラインは、ボーリング調査をもとに、地質屋さんが予想したものです（固さも）。真ん中辺りにある半円状の領域は、この山に掘削予定のトンネル範囲です。
トンネルを掘削したら、山の応力状態（内部の力）がどうなるか？を確認するために、有限要素法を行います。つまり落盤などが心配なのです。応力状態は変形量からバネモデルのバネの伸び縮みを計算し、伸び縮み量にバネ定数をかければ得られるので、結局やる事はさっきといっしょです。
ただし今度は、解析領域の形がやや複雑な事、材料定数（バネ定数）が領域ごとに異なる事から、最低限Ａ，Ｂ，Ｃの領域を、解析ソフトに入力する必要があります。トンネル範囲がたくさんの領域に分かれているのは、地盤改良などの工法を併用し、地山の性質（バネ定数）が変わってくるからで、これらも入力する必要があります。その他のラインは、綺麗なメッシュを得るための補助ラインですが、それらによってさらに領域が増えます。

解析ソフトに領域を認識させた後は、領域ごとに、通常はAdvancing Front法を用いて半自動のメッシュ生成を行います。得られたメッシュは、図-4ですが、ここまで来るのが大変でした。
図-3の丸は幾何学的な点を表し、図-3の形状は点をつなぐラインの集積です。領域を定義するとは、そのラインの集積に位相情報（点のつながり方）を与えるという事です。基本は、領域境界の全ての点を、境界に沿う順番で手入力です。そうすると例えば、右側のＢの領域の認識などは、けっこうな面倒地道な作業になります。しかも使用した馴染みのソフトは、点の認識が甘いという欠点を持ち、領域に目に見えない隙間ができる事も、しばしばでした。その上、馴染みのこの古いソフトは、既にWindowsの古典標準と化している「Undo機能」を持たない始末です。つまり隙間をなくする作業は、全て「最初からやり直し」なのです。
昔の領域認識作業は、紙の図面上でデシタイザーという機械を用い、虫眼鏡のターゲット内に一点一点を捉えて位置を記録し、位置の数値データをPCに食わせるという、今思えば、とても人間がやるような作業ではありませんでした。それを画面上のクリック一発で処理する、このソフトが現れた時には画期的で、不便さに気づいていなかったのです。

で、領域の自動認識はできないか？、と考えました。ダイクストラ法の出番です。図-5に示すように、点Ａから出発して、点Ａに戻るような最短経路を探せばいいじゃないか？、という発想です。それはダイクストラ法で、出発点と目的点を一致させれば可能です。
ただこの方法には、2つの欠点があるのは、最初からわかっていました。一つは島が残る事です。
各領域の分岐点（三叉路以上になっている交差点）から出発点に戻る探索を順次行うとすると、外周が短い領域から順番に認識されます。またある分岐点の分岐辺が属する領域が確定したなら、その点から出発する探索においては、確定辺は二度と探索するべきでありません。同じ領域が二重に生じるからです。
そうすると領域の大きさによっては、その周囲の領域が全て先に認識され、図-6の「島」と図示した領域については、探索が行われない事になります。
これへの解決策は、たんに無視すればOKだと気づきました。初期の状態でやれるところまでダイクストラをやった後で、2つの領域に属するラインと点を全て削除すれば、孤立した島だけが残ります。孤立した島では探索不要なので、それらの認識は一瞬です。
もう一つの欠点は、複合領域が出来うる事です。出発点に戻る最短経路が、必ずしも単一領域に対応しない事です。図-7はその例です。赤のルートは青のルートより短い、という結果になります。それで初回は赤のルートが認識され、複合領域になります。

可能性としては、赤と青が正確に等長という場合もあり得ます。その場合は、ダイクストラ法の最小値判定で、偶然に赤が選ばれた訳です。
複合領域になったら、後述の繰り返しに持ち込んで処理できるので、まずは無視してOKとわかりました。ただし複合領域判定は必要です。それには生成された位相情報を利用できます。
単一でも複合でも、とにかく領域として認識されたら、領域境界上の点の順序集合が得られます（右周りか左周り）。その中から三叉路以上になっている点を選び、分岐が境界上にない事を確認して、分岐に沿って1ライン進みます。図-7のC1とC2です。全ての三叉路以上の交差点で、C1のような内部点がなければ、単一領域です。
点の領域内部／外部／境界上の判定にも、位相情報を利用できます。点と境界点を結ぶ直線の角度を、順序集合に従って右回りか左回りかで順次計算し、角度増分を積算すれば結果は、外部：0，内部：±2π，境界：±πです。これは、順序集合に従って自分が境界点を目視していった場合、自分の首がどのように回るか想像すれば、すぐわかります。0，π，2πの分離は非常に明確なので、コンピューターの数値誤差を考慮しても、まず間違う事はありません。
以上の処理を完了すれば、残るのは孤立した単一島か複合島です。複合島があれば単一島は全て無視し、複合島の内部探索を行います。

複合島の外周は、その外側の（消してしまった）領域の確定に伴って、確定辺になっています。一方内部の辺は、まわりの領域が確定していないので、全て未確定辺です。よって今までのアルゴリズムをたんに繰り返し適用するだけで、図-8の赤で示したような探索が自動で始まり、複合領域は単一領域にいつかは分解されて最後に残るのは、もともとの孤立島と、複合島から発生した孤立島のみです（図-9）。
図-9の状態まで来たら、後は孤立島を全て認識させて終了・・・と思ったのですが、最後の問題がありました。

島の定義はこうです。
　　「境界上の全ての点が、二分木である事」。
これは全ての点が交差点でない事を言っていますが、この条件に、解析領域の最外周が該当してしまいます。
解析領域の最外周を領域として認識させると、内部の全ての単一領域と重複が生じて、領域ごとにメッシュを切るという作業ができなくなります。
現在は削除したラインと点を復旧した後、全ての点を含む領域は不可、という方法で対処していますが、もっと完全な方法が必要だと、その後思いました。
というのは、ここまでの作業を実行するプログラムは、既にApplicationの態を成しています。このまでで終わらせるのは勿体ないので、領域ごとの完全自動メッシュまでやろうと思っています。その価値は、図-3にようにやたらと領域を増やさずに、必要最小限の内部領域だけで解析領域を定義し、何の手もかけずに自動メッシュできる事です。
図-3のような「図面」は、いくら製図ソフトが発達したとは言え、結局は手で描きます。また領域入力後に図-4のメッシュを発生させるため、領域境界の分割数を手動で調整するという、もう一山作業もありました。
完全自動の領域認識−メッシュ生成が遭遇する典型的状況は、図-10のようなものだと思います。つまり孤立領域が基本です。という事は、今まで外周しか持たない単一領域で考えて来ましたが、内周を持つ領域への対処が是非必要です。孤立が基本だからです。それには、それぞれを外周しか持たない単一領域として認識させ、重複を起こさせないために領域の包含関係を整理して、重複を持つ領域に内周を追加する必要があります。つまり領域の包含順序の階層化です。位相情報を利用して、領域間の含む／含まない関係はすぐ判定できますが、単純にやると多分木を使った経路探索のような目に合うような気がして、現在は方針を検討中です。

とは言え、ここまでの結果にはいちおう満足しています。最後に領域認識過程の図をあげます。

NO.1893　　　　　ｅ　と π について　　 　　　　2010.12.25.　　水の流れ

第２５１回数学的な応募問題

先日同僚との会合で次のような問題を頂きました。 その場ですぐに発想が出てきませんでしたが、・・・

ｅ（超越数）とπ（円周率）について

　　ｅ^πとπ^ｅ　の大小を調べよ。

＜水の流れ：ある関数のグラフを利用して解きましたが・・・他の解法があれば・・・＞

注：この記事に関する投稿の掲載は、2011年１月10日以降とします。

戻る