为什么这个 matlab 和 C++ 代码会产生不同的结果

Question

+我正在尝试在C++中实现 matlab 算法。

这是 matlab 代码：

p = 3;
K = [3 4 5; 4 5 6; 7 8 9];
e = ones(p,1);
K2 = K - (1/p)*K*ones(p) - 1/p + (1/p^2)*sum(K(:))
[V_K,D_K] = eig(K2);

虽然这是使用 OpenCV 的类似C++代码：

float data[] = {3, 4, 5,
                4, 5, 6,
                7, 8, 9};
cv::Mat m(3, 3, CV_32F, data);
float p = K.rows;
cv::Mat CK = K - (1/p)*K*cv::Mat::ones(p,p,CV_32F) - 1/p + (1/std::pow(p,2))*cv::sum(K)[0];
cv::Mat eigenvalues(1,p,CK.type()), eigenvectors(p,p,CK.type());
cv::eigen(CK,eigenvalues,eigenvectors);

矩阵代码打印：

CK =
4.3333    5.3333    6.3333
4.3333    5.3333    6.3333
4.3333    5.3333    6.3333
0.5774    0.6100   -0.1960
0.5774   -0.7604   -0.6799
0.5774    0.2230    0.7066
16.0000         0         0
      0   -0.0000         0
      0         0    0.0000

虽然C++代码产生：

CK=[4.3333335, 5.3333335, 6.3333335;
 4.3333335, 5.3333335, 6.3333335;
 4.333333, 5.333333, 6.333333]
eigenvectors=[0.53452265, 0.56521076, 0.62834883;
 -0.41672006, 0.8230716, -0.38587254;
 0.73527533, 0.05558794, -0.67548501]
eigenvalues=[17.417906;
 -0.33612049;
 -1.0817847]

如您所见，这些值完全不同（即使是CK的值！为什么会发生这种情况，我该如何避免？

请注意，我不完全确定我的C++实现是否正确！

我发现这个问题和这个问题相关，但它们似乎与细微的差异有关，而这里的错误是巨大的！

更新：

我试图遵循评论和答案中的建议。不幸的是，提出的解决方案都没有解决问题。首先，我尝试float精确地使用Eigen库。这是使用 Eigen::Map 结构的代码，如下所述：

//in order to map OpenCV matrices to Eigen, they need to be continous
assert(CK.isContinuous());
Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> CKEigenMapped (CK.ptr<float>(), CK.rows, CK.cols);
Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> CKEigen = CKEigenMapped;
Eigen::EigenSolver<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> es (CKEigen,true);
std::cout<<"Eigenvalues:"<<std::endl<< es.eigenvalues() << std::endl;
std::cout<<"Eigenvectors:"<<std::endl<< es.eigenvectors() << std::endl;

然后我尝试通过CK.convertTo(CK, CV_64F)从float转换为double：

//Double precision
CK.convertTo(CK, CV_64F);
Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> CKEigenMappedD (CK.ptr<double>(), CK.rows, CK.cols);
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> CKEigenD = CKEigenMappedD;
Eigen::EigenSolver<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> esD (CKEigenD,true);
std::cout<<"Eigenvalues:"<<std::endl<< esD.eigenvalues() << std::endl;
std::cout<<"Eigenvectors:"<<std::endl<< esD.eigenvectors() << std::endl;

最后，我尝试使用cv2eigen函数（我认为Eigen::Map可能是错误的），如下所述：

//Double precision, cv2eigen
Eigen::MatrixXd X=Eigen::MatrixXd(CK.rows,CK.cols);
cv2eigen(CK,X);
Eigen::EigenSolver<Eigen::MatrixXd> esDD (X,true);
std::cout<<"Eigenvalues:"<<std::endl<< esDD.eigenvalues() << std::endl;
std::cout<<"Eigenvectors:"<<std::endl<< esDD.eigenvectors() << std::endl;

这些是对应于前 3 个解决方案的结果：

Eigenvalues:
(-4.17233e-07,0)
          (16,0)
(-3.37175e-07,0)
Eigenvectors:
(-0.885296,0)   (0.57735,0)  (-0.88566,0)
 (0.328824,0)   (0.57735,0)  (0.277518,0)
 (0.328824,0)   (0.57735,0)  (0.372278,0)
Eigenvalues:
          (16,0)
  (8.9407e-08,0)
(-1.88417e-16,0)
Eigenvectors:
  (0.57735,0)  (0.480589,0)  (0.408248,0)
  (0.57735,0)  (0.480589,0) (-0.816497,0)
  (0.57735,0) (-0.733531,0)  (0.408248,0)
Eigenvalues:
          (16,0)
  (8.9407e-08,0)
(-1.88417e-16,0)
Eigenvectors:
  (0.57735,0)  (0.480589,0)  (0.408248,0)
  (0.57735,0)  (0.480589,0) (-0.816497,0)
  (0.57735,0) (-0.733531,0)  (0.408248,0)

正如您所注意到的：

1. 它们都与 Matlab 结果都不对应 （ ：'（ ）
2. 使用 double 和 float 是有区别的
3. 使用 Eigen::Map 和 cv2eigen 之间没有区别

请注意，我不是Eigen专家，我可能以错误的方式使用Eigen::EigenSolver。

更新 2：

这开始一团糟了！这是使用Amradillo的代码。请注意，A具有相同的 K2 值（C++ 中CK）：

arma::mat A(3,3);
A   << 4.333333333333333 << 5.333333333333333 << 6.333333333333333 <<arma::endr
    << 4.333333333333333 << 5.333333333333333 << 6.333333333333333 <<arma::endr
    << 4.333333333333333 << 5.333333333333333 << 6.333333333333333 <<arma::endr;
arma::cx_vec eigval;
arma::cx_mat eigvec;
eig_gen(eigval,eigvec,A);
std::cout<<"eigval="<<std::endl<<eigval<<std::endl<<"eigvec="<<std::endl<<eigvec<<std::endl;

这些是打印值：

eigval=
    (+1.600e+01,+0.000e+00)
    (-4.010e-17,+3.435e-16)
    (-4.010e-17,-3.435e-16)
eigvec=
    (+5.774e-01,+0.000e+00)    (-5.836e-02,+3.338e-01)    (-5.836e-02,-3.338e-01)
    (+5.774e-01,+0.000e+00)    (+7.174e-01,+0.000e+00)    (+7.174e-01,-0.000e+00)
    (+5.774e-01,+0.000e+00)    (-5.642e-01,-2.284e-01)    (-5.642e-01,+2.284e-01)

说真的，所有这些库有什么问题？他们甚至彼此不完全同意！

Answer 1

cv::eigen假设输入矩阵是对称的，而你的不是。 这就是为什么存在差异的原因。

我相信 openCV 不支持非对称矩阵的特征向量，您可能需要使用另一个库。

更新：PCA（主成分分析）是一种特征向量分解，所以你可以这样做，但最好是使用一些特定的数学库，如 EIGEN 或 ARMADILLO。