はじめに
本記事はOpenCV Advent Calendar 2021 の10日目の記事です．
この記事では，「逆行列は計算してはいけない」をOpenCVでどうやって実現するのか，一番簡単な最小二乗法を例にして，特異値分解を使って行う方法を説明します．
特異値分解 (singular value decomposition: SVD) とは，行列$A$を下記の3つの行列の積に分解する方法です．
A = U\Sigma V^t\\
\Sigma = \rm{diag} \{ \sigma_1, \sigma_2, \cdots, \sigma_n\}\\
\sigma_n は特異値
本記事は，この行列分解を使って最小二乗を行う方法を説明します．
最小二乗と一般逆行列
行列Aの本数が多い場合下記行列の解は最小二乗で解かれることが多いです．
b = Ax\\
x = \arg\min_x\|Ax-b\|_2^2
一般逆行列を以下とすると解は
A^+ = (A^tA)^{-1}A^t\\
x = A^+b
で求めることができます．
これを直接計算するようなOpenCVのコードは以下です．
Mat x = (A*A.t()).inv()*b;
ただし，この計算は，逆行列の計算$（A^tA)^{-1}$が途中に入るため，数値計算精度の問題が生じやすく，実際には特異値分解を使ってこの解を求めたほうが精度が高くなります．
一般に数値計算をする際，逆行列計算はだいたい省略できるので逆行列を求めることになったら適切な計算方法がないかいったん立ち止まって考えるとよいです．
（回転行列，ユニタリ行列などの転置するだけで逆行列になる場合などを除く．）
精度を上げて計算しようとすると簡単にランク落ちを始めます．
詳細はリンク先の一般逆行列に書いてあることに任せますが，ざっくり一般逆行列は特異値分解の成分を使って以下で表せます．
A^+ = V\Sigma^{-1}U^t \\
\Sigma^{-1} = \rm{diag} \{ 1/\sigma_1, 1/\sigma_2, 1/\cdots, 1/\sigma_n \}
つまり，特異値分解を使って以下で$x$が表せます．
x = V\Sigma^{-1}U^t b
OpenCVによる特異値分解と最小二乗計算
OpenCVによるSVDは，呼び出し方がいくつかありますが，例えば以下のように使えます．
Mat w, u, vt;
cv::SVDecomp(A, w, u, vt, SVD::FULL_UV);
ここで，uはU，vtはVの転置した行列で，wはΣの対角要素を1次元のベクトルとして保持したものです．
最小二乗を計算するには，Σの逆数だったり，Uの転置，Vtの転置の転置等を計算すれば最小二乗の計算ができますが，もっと便利が関数backSubstが用意されています．
SVD::backSubst(w, u, vt, b, x);
この関数により，最小二乗をするために，以下の式（再掲）を求めるために，SVDの計算結果を適切に適用してくれます．
この関数は，w,u,vt,bを入力すれば，dest相当のxに計算結果を出力してくれます．
x = V\Sigma^{-1}U^t b
線形フィッティングの例
ax+bの直線をSVDを使って最小二乗で計算する例を下記に示します．
実際には，観測したx,yの値を入れますが，シミュレーションとして下記を入れます．
- a=2.0, b=0.5として，観測結果は標準偏差0.1のノイズが乗ったもの．
- 観測位置ｘも乱数で与える．
- ax+bの式なので，A(n, 1)は常時1です．
//y = c0*x+c1
double c0 = 2.0;
double c1 = 0.5;
const int n = 1000;
Mat A(n, 2, CV_64F);
Mat b(n, 1, CV_64F);
RNG rng;
RNG noise;
double sigma = 0.1;
for (int i = 0; i < n; i++)
{
double x = rng.uniform(0.0, 100.0);
double y = c0 * x + c1 + noise.gaussian(sigma);
A.at<double>(i, 0) = x;
A.at<double>(i, 1) = 1.0;
b.at<double>(i) = y;
}
Mat w, u, vt, x;
cv::SVDecomp(A, w, u, vt, SVD::FULL_UV);
SVD::backSubst(w, u, vt, b, x);
cout << x << endl;
return 0;
出力
[1.999813248601801;
0.5132182461533308]
まとめ
OpenCVによる特異値分解と最小二乗計算について説明しました．
ストレートに考えると逆行列を呼び出したくなりますが，
「逆行列を計算するのはやめましょう．」
遅いのに加えて，数値計算精度が危険です．
ここまで説明してなんですが，OpenCVの数値計算精度等は，よく洗練された行列ライブラリと比べて，はてなマークがつく箇所がいくつもあります．
バグではありませんが，より正確な結果が必要な時はBLASやBLAS，やっぱりBLASを使いましょう．
BLASの命名規則に心が折れた人は，Eigenなどのもっとユーザーフレンドリーなものを使いましょう．
OpenCVでの行列演算は，正直あまりお勧めできるものではありません．
フルランクに近い，普通の計算は何やっても変わりませんが，ギリギリランク落ちするような計算は，結構不安な挙動をします．
カメラキャリブレーション周りとか，ちゃんとしたライブラリで書き直したら精度上がる気がします．
あと，キャリブレーションパターンを何百枚も入れて最適化すると，数時間単位で計算時間がかかるの絶対にコードの作り方が間違ってます．．．
おまけ：重み付き最小二乗法
重み付き最小二乗法は重みを表す対角行列$W$を用いて下記で表せます．
サンプル点がそれぞれ，どれくらいの重みをもつかを$W$の要素で決めます．（中で二乗されるのでsqrtを取っておく．）
x = (A^tW^2A)^{-1}AW^2b
これはOpenCVのSVDでは，こんな感じで呼び出せば使えます．証明は頑張って数式を展開してみてください．
（OpenCVの挙動的にこれであっていると思いますが，間違ってたら誰か教えてください．）
cv::SVDecomp(W * A, w, u, vt, SVD::FULL_UV);
SVD::backSubst(w, u, vt, W * b, x);
Eigenを使った場合はこんな感じです．
backSubstの代わりにsolveメソッドがあります．
下記以外にもSVDの計算にJacobi法も指定できます．
Eigen::BDCSVD<Eigen::MatrixXd> svd(W*A, Eigen::ComputeThinU | Eigen::ComputeThinV);
Eigen::MatrixXd x = svd.solve(W*b);
また，こんな感じでEigenとOpenCVは行き来できます．
cv::Mat x_cv
cv::eigen2cv(x, x_cv);