OpenCV中LU分解实现矩阵求逆invert(DECOMP_LU)

一时失言乱红尘 2022-08-21 11:57 279阅读 0赞

OpenCV3.0中实现矩阵求逆有四种方法（LU、cholesky、eig以及SVD），使用187\*187矩阵测试，100轮计算耗时如下（cholesky方法并没有得到结果）：

![Center][]

179\*179矩阵测试，100轮耗时如下：

![Center 1][]

Cholesky方法对输入矩阵限制严格（共轭对称和正定矩阵）LU分解耗时较短且对输入矩阵没有限制，因此重点对LU分解进行分析。其实现源码如下

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

OpenCV实现的LU分解并没有按照原始的LU分解求解L和U矩阵的逆然后再进行矩阵乘法。

![Center 2][]

其思想是对增广矩阵做初等行变换，为了将除法转换为乘法，将对角线矩阵设置为其倒数，后面只计算乘法即可。由于工作需求，需要将其进行定点化在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;  
\}

手动计算了一个矩阵求逆

![Center 3][]

[Center]: /images/20220724/de98a1a2dadf424a9feb4ebf1fbb8138.png
[Center 1]: /images/20220724/2119d0ce28f2419b9a1179d15dee17d7.png
[Center 2]: /images/20220724/9ce8133e89d04571b9829ba045beaea7.png
[Center 3]: /images/20220724/d5d649f98e534e7fa47290f9582b873f.png