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

我正在使用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 col具有145k nnz元素而B具有
36k行3 cols和110k nnz元素(注意B实际上是密集的)需要我的桌面74s执行x = solve(A，B)。

这些时间正常吗？我猜答案取决于机器
该代码正在运行。我有AMD FX 6300，所以我想
硬件不是主要问题。

如果确实很慢，可能是什么问题？ 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循环的一次迭代。
要运行以上代码，您将需要:

本征

C++ 11(如果不使用typedef代替)

我上载了here

的矩阵A和B



编辑

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

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

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

#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 s

LSCG使用时间矩阵:52.7668 s

使用自动的SimplicialLDLT时间:0.216593 s

使用矩阵的SimplicialLDLT时间:0.239976 s

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

当我使用auto时，求解器实际上是在求解系统吗？

是否将规划求解对象隐式转换为矩阵对象
当我使用矩阵说明符时？由于仅当使用LSCG时
前两种情况的变化是使用auto和Matrix
分别将此隐式转换为矩阵
52.7668-0.000415秒？

是否有更快的方法可以从中提取解矩阵
解决对象？

最佳答案

给定矩阵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);

将LSCG的公差设置为1E-14之后，LSCG的时间为0.15s，而LSCG的时间为6s。

关于c++ - Eigen LSCG求解器性能问题，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/46014719/