2013年02月16日

LU分解で行列式と逆行列の計算 その2

LU分解のメモです.一年ぐらい前の記事に書いていましたが,

AX = B
A = LU
LUX = B
X = U-1L-1B

となりますので,逆行列を求めたい場合は,

X = A-1,B = E

と考えれば解けます.

行列式が0ではないが0に近い値をとる場合,LU分解中に桁落ちなどで0除算が発生して解けない場合があります.その対策として,一般的にはpivot選択が行われているようです[1].

pivot選択は,大きい値を探して順番を入れ替えるという処理が必要であるため,入れ替えた順番を保存するためのメモリが必要であったり入替えが面倒ということがあるかもしれません.

LU分解で検索したときに,対角成分で割ると数値誤差が減るという情報[2]がありましたので,ちょっと試してみました.

数式で考えると次のようになると思います.

AA-1 = E … (1)

左から対角成分を1にするための行列Mを掛けます.
M_001.png
ここで,a11, ... , annは行列Aの対角成分で,その逆数をもつ行列がMです.

MAA-1 = M

このMAをLU分解します.

MA = LU
LUA-1 = M

元の(1)式と比べると,EMに変わっただけです.

まとめると,通常のLU分解との違いは以下のようになります.
  1. 逆行列を求めたい行列Aの対角成分の逆数を持つ行列Mを作成
  2. AMを掛けてMAを作成(対角成分が1になる)
  3. MAをLU分解する(pivot選択しなくても良い)
  4. LUX = MとしてXについて解く(X = A-1となって逆行列が求まる)
やっていることはpivot選択と殆ど変わらない(大きい値を優先して選ぶか,予め対角成分の値で正規化?しておくかの違い)と思います。

実験として,4×100の適当な正規分布に従うデータを生成して,その共分散行列(4x4)の逆行列をLU分解で計算してみました.データ生成の際に,行列式が小さくなるように正規分布の分散を調整しました.

行列式が0.0001辺りの行列は,pivot選択しない場合は解けませんでしたが,今回のやり方でpivot選択することなく解けるようになりました.

ただ,やはり限界はあるようで,0.00001辺りになるとpivot選択しても今回のやり方でも解けませんでした.

精度的にも特に変わらないので,結論としては既にPivot選択でやっているのなら無理に今回の方法にする必要もなさそうです.

[1] http://akita-nct.jp/yamamoto/lecture/2004/5E/linear_equations/text/html/node4.html
[2] http://okwave.jp/qa/q2363135.html
web拍手 by FC2
posted by シンドラー at 02:19 | Comment(0) | TrackBack(0) | OpenGL & Metasequoia | このブログの読者になる | 更新情報をチェックする

2013年02月04日

点広がり関数とフーリエ変換について その2

続きです。せっかくなのでWienerフィルタを使用したボケ除去まで試してみます。

点広がり関数が分かっていて、ノイズが入っていない理想的な状態ですが、やはりぼかして情報量が減っているのか完全な復元はできないみたいですね。

本のサンプルページを参考にしているのは良くない気もしますね…。

参考サイト
[1] http://pub.maruzen.co.jp/realize/search/sample/RLB220496.pdf

続きを読む
web拍手 by FC2
posted by シンドラー at 23:51 | Comment(0) | TrackBack(0) | Image Processing | このブログの読者になる | 更新情報をチェックする

広告


この広告は60日以上更新がないブログに表示がされております。

以下のいずれかの方法で非表示にすることが可能です。

・記事の投稿、編集をおこなう
・マイブログの【設定】 > 【広告設定】 より、「60日間更新が無い場合」 の 「広告を表示しない」にチェックを入れて保存する。