EIGEN 3.3 SVD.SOLVE返回错误的值
Eigen 3.3 SVD.solve returning wrong values
从特征3.2.7更新到3.3.4后,我遇到了jacobisvd.solve的问题,在返回一个非常错误的结果。BDCSVD产生相同的结果。
问题可以使用以下代码复制:
#include <Eigen/Eigen>
#include <Eigen/SVD>
int main(int argc, const char * argv[]) {
// M is calculated beforehand and does not change with different Eigen versions
Eigen::Matrix<double, 12, 12> M = Eigen::Matrix<double, 12, 12>::Zero();
M(0,0) = 27; M(0,2) = 4.3625039999999995; M(0,3) = -9; M(0,5) = -2.1812519999999997; M(0,6) = -9;
M(1,1) = 27; M(1,2) = 3.2718720000000001; M(1,4) = -9; M(1,7) = -9; M(1,8) = -1.6359360000000001;
M(2,0) = 4.3625039999999995; M(2,1) = 3.2718720000000001; M(2,2) = 2.4780489612000003;
M(2,3) = -2.1812519999999997; M(2,5) = -0.82601632039999995; M(2,7) = -1.6359360000000001;
M(2,8) = -0.82601632039999995; M(3,0) = -9; M(3,2) = -2.1812519999999997; M(3,3) = 9;
M(4,1) = -9; M(4,4) = 9; M(5,0) = -2.1812519999999997; M(5,2) = -0.82601632039999995;
M(5,5) = 0.82601632039999995; M(6,0) = -9; M(6,6) = 9; M(7,1) = -9; M(7,2) = -1.6359360000000001;
M(7,7) = 9; M(8,1) = -1.6359360000000001; M(8,2) = -0.82601632039999995; M(8,8) = 0.82601632039999995;
Eigen::JacobiSVD<Eigen::MatrixXd> svd(M, Eigen::ComputeFullU);
Eigen::Matrix<double, 12, 12> ut = svd.matrixU();
Eigen::Matrix<double, 6, 10> l_6x10;
Eigen::Matrix<double, 4, 6> dv0;
Eigen::Matrix<double, 4, 6> dv1;
Eigen::Matrix<double, 4, 6> dv2;
for(int i = 0; i < 4; i++) {
int a = 0, b = 1;
for(int j = 0; j < 6; j++) {
dv0(i, j) = ut(3 * a + 0, 11 - i) - ut(3 * b + 0, 11 - i);
dv1(i, j) = ut(3 * a + 1, 11 - i) - ut(3 * b + 1, 11 - i);
dv2(i, j) = ut(3 * a + 2, 11 - i) - ut(3 * b + 2, 11 - i);
b++;
if (b > 3) {
a++;
b = a + 1;
}
}
}
for(int i = 0; i < 6; i++) {
l_6x10(i,0) = (dv0(0, i) * dv0(0, i) + dv1(0, i) * dv1(0, i) + dv2(0, i) * dv2(0, i));
l_6x10(i,1) = 2.0f * (dv0(0, i) * dv0(1, i) + dv1(0, i) * dv1(1, i) + dv2(0, i) * dv2(1, i));
l_6x10(i,2) = (dv0(1, i) * dv0(1, i) + dv1(1, i) * dv1(1, i) + dv2(1, i) * dv2(1, i));
l_6x10(i,3) = 2.0f * (dv0(0, i) * dv0(2, i) + dv1(0, i) * dv1(2, i) + dv2(0, i) * dv2(2, i));
l_6x10(i,4) = 2.0f * (dv0(1, i) * dv0(2, i) + dv1(1, i) * dv1(2, i) + dv2(1, i) * dv2(2, i));
l_6x10(i,5) = (dv0(2, i) * dv0(2, i) + dv1(2, i) * dv1(2, i) + dv2(2, i) * dv2(2, i));
l_6x10(i,6) = 2.0f * (dv0(0, i) * dv0(3, i) + dv1(0, i) * dv1(3, i) + dv2(0, i) * dv2(3, i));
l_6x10(i,7) = 2.0f * (dv0(1, i) * dv0(3, i) + dv1(1, i) * dv1(3, i) + dv2(1, i) * dv2(3, i));
l_6x10(i,8) = 2.0f * (dv0(2, i) * dv0(3, i) + dv1(2, i) * dv1(3, i) + dv2(2, i) * dv2(3, i));
l_6x10(i,9) = (dv0(3, i) * dv0(3, i) + dv1(3, i) * dv1(3, i) + dv2(3, i) * dv2(3, i));
}
Eigen::Matrix<double, 6, 4> L_6x4;
for(int i = 0; i < 6; i++) {
L_6x4(i, 0) = l_6x10(i, 0);
L_6x4(i, 1) = l_6x10(i, 1);
L_6x4(i, 2) = l_6x10(i, 3);
L_6x4(i, 3) = l_6x10(i, 6);
}
// L_6x4 has the following values on 3.2.7 (everything else is 0):
//
// L_6x4(0,2) = 1;
// L_6x4(1,0) = 1;
// L_6x4(1,1) = 1;
// L_6x4(5,0) = -1.137432760287006;
// L_6x4(5,2) = -1.1374327602870071;
// L_6x4(5,3) = -1.1374327602870049;
//
//
// on 3.3.4 it has the following slightly different values:
//
// L_6x4(0,2) = 1;
// L_6x4(1,0) = 1;
// L_6x4(1,1) = 1;
// L_6x4(5,0) = -1.1374327602869998;
// L_6x4(5,2) = -1.1374327602870271;
// L_6x4(5,3) = -1.1374327602869889;
// Rho is calculated beforehand and does not change with different Eigen versions
Eigen::Matrix<double, 6, 1> Rho;
Rho << 0.25, 0.140625, 0, 0.390625, 0.25, 0.140625;
Eigen::JacobiSVD<Eigen::MatrixXd> l_6x4(L_6x4, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::Vector4d B4 = l_6x4.solve(Rho);
// The slight difference in L_6x4 apparently causes an insane difference in the result of solve.
//
// Eigen 3.2.7: B4={0.056766494562302421, 0, 0, -0.064568070601816963}
// Eigen 3.3.4: B4={-4773392957911.6992, 0, 0, -4196637484493.9165}
return 0;
}
用特征egen 3.3.4计算的B4值使我认为这可能是一些奇怪的内存对齐问题,SDK中的错误或希望在此程序中有一个错误。
我尝试用-deigen_dont_align_statitational或-deigen_dont_vectorize和-deigen_disable_unaligned_array_array_assert对其进行编译。我还尝试了这些矩阵上的eigen :: dontalign。这些都不会对结果产生任何影响。
用NDK-R14B Clang C _静态测试在Android(4.4和5.1)上测试和Mac(MacOS Sierra 10.12.6),带有:" clang -std = C 14 -G -O0 -I/eigen_path main.cpp -o main"
这两个答案实际上非常相似,并且产生相似的残留物。这是因为您的矩阵具有两个零值,而第三个非常接近机器精度。因此,根据圆形错误,它被认为是0,在这种情况下,您可以在解决方案的其余3D子空间内获得最小的范围解决方案:
{0.056766494562302421, 0, 0, -0.064568070601816963}
和0.91的残差。否则,它被认为是有意义的,然后在较小的2D子空间(对应于两个零单数值)中获得最小的范围解决方案:
{-4773392957911.6992, 0, 0, -4196637484493.9165}
较小的残差为0.89。因此,第二个解决方案在某种意义上是更准确的,但是如果您更喜欢具有较小规范但剩余较高的一个解决方案,那么您可以调整阈值,以使第三个单数值被认为为零。这是通过SetThreshold完成的:
Eigen::JacobiSVD<Eigen::MatrixXd> l_6x4;
l_6x4.setThreshold(1e-14);
l_6x4.compute(L_6x4, Eigen::ComputeFullU | Eigen::ComputeFullV);
Eigen::Vector4d B4 = l_6x4.solve(Rho);
- 警告处理为错误这里有什么问题
- "error: no matching function for call to"构造函数错误
- boost::进程间消息队列引发错误
- C++,OpenCV,尝试显示图像时"OpenCV(4.3.0) Error: Assertion failed (size.width>0 && size.height>0)"此错误
- 有关插入适配器的错误。[错误]请求从 'back_insert_iterator<vector<>>' 类型转换为非标量类型
- QT在错误的班级中寻找空位
- vector.resize()中的分配错误
- 代码在main()中运行,但在函数中出现错误
- 释放错误后堆使用
- (C++)分析树以计算返回错误值的简单算术表达式
- Project Euler问题4的错误解决方案
- 我的字符计数代码计算错误.为什么
- 从"int*"强制转换为"unsigned int"会丢失精度错误
- 尝试导入pybind-opencv模块时出现libgtk错误
- CMake项目Boost库错误:Boost/config/compiler/gcc.hpp:165:10:致命错误:cs
- 在某些循环内使用vector.push_back时出现分段错误
- MSVC多行宏编译器错误
- 静态数据成员的问题-修复链接错误会导致编译器错误
- 为什么在运行时没有向我们提供有关分段错误的更多信息?
- EIGEN 3.3 SVD.SOLVE返回错误的值