使用特征3线性系统解算器得出错误结果

Question

我正在用Eigen 3将一个MATLAB代码移植到C，我决定构造我的线性解算器，而不是从矩阵对象调用它，这样我就可以重用它了。不幸的是，它没有产生预期的结果。经过几次测试后，我将问题追溯到了似乎是线性解算器对象的地方，如以下相关代码所示：

MatrixXcd M(6,6), Y(6,6), dvalor(6,6);
// Initialization of the matrices...

FullPivHouseholderQR<MatrixXcd> mdivide;

auto result = mldivide.compute(Y).solve(M).eval() *
    dvalor * mldivide.compute(M).solve(Y).eval();

std::cout << "result:" << std::endl
    << result << std::endl << std::endl;

std::cout << "using fullPivHouseholderQr:" << std::endl
    << Y.fullPivHouseholderQr().solve(M) * dvalor *
    M.fullPivHouseholderQr().solve(Y) << std::endl << std::endl;

std::cout << "using jacobi:" << std::endl
    << Y.jacobiSvd(ComputeThinU | ComputeThinV).solve(M) * dvalor *
    M.jacobiSvd(ComputeThinU | ComputeThinV).solve(Y) << std::endl << std::endl;

std::cout << "using inverse:" << std::endl
    << Y.inverse() * M * dvalor * M.inverse() * Y << std::endl << std::endl;

结果是：

result:
(0.564196,0.00606298)       (-0.15282,-0.00179293)      (-0.16564,-0.000220726)      (-0.16564,-0.000220726)       (-0.15282,-0.00179293)       (-0.179235,0.00541725)
(2.73563e+184,-8.89017e+185)              (1,1.12229e-10)   (-6.04297e-11,5.62259e-12)   (-6.02256e-11,2.60221e-12)  (-5.33448e-11,-3.05694e-11)   (-6.55427e-11,8.17723e-12)
(2.11057e+184,-8.81209e+185)   (-5.53615e-11,4.54731e-12)              (1,1.11384e-10)   (-5.99088e-11,9.29495e-12)   (-5.51744e-11,1.52693e-12)  (-6.21597e-11,-3.11176e-11)
(2.11057e+184,-8.81209e+185)      (-0.171319,-0.00113681)      (0.134093,-2.05541e-05)       (0.865907,2.05542e-05)        (0.171319,0.00113681)     (-0.140389,-0.000216684)
(2.73563e+184,-8.89017e+185)        (0.218904,0.00143522)      (-0.171337,3.98397e-05)      (0.171337,-3.98397e-05)       (0.781096,-0.00143522)       (0.179383,0.000262654)
(1.4668e+184,-8.78677e+185)   (-5.52732e-11,5.99039e-12)   (-5.69133e-11,-3.2223e-11)   (-5.97328e-11,9.63533e-12)   (-5.51767e-11,4.01574e-12)               (1,1.1063e-10)

using fullPivHouseholderQr:
(1,1.1063e-10)  (-5.51766e-11,4.01578e-12)  (-5.97331e-11,9.63535e-12)  (-5.6913e-11,-3.22229e-11)  (-5.52733e-11,5.99029e-12)   (-6.46967e-11,1.0973e-11)
(-6.54343e-11,6.20279e-12)             (1,1.12229e-10)  (-6.04297e-11,5.62259e-12)  (-6.02256e-11,2.60221e-12) (-5.33448e-11,-3.05694e-11)  (-6.55427e-11,8.17723e-12)
(-6.48883e-11,1.07408e-11)  (-5.53615e-11,4.54731e-12)             (1,1.11384e-10)  (-5.99088e-11,9.29495e-12)  (-5.51744e-11,1.52693e-12) (-6.21597e-11,-3.11176e-11)
(-6.21592e-11,-3.11175e-11)  (-5.51743e-11,1.52697e-12)  (-5.99093e-11,9.29482e-12)             (1,1.11385e-10)  (-5.53614e-11,4.54728e-12)  (-6.48881e-11,1.07408e-11)
(-6.55427e-11,8.1773e-12) (-5.33444e-11,-3.05693e-11)  (-6.02261e-11,2.60223e-12)  (-6.04293e-11,5.62257e-12)             (1,1.12229e-10)  (-6.54344e-11,6.20272e-12)
(-6.46967e-11,1.0973e-11)  (-5.52732e-11,5.99039e-12)  (-5.69133e-11,-3.2223e-11)  (-5.97328e-11,9.63533e-12)  (-5.51767e-11,4.01574e-12)              (1,1.1063e-10)

using jacobi:
(1,1.1063e-10)  (-5.51761e-11,4.01562e-12)   (-5.9733e-11,9.63488e-12) (-5.69134e-11,-3.22232e-11)    (-5.5273e-11,5.9903e-12)  (-6.46969e-11,1.09735e-11)
(-6.54346e-11,6.20257e-12)             (1,1.12228e-10)  (-6.04303e-11,5.62232e-12)   (-6.0225e-11,2.60224e-12) (-5.33448e-11,-3.05694e-11)  (-6.55414e-11,8.17752e-12)
(-6.48874e-11,1.07408e-11)  (-5.53617e-11,4.54725e-12)             (1,1.11384e-10)  (-5.99097e-11,9.29465e-12)  (-5.51742e-11,1.52688e-12) (-6.21596e-11,-3.11175e-11)
(-6.21608e-11,-3.11175e-11)   (-5.5174e-11,1.52688e-12)  (-5.99104e-11,9.29433e-12)             (1,1.11384e-10)  (-5.53612e-11,4.54728e-12)  (-6.48873e-11,1.07413e-11)
(-6.55423e-11,8.17718e-12) (-5.33459e-11,-3.05695e-11)  (-6.02258e-11,2.60206e-12)   (-6.04294e-11,5.6224e-12)             (1,1.12229e-10)   (-6.5434e-11,6.20302e-12)
(-6.46966e-11,1.09729e-11)   (-5.5273e-11,5.99037e-12) (-5.69134e-11,-3.22233e-11)  (-5.97333e-11,9.63543e-12)  (-5.51767e-11,4.01582e-12)              (1,1.1063e-10)

using inverse:
(1,1.1063e-10)  (-5.51768e-11,4.01574e-12)  (-5.97329e-11,9.63535e-12)   (-5.6914e-11,-3.2223e-11)  (-5.52733e-11,5.99031e-12)   (-6.46967e-11,1.0973e-11)
(-6.54338e-11,6.20274e-12)             (1,1.12229e-10)   (-6.04293e-11,5.6225e-12)   (-6.02264e-11,2.6022e-12) (-5.33451e-11,-3.05693e-11)   (-6.55425e-11,8.1773e-12)
(-6.48877e-11,1.07408e-11)  (-5.53615e-11,4.54726e-12)             (1,1.11384e-10)  (-5.99097e-11,9.29488e-12)  (-5.51743e-11,1.52697e-12) (-6.21592e-11,-3.11176e-11)
(-6.21587e-11,-3.11176e-11)  (-5.51744e-11,1.52696e-12)  (-5.99093e-11,9.29489e-12)             (1,1.11384e-10)  (-5.53616e-11,4.54727e-12)   (-6.4888e-11,1.07408e-11)
(-6.55419e-11,8.17731e-12) (-5.33448e-11,-3.05693e-11)   (-6.02257e-11,2.6022e-12)   (-6.04301e-11,5.6225e-12)             (1,1.12229e-10)  (-6.54343e-11,6.20273e-12)
(-6.46964e-11,1.0973e-11)  (-5.52732e-11,5.99031e-12)  (-5.69131e-11,-3.2223e-11)  (-5.97336e-11,9.63535e-12)  (-5.51767e-11,4.01574e-12)              (1,1.1063e-10)

该算法在MATLAB中的第一次迭代产生了一个近似单位矩阵，可以在C中通过反转矩阵并从矩阵对象调用解算器来观察该矩阵。

但使用解算器对象时，结果完全错误。我还尝试了eigen文档页面中的其他解算器，得到了相同的结果：https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html.

我是否遗漏了一个步骤，或者从矩阵调用解算器会以不同的方式初始化解算器？

本征版本是3.3。我还没有使用任何编译器优化，只是C 11：

g++ -I. -std=c++11 main.cpp

矩阵的大小是恒定的，它将在一个大循环中执行。因此，通过拉出解算器可以避免的内存分配越多越好。

谢谢。

Answer 1

在第一个示例中，您在同一个表达式中使用了两次mdivide。。compute（）更改内部状态。

使用两个mDM对象，一个用于Y，一个用于M，可能会得到更好的结果。

Answer 2

小心使用auto，查看常见陷阱，基本上在您的情况下结果是产品