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

问题描述

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

我选择了带有 LLT 分解的 Eigen3 库来求解线性系统，与其他方法（如 LU）相比，它要快得多。该系统解决方案在配备 4.8 GHz 处理器的英特尔 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);
    }

解决方法

您有一个稀疏矩阵，但您在 Eigen 中将其表示为一个稠密矩阵。您拥有的矩阵文件也很密集，如果以稀疏形式（例如市场格式）存储使用起来会更方便。

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

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

然后在我的 PC 上（应该比你的慢，它只有一个旧的 4770K）实际因子/求解只需要 0.01 秒。

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