线性代数方程 Ax = By 在 Eigen C++ 中的最小二乘解答案

【问题标题】：Least Squares Solution of Linear Algerbraic Equation Ax = By in Eigen C++线性代数方程 Ax = By 在 Eigen C++ 中的最小二乘解
【发布时间】：2014-02-28 17:11:54
【问题描述】：

我有一组矩阵形式的线性代数方程，Ax=By。其中 A 是 36x20 的矩阵，x 是大小为 20 的向量，B 是 36x13 和 y 是13x1。 排名(A)=20。因为系统是超定的（方程数量多于变量），所以最小二乘解是可能的，即； x = (A^TA)^-1A^TBy。我想要解决方案，以便最小化残差 e = Ax-By。

使用 C++ 的Eigen/Dense 库，我已经制定了所有矩阵等。我尝试了此页面上描述的方法 Eigen Tutorial!

我猜这个页面中描述的方法只适用于方阵。因为当它尝试运行它时会出错。

 x = A.jacobiSvd( ComputeThinU | ComputeThinV ).solve(B*y);

错误

 /usr/include/eigen3/Eigen/src/SVD/JacobiSVD.h: In member function 'const    
 Eigen::internal::solve_retval<Eigen::JacobiSVD<MatrixType, QRPreconditioner>, Rhs> 
 Eigen::JacobiSVD<MatrixType, QRPreconditioner>::solve(const 
 Eigen::MatrixBase<OtherDerived>&) const [with Rhs = 
 Eigen::GeneralProduct<Eigen::Matrix<float, 36, 13>, Eigen::Matrix<double, -1, 1>, 4>; 
 _MatrixType = Eigen::Matrix<float, 36, 20>; int QRPreconditioner = 2]':
 /usr/include/eigen3/Eigen/src/SVD/JacobiSVD.h:658:5: warning: control reaches end of  
 non-void function [-Wreturn-type]
 make[2]: *** [src/CMakeFiles/spacebot_actuationKinematics.dir 
 /ActuationKinematics.cpp.o] Error 1
 make[1]: *** [src/CMakeFiles/spacebot_actuationKinematics.dir/all] Error 2
 make: *** [all] Error 2

【问题讨论】：

页面支持二维码分解。你试过了吗？它也适用于最小二乘的矩形矩阵 A，以及 SVD 方法。您的错误来自报告函数的某些特殊情况下缺少返回语句，请检查 Eigen 的错误报告。
y 是什么意思？ LS 方法是找到 x 的线性组合的参数，该参数最适合 x 观察到的 y 观察。因此，您将 x 值和 y 值替换为方程。你有 36 个观察值吗？
@LutzL 是的，我试过了。这是输出。断言失败，因为 rows == cols。 ActuationKinematicsTest: /usr/include/eigen3/Eigen/src/LU/Inverse.h:334: const Eigen::internal::inverse_impl Eigen::MatrixBase::inverse() const [with Derived = Eigen ::Matrix]: 断言 `rows() == cols()' 失败。中止（核心转储）
是的，当然，如果您尝试使用 LU 分解来计算逆矩阵，这是可以预料的。这仅适用于常规方阵。使用 QR 和 SVD 一样快但更稳定的因式分解，不计算逆矩阵，只求解线性系统。
这是一个相当强烈的声明。如果 A 具有满秩和合理的条件数，则 QR 在数值上与 SVD 一样稳定。 LU 只有在使用完整的、昂贵的旋转时才能接近。如果矩阵 A 秩不足或几乎如此，您将需要 SVD 来确定有效秩。但是，如果您决定奇异值是很小还是实际上为零，那么线性系统的任何解决方案都将极大地依赖于阈值。

标签： c++ linear-algebra eigen

【解决方案1】：

您的求解（B*y）中的矩阵乘法似乎有问题；部分。尝试单独做 B*y 并使用 solve(result);而是。

Eigen::GeneralProduct<Eigen::Matrix<float, 36, 13>, Eigen::Matrix<double, -1, 1>, 4>

那句话让我产生了这种怀疑。它说 y 变量的大小为 -1x1，因此您的程序无论如何都不会运行，因为它不能与矩阵相乘。

另外，教程说...

A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << endl;

我不知道 Eigen 究竟是如何工作的，但这似乎是问题所在。

【讨论】：

如果矩阵 B 的顺序是 m x n 并且向量 y 的顺序是 p x q 那么对于乘法 p == n 并且在我的情况下它总是正确的。如果我们在参数容器之外或在参数容器中进行乘法，有什么区别？这条线没有更多的问题。 QR 分解适用于我的情况，但 x 的解决方案并不精确。因为我的 A,B 任何 y 都会随着时间而变化。
@wafeeq 取决于方法重载，它如何接受它的参数。在括号内完成时，它可能无法解决乘法问题。我知道矩阵是正确的，从程序员的角度来看，似乎第二个 arg 没有被正确处理，只需查看 gcc 的调试信息。我确定的唯一方法是更深入地调查 eigen .

【解决方案2】：

如文档中所述，ComputeThin* 选项仅适用于Dynamic 大小的矩阵。对于固定大小，您必须使用ComputeFull*。不过，在您的情况下，最好使用Dynamic 大小矩阵，即MatrixXf。固定大小的矩阵只对非常小的矩阵有意义。

最后，ColPivHouseholderQR 可能是最小二乘求解的更好选择。 SVD 有点矫枉过正。

【讨论】：

我想这些矩阵是否是动态分配的内存与它们无关。我已经在我的代码中使用了 MatrixXd。问题是最小二乘解的误差是不可接受的。它很大。 error = 通过最小二乘法找到 x 的值后的 Ax-By。