Mat.inv() yielding all zeroes in opencv

Question

I have the following code :

    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;

The Transpose works correctly, so does the matrix multiplication. Thus the Mat temp1 has a size of 960x960. However, when I do temp2 =temp1.inv(), I recieve all zeroes in temp2. I mean zeroes is all of the 960x960 cells. Also, R is of type CV_32FC1 only. So it is probably not a datatype issue. I cannot understand the issue here. I googled so much. Can you please help.

EDIT

I am copying below the gdb output for the Mat::inv() function. I am having a hard time figuring it all out, but if someone is more familiar with OpenCV, maybe it will be of help :)

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

Most likely, the determinant is zero.

From Wikipedia:

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.

You can display the determinant like so...

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

From the documentation for Mat::inv(), there are three methods to choose from:

• DECOMP_LU (default) is the LU decomposition. The matrix must be non-singular.
• DECOMP_CHOLESKY is the Cholesky decomposition for symmetrical positively defined matrices only. This type is about twice faster than LU on big matrices.
• DECOMP_SVD is the SVD decomposition. If the matrix is singular or even non-square, the pseudo inversion is computed.

From the documentation for invert(), which is presumably used internally by 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.

This agrees with the results that you are seeing.

---

notes about the math

I'm no mathematician, but I get the impression that inverting a matrix can be a messy business -- all the more so if your matrix is very large. In fact, it may be true that these inverses exist in principle, but are practically impossible to calculate with any accuracy. In running some experiments with your code, I found that in many cases I would get determinants that were not exactly zero, but were very close to zero -- perhaps indicating that numerical precision may be the limiting factor. I tried specifying the matrices using 64-bit values instead of 32, and got different, but not necessarily better answers.

It may be useful to recognize that, based on the way you are calculating the temp1 matrix, it will always be symmetric. The DECOMP_CHOLESKY method is specifically designed to work on symmetric positive definite matrices, so using that might provide some advantages.

Experimentally, I found that normalizing (as suggested by @cedrou) makes it more likely that the inverse function returns a non-zero matrix (with DECOMP_LU but not with DECOMP_CHOLESKY). However, based on my guesses of how you might be initializing the R matrix, the resulting matrices never seemed to satisfy the definition of an inverse: A*inverse(A)=Identity. But you don't necessarily care about that -- which is perhaps why the SVD method computes a pseudo inverse.

Finally, it seems that this deeper question of why inversion is failing might be a math question rather than a programming question. Based on that I did some searching on the math site, and it turns out that someone has already asked how to do this very thing: https://math.stackexchange.com/questions/182662

---

notes about debugging

Based on your debug trace, I am inclined to think that the part that you are interested in was compiled into a non-traceable library and skipped over when you ran step. In other words, that mysterious blank line after your first step represents the part where it actually ran the inv() function. After that it is assigning the result to temp2 and destructing temporary objects. So your debug trace doesn't tell us anything about what is happening inside of 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);

I ran a debugger on this myself and was able to trace through the inner call to invert() and watch it decide to fail based on an internal analysis of the matrix (determining that it was not invertible) -- and therefore return a matrix filled with zeros, matching what you have reported.

The invert() function is defined in cxlapack.cpp, in case you are interested in taking a look at the source code.