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;
}
手動計算了一個矩陣求逆