我将Armadillo 5.200与Visual Studio 2013结合使用来解决线性方程组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给我的解决方案，并且对我正在执行的计算具有逻辑意义)。为什么 Armadillo 给我不正确的解决方案？
最佳答案
解决方案Armadillo是正确的，因为您解决了预定过多的系统。
因此，通过不同方法获得的解决方案可能会有所不同。
下面的代码演示了如何使用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;
}
与往常一样，在编译时指示大小矩阵可以减少内存分配的成本，因此可以使用它来加快计算速度。
关于c++ - Armadillo的Solve(A，b)返回与Matlab，Eigen不同的答案，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/31705168/