- Java 双重比较
- java - 比较器与 Apache BeanComparator
- Objective-C 完成 block 导致额外的方法调用?
- database - RESTful URI 是否应该公开数据库主键?
我使用 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
这里的 Eigen 解是正确的(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/
我正在使用 spyder - python。我想测试我的代码。我已关注 pip install spyder-unittest和 pip install pytest .我已经重新启动了内核并重新启动
我编写了一个简单的代码来匹配和选择/丢弃两个文件中的元素。此代码适用于包含 Id 对的 $file_in: Id1 Id2 Id1 Id3 Id1 Id4 Id3 Id4 Id3 Id5 Id3 Id
问题陈述:当 for 循环第二次执行时,我收到 Stale Element Exception。 描述: 我正在使用 for 循环来处理表元素。在第一次迭代中,它将在页面上搜索所需的元素。如果该元素在
我编写了一个简单的代码来匹配和选择/丢弃两个文件中的元素。此代码适用于包含 Id 对的 $file_in: Id1 Id2 Id1 Id3 Id1 Id4 Id3 Id4 Id3 Id5 Id3 Id
我正在尝试用 Java 创建康威的生命游戏,但似乎遇到了障碍。我创建了一个名为“Cell”的类,它包含一个 boolean 变量,该变量本质上确定细胞是活的还是死的,以及在需要时杀死或创建细胞的方法。
我有一张 table userid points expires我只想对所有匹配特定 userid 的点求和和 CURTIME() < expires . 所以我做到了: select *, su
题目地址:https://leetcode.com/problems/solve-the-equation/description/ 题目描述: Solve a given equation an
自从工作流构建问世以来,我一直在使用 TFS 构建。我从未见过使用以下任何一项的需要或愿望: Copy Build output to the server 什么服务器?什么目录? `将构建输出复制到
我有这个矩阵 a = {{2, -2, -4}, {-2, 5, -2}, {-4, -2, 2}} 然后我解出了一个缺少一项的方程。方程的形式为逆[p].a.p == q其中 p 是缺少条目 (x5
下面的代码解决了 hanoi 使用预定义函数 moveLOD、swapLOI 和 swapLID 返回 Action 列表的问题。 MoveLOD:将 1 个圆盘从第一个位置移动到三元组第三个位置中的
基本上,我通过从文本文件中读取一组整数来实现 AVL 树,然后使用 add() 方法填充树。另外,程序应该按顺序打印整数集。 当我运行该程序时,会弹出 StackOverflowError。我认为这个
我在 html 中有一个 span 标签,这个 可以在页面上多次出现。 我只想选择第一个跨度标签。我就是这样做的 var $myvar= $(".mydiv .indiv .myspan:first"
来自 sympy solve() 的解决方案是否以某种方式排序?它是从最小到最大的解决方案吗? 如何强制解决方案的非负性? 在我的问题中,我需要独特的最小正解决方案。我感谢所有的帮助 最佳答案 它们不
我正在尝试打开一个文件,但我收到了: 该进程无法访问该文件,因为它正被另一个进程使用。该文件是一个 XML 文档。谁能帮忙? string activeDirectory = @"X:\SubGr
如何解决此内存泄漏...我什至在最后释放它,如图片所示,但它仍然存在。在 if 语句几乎 10-15 条件下,它像给定的代码一样使用......但最后我发布了它。 LoginResponse *res
我为不同的日子编写了以下代餐代码,但我每天都吃同样的餐。我想隔天吃“肉”和“素食”食物。 my dataframe is as follows: id name
这会创建一个点列表,sympy solve() 方法应该返回 x 的值。相反,它返回另一个方程式,我不确定为什么。 ogrid() 和 ravel() 正在绘图中创建一个点列表,这是在 Matplot
只是好奇它是否使用高斯消去法或其他等效方法? 最佳答案 来自numpy docs : solve is a wrapper for the LAPACK routines dgesv and zges
首先, solve_poly_system( seq, *gens, **args), 有人确切知道solve_poly_system的参数是什么意思吗? 我有系统, rd = λk ua = λk
我正在尝试使用 sympy 的 solve 命令求解方程,但我的输出是一个空列表 [ ]。我认为这可能发生的唯一原因是没有解决方案,但我怀疑这就是原因。任何人都知道为什么我没有得到答案?谢谢! fro
我是一名优秀的程序员,十分优秀!