使用 Eigen 求解线性方程组

Posted 2023-02-22

【中文标题】使用 Eigen 求解线性方程组

【英文标题】：Solving Systems of Linear Equations using Eigen

【发布时间】：2018-09-04 02:55:52

【问题描述】：

我目前正在使用 C++ 进行流体模拟，该算法的一部分是求解线性方程组的稀疏系统。人们建议为此使用库 Eigen。我决定使用我编写的这个简短程序对其进行测试：

#include <Eigen/SparseCholesky>

#include <vector>
#include <iostream>

int main() 
    std::vector<Eigen::Triplet<double>> triplets;
    triplets.push_back(Eigen::Triplet<double>(0, 0, 1));
    triplets.push_back(Eigen::Triplet<double>(0, 1, -2));
    triplets.push_back(Eigen::Triplet<double>(1, 0, 3));
    triplets.push_back(Eigen::Triplet<double>(1, 1, -2));

    Eigen::SparseMatrix<double> A(2, 2);
    A.setFromTriplets(triplets.begin(), triplets.end());

    Eigen::VectorXd b(2);
    b[0] = -2;
    b[1] = 2;

    Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> chol(A);
    Eigen::VectorXd x = chol.solve(b);

    std::cout << x[0] << ' ' << x[1] << std::endl;

    system("pause");

它给出了这两个方程：

x - 2y = -2

3x - 2y = 2

正确的解决方法是：

x = 2

y = 2

但问题是程序运行时会输出： 0.181818 -0.727273

这是完全错误的！我已经调试了几个小时，但这是一个非常短的程序，我完全按照 Eigen 网站上的教程进行操作。有谁知道是什么导致了这个问题？

附：我知道我使用的类是用于稀疏矩阵的，但它们与普通 Matrix 类之间的唯一区别是元素的存储方式。

【问题讨论】：

你已经建立了一个转置矩阵

@n.m.我以为可能是这样，但即使我转置它，我仍然得到错误的答案：/

@n.m.你在说什么？即使我转置系数矩阵，这些答案仍然不起作用

嗯，不，对不起，我搞砸了我的数学。这不是针对转置矩阵的解决方案，而是针对 (1,3,3,-2) 的解决方案。

这只是合乎逻辑的，因为单纯的 Cholesky 仅适用于自伴矩阵！它在 Eigen 中也已被弃用。

【参考方案1】：

SimplicialCholesky 用于对称正定 (SPD) 矩阵，您组装的矩阵甚至不是对称的。默认情况下它只读取下三角形部分中的条目而忽略其他条目，因此它解决了：

x + 3y = -2
3x -2y = 2

如您所见，对于非对称平方问题，您需要使用基于 LU 的直接求解器或迭代求解器领域中的 BICGSTAB。这一切都在doc中进行了总结。

【讨论】：

并由 eigen 的主要开发者回答！感谢您在 eigen 方面的出色工作。

【参考方案2】：

您应该使用能够处理非对称稀疏矩阵的求解器。另一种可能的方法是寻求不是原始系统 [A]x=b 而是 [A]T*[A]x=[A]T*b 的解决方案，其中 [A]T 代表 [A] 转置.后一个系统的矩阵是对称的并且是正定的（只要 [A] 是非奇异的）。唯一的缺点是如果原来的 [A] 在这个意义上不是“好”的，那么 [A]T[A] 可能是相当病态的。只是一个旨在解决此类问题的软件示例： http://members.ozemail.com.au/~comecau/CMA_LS_Sparse.htm