7

我有以下代码：

    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;

转置工作正常，矩阵乘法也是如此。因此，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
4

3 回答 3

10

最有可能的是，行列式为零。

来自维基百科：

不可逆的方阵称为奇异或退化。方阵是奇异的当且仅当其行列式为 0。

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

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

从Mat::inv() 的文档中，有三种方法可供选择：

  • DECOMP_LU（默认）是 LU 分解。矩阵必须是非奇异的。
  • DECOMP_CHOLESKY 是对称正定义矩阵的 Cholesky 分解。这种类型在大矩阵上比 LU 快大约两倍。
  • DECOMP_SVD 是 SVD 分解。如果矩阵是奇异的或什至是非方形的，则计算伪逆。

从invert() 的文档中，大概是 Mat::inv() 内部使用的：

在 DECOMP_LU 方法的情况下，函数返回 src 行列式（ src 必须是正方形）。如果为 0，则矩阵不反转，dst 用零填充。

这与您看到的结果一致。

关于数学的笔记

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

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

在实验中，我发现归一化（如@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中定义，以防您有兴趣查看源代码。

于 2013-03-03T09:27:09.097 回答
3

对于随机矩阵R乘积R^T*R可能是奇异的。因此任何类型的 LU 分解都会过早停止，导致输出为零。

为了克服这一点，可以反转矩阵R^T*R+alpha*I。这I是单位矩阵，alpha- 一些正数。如果 alpha 接近于零且R^T*R不是奇异的，R^T*R+alpha*I则 的倒数接近 的倒数R^T*R。有关详细信息，请参阅Tikhonov 正则化

另一种情况是矩阵R^T*R不是奇异的而是病态的。大型非结构化矩阵的条件数可能很大，导致矩阵求逆的奇怪行为（LU 分解仅适用于相对较小的条件数）。

关于规范化

归一化矩阵将改善反演行为，因为它减少了条件数。

于 2013-03-05T10:26:52.223 回答
2

我得出的结论与@berak 相同。希望下面的实验对你有所帮助。

我用一个填充了一些随机值的矩阵尝试了你的代码（正态分布以 0 为中心，标准偏差为sigma。随着sigma增加，最终矩阵的值减少。

这是我使用的代码：

cv::Mat R(960,960, CV_32FC1);

double sigma[] = { 1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0, 1000000.0, 10000000.0 };
cv::Scalar mean[_countof(sigma)] = {0};
cv::Scalar stdv[_countof(sigma)] = {0};
for (int i = 0; i < _countof(sigma); i++)
{
    cv::randn(R, cv::Scalar::all(0.0), cv::Scalar::all(sigma[i]));
    cv::Mat temp2 = (R * R.t()).inv();

    cv::meanStdDev(temp2, mean[i], stdv[i]);
}

以下是增加 sigmas 的输出矩阵的均值和标准差：

sigma       mean        stddev
1.0         3.94e-004   1.32
10.0        1.25e-004   3.82e-002
100.0       3.32e-007   1.09e-004
1000.0      2.40e-009   2.23e-006
10000.0     9.82e-012   1.05e-008
100000.0    2.23e-013   1.73e-010
1000000.0   1.44e-015   2.88e-012
10000000.0  9.61e-017   2.77e-014

因此，您的解决方案是标准化您的输入矩阵，以便所有值都适合 [0;1] 或 [-1;1]。

于 2013-02-28T15:06:25.277 回答