1. 程式人生 > >OpenCV中LU分解實現矩陣求逆invert(DECOMP_LU)

OpenCV中LU分解實現矩陣求逆invert(DECOMP_LU)

OpenCV3.0中實現矩陣求逆有四種方法(LU、cholesky、eig以及SVD),使用187*187矩陣測試,100輪計算耗時如下(cholesky方法並沒有得到結果):


179*179矩陣測試,100輪耗時如下:


Cholesky方法對輸入矩陣限制嚴格(共軛對稱和正定矩陣)LU分解耗時較短且對輸入矩陣沒有限制,因此重點對LU分解進行分析。其實現原始碼如下

LU分解比列主元消去法耗時短很多。

OpenCV實現的LU分解並沒有按照原始的LU分解求解L和U矩陣的逆然後再進行矩陣乘法。


其思想是對增廣矩陣做初等行變換,為了將除法轉換為乘法,將對角線矩陣設定為其倒數,後面只計算乘法即可。由於工作需求,需要將其進行定點化在DSP上實現,下一篇部落格分享定點化程式碼。

// LU decompose method to solve inverse matrix
BOOL MatrixInv(F32 B[3][3], const S32 A[3][3], const U32 n)
{
    S32 ii;
    U32 i, j, k;
    F32 temp, d, s, alpha;
    F32 t[3][3];


    for (i = 0; i < n; i++)
    {
        for (j = 0; j < n; j++)
        {
            t[i][j] = (F32)(A[i][j]);
            B[i][j] = 0.0f;
        }
        B[i][i] = 1.0f;
    }


    for (i = 0; i < n; i++)
    {
        // find principal element
        k = i;
        for (j = i + 1; j < n; j++)
            if (fabs(t[j][i]) > fabs(t[i][i]))
                k = j;


        // principal element is zero so A is Not full rank matrix and has no inverse matrix
        if (fabs(t[k][i]) < FLT_EPSILON)
        {
            return FALSE;
        }
        // if principal element is not in row i,switch row i and row k
        if (k != i)
        {
            for (j = 0; j < n; j++)
            {
                temp = t[i][j];
                t[i][j] = t[k][j];
                t[k][j] = temp;
                temp = B[i][j];
                B[i][j] = B[k][j];
                B[k][j] = temp;
            }
        }


        temp = -1*(t[i][i]);//d = -1.0f / t[i][i];
        d = 1 / temp;
        for (j = i + 1; j < n; j++)
        {
            alpha = t[j][i] * d;
            for (k = i + 1; k < n; k++)
                t[j][k] += alpha*t[i][k];
            for (k = 0; k < n; k++)
                B[j][k] += alpha*B[i][k];
        }


        t[i][i] = -d;
    }


    for (ii = n - 1; ii >= 0; ii--)
        for (j = 0; j < n; j++)
        {
            s = B[ii][j];
            for (k = ii + 1; k < n; k++)
                s -= t[ii][k] * B[k][j];
            B[ii][j] = s*t[ii][ii];
        }


    return TRUE;
}

手動計算了一個矩陣求逆