'Implementation of Matlab matrix inverse function in C++ with Eigen

So I need to re-write matrix right-handed division from Matlab to C++:

At = (xPow*yPow')/(yPow*yPow');

I mocked some matrices:

>> xPow*yPow'

ans =

0.0004    0.0040    0.0004    0.0004
0.0014    0.0263    0.0014    0.0014
0.0004    0.0012    0.0004    0.0004
0.0012    0.0053    0.0012    0.0012

and

>> yPow*yPow'

ans =

0.0001    0.0004    0.0001    0.0001
0.0004    0.0256    0.0004    0.0004
0.0001    0.0004    0.0001    0.0001
0.0001    0.0004    0.0001    0.0001

Matlab returns same result for (xPow*yPow')/(yPow*yPow') as well as xPow*yPow' * inv(yPow*yPow'). 

>> xPow*yPow' * inv(yPow*yPow')

ans =

36.1259    0.1127  -30.3163   -2.6999
40.6472    0.8810  -19.7529  -11.8430
-1.5578   -0.0182   12.1397   -7.0087
124.4466    0.0594 -130.0163   16.6710

>> At = (xPow*yPow')/(yPow*yPow')

At =

 36.1259    0.1127  -30.3163   -2.6999
 40.6472    0.8810  -19.7529  -11.8430
 -1.5578   -0.0182   12.1397   -7.0087
124.4466    0.0594 -130.0163   16.6710

Eigen library has function .inverse() so I thought I might use it to implement division of those matrices:

xyPowMult(0,0) = 0.0004;
xyPowMult(0,1) = 0.0040;
xyPowMult(0,2) = 0.0004;
xyPowMult(0,3) = 0.0004;
xyPowMult(1,0) = 0.0014;
xyPowMult(1,1) = 0.0263;
xyPowMult(1,2) = 0.0014;
xyPowMult(1,3) = 0.0014;
xyPowMult(2,0) = 0.0004;
xyPowMult(2,1) = 0.0012;
xyPowMult(2,2) = 0.0004;
xyPowMult(2,3) = 0.0004;
xyPowMult(3,0) = 0.0012;
xyPowMult(3,1) = 0.0053;
xyPowMult(3,2) = 0.0012;
xyPowMult(3,3) = 0.0012;

yyPowMult(0,0) = 0.0001;
yyPowMult(0,1) = 0.0004;
yyPowMult(0,2) = 0.0001;
yyPowMult(0,3) = 0.0001;
yyPowMult(1,0) = 0.0004;
yyPowMult(1,1) = 0.0256;
yyPowMult(1,2) = 0.0004;
yyPowMult(1,3) = 0.0004;
yyPowMult(2,0) = 0.0001;
yyPowMult(2,1) = 0.0004;
yyPowMult(2,2) = 0.0001;
yyPowMult(2,3) = 0.0001;
yyPowMult(3,0) = 0.0001;
yyPowMult(3,1) = 0.0004;
yyPowMult(3,2) = 0.0001;
yyPowMult(3,3) = 0.0001;

AtTemp = xyPowMult * yyPowMult.inverse();

cout << " x*y' " << endl;
cout << xyPowMult << endl << endl;

cout << " y*y' " << endl;
cout << yyPowMult << endl << endl;

cout << " x*y' * inv(y*y') " << endl;
cout << xyPowMult * yyPowMult.inverse() << endl << endl;

And the console outcome shows different result:

 x*y'
0.0004  0.004 0.0004 0.0004
0.0014 0.0263 0.0014 0.0014
0.0004 0.0012 0.0004 0.0004
0.0012 0.0053 0.0012 0.0012

 y*y'
0.0001 0.0004 0.0001 0.0001
0.0004 0.0256 0.0004 0.0004
0.0001 0.0004 0.0001 0.0001
0.0001 0.0004 0.0001 0.0001

 x*y' * inv(y*y')
-1.#IND -1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND -1.#IND
-1.#IND -1.#IND -1.#IND -1.#IND

So my questions are:

  1. Why the result is different?
  2. How can I implement matrix "slash" division with usage of Eigen?
  3. If not Eigen, then could you suggest any simple way?
c++matlabmatrixeigen


Solution 1:[1]

You have two issues. First, as Ilya Popov pointed out, y*y' is singular. Actually, so is x*y'. Second, matlab's \ operator actually solves a system of linear equations (Ax=B --> solves for x). Solving singular (or near singular) matrices using naive methods (e.g. inverse()) is not usually a good idea.

So to do the same with Eigen, you would set up the equations to solve and use the solution (intro). However, you would choose a specific algorithm based on prior knowledge. For example, 

A.fullPivLu().solve(b);

would give you A\b using LU.

Solution 2:[2]

As you can see, the last two rows of y*y' are identical. The matrix is singular (as is any matrix like v*v'). Such a matrix does not have an inverse. So you cannot expect any meaningful result here. Strange that Matlab produces anything. Are you sure your formula is correct?

Solution 3:[3]

Here:

#include <Eigen/Dense>

Eigen::MatrixXd inv_eigen (Eigen::MatrixXd m)
{
    // Matlab: inv(X) = X^(-1) = X\speye(size(X))
    return m.fullPivLu().solve(Eigen::MatrixXd::Identity(m.rows(), m.cols()));
}

Result:

X = 3×3

  1     0     2
 -1     5     0
  0     3    -9

Matlab inv(X)

  0.8824   -0.1176    0.1961
  0.1765    0.1765    0.0392
  0.0588    0.0588   -0.0980

Eigen inv(X)

  0.882353  -0.117647   0.196078
  0.176471   0.176471   0.0392157
  0.0588235  0.0588235 -0.0980392

Sources

This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.

Source: Stack Overflow

Solution Source
Solution 1 Avi Ginsburg
Solution 2 Ilya Popov
Solution 3 Cris Luengo

Related Questions