c++ - Mat.inv() 在 opencv 中产生全零

Question

我有以下代码:

    cv::Mat temp0 = R.t();
    cv::Mat temp1 = R * temp0;
    cv::Mat temp2(960, 960, CV_32FC1);
    temp2 = temp1.inv();

    cv::Size s = temp1.size();
    std::cout<<s.height<<" "<<s.width<<std::endl;

    std::cout<<cv::format(temp1, "numpy" ) << std::endl;
    std::cout<<cv::format(temp2, "numpy" ) << std::endl;

Transpose 工作正常，矩阵乘法也是如此。因此 Mat temp1 的大小为 960x960。但是，当我执行 temp2 =temp1.inv() 时，我在 temp2 中收到所有零。我的意思是零是所有 960x960 单元格。此外，R 仅属于 CV_32FC1 类型。所以这可能不是数据类型问题。我无法理解这里的问题。我用谷歌搜索了很多。你能帮忙吗？

编辑

我正在复制 Mat::inv() 函数的 gdb 输出。我很难弄清楚，但如果有人更熟悉 OpenCV，也许会有帮助:)

Breakpoint 1, CreateShares::ConstructShares (this=0x80556d0, channel=..., k=2, n=4) at CreateShares.cpp:165
165     temp2 = temp1.inv();
(gdb) step

cv::Mat::operator= (this=0xbffff294, e=...) at /usr/include/opencv2/core/mat.hpp:1373
1373        e.op->assign(e, *this);
(gdb) 
1374        return *this;
(gdb) step
1375    }    
(gdb) step
cv::MatExpr::~MatExpr (this=0xbfffef64, __in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:1167
1167    class CV_EXPORTS MatExpr
(gdb) step
cv::Mat::~Mat (this=0xbfffefdc, __in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:295
295     release();
(gdb) step
cv::Mat::release (this=0xbfffefdc) at /usr/include/opencv2/core/mat.hpp:381
381     if( refcount && CV_XADD(refcount, -1) == 1 )
(gdb) step
383     data = datastart = dataend = datalimit = 0;
(gdb) step
384     size.p[0] = 0;
(gdb) step
385     refcount = 0;
(gdb) step
386 }
(gdb) step
cv::Mat::~Mat (this=0xbfffefdc, __in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:296
296     if( step.p != step.buf )
(gdb) step
298 }
(gdb) step
cv::Mat::~Mat (this=0xbfffefa4, __in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:295
295     release();
(gdb) step
cv::Mat::release (this=0xbfffefa4) at /usr/include/opencv2/core/mat.hpp:381
381     if( refcount && CV_XADD(refcount, -1) == 1 )
(gdb) step
383     data = datastart = dataend = datalimit = 0;
(gdb) step
384     size.p[0] = 0;
(gdb) step
385     refcount = 0;
(gdb) step
386 }
(gdb) step
cv::Mat::~Mat (this=0xbfffefa4, __in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:296
296     if( step.p != step.buf )
(gdb) step
298 }
(gdb) step
cv::Mat::~Mat (this=0xbfffef6c, __in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:295
295     release();
(gdb) step
cv::Mat::release (this=0xbfffef6c) at /usr/include/opencv2/core/mat.hpp:381
381     if( refcount && CV_XADD(refcount, -1) == 1 )
(gdb) step
383     data = datastart = dataend = datalimit = 0;
(gdb) step
384     size.p[0] = 0;
(gdb) step
385     refcount = 0;
(gdb) step
386 }
(gdb) step
cv::Mat::~Mat (this=0xbfffef6c, __in_chrg=<optimized out>) at /usr/include/opencv2/core/mat.hpp:296
296     if( step.p != step.buf )
(gdb) step
298 }
(gdb) step
CreateShares::ConstructShares (this=0x80556d0, channel=..., k=2, n=4) at CreateShares.cpp:167
167     cv::Size s = temp1.size();
(gdb) step
cv::Mat::MSize::operator() (this=0xbffff284) at /usr/include/opencv2/core/mat.hpp:705
705     return Size(p[1], p[0]);
(gdb) step
cv::Size_<int>::Size_ (this=0xbffff2f8, _width=960, _height=960) at /usr/include/opencv2/core/operations.hpp:1624
1624        : width(_width), height(_height) {}
(gdb) step
cv::Mat::MSize::operator() (this=0xbffff284) at /usr/include/opencv2/core/mat.hpp:706
706 }
(gdb) step

Answer 1

行列式很可能为零。

来自维基百科:

A square matrix that is not invertible is called singular or degenerate. A square matrix is singular if and only if its determinant is 0.

您可以像这样显示行列式...

std::cout<<"determinant(temp1)="<<cv::determinant(temp1)<<"\n";

来自 the documentation for Mat::inv() ，共有三种方法可供选择:

• DECOMP_LU(默认)是 LU 分解。 矩阵必须是非奇异的。
• DECOMP_CHOLESKY 是仅用于对称正定义矩阵的 Cholesky 分解。这种类型在大矩阵上比 LU 快两倍。
• DECOMP_SVD 是 SVD 分解。如果矩阵是奇异矩阵甚至是非方阵，则计算伪逆。

来自 the documentation for invert() , 大概由 Mat::inv() 内部使用:

In the case of DECOMP_LU method, the function returns the src determinant ( src must be square). If it is 0, the matrix is not inverted and dst is filled with zeros.

这与您看到的结果一致。

---

数学笔记

我不是数学家，但我的印象是求逆矩阵可能是一件很麻烦的事情——如果您的矩阵非常大，则更是如此。事实上，这些逆函数原则上可能确实存在，但实际上不可能精确计算。在对您的代码进行一些实验时，我发现在许多情况下我会得到不完全为零但非常接近于零的行列式——这可能表明数值精度可能是限制因素。我尝试使用 64 位值而不是 32 位值指定矩阵，得到了不同但不一定更好的答案。

根据您计算 temp1 矩阵的方式，认识到它总是对称可能会很有用。 DECOMP_CHOLESKY 方法专门设计用于处理对称 positive definite。矩阵，因此使用它可能会提供一些优势。

在实验上，我发现归一化(如@cedrou 所建议的那样)使逆函数更有可能返回一个非零矩阵(使用 DECOMP_LU 而不是使用 DECOMP_CHOLESKY).但是，根据我对您可能如何初始化 R 矩阵的猜测，生成的矩阵似乎永远不会满足逆矩阵的定义:A*inverse(A)=Identity。但您不一定关心这一点——这也许就是 SVD 方法计算伪逆的原因。

最后，为什么反演失败这个更深层次的问题似乎是一个数学问题，而不是一个编程问题。基于此，我在数学网站上进行了一些搜索，结果发现有人已经问过如何做这件事:https://math.stackexchange.com/questions/182662

---

调试注意事项

根据您的调试跟踪，我倾向于认为您感兴趣的部分被编译到一个不可跟踪的库中，并在您运行 step 时跳过。换句话说，第一个 step 之后的神秘空白行代表它实际运行 inv() 函数的部分。之后，它将结果分配给 temp2 并销毁临时对象。因此，您的调试跟踪不会告诉我们有关 inv() 内部发生的事情的任何信息。

165     temp2 = temp1.inv();
(gdb) step

cv::Mat::operator= (this=0xbffff294, e=...) at /usr/include/opencv2/core/mat.hpp:1373
1373        e.op->assign(e, *this);

我自己对此运行了一个调试器，并且能够通过对 invert() 的内部调用进行跟踪，并观察它根据对矩阵的内部分析决定失败(确定它不是可逆)——因此返回一个用零填充的矩阵，与您报告的内容相匹配。

invert() 函数在 cxlapack.cpp 中定义，如果您有兴趣查看源代码。