Матричная обратная функция библиотеки Eigen возвращает NaN

Question

Это новый созданный пример с той же проблемой. 2 матрицы "firstMultiplyMat", "a" имеют одинаковые значения, но указаны разными способами. Тем не менее, инверсия каждого из них различна. Невозможно получить обратное значение для «firstMultiplyMat», тогда как «a» имеет правильную рассчитанную обратную матрицу.

#include<iostream> 
#include<stdio.h> 
#include <Eigen/Core> 
#include <Eigen/Dense> 

using namespace std; 
using namespace Eigen; 
using Eigen::MatrixXd; 

struct MemberLocation{ 
int id; 
double x; 
double y; 
}; 

int main(){ 

//Attributes 
int M = 2 ;//anchor 
int N = 2 ;//nonanchor 
int NA = 4; 
int p = 0; 
int numOfIterations =2; 
MemberLocation *member1 = new MemberLocation[M]; 
MemberLocation *member2 = new MemberLocation[N]; 
MatrixXd neighbourdistanceEst(M,N); 
MatrixXd neighbourdistanceDeriv(M,NA); 
MatrixXd distanceDerivTranspose(NA,M); 
MatrixXd firstMultiplyMat(NA,NA); 
MatrixXd firstMultiplyMatInverse(NA,NA); 

//Structs initilaization 
member1[0].id = 5; 
member1[0].x = 9.31301; 
member1[0].y = 19.3955; 

member1[1].id = 2; 
member1[1].x = 46.6279; 
member1[1].y = 0.00571905; 

member2[0].id = 4; 
member2[0].x = 11.7718; 
member2[0].y = 7.99507; 

member2[1].id = 6; 
member2[1].x = 23.6158; 
member2[1].y = 3.80408; 

// Filling "neighbourdistanceDeriv" matrix 
for (int i = 1 ; i < numOfIterations ; i++){ 

for (int j = 0 ; j < M ; j++){ 
for (int k = 0 ; k < N ; k++){ 

    int id1 = member2[k].id; 
    int id2 = member1[j].id; 
    double xDiff = member2[k].x - member1[j].x; 
    double yDiff = member2[k].y - member1[j].y; 
    double distance = pow(xDiff,2) + pow(yDiff,2);   

    neighbourdistanceEst(j,k) = sqrt(distance); 

    neighbourdistanceDeriv(j,p) = xDiff/neighbourdistanceEst(j,k); 
    neighbourdistanceDeriv(j,p+1) = yDiff/neighbourdistanceEst(j,k); 

    p+=2; 

} 
p = 0; 
} 
} 

// operations on "neighbourdistanceDeriv" matrix 
distanceDerivTranspose = neighbourdistanceDeriv.transpose(); 
firstMultiplyMat = distanceDerivTranspose*neighbourdistanceDeriv; 
firstMultiplyMatInverse = firstMultiplyMat.inverse(); 

// printing "firstMultiplyMat" matrix and its inverse 
std::cout << "firstMultiplyMat:\n" << firstMultiplyMat << std::endl; 
std::cout << "inverse of firstMultiplyMat:\n" << firstMultiplyMatInverse << std::endl<< std::endl; 

// fixed array with same values obtained in "neighbourdistanceDeriv" matrix at runtime 
MatrixXd a(4,4); 
MatrixXd b(4,4); 
a(0,0) = 0.994532; 
a(0,1) = -0.423853; 
a(0,2)= 1.10423 ; 
a(0,3) = -0.314096; 
a(1,0) = -0.423853 ; 
a(1,1) = 1.00547 ; 
a(1,2) = -0.881237 ; 
a(1,3) = 0.756726 ; 
a(2,0) = 1.10423 ; 
a(2,1) = -0.881237; 
a(2,2) = 1.43045; 
a(2,3) = -0.658827; 
a(3,0) = -0.314096; 
a(3,1) = 0.756726 ; 
a(3,2) = -0.658827; 
a(3,3) = 0.569549 ; 
b = a.inverse(); 
std::cout << "matrix a :\n" << a <<std::endl; 
std::cout << "inverse of matrix a :\n" << b << std::endl; 
//std::cout<<"hi \n"; 
return 0; 
} 

---

Выход:

firstMultiplyMat: 
0.994534 -0.423857 1.10423 -0.314099 
-0.423857 1.00547 -0.881237 0.756725 
1.10423 -0.881237 1.43045 -0.658827 
-0.314099 0.756725 -0.658827 0.569548 
inverse of firstMultiplyMat: 
-nan -nan -nan -nan 
-nan -nan -nan -nan 
-nan -nan -nan -nan 
-nan -inf -inf inf 

matrix a : 
0.994532 -0.423853 1.10423 -0.314096 
-0.423853 1.00547 -0.881237 0.756726 
1.10423 -0.881237 1.43045 -0.658827 
-0.314096 0.756726 -0.658827 0.569549 
inverse of matrix a : 
-119414  -1681.81  132361  89489.1 
-1681.81  942762  10176.4 -1.24175e+06 
132361  10176.4  -146638  -110149 
89489.1 -1.24175e+06  -110149 1.57177e+06 

Answer 1

Матрица firstMultiplyMat само по себе является тонкой памяти стрелке; если вы введете значения в них, которые вы положили в a, он будет инвертироваться правильно. Ваша проблема заключается в значениях, которые находятся внутри firstMultiplyMat, которые не являются точно теми же значениями, что и печатаются, потому что напечатанное усечено.

По существу, ваша матрица не может быть обратимой. Вы можете обойти это, просто не вычисляя обратное и вместо того, чтобы с помощью QR-решатель, например:

HouseholderQR<MatrixXd> qr(A); 
x = qr.solve(b); // computes A^-1 * b 

Более подробная информация о HouseholderQR доступна here