スポンサーサイト

--年--月--日 --:--

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Occlusion Query

2008年12月25日 00:29

リアルタイム 3D において Occlusion Culling は長年の課題でした。
実際に描画されるかどうかは Depth Test や Alpha Test が行われた後にしかわからないわけですから、Depth Buffer を参照する術を持たない(あっても物凄く非効率だった)Shader 技術が普及する以前の GPU では現実的な話ではありませんでした。

もっとも、頂点シェーダーが重い!というのは最近ではあまり考えられなくなりましたが、それでも描画がされないのならモーフィング処理などを行わなくても良くなるので、効果がないわけではありません。
さらに、最近の GPU は unified shader に対応しているのが当たり前になってるので、頂点シェーダーを通さなければそれだけ Pixel Shader に時間を費やせます。
また、オブジェクト描画によるフラグメントの更新が行われないのがわかるのなら、そのオブジェクトが使用する Texture を転送しなくてもいいワケなので、バス幅の節約にもなります。

さて、こんな Occlusion Culling、フラグメントまでいったかどうかの判定ですむわけですから、最近のハードウェアなら簡単に実装できそうな気もします。

実は、OpenGL では 1.5 (2003年7月にリリース)から採用されている機能で、DirectX には 10 から Occlusion Predicate Query という名前で実装されています。
二つとも、Bounding Sphere や Axis Aligned Bounding Box、Oriented Bounding Box などの簡易オブジェクトによる本描画前の遮断問い合わせを実現するための機能です。
従来は難しかったフラグメント単位での問い合わせが出来る点でなかなかに強力です。
OpenGL では GLSL の本格導入前(GLSL の本格導入は 2.0 から)から実装されているのが面白いですね。


そんな Occlusion Culling、GL での実装は簡単です。

1:まずオブジェクト ID を生成する(別になくても OK)
2:glBeginQuerry で遮断問い合わせ開始を宣言
3:実際に簡易的なオブジェクトの描画を行う
4:glEndQuerry で遮断問い合わせ完了を宣言
5:Depth Test をパスしたサンプル数を取得する
6:生成したオブジェクト ID を削除する。

です。


ちょっとここで注意。
windows で OpenGL をプログラムする場合、以下の点に注意しないといけません。

OpenGL を管轄する khronos グループでは OpenGL SDK という名前の SDK は配布していません!
OpenGL 1.1 までなら Windows Platform SDK に付属している OpenGL32.dll で使用できます。
  
OpenGL 1.2 以上 (1.5 や 2.0 など)の機能は Windows 付属のライブラリ(OpenGL32.dll)では扱えません。
使う場合は、OpenGL の HP から glext.h を落としてきてインクルードし、wglGetProcAddress 関数などで取得します。
  
 例:
  glGenFramebuffersEXT = (PFNGLGENFRAMEBUFFERSEXTPROC)wglGetProcAddress( "glGenFramebuffersEXT" );

また、glew というライブラリを使うことで、wglGetProcAddress を使用しなくても、いきなり glGenFramebuffersEXT(...) という感じで意識せずに使うことができます(オススメ)



wglGetProcAddress、glewInit などを使用する場合、必ず一つはコンテキストがないと失敗します。
コンテキストは wglCreateContext 関数で作成できます。



特に 1.2 以降にも対応している最新の OpenGL32.dll を探そうとするとハマるので注意してください。。
そんな gl.h や OpenGL32.dll は存在しません。
1.2 以降の関数には wglGetProcAddress や glewInit を通じてのみ、アクセス可能です。



閑話休題。
では、遮断問い合わせ処理に必要な関数をそれぞれ上げてみましょう。


void glGenQueries( GLsizei n, GLuint* ids )

ids に n 個の問い合わせオブジェクトの識別 ID を返します。




void glBeginQuery( GLenum target, GLuint id )

クエリー用の描画開始を宣言します。
target には GL_SAMPLES_PASSED だけを指定できます。
id には glGenQueries で生成した ID を渡してください。
フラグメント処理が通ったかどうかの判定だけなのですが、レンダーステートは反映されます。
ですので、レンダリング時間を増やすようなモードは無効にしたほうがいいでしょう。
(例えば Stencil 処理など)




void glEndQuery( GLenum target )

クエリー用の描画終了を宣言します。
target には同じく GL_SAMPLES_PASSED だけ指定可能です。




void glGetQueryObjectiv( GLuint id, GLenum pname, GLint* params )
void glGetQueryObjectuiv( GLuint id, GLenum pname, GLint* params )

二つとも、まったく同じ内容の関数です。
指定したクエリー id のフラグメント処理が行われたかどうかを問い合わせます。
id には glGenQueries で生成し、glBeginQuery で描画が行われた ID を渡してください。
pname には GL_QUERRY_RESULT、GL_QUERY_RESULT_AVAILABLE が指定できます。
GL_QUERRY_RESULT の場合は、params にはパスしたフラグメントかサンプル(MRT が有効な場合)数が帰ります。
これで params が 0 だった場合は、オブジェクトが全て遮断されていることになります。
GL_QUERY_RESULT_AVAILABLE を指定した場合、params には GL_TRUE か GL_FALSE が入ります。
GL_TRUE の場合は、渡した id での問い合わせが可能なことを意味し、GL_FALSE の場合は不可能なことを意味します。




void glDeleteQueries( GLsizei n, const GLuint* ids )

指定した n 個の ids を削除します。
ids には glGenQueries で生成した ID を渡してください。




以上の関数を使うことで occlusion culling が出来ます。
今回は GL の説明しかしませんでしたが、DirectX の場合も同じような機能ですので、そちらを参考にして下さい。

関数型プログラム

2008年11月19日 23:53

関数型プログラムを実験してみました。

実装した内容は、カリー化なんですが、思いのほか面倒でした。。。
けど、boost の bind2nd とかよりは格段に使いやすいですよ。
なかなか面白いものができたので公開してみます。

#include
#include "macro.h"
#include "typelist.h"
#include "service.h"
#include "functor.h"


double FunctorTest( float arg1, int arg2, char *arg3 )
{
  printf( "%-8s : (%.2f + 1) * %2d = ", arg3, arg1, arg2 );
  return double((arg1+1) * arg2);
}

template< typename Func >
double FunctorCallTest( Func f )
{  return f(5.0f, "call");  }

void main()
{
  Functor functorFactory;
  char str[] = "hoge";

  // bind テスト
  printf( "%6.2f\n", functorFactory(FunctorTest, Binder::unbind, 10, str )( 1.23f ) );
  printf( "%6.2f\n", functorFactory(FunctorTest, 2.46f, Binder::unbind, str )( 5 ) );
  printf( "%6.2f\n", functorFactory(FunctorTest, 3.69f, 3, Binder::unbind)( "hogehoge" ) );

  // bind することで違う関数を生成しているのを試すテスト
  printf( "%6.2f\n", FunctorCallTest( functorFactory(FunctorTest, Binder::unbind, 20, Binder::unbind) ) );
}



こんなプログラムを走らせると、、、

hoge : (1.23 + 1) * 10 = 22.30
hoge : (2.46 + 1) * 5 = 17.30
hogehoge : (3.69 + 1) * 3 = 14.07
call : (5.00 + 1) * 20 = 120.00



こんな結果が返ってきます。

バインドしたくない部分は Binder::unbind で置き換えてるわけですね。



ソースは以下。

ヘッダー
macro_h.txt
service_h.txt
typelist_h.txt
functor_h.txt

コード
sample_cpp.txt

macroとservice、typelist はTypeTraisとかTuple、Typelistを定義しているだけです。
Functor が今回最大の目玉。
特にカリー化している最中の処理は面白いです。
Tupleに収まっている引数をマージして新しいTupleを作成、それを引数として関数を呼んでいます。
かなり難解なので、良かったら読みほどいてくださいw

ちなみに、汎用化のためにマクロをこれでもか!というぐらい使用していますので、物凄く見づらいソースになっていますw

PSM(Perspective Shadow Maps)とLiSPSM(Light Space Perspective Shadow Maps)のアルゴリズム

2008年11月04日 00:21

PSM、LSPSM のアルゴリズムを自分なりに解釈しまとめてみたのでアップしてみます。
参考になれば幸いです。

尚、絵は後で追加します。多分。

PSM のアルゴリズム



ワールド座標上のディレクショナルライトの方向をカメラローカルの座標に変換する
この表示用のローカル座標系を保持しているカメラを表示カメラとする。




1で変換した表示カメラローカル座標でのディレクショナルライトを同次座標上で歪ませる。
ライトは、無限遠平面上に位置させます(通常は 1 の同次座標値を 0 と置くことで表現できる)
このとき、ライトを実際に「置く」のがポイント。
ディレクショナルライトを影響距離が無制限の点光源を無限遠平面に置いたと見なしている。
また、行列変換はあくまで遮断を行うための行列です。カメラ位置を「下げる」平行移動は行いません。
(透視投影行列の [4][x] 成分はカメラ位置が原点なので近平面を原点にするための移動にあたります。これを無視するということです。また、[3][4] 成分は透視投影による歪み(fovy による歪み)に当たります。 この[3][4]成分が重要で、これにより透視投影が実現されています。詳しくは依然書いた透視投影の記事を参考にしてください)




同次座標上で歪んだ2のベクトルは歪んだことにより、 (視線に対して垂直の場合以外は)無限遠平面から外れるので、3 次元座標上の点になる。
もし、視線に対して垂直の場合(要するに光源位置は相変わらず無限平面にある場合)はライトの位置とカメラの位置での歪みは発生しないことになるので通常のシャドウマップと同じ意味になる。
よって、その場合は通常のシャドウマップを適応する。




3で求めた光源位置から透視投影変換後の空間(表示空間と仮に呼びます)の中心(0,0,0.5)を見るようなカメラを作成する。
これを通常のカメラと区別するためライトカメラと呼ぶ。
表示空間は -1 <= x,y <= 1、 0 <= z <= 1 になるので、表示空間の中心は、(0, 0, 0.5) ということになる。
このライトカメラから表示空間を投影することで、表示空間が歪む。




4で求めたライトカメラからの映像を平面に投影する。
このとき、表示カメラの近平面が画面いっぱいに上手く収まるように
視野角や近平面、遠平面を調整する。
  
具体的には、、、
  
視野角は表示空間の中心から一番遠い頂点までの距離を利用し逆算します。
  rad : 表示空間の中心から一番遠い(表示空間を構成する)頂点までの距離
  dist : ライトカメラから表示空間の中心(注視点)までの距離
とすると、
  fovy = 2 * atan( rad / dist )
になります。
  
近平面までの距離は、中心から z 軸のマイナス方向に一番遠い点(カメラに一番近い点)を取ります。
よって、
  pt_near : ライトカメラローカル座標上で一番 z 軸が小さい頂点の z 値
  center : ライトカメラローカル座標上での表示空間の中心の z 値
  dist : ライトカメラから近平面の中心(注視点)までの距離
とすると、
  plane_near = dist + (pt_near - center)
になる。
  
同様に遠平面までの距離は
  rad : ライトカメラローカル座標上で一番 z 軸が大きい頂点の z 値
  center : ライトカメラローカル座標上での近平面の中心の z 値
  dist : ライトカメラから近平面の中心(注視点)までの距離
とすると、
  plane_far = dist + (pt_far - center)
になる。

以上、これらの値を使って透視投影行列を作成する。




4と5で出来た行列は表示カメラにより投影された後の処理なので、 通常の描画処理の後に4→5と適応する。
これで PSM 用のシャドウマップが出来る。



補足


3の時点で、同次座標系の値が - の場合、無限遠平面にあったカメラ位置が
カメラの近平面より後ろ側に透視投影されることを意味します。
これだとライトカメラからみた近平面が裏になるので通常の透視投影ではカリングされ表示されません。
よって、この場合は、5で使用するライトカメラの投影行列の近平面を - にした値で対処します。
要するに、さらに裏返しにするわけです。
こうすることで、カリングにより除外されることがなくなります。





LiSPSM のアルゴリズム

PSM の欠点を補うために新しく考案されたアルゴリズムです。
PSM のアルゴリズムで気付かれたかもしれませんが、シャドウマップ関連の技法というのは
ある3次元の点を2次元上の決まった位置に拘束する関係式さえあれば、どんな関係式でもシャドウマップを生成することが出来、
また、それを参照することで影を落とせます。

かなり強力な技法ですね。
シャドウマップ技法がメインストリームになったのも頷けます。

LiSPSM では PSM の技法に一捻り加えます。



表示カメラの視線(-Z方向)ベクトル(表示視線ベクトル)とライトカメラの視線ベクトル(ライトベクトル)の外積を求めます。
両視線が形成する平面に対して垂直なベクトルを求めるわけです。
この時点で、表示視線ベクトルとライトベクトルの成す角度が 0 or 180 (すなわち、外積が縮退してしまう)ならライト位置による歪みがないので通常のシャドウマップ処理にします。

  
  

1で求めたベクトルとライトベクトルの外積を求めます。
両視線ベクトルが形成する平面上且つライトベクトルに垂直なベクトルを求めているわけです。

  
  

ライトベクトル、1のベクトル、2のベクトルで新しい正規直交座標系が出来ました。
この座標系を involve カメラ座標系とします。
(名称は適当に決めています)
このとき、
  ライトベクトル - x
  1 のベクトル - y
  2 のベクトル - z
にあたります。
影を落とす全てのオブジェクト(Caster)と表示カメラの視錐台をこの座標系で変換します。




involve カメラ座標系で全てを内包する Axis Align Bounding Box を作成します。
この AABB の z の小さい面が近平面の距離になります。

  


表示視線カメラの近平面までのベクトルを2のベクトルに射影し、値を計算します。
  
単純に
  Near = dot( 2のベクトル, 表示視線ベクトル) * 表示カメラの近平面までの距離
          … dot は内積を求める関数
でOKです。
  
ここまでの計算は表示カメラローカル上で行っているので、結果的に
  Near = 2のベクトルの z 値 * 表示カメラの近平面までの距離
になります。




4で求めた近平面(AABB の z が一番小さい面)の中心から、5で求めた距離を引いてやった位置が
involve カメラの位置になります。
ちなみに、5で求めた値は大まかな位置で、もっと大きい値にしてもいいですし、正の範囲ならもっと小さい値にしても問題ありません。
しかし、あまりにも小さい値にしすぎると、involve カメラの遠平面に近い点が、2次元に投影した時、極端に中心に寄ってしまいます。
少し考えればわかると思いますが、小さい値にする=視野角が広がることを意味するので、パースが極端についてしまうわけです。
逆に極端に大きくしてしまうと、今度はパースがほとんどつかなくなってしまうので、表示カメラの視錐台の範囲外の部分が多く描画され
テクスチャ領域を無駄にしてしまいます。
5の値を目安に、いい感じの値を入れてみてください。

  
  

カメラ位置が求まり、座標系も求まりました。
あとは透視投影変換用の行列を作ることが出来れば、2 次元に落とすことが出来ます。
  
透視投影行列で使われる近平面までの距離は、6で求めた値を使用します。
遠平面までの距離は AABB の各頂点の z が一番大きい値になります。
後は、近平面の範囲ですが、、、
まず、5で求めた各頂点と視点を結ぶベクトルを求め、これを involve カメラ座標系の z = 1 の値に写像します。
  
  vec = involve カメラ座標系での Caster ポジション - 5のポジション
  vec /= vec.z
  
ということですね。
これで、vec.z = 1 のときのベクトルが求まります。
  
このベクトルを近平面まで伸ばした時の最大 x と最大 y が近平面での高さになります。
(実際は、中心から x までの距離なので、×2すると矩形の幅と高さになる)
よって、、、
  
  vec *= near
  near_w = 2 * max( near_w, vec.x );
  near_y = 2 * max( near_y, vec.y );
  
ということです。

  
  

ここまでで、(ライトに垂直な)involve カメラによる2次元座標への投影が終わりました。
しかし、このままでは通常のライト座標系の z 軸ではないので、z テストが成功しません。
よって、y 軸を中心に 90 度回転してやり、ライトと z 軸の方向をそろえます。
ここで、involve カメラ座標系では
  x,y -> -1 <= x,y <= 1
   z -> 0 <= z <= 1
の範囲になっています。
 
これを、y 軸周りに 90 度回転させているので、、
   y,z -> -1 <= y,z <= 1
   x -> 0 <= x <= 1
となってしまっています。
 
これでは、ライトカメラの近平面で x 軸方向に半分までしか描画されませんし、
z テストも不正な値で処理されてしまいます。
 
これを回避するため、まずはスケール変換から行います。
 
x は 0 <= x <= 1 となっています。
これを -1 <= x <= 1 と拡張したいので2倍にします。
  x -> 0 <= x <= 2
  
逆に z は -1 <= z <= 1 となっているので、これを 0 <= z <= 1 と収まるように半分にします。
  z -> -0.5 <= z <= 0.5

y は通常のままで ok なので変化ナシです。
 
次に、移動処理。
 
拡張の結果、x は 0 <= x <= 2 となっています。
これを -1 <= x <= 1 としたいので、-1 します。
  x -> -1 <= x <= 1
  
次に z は -0.5 <= z <= 0.5 となっていますので、これを 0 <= z <= 1 としたいので、+ 0.5 します。
  z -> 0 <= z <= 1

  
  


これで、全ての計算が終わりました。
ここまで求めた行列は、キャスターが全て納まり、尚且つ表示カメラの視錐台の範囲で歪むように作られています。
  
この行列を用いて、キャスターの位置を歪ませてやれば、視錐台に影響するキャスターを逃すことがなく、且つ都合のいい座標変換を行うことが出来ます。
 
最終的な行列の計算順序は
  カメラ行列 -> LiSPSM 行列 → カメラ透視投影行列
ということになります。
 
今までの処理はカメラのローカル座標系上での変換処理だったことを思い出してください。




LiSPSM も PSM も考え方の出発点はカメラ空間がライト方向からどう見えるか?ということと、同時座標系の根本的な考え方を利用して作られています。
これをさらに発展させた考え方に TSM がありますが、また別の機会に紹介します。

ものすっごい久しぶりの株

2008年10月21日 00:08

最近の株安を受けて、株を買う決意をしました。

善は急げ。
ってことで、商船三井と住友物産を購入。
二つとも配当利回りがいい!
住友物産は+円高による資源の買いあさりを期待。


海運業は、これから不況になるので物流が減るかなと思ったんですが、
輸出国の日本ではなくてはならない業種なので楽観視してます。
物を作るときの部品や出来上がった製品、すべて海運だしね。


ってことで、久々の株。
今度は前回のようにスイングで運用じゃなくてホールドで運用するつもりなので配当利回りを優先してみました。
・・・でもキャピタルゲインもほしいなぁ(・∀・)
機が熟せば売ります( ・`ω・´)

ってことで、あがるの期待してます!
がんばれニッポン!

線形相補性問題

2008年10月08日 23:45

今まで線形計画問題について利用する事例や解法、グラフ理論との関係などを説明してきました。
以上を踏まえて、線形相補性問題(LCP)に入ります。
線形相補性問題とは目的関数が2次関数で制約条件が1次関数の非線形計画問題のことで、線形制約凸2次計画と同じ意味になります。

実際、ゲーム制作のどんな場面で使用されているか?というと・・・主に物理演算で使用します。
昨今、ミドルウェアとして出回っている物理エンジン(Havok や Open Dynamics Engine など)は制約式を速度や角速度、目的関数を距離などと置き、LCP を解くことで物体間の衝突解決やラグドール処理などを実現しています。
例えば2つの凸体間の最短距離を求めたい場合を考えてみます。
物体上の点をP1・P2とすると平方距離を |P1-P2|^2 と表現できます。これに物体上の点であるという制約式を加えて、最小値を求めるような式を作ればいいわけですね。
また、元になる幾何的な理論はエピポーラ幾何(ステレオビジョンの幾何学)などにも応用できます。
ゲーム分野でも応用範囲の広い分野です。


そもそも非線形計画問題(NLP)とは何でしょう?
線形計画問題が線形に変化する(1次式)問題の最適解を求めていたわけですから、非線形計画問題とは非線形(例えば2次関数とか)に変化する問題の最適解を求める問題なのかな?と推測できます。
その通りですw

非線形計画問題とは2次以上の目的関数および制約条件式の最適解を求める問題全般を指します。
この問題は線形計画問題のように単純に解が求まるか?というと高次関数についてはほぼムリです。
なぜなら4次などの高次関数は部分的な極限を複数含むので、解の探索において部分的な極限に行き着いてしまうとそこが最適解と判断してしまい探索をやめてしまうからです。
ですから、ある一定範囲内での部分的な最適解は求まりますが全体的な最適解を求めるのは極端に難しいわけです。
ただ、3次以上の問題は2次の問題を包括的に含むので、適当な範囲内で区切ってやれば2次問題の組み合わせとして帰結させることが出来ます。
(詳しい人の突っ込み待ってます)


上記の理由から、2次計画問題は非線形計画問題の中でも特別な問題として処理されます。
線形相補性問題はそんな2次計画問題の中で以下の制約で成り立った問題です。
 ・目的関数が2次式
 ・制約条件が1次式
 ・変数(目的関数や制約条件に使用されている未知数)が非負


実際に2次計画問題を解く方法として Lemke 法や内点法というのがあります。
Lemke 法というのは、シンプレックス法を拡張した様な手法で、内点法というのは幾何的な手法で徐々に解に近づけていく方法です。
大規模な(制約式が20以上あるとか)問題に対しては内点法の方が処理が早いですが、小規模な問題に対しては並列化が出来る Lemke 法の方が処理が早いと思われます。
また、内点法は反復的に最適解を求めていくものなので、精度は反復回数と初期値に依存してしまいます。
内点法についてはもう少し詳しく説明できるようになったら補足します。。(内点法自体、概念しか勉強していないんです;)



そんなわけで、Lemke 法に絞って説明します。

Lemke 法。正確には Lemke の相補吐き出し法というらしいです。
やってることはシンプレックス法に非常に良く似ています。
シンプレックス法ではスラック変数を導入し基底を入れ替えていくだけでは解けない問題があるのですが、
そういう場合、人為変数というものを導入して解いていきます(2段階単体法っていうのかな?)
Lemke 法では最初からこの人為変数を導入し、解いていくことになります。
また、相補条件というのが追加されていることと双対条件を利用してること以外は基本的にシンプレックス法の拡張です。



以下説明。
かなり長いかも。


まず適当な凸二次計画法を定義します。
目的関数を T、 制約式を C と記述していきます。

 T:
s_①-1

 C :
s_①-2


まず、目的関数のヘッセ行列(偏微分の係数を対象になるように配置した行列)を作ります。

s_②


ヘッセ行列を利用すれば、2次項の係数を行列として簡単に表せることが出来ます。
よって元の式は、

s_③


と書けます。
ここまでが前段階です。

ちなみに、以前少し書きましたが、ヘッセ行列の特徴として正値対称行列ならばその関数は極限:最小値を持ち、
負値対称行列なら極限:最大値を持ちます。
行列の固有値が全て正(正値対称行列)なのか負(負値対称行列)なのかで最小を持つか最大を持つかが異なるわけですね。

なんでこう言えるのか?ってことを簡単に説明すると、
1:2次式は平行移動で2次項と定数項だけに変換できる(※)
2:定数項を無視すると、1の式はヘッセ行列と各軸の掛け算で表せる。
s_④-1
3:ヘッセ行列は対称行列になるのでハウスホルダー法などで主軸変換が必ず可能、尚且つかける行列は正規直行行列(要するにただの座標変換)
s_④-2
4:正規直行行列が掛かってるだけなので極限の有無には影響しないことを理由に、元の問題は単純に s_④-3 と書ける。
5:ここで、s_④-4 には主軸変換した後のヘッセ行列の固有値が入り、[画像④-5]は x がどんな実数でも必ず正なので、結果的にs_④-4の値に依存する。
  また、極値を持つのは s_④-4が全て同じ符号の時のみ。
  (s_④-4で正負両方が存在する場合は鞍点型になるし、どれか 0 が入ると 0 が入った軸上では一定の値を持つことになるので極値を取る点は無限に存在する)


 
※間単に照明すると、、、
 2次式の一般系s_*1はヘッセ行列s_*2-1を用いて書くと、s_*3になります。(s_*2-2 a s_*2-2 x
 
 これに平行移動量 p を追加すると、以下の様な式になります。
 s_*4-1
 
 ここで①の変形には、s_*4-2s_*4-3 を用いてます。
 
 ①の式にs_*5が成り立つように p を選べば、
 s_*6-1
 のように変形できます。
 ここで、s_*6-2を利用しています。
 
 よって、どんな2次式でも p に適切な値を入れてやれば1次項が消えますので、2次項と定数項のみで表現できるわけです。


 

実際は、固有値を求める必要は必ずしもなくて、シルベスタの定理というのを利用すれば主小行列式を順番に見ていくだけで解ります。
小行列式は一つ大きい小行列式の中でも使われるので、再帰的に行列式を求める割には計算量も掛かりません。
(公式などを使わずに行列のランクだけみて行列式を求める場合、小行列を再帰的に作って行列式を求めていきますから、コスト的にはまったく一緒ですね)

実際、固有値を求めるコストの方が圧倒的に高いです(特異値分解(SVD) だったり QR 法だったり)



閑話休題。
それでは Lemke 法の処理手順を書いていきます。

っと、その前に。
実は制約式も大きな行列に出来ます。
挙げた例の制約条件をよく見てみましょう。
ただの一次式ですよね?
等号が同じ方向なら纏めてしまえそうです。
っていうことで、、、

s_①-2

このような制約式は

s_⑤

と書くことで、

s_⑥

このような行列を利用した形で記述できます。


ここまでで前提条件が揃いました。
では纏めて行きます。
例で挙げた式を解いていきます。

1:目的関数と制約式を洗い出す。

2:目的関数と制約式を纏めた大きな行列を作る。
  s_⑧
  
3:2の行列との関係式を作る。
  s_⑨-1 alpha
  s_⑨-1 beta
  s_⑨-1 c
  s_⑨-1 delta
  s_⑨-1 gamma
  s_⑨-1 tau
  と置き、
  s_⑨-2
  を作る。
  
  αは元の問題の最適解、βは双対問題の最適解で、τは元の問題のスラック変数、σは双対問題のスラック変数になります。
  α、β、σ、τはまだ解らない値なので、不定数に置き換えます。
  
  さらに、これらを以下のように置き換えると、
  s_⑩-1
  s_⑩-2
  s_⑩-3 delta
  s_⑩-3 sdelta
  
  こんな形の一次式の形になります。
  s_⑩-4
  
4:この時点で、s_qi)0.pngが成り立つなら、s_δi=qis_δi と置けば、全ての条件が成り立ちます。
  よって、その場合はここで終了です。
  q の要素に一つでも 0 より小さい値が入ってる場合は 5 へ。
  
5:3の行列を解く為、人為変数[画像(R)]を導入します。
  人為変数を十分大きな値にしておけば、s_artifact.pngs_δi となり、式が成り立ちます。
  なんで、こんな変数を導入できるか?というと、、、一番最後、制約式と目的関数を満たすものが求まったときに、
  人為変数の係数が 0 になっていれば問題ないので OK なんです(結構強引)
  例 :
    s_□1-1
    s_□1-2
    ↓
    s_□2-1
    s_□2-2


6:ここからが本番。
  まず、定数項が一番小さい値を人為変数に置き換えます。
  人為変数は任意の数で、且つ式を成り立たせるためだけに導入した変数です。
  s_qi)0.png(この場合 s1~s6)を成り立たせれば人為変数の目的としては OK なので、
  式を成り立たせるためには一番小さな値を人為変数の値として代入しておけば、s_qi)0.png(この場合x1~x6)でも成り立ちますし、全ての条件を満たします。
  
  よって、一番小さい定数項の式を探し出し、基底を人為変数に置き換えます。
  例 :
    s_□2-1
    s_□2-2
          ↓
    s_□3-1
    s_□3-2

  
7:ここで、人為変数が基底になっているので、基底になっていない変数の対(これを非基底対と称します)が必ず1つだけ存在します。
  逆にこの時点で人為変数が基底になっていなく、人為変数が 0 でも制約式を満たすのならば終了です。

  また、非基底対の片方は直前まで基底になっていた変数です。
  これを基底にしても直前の状態に戻るだけですので、非基底対で直前まで基底になっていなかった変数を新しい基底に選びます。

  次に、この新しい基底を増やしていくと一番早く定数項と等しくなる式を選びます。
  要するに各基底式において、新しく選んだ基底と定数項以外を 0 と見なして、定数項と同じ値になるように新しく選んだ基底に変数を代入して行き、最も小さな値で定数項と同じになる式を選び基底を取り替えるわけです。
  ・・・書いてて意味が解らなくなってるので、例で詳しく書きます。
  
  例 :
    ここで、非基底対は s5, x5。s5 は前回の基底なので x5 を新たな基底に選ぶとします。
    そうすると、基底の部分と定数項以外はとりあえず関係ないので 0 と見なすと、、、
    s_□4-1

    ここで、実際に関係あるのは基底 s1 と基底 s3 の式だけですね。
    よって、この二つの式で x5 を徐々に増やしていくと一番早く定数項と等しくなる式を見つけます。
    まぁ、考えればすぐわかりますが、 s3 の式の方が早いですね。
    x5 に 21 / 5 を代入すれば定数項と同じ 21 になります。

    そんなわけで、今度は s3 の基底式で x5 と s3 を入れ替えます。
    s_□4-2
    s_□4-3
    
    新しい基底式が出来上がったので、これを他の基底式に代入していきます。
    っていっても、今までの流れの中からでもわかりますが関係あるのは s1 の基底式だけですね。

  
8:基底をいくら増やしても基底変数が負にならないときは終了。それ以外は 7を繰り返し行います。

例で挙げた式を繰り返しやっていくのは面倒なので、答えだけ記述します。

答えは
x = 3.8
y = 1.42
z = 2.64
となります。



これで、LCP の問題が解けます。
実際は堂々巡り(退化)したときの処理を追加したりしないといけません。たしかw




2~3が成立する理由ですが、、、LCPの最適解はヘッセ行列を制約式と見立てたLPの最適解と解が
等しいということから来ています。

上が成立する理由は最適解じゃないと仮定して式を変形していくと矛盾した結果になるので解るのが一つと、
LP の場合双対問題が成立するのを利用してさらに式を変形していけば求まります。

んが!かなり長いので割愛します。
時間があれば次回乗っけようと思います。


んで、サンプルソース。

  LCP.txt
  
かなり適当に組んでいるので、まともに動く保障はありません。
それとソースがなかなかに汚いので、参考にすらならないかもしれません。

概念がわかれば、シンプレックス法の延長ですからプログラムに起こすのは簡単だと思います。
例外処理が面倒でしたけど。。
これより解りやすくて見やすくて堅牢なのが出来たら、ぜひご教授くださいw
(僕が作ったソースじゃ遅いと思うし)


最近の記事


上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。