0

In my algorithm I employ a sparse matrix inverse operation and I solve it by using the A*x=b method using the QR decomposition method. On Matlab the QR operation runs fine.

However, when I tried to convert the code to C++ using the Eigen library, I didn't get the same answer. In some cases, there is a shift of some value along each element of vector x compared to the results in Matlab. This value which causes the shift is however constant through all the elements in the vector.

A glimpse of what I do:

Eigen::SparseMatrix<float> A(m, n);
Eigen::VectorXf b;
Eigen::SparseQR<Eigen::SparseMatrix<float>, Eigen::COLAMDOrdering<int>> solver; 
solver.compute(A);
Eigen::VectorXf x = solver.solve(b);

x is my final vector which contains the result of A.inverse()*b isn't it?

Additionally, I tried to solve it as a full matrix but still produced different answers on C++ compared to Matlab.

Did anyone here face similar problem ? If yes, any help or pointers are welcome. On the other hand, if there is something wrong with my understanding, any correction is also appreciated.

Thanks.

user1389578
  • 21
  • 1
  • 3
  • Please, use double not float if you want accuracy, and you can check the accuracy of the returned solution by printing the residual `(A*x-b).norm()/b.norm()`. If it is small enough, then you got a solution. Note that if the problem is not full rank, then the solution of your problem is not unique. – ggael Mar 27 '15 at 09:06
  • Hi, thanks for the hint. Without setting the Pivot threshold, it was not a full rank but I manually set a low threshold value and now I get a full rank matrix. But again, I would expect certain level of accuracy on the first 2-3 decimal places using float. So you say I need double to get perfect accuracy wrt Matlab ? Because I tried and it gives different results and not accurate. – user1389578 Mar 30 '15 at 08:43
  • What is the value of `(A*x-b).norm()/b.norm()` when using double? – ggael Apr 01 '15 at 07:17
  • Hi, I am sorry I forgot to update ! I calculated the residul error with float itself and I got desired results. I used the value as pivoting threshold and got full rank and then I get desired results. Thanks a lot for the pointer. – user1389578 Apr 02 '15 at 08:05

0 Answers0