#連立一次方程式
AX=B
を満たすxを求める問題は非常に良く使う数値計算です.解法アルゴリズムとしては,様々なものがあります.
- シンプルに逆行列A-1を求める(数値計算に詳しい研究室だとまずはじめに教えられること.「逆行列の計算はやってはいけない.遅すぎる.」O(n3).)
- ガウスの消去法 (O(n3).逆行列を求めるのと同じ.遅すぎる.講義で習った方法.)
- LU分解(O(n3))
- コレスキー分解(Cholesky)
- 固有値分解
最小二乗問題を解く場合.Aが正則な行列ではない場合
- 特異値分解
- QR分解
#OpenCVの実装
OpenCVでは,これをcv::solve
で提供しています.下記関数は,線形システム(入力行列Aが正則,非特異である場合)か,最小二乗問題を解く関数.
bool cv::solve(InputArray A, InputArray B, OutputArray X int flags=DECOMP_LU)
flags:
- DECOMP_LU: LU分解
- DECOMP_CHOLESKY: コレスキー分解. 行列Aは正定値対称行列であること.LU分解より高速.
- DECOMP_EIG: 固有値分解.行列Aは対称行列であること.
- DECOMP_SVD: 特異値分解(singular value decomposition: SVD)
- DECOMP_QR: QR分解(現在,SVDをかわりに呼び出す実装になっているためSVDと答えが同じです.いいのか!?).
- DECOMP_NORMAL: 上記以外だった場合,下記システムをかわりに解く.
A^TAX=A^TB
実装では,LU分解やコレスキー分解は,halモジュールを読んでいます.固有値分解ではヤコビ法を自前で実装して解いています.
特異値分解では,JacobiSVD法を自前で実装して解いています.QR分解はSVDと同じとして,SVDで解いてます.条件が揃えば,SVDは固有値分解と同じになるため演算が速い固有値分解をコールするように作られています.
詳細はopencv/core/lanpack.cppのコードを参照してください.
しかしながら,精度や速度を考えるとより適切なアルゴリズムが多く提案されています.Eigenは以下の演算がモジュールとして用意されるほど数値計算方法が充実しています.もし,大きな行列の線形解法が必要なときは,速度・精度の面からOpenCVのソルバを使うのが適切か一度検討してみましょう.
- LU
- Cholesky
- Householder
- SVD
- QR
- Eigenvalues
#まとめ
OpenCVの線形ソルバについてまとめました.
線形ソルバとしてはOpenCVは小回りを利かせるためには,必要最低限のものはそろえていますが,充実しているとはいえません.
より深く勉強するためには,Eigenのチュートリアルなどが参考になります.
また,このスライドやこの資料も良くまとまっています.
solveするにもいろいろあるということが分かると思います.