OpenCVをつかってY=AXを満たすAを最小二乗法で求める

[前回](http://qiita.com/tkymx/items/b1dd4ed1693d288f59dc)に続き今回はY=AXのAの部分を求めていきたいと思います．

#解きたい問題

$Y=AX$

今回はYとXが固定されているときのAを求めてみたいと思います．

$A = YX^{-1}$

と，いきたいところですが$X$が正方行列でなくてはいけないそうなので，以下のように変形して$A$を求めていきたいと思います．

$
\begin{align*}
Y & = AX \\\\
YX^T  & = AXX^T \\\\
YX^T(XX^T)^{-1} & = A 
\end{align*}
$

これで$A$が求まりました．それではOpenCVのソースコードを見てみましょう．

```cpp

#define _CRT_SECURE_NO_WARNINGS

#include <stdio.h>
#include <opencv\cv.h>

#define SIZE 4
#define COUNT 6

void show(CvMat* mat, const char* str)
{
	printf("%s\n",str);
	for (int i = 0; i < mat->rows; i++)
	{
		for (int j = 0; j < mat->cols; j++)
		{
			printf("%f ", cvmGet(mat, i, j));
		}
		printf("\n");
	}
}

int main(void)
{
	CvMat *mA, *mX , *mAX, *mX_T, *mXX_T, *mXX_T_inv, *mY, *mYX_T;

	mA = cvCreateMat(SIZE, SIZE, CV_32F);
	mX = cvCreateMat(SIZE, COUNT, CV_32F);
	mX_T = cvCreateMat(COUNT, SIZE, CV_32F);
	mXX_T = cvCreateMat(SIZE, SIZE, CV_32F);
	mXX_T_inv = cvCreateMat(SIZE, SIZE, CV_32F);
	mY = cvCreateMat(SIZE, COUNT, CV_32F);
	mAX = cvCreateMat(SIZE, COUNT, CV_32F);
	mYX_T = cvCreateMat(SIZE, SIZE, CV_32F);

	cvSet(mA, cvScalar(0.0f));
	for (int i = 0; i < SIZE; i++)
	{
		for (int j = 0; j < COUNT; j++)
		{
			cvmSet(mX, i, j, rand()%10);
			cvmSet(mY, i, j, rand() % 10);
		}
	}


	// 逆行列を求める。
	cvmTranspose(mX, mX_T);

	show(mX_T, "mX_T");

	cvmMul(mX, mX_T , mXX_T );

	show(mXX_T, "mXX_T");

	cvmInvert(mXX_T, mXX_T_inv);

	show(mXX_T_inv, "mXX_T_inv");

	// ベクトルAを求める。
	cvmMul( mY, mX_T, mYX_T );

	show(mYX_T, "mYX_T");

	cvmMul(mYX_T, mXX_T_inv, mA);

	show(mA, "mA");

	cvmMul(mA, mX, mAX);

	show(mAX, "mAX");

	printf("mY=mAmX\n");

	// 出力
	for (int i = 0; i < SIZE; i++)
	{
		printf("(");

		for (int j = 0; j < COUNT; j++)
		{
			printf("%f ", cvmGet(mY, i, j));
		}

		printf(") = (");

		for (int j = 0; j < SIZE; j++)
		{
			printf("%f ", cvmGet(mA, i, j));
		}

		printf(") (");

		for (int j = 0; j < COUNT; j++)
		{
			printf("%f ", cvmGet(mX, i, j));
		}
		printf(")\n");
	}

	cvReleaseMat(&mA);
	cvReleaseMat(&mX);
	cvReleaseMat(&mX_T);
	cvReleaseMat(&mXX_T);
	cvReleaseMat(&mXX_T_inv);
	cvReleaseMat(&mY);
	cvReleaseMat(&mAX);
	cvReleaseMat(&mYX_T);

	return 0;
}
```
**実験結果**

```
mX_T
1.000000 1.000000 2.000000 7.000000
4.000000 1.000000 1.000000 5.000000
9.000000 5.000000 8.000000 3.000000
8.000000 7.000000 7.000000 2.000000
2.000000 1.000000 1.000000 3.000000
5.000000 2.000000 9.000000 1.000000
mXX_T
191.000000 118.000000 181.000000 81.000000
118.000000 81.000000 111.000000 46.000000
181.000000 111.000000 200.000000 69.000000
81.000000 46.000000 69.000000 97.000000
mXX_T_inv
0.106271 -0.089161 -0.040635 -0.017554
-0.089161 0.128615 0.006182 0.009063
-0.040635 0.006182 0.036640 0.004937
-0.017554 0.009063 0.004937 0.017158
mYX_T
140.000000 97.000000 151.000000 94.000000
100.000000 70.000000 104.000000 87.000000
145.000000 87.000000 118.000000 97.000000
71.000000 45.000000 64.000000 105.000000
mA
-1.556606 1.778637 0.907525 0.779881
-1.367388 1.518452 0.609352 0.885200
1.154647 -0.130138 -0.551793 0.490036
-0.910809 0.804589 0.256460 1.279059
mAX
7.496251 0.359145 4.483577 7.910050 1.912594 4.721852
7.566167 1.084252 2.816179 5.725918 2.048627 2.569327
3.351180 6.386842 6.796908 5.443742 3.097474 1.036865
9.360111 3.813107 1.714520 2.698987 3.076607 0.642331
mY=mAmX
(7.000000 0.000000 4.000000 8.000000 4.000000 5.000000 ) = (-1.556606 1.778637 0.907525 0.779881 ) (1.000000 4.000000 9.000000 8.000000 2.000000 5.000000 )
(7.000000 1.000000 2.000000 6.000000 4.000000 3.000000 ) = (-1.367388 1.518452 0.609352 0.885200 ) (1.000000 1.000000 5.000000 7.000000 1.000000 2.000000 )
(2.000000 6.000000 5.000000 6.000000 8.000000 2.000000 ) = (1.154647 -0.130138 -0.551793 0.490036 ) (2.000000 1.000000 8.000000 7.000000 1.000000 9.000000 )
(9.000000 4.000000 1.000000 3.000000 4.000000 1.000000 ) = (-0.910809 0.804589 0.256460 1.279059 ) (7.000000 5.000000 3.000000 2.000000 3.000000 1.000000 )
```

#考察
実験結果で表示されているmAXが演算によって求めたmYの値になります．mY=mAmXで示されているmYの値と比べてみてください．だいたいは同じですが多少異なっていることが分かります．これは，演算によって情報が抜けてしまったからです．今回最小二乗法に用いたYとXは4*6の行列です．しかし，Aは4*4しかありません．よってYとXの変換を完璧に行うにはAの表現力が足りず，結果として情報が抜けてしまい，このような結果になってしまいました．