使用 Eigen3 快速求解稀疏正定线性系统

Question

我需要对具有稀疏对称正定矩阵 (2630x2630) 的线性系统求解数百万次。我已经在 Mathematica 中绘制了矩阵，如下所示。

我选择了具有 LLT 分解的 Eigen3 库来求解线性系统，与 LU 等其他方法相比，它要快得多。该系统解决方案在配备 4.8 GHz 处理器的 intel 10700 中耗时 0.385894 秒。代码：

#include <Eigen/Core>
using namespace Eigen;
using namespace std;


int main()
{
    ifstream datamatrix("matrix.txt");
    ifstream datavec("vector.txt");
    int m=2630;
    MatrixXd A(m,m);
    VectorXd b(m);
    for (int i = 0; i < m; i++)
    {
     for (int j = 0; j < m; j++)
     {
         datamatrix >> A(i,j);
     }
     datavec >> b(i);
    }
    chrono::steady_clock sc;
    auto start = sc.now();
    VectorXd x = A.llt().solve(b);
    auto end = sc.now();
    // measure time span between start & end
    auto time_span = static_cast<chrono::duration<double>>(end - start);
    cout << "Operation took: " << time_span.count() << " seconds.";
}

是否可以使用 Eigen3 或 MKL 来加速它？

矩阵和向量文件可以在这里下载：

矩阵：https://www.dropbox.com/s/k521t91cd8u7t5h/matrix.txt?dl=0

向量：https://www.dropbox.com/s/ldajnzl2qj3a7zh/vector.txt?dl=0

编辑

我发现大部分时间花在填充稀疏矩阵上。使用下面的代码使用 Triplet 来填充稀疏矩阵，速度要快得多！观察。 MatDoub 是来自数字食谱代码的全稠密矩阵 - nr3.h

    void slopeproject::SolveEigenSparse(MatDoub A, MatDoub b, MatDoub& x)
    {
    
        typedef Eigen::Triplet<double> T;
        std::vector<T> tripletList;
        int sz=A.nrows();
    
//approximated number of non-zero entries in the matrix
        tripletList.reserve(sz*50);
       // tripletList.reserve(80000);
    
        x.assign(sz, 1, 0.);
        SparseMatrix<double> AA(sz, sz);
        VectorXd bbb(sz);
        for (int i = 0; i < sz; i++)
        {
            for (int j = 0; j < sz; j++)
            {
//checking if the value of the dense matrix is zero. If not it is appended to the tripletlist

                if(fabs(A[i][j])>1.e-12)
                {
                    tripletList.push_back(T(i,j,A[i][j]));
                }
            }
            bbb(i) = b[i][0];
        }
    //transfer from the tripletlist to the sparse matrix very fast

        AA.setFromTriplets(tripletList.begin(), tripletList.end());
    //optional. I dont know wath this function do.
        AA.makeCompressed();
        SimplicialLLT< SparseMatrix<double> > solver;
//solve the system
        VectorXd xx = solver.compute(AA).solve(bbb);
        for(int i=0;i<sz;i++)x[i][0]=xx(i);
    }

Answer 1

您有一个稀疏矩阵，但您在 Eigen 中将其表示为稠密矩阵。你手头的矩阵文件也是密集的，如果是稀疏的形式使用起来会更方便，比如Market格式。

如果我将矩阵更改为稀疏矩阵，并使用

SimplicialLLT< SparseMatrix<double> > solver;
VectorXd x = solver.compute(A).solve(b);

然后在我的电脑上（应该比你的慢，它只有一个旧的4770K）实际factor/solve只需要0.01秒。

顺便说一句，即使是密集求解也应该比您看到的更快，所以我猜您没有在编译器设置中启用 SIMD。

使用 Eigen3 快速求解稀疏正定线性系统

Fast solve a sparse positive definite linear system with Eigen3

c++

linear-algebra

solver

intel-mkl

eigen