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

前回に続き今回は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のソースコードを見てみましょう．

#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の表現力が足りず，結果として情報が抜けてしまい，このような結果になってしまいました．