犰狳的求解（A，b）从 Matlab，Eigen 返回不同的答案
Posted
技术标签:
【中文标题】犰狳的求解（A，b）从 Matlab，Eigen 返回不同的答案【英文标题】：Armadillo's solve(A, b) returning different answer from Matlab, Eigen 【发布时间】：2015-07-29 15:44:05 【问题描述】：我正在使用带有 Visual Studio 2013 的 Armadillo 5.200 来求解线性方程组 xA=b。我的代码通过求解 x=(A'\b')' 来计算这个，我正在使用 Armadillo 的 solve() 函数来求解线性方程组。
我的问题是返回给我的解决方案不正确 - 它与相同问题的 Eigen 和 Matlab 解决方案不同，它们匹配。我不只是首先使用 Eigen 的原因是它执行得太慢（我必须进行数百万次计算），我想看看使用 Armadillo 是否能加快它的速度。我现在也很好奇为什么会发生这个错误。
这是一个说明问题的示例代码：
#include "stdafx.h"
#include <iostream>
#include <armadillo>
#include <Eigen/Dense>
// Matrix operation to find world coordinates.
arma::mat mrdivide(arma::mat A, arma::mat B)
arma::mat A_t = A.t();
arma::mat B_t = B.t();
return solve(A_t, B_t).t();
Eigen::MatrixXd mrdivide(Eigen::MatrixXd A, Eigen::MatrixXd B)
Eigen::MatrixXd A_t = A.transpose();
Eigen::MatrixXd B_t = B.transpose();
return ((A_t).colPivHouseholderQr().solve(B_t)).transpose();
int main(int argc, char* argv[])
Eigen::Matrix<double, 4, 3> M_eig;
M_eig << 761.544, 0, 0,
0, 761.544, 0,
639.5, 399.5, 1.0,
3.762513283904080e+06, 1.824431013104484e+06, 9.837714402800992e+03;
arma::mat M_arma;
M_arma << M_eig(0, 0) << M_eig(0, 1) << M_eig(0, 2) << arma::endr
<< M_eig(1, 0) << M_eig(1, 1) << M_eig(1, 2) << arma::endr
<< M_eig(2, 0) << M_eig(2, 1) << M_eig(2, 2) << arma::endr
<< M_eig(3, 0) << M_eig(3, 1) << M_eig(3, 2) << arma::endr;
Eigen::Matrix<double, 1, 3> pixelCoords_eig;
pixelCoords_eig << 457, 520, 1;
arma::mat pixelCoords_arma;
pixelCoords_arma << 457 << 520 << 1 << arma::endr;
Eigen::Matrix<double, 1, 4> worldCoords_eig;
arma::mat worldCoords_arma;
worldCoords_arma = mrdivide(M_arma, pixelCoords_arma);
worldCoords_eig = mrdivide(M_eig, pixelCoords_eig);
std::cout << "world coords arma: " << worldCoords_arma << std::endl;
std::cout << "world coords eig: " << worldCoords_eig << std::endl;
我的输出是：
world coords arma: 5.3599e-002 4.242e-001 1.3120e-001 8.8313e-005
world coors eig: .0978826 .439301 0 0.00010165
这里的特征解是正确的（Matlab 给我的，并且对我正在执行的计算具有逻辑意义）。为什么犰狳会给我错误的解决方案？
【问题讨论】：
【参考方案1】：解决方案犰狳是正确的，因为您解决了超定系统。 因此不同方法得到的解可能不同。 下面的代码演示了如何使用 Armadilio 库，得到类似 Matlab 的结果。 为简单起见，我没有完全实现算法，也没有检查参数函数的正确性，因此交换行是人为的（并不是真正必要的）。
# include <iostream>
# include <armadillo>
# include <algorithm>
# include <array>
arma::mat qr_solve2(const arma::mat &A,const arma::mat &B);
arma::mat mrdivide(arma::mat A, arma::mat B)
return solve(A.t(), B.t());
int main()
std::array<double,12> arr =
761.544 , 0 ,639.5,
3.762513283904080e+06,0 , 761.544,
399.5,1.824431013104484e+06,0,
0, 1.0, 9.837714402800992e+03
;
arma::mat A(4,3);
std::copy(arr.begin(),arr.end(),A.begin());
arma::mat B;
B << 457 << 520 << 1 << arma::endr;
arma::mat V = mrdivide(A,B);
std::cout<<V<<std::endl;
std::cout<<A.t()*V<<std::endl;
arma::mat Z = qr_solve2(A.t(),B.t());
std::cout<<Z<<std::endl;
std::cout<<A.t()*Z<<std::endl;
return 0;
arma::mat qr_solve2(const arma::mat &A,const arma::mat &B)
arma::mat Q, R;
arma::qr(Q,R,A);
unsigned int s = R.n_rows-1;
arma::mat R_ = R( arma::span(0,s), arma::span(0,s-1) ) ;
R_ = arma::join_horiz(R_,R.col(s+1));
arma::mat newB = Q.t()*B;
arma::mat X(s+1,1);
for (int i = s; i >= 0; i--)
X[i] = newB[i];
for (int j = s; j != i; j--)
X[i] = X[i] - R_(i,j) * X[j];
X[i] = X[i]/R_(i,i);
arma::mat res = X( arma::span(0,s-1), arma::span(0,0) ) ;
res = arma::join_vert(res,arma::mat(1,1, arma::fill::zeros));
res = arma::join_vert(res,X.row(s));
return res;
与往常一样，在编译时指示大小矩阵会降低内存分配的成本，因此可以使用它来加快计算速度。
【讨论】：
【参考方案2】：关于 Eigen 的速度，您可以通过以这种方式移除所有堆分配来显着提升（几乎快一个数量级）：
Matrix<double, 1, 4> mrdivide(const Matrix<double, 4, 3> &A,
const Matrix<double, 1, 3> &B)
return A.transpose().colPivHouseholderQr().solve(B.transpose());
您还可以通过对 M_eig 使用行优先矩阵使 A.transpose() 是列优先矩阵来多节省 10%：
Matrix<double, 4, 3,RowMajor> M_eig;
Matrix<double, 1, 4> mrdivide(const Matrix<double, 4, 3, RowMajor> &A,
const Matrix<double, 1, 3> &B);
最后，由于您的问题在数值上条件良好，您还可以使用基于 Cholesky 的求解器来获得额外的 x2 加速（在这种情况下，保留 M_eig 的默认存储）：
Matrix<double, 1, 4> mrdivide(const Matrix<double, 4, 3> &A, const Matrix<double, 1, 3> &B)
return (A*A.transpose()).ldlt().solve(A*B.transpose());
为了完整起见，这里是一个实现 QR 和 LDLT 解决方案的独立示例：
#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
Matrix<double, 1, 4> mrdivide_qr(const Matrix<double, 4, 3> &A,
const Matrix<double, 1, 3> &B)
return A.transpose().colPivHouseholderQr().solve(B.transpose());
Matrix<double, 1, 4> mrdivide_ldlt(const Matrix<double, 4, 3> &A, const Matrix<double, 1, 3> &B)
return (A*A.transpose()).ldlt().solve(A*B.transpose());
int main(int argc, char* argv[])
Matrix<double, 4, 3> M_eig;
M_eig << 761.544, 0, 0,
0, 761.544, 0,
639.5, 399.5, 1.0,
3.762513283904080e+06, 1.824431013104484e+06, 9.837714402800992e+03;
Matrix<double, 1, 3> pixelCoords_eig;
pixelCoords_eig << 457, 520, 1;
Matrix<double, 1, 4> worldCoords_eig;
worldCoords_eig = mrdivide_qr(M_eig, pixelCoords_eig);
std::cout << "world coords using QR: " << worldCoords_eig << std::endl;
worldCoords_eig = mrdivide_ldlt(M_eig, pixelCoords_eig);
std::cout << "world coords using LDLT: " << worldCoords_eig << std::endl;
mrdivide_ldlt 与mrdivide_qr 一样快，后者本身比问题的mrdivide 函数快得多。
【讨论】：
使用你的第一个和第二个建议（M_eig 的换位不会保持上三角，所以我不认为我可以使用 ldlt() 不幸的是）我得到一个“YOU_MIXED_MATRICES_OF_DIFFERENT_SIZES”编译错误。回到使用MatrixXd 而不是 Matrix<double, 4, 3> 和 Matrix<double, 1, 3> 可以消除这个错误，但这让我回到了原来的速度问题。
想通了——返回行应该是return A.transpose().colPivHouseholderQr().solve(B.transpose()).transpose();。这确实加快了速度，谢谢！
hm，最后一个转置应该不需要，因为在赋值中转置向量是隐式的。另外，对于 ldlt，不是你通过 A*A.transpose() 和 A*B.transpose() 来获得最小范数解，而不是原始矩阵！
A.transpose().colPivHouseholderQr().solve(B.transpose()) 的解决方案是一个 4x1 矩阵，但是，由于我正在为 1x4 矩阵求解 (A'\b')'，我想我最终确实需要最终的转置。 LDLT 也不适合我——我认为我不能使用它，因为我将用不是上三角或下三角的矩阵替换 M_eig——但就目前而言，我不断收到缓冲区溢出错误。
嗯，你一定在做一些奇怪的事情。我用一个独立的示例编辑了答案。以上是关于犰狳的求解（A，b）从 Matlab，Eigen 返回不同的答案的主要内容，如果未能解决你的问题，请参考以下文章