EIGEN LSCG求解器性能问题

Eigen LSCG solver performance issue

本文关键字:性能 问题 LSCG EIGEN      更新时间:2023-10-16

我正在使用eigen的LSCG迭代求解一个我使用稀疏矩阵表示的过度确定系统,我相信太慢了。通过迭代,我的意思是:

    //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行和145K NNZ元素的18K COLS,B具有36K行3 col和110k nnz元素(请注意B实际上是密集的),我的桌面74s才能执行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;
}

上面是上述伪代码中的WHIL循环的一种迭代。要运行上述代码,您需要:

  1. eigen
  2. C 11(如果不用Typedefs替换)
  3. 我在这里上传的矩阵A和B


编辑

ggael提议使用SimplicialDlt声称与我的问题相比,它的性能更好。

为了将特定的求解器的性能与特定问题进行比较,eigen提供了基准测试。不幸的是,我无法使用它,因为只能与它一起使用方形矩阵。

我编辑了比较LSCG和SimplicialDlt的代码(请不要对代码的长度灰心,因为它由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

使用自动的SimplicialDlt时间:0.216593 S

使用矩阵的SimplicialDlt时间:0.239976 S

作为结果,当我使用eigen :: matrixxd而不是自动时,我得到了GGAEL在他的一篇评论中所描述的情况,即LSCG在没有设置更高容忍度和SimpleicialDlt的情况下无法实现解决方案。

  1. 当我使用自动时,求解器实际上是在解决系统吗?
  2. 是被隐式转换为矩阵对象的求解器对象当我使用矩阵指定符时?因为使用LSCG唯一在前两种情况下,更改是使用自动和矩阵分别对矩阵的隐式转换为52.7668-0.000415秒?
  3. 是否有一种更快的方法来提取溶液矩阵从解决对象?

给定矩阵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.15s,与将LSCG耐受性的LSCG的6s相比。