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!

Henri Menke · Answer 1 · 2017-06-13T01:08:08.003

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

Thanks for replying so quickly! 1. I don't know what MKL is. 2. Isn't it strange that I'm not even getting errors? Could it be a linker/compiler issue? Or maybe I'm missing a library? 3. Unfortunately, I do need inverse, and I'm aware that inverse is usually a bottleneck, thank you! — Nick Mason Kroeger, Jun 13 '17 at 01:04
Please use the [voting system](https://stackoverflow.com/help/why-vote) instead of posting "thanks". But anyway, you're welcome :) — Henri Menke, Jun 13 '17 at 01:07
MKL is Intel's implementation of the LAPACK library with special optimisations for their processors. This can give you quite some speedup. Also a lot of algorithms are implemented in parallel using OpenMP which gives you additional speed. — Henri Menke, Jun 13 '17 at 01:09

score 2 · Answer 2 · answered Jun 13 '17 at 00:44

"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.

That doesn't matter. Armadillo code is Armadillo code, I just called it from R, and they both use the same LAPACK backend. So I gave you an existence proof: no issue with Armadillo, likely local issue at your end you need to resolve. — Dirk Eddelbuettel, Jun 13 '17 at 01:03
Thanks, I'm aware it's a local issue, I'm sure it's not an Armadillo issue, but I'd appreciate any insight to what the issue might be, thanks! — Nick Mason Kroeger, Jun 13 '17 at 01:07
Start with a _complete and reproducible yet minimal example_. We don't have your code. You did not show how you compiled and linked etc pp. — Dirk Eddelbuettel, Jun 13 '17 at 01:07

Armadillo C++ doesn't find matrix inverse

2 Answers2