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;
}
手动计算了一个矩阵求逆