Armadillo C++ doesn't find matrix inverse

Question

I'm using Armadillo & C++ and I'm trying to find the inverse of a Matrix, however, the inverse just returns the matrix itself.

It seems to me that there isn't any computation. Also, there are no errors thrown.

I am using the following header:

#include <armadillo>
using namespace std;
using namespace arma;

and I have been using Armadillo for a couple days and ran through several matrix manipulations that work properly.

Input:

mat A = randu<mat>(5,5);
A.print("A: ");
mat B = inv(A);
B.print("inv(A): ");

Output:

A: 
   0.0013   0.1741   0.9885   0.1662   0.8760
   0.1933   0.7105   0.1191   0.4508   0.9559
   0.5850   0.3040   0.0089   0.0571   0.5393
   0.3503   0.0914   0.5317   0.7833   0.4621
   0.8228   0.1473   0.6018   0.5199   0.8622
inv(A): 
   0.0013   0.1741   0.9885   0.1662   0.8760
   0.1933   0.7105   0.1191   0.4508   0.9559
   0.5850   0.3040   0.0089   0.0571   0.5393
   0.3503   0.0914   0.5317   0.7833   0.4621
   0.8228   0.1473   0.6018   0.5199   0.8622
Process finished with exit code 0

Question:

Why isn't inv(ofAMatrix) working, any hints or ideas? Thanks!

Answer 1

This works just fine with Armadillo 7.900.1 with Intel (R) MKL backend and Clang 5.0.

You should never take the inverse of a matrix unless it is absolutely necessary. Also you have to make sure that the inverse actually exists otherwise the algorithm will happily output garbage. If you want to compute the inverse of A to find x as in

x = A^-1 b

it is better to solve the linear system

A x = b

instead. These solvers are much faster and have way better convergence.

#include <armadillo>

int main()
{
  arma::mat A = { { 0.0013 , 0.1741 , 0.9885 , 0.1662 , 0.8760 } ,
                  { 0.1933 , 0.7105 , 0.1191 , 0.4508 , 0.9559 } ,
                  { 0.5850 , 0.3040 , 0.0089 , 0.0571 , 0.5393 } ,
                  { 0.3503 , 0.0914 , 0.5317 , 0.7833 , 0.4621 } ,
                  { 0.8228 , 0.1473 , 0.6018 , 0.5199 , 0.8622 } };
  A.print("A: ");
  arma::mat B = arma::inv(A);
  B.print("inv(A): ");
  arma::mat I = A*B;
  I.print("I: ");
}

Output:

A: 
   0.0013   0.1741   0.9885   0.1662   0.8760
   0.1933   0.7105   0.1191   0.4508   0.9559
   0.5850   0.3040   0.0089   0.0571   0.5393
   0.3503   0.0914   0.5317   0.7833   0.4621
   0.8228   0.1473   0.6018   0.5199   0.8622
inv(A): 
    0.4736   -1.7906    4.4377    2.2515   -2.4784
    2.9108   -3.1697   12.1159    7.7356  -11.1675
    2.5212   -2.8557    6.8074    4.7142   -6.1801
   -1.0317    0.9400   -2.3230    0.2413    1.3297
   -2.0869    3.6766   -9.6555   -6.9062    8.9447
I: 
   1.0000e+00   1.1340e-16  -1.8134e-15  -6.4918e-16  -4.8899e-17
   7.6334e-17   1.0000e+00  -9.1810e-16  -9.4668e-16   8.7907e-16
   2.5424e-16  -4.3981e-16   1.0000e+00   9.2981e-16  -2.0864e-15
   9.3036e-17  -2.6745e-17   7.5137e-16   1.0000e+00  -8.1372e-16
   4.3422e-16  -4.2293e-16   1.1321e-15   1.0687e-15   1.0000e+00

Answer 2

"Works for me" as they. Driving this from R and RcppArmadillo:

First, we read the matrix and use the generalized inverse from the MASS package:

R> M <- as.matrix(read.table(text="0.0013   0.1741   0.9885   0.1662   0.8760
   0.1933   0.7105   0.1191   0.4508   0.9559
   0.5850   0.3040   0.0089   0.0571   0.5393
   0.3503   0.0914   0.5317   0.7833   0.4621
   0.8228   0.1473   0.6018   0.5199   0.8622"))
M <- as.matrix(read.table(text="0.0013   0.1741   0.9885   0.1662   0.8760
+    0.1933   0.7105   0.1191   0.4508   0.9559
+    0.5850   0.3040   0.0089   0.0571   0.5393
+    0.3503   0.0914   0.5317   0.7833   0.4621
+    0.8228   0.1473   0.6018   0.5199   0.8622"))
R> M
         V1     V2     V3     V4     V5
[1,] 0.0013 0.1741 0.9885 0.1662 0.8760
[2,] 0.1933 0.7105 0.1191 0.4508 0.9559
[3,] 0.5850 0.3040 0.0089 0.0571 0.5393
[4,] 0.3503 0.0914 0.5317 0.7833 0.4621
[5,] 0.8228 0.1473 0.6018 0.5199 0.8622
R> MASS::ginv(M)
          [,1]      [,2]     [,3]      [,4]      [,5]
[1,]  0.473579 -1.790599  4.43767  2.251542  -2.47842
[2,]  2.910752 -3.169657 12.11587  7.735612 -11.16755
[3,]  2.521167 -2.855651  6.80743  4.714239  -6.18015
[4,] -1.031667  0.940028 -2.32302  0.241345   1.32967
[5,] -2.086858  3.676647 -9.65548 -6.906203   8.94472
R> 

The we use RcppArmadillo:

R> Rcpp::cppFunction("arma::mat armaInv(arma::mat x) { return arma::inv(x); }", depends="RcppArmadillo")
R> armaInv(M)
          [,1]      [,2]     [,3]      [,4]      [,5]
[1,]  0.473579 -1.790599  4.43767  2.251542  -2.47842
[2,]  2.910752 -3.169657 12.11587  7.735612 -11.16755
[3,]  2.521167 -2.855651  6.80743  4.714239  -6.18015
[4,] -1.031667  0.940028 -2.32302  0.241345   1.32967
[5,] -2.086858  3.676647 -9.65548 -6.906203   8.94472
R> 

Same answer both ways.