スポンサーサイト

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

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

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

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ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。