c++ - Eigen LSCG 求解器性能问题

Question

我正在使用 Eigen 的 LSCG 迭代求解一个我用稀疏矩阵表示的超定系统，我相信也是如此 slow .迭代我的意思是这样的:

    //The following is pseudo code
    main() {
    //compute A
    A=..
    //compute B
    B=..
    while(someCondition)
    {   
        x=solve(A,B)
        //based on x alter A and B  
        A=..
        B=..
    }       
}

例如，当 A 有 36k 行和 18k 列以及 145k nnz 元素时，B 有 36k 行 3 列和 110k nnz 元素(注意 B 实际上是密集的)我的桌面需要 74 秒才能执行 x=solve(A,B)。

1. 这些时间正常吗？我想答案取决于机器 代码正在运行。我有一个 AMD FX 6300，所以我想 硬件不是主要问题。
2. 如果确实很慢，可能是什么问题？ B 矩阵实际上不是稠密的这一事实会减慢求解步骤吗？

为了在你的机器上测试时间，我写了一些简单的测试代码:

#include <Eigen/Sparse> //system solving and Eigen::SparseMatrix
#include <ctime> //measure time to execute
#include <unsupported/Eigen/SparseExtra> //loadMarket

using SpMatrix = Eigen::SparseMatrix<double>;
using Matrix = Eigen::MatrixXd;
int main() {

  //load A and B
  SpMatrix A, B;
  Eigen::loadMarket(A, "/AMatrixDirectory/A.mtx");
  Eigen::loadMarket(B, "/BMatrixDirectory/B.mtx");

  std::clock_t start;

  start = std::clock();

  Eigen::LeastSquaresConjugateGradient<SpMatrix> solver;
  solver.compute(A);
  Matrix x = solver.solve(B);

  std::cout << "Time: " << (std::clock() - start) / (double)(CLOCKS_PER_SEC)
            << " s" << std::endl;

  return 0;
}

以上是上述伪代码中 while 循环的一次迭代。 要运行上述代码，您需要:

1. Eigen
2. C++11(如果不使用 typedef 替换 using)
3. 我上传的矩阵A和B here

---

---

编辑

ggael 建议使用 SimplicialLDLT在我的问题中声称它比 LSCG 具有更好的性能。

为了将 Eigen 的求解器性能与特定问题进行比较，Eigen 提供了 BenchmarkRoutine .不幸的是我无法使用它，因为它只能使用方阵。

我编辑了比较 LSCG 和 SimplicialLDLT 的代码(请不要因为代码的长度而气馁，因为它由 4 个彼此相似的 block 组成，只有一些小的变化):

#include <Eigen/Sparse> //system solving and Eigen::SparseMatrix
#include <ctime>        //measure time to execute
#include <unsupported/Eigen/SparseExtra> //loadMarket

// Use typedefs instead of using if c++11 is not supported by your compiler
using SpMatrix = Eigen::SparseMatrix<double>;
using Matrix = Eigen::MatrixXd;

int main() {
  // Use multi-threading. If you don't have OpenMP on your system
  // it will simply use 1 thread and it will not crash. So do not worry about
  // that.
  Eigen::initParallel();
  Eigen::setNbThreads(6);

  // Load system matrices
  SpMatrix A, B;
  Eigen::loadMarket(A, "/home/iason/Desktop/Thesis/build/A.mtx");
  Eigen::loadMarket(B, "/home/iason/Desktop/Thesis/build/B.mtx");

  // Initialize the clock which will measure the solvers
  std::clock_t start;

  {
    // Reset clock
    start = std::clock();
    // Solve system using LSCG
    Eigen::LeastSquaresConjugateGradient<SpMatrix> LSCG_solver;
    LSCG_solver.setTolerance(pow(10, -10));
    LSCG_solver.compute(A);
    // Use auto specifier to hold the return value of solve. Eigen::Solve object
    // is being returned
    auto LSCG_solution = LSCG_solver.solve(B);
    std::cout << "LSCG Time Using auto: "
              << (std::clock() - start) / (double)(CLOCKS_PER_SEC) << " s"
              << std::endl;
  }
  {
    // Reset clock
    start = std::clock();
    // Solve system using LSCG
    Eigen::LeastSquaresConjugateGradient<SpMatrix> LSCG_solver;
    LSCG_solver.setTolerance(pow(10, -10));
    LSCG_solver.compute(A);
    // Use Matrix specifier instead of auto. Implicit conversion taking place(?)
    Matrix LSCG_solution = LSCG_solver.solve(B);
    std::cout << "LSCG Time Using Matrix: "
              << (std::clock() - start) / (double)(CLOCKS_PER_SEC) << " s"
              << std::endl;
  }
  {
    // Reset clock
    start = std::clock();
    // Solve system using SimplicialLDLT
    Eigen::SimplicialLDLT<SpMatrix> SLDLT_solver;
    SLDLT_solver.compute(A.transpose() * A);
    auto SLDLT_solution = SLDLT_solver.solve(A.transpose() * B);

    std::cout << "SimplicialLDLT Time Using auto: "
              << (std::clock() - start) / (double)(CLOCKS_PER_SEC) << " s"
              << std::endl;
  }
  {
    // Reset clock
    start = std::clock();
    // Solve system using SimplicialLDLT
    Eigen::SimplicialLDLT<SpMatrix> SLDLT_solver;
    SLDLT_solver.compute(A.transpose() * A);
    // Use Matrix specifier instead of auto. Implicit conversion taking place(?)
    Matrix SLDLT_solution = SLDLT_solver.solve(A.transpose() * B);

    std::cout << "SimplicialLDLT Time Using Matrix: "
              << (std::clock() - start) / (double)(CLOCKS_PER_SEC) << " s"
              << std::endl;
  }
  return 0;

以上代码产生如下结果:

LSCG 时间使用自动:0.000415 秒

使用矩阵的 LSCG 时间:52.7668 秒

使用自动的 SimplicialLDLT 时间:0.216593 秒

使用矩阵的 SimplicialLDLT 时间:0.239976 秒

结果证明，当我使用 Eigen::MatrixXd 而不是 auto 时，我得到了 ggael 在他的评论之一中描述的情况，即 LSCG 在没有设置更高容差的情况下无法实现解决方案，而 SimplicialLDLT 速度要快得多。

1. 当我使用 auto 时，求解器是否真的在求解系统？
2. Solver 对象是否被隐式转换为 Matrix 对象 当我使用 Matrix 说明符时？因为当使用 LSCG 时唯一 前两种情况的变化是使用 auto 和 Matrix 分别，这个隐式转换到 Matrix 需要 52.7668-0.000415 秒？
3. 有没有更快的方法从 解决对象？

Answer 1

给定矩阵 A 的稀疏模式，明确地形成正规方程并使用直接求解器(如 SimplicialLDLT)会更快。最好对 B 使用密集类型，因为它无论如何都是密集的，并且在内部，所有求解器都会将稀疏的右侧转换为密集列:

Eigen::MatrixXd dB = B; // or directly fill dB
Eigen::SimplicialLDLT<SpMatrix> solver;
solver.compute(A.transpose()*A);
MatrixXd x = solver.solve(A.transpose()*dB);

这需要 0.15 秒，而将 LSCG 的公差设置为 1E-14 后，LSCG 需要 6 秒。