在 C++ 中解决稀疏线性系统的最佳方法 - GPU 可能吗？

Question

我目前正在做一个我们需要解决的项目

|Ax - b|^2.

在这种情况下，A是一个非常稀疏的矩阵A'A，每行最多有 5 个非零元素。

我们正在处理图像，其中的维度A'A是NxNN 是像素数。在这种情况下N = 76800。我们计划去RGB然后维度将是3Nx3N。

在 matlab 中(A'A)\(A'b)，使用双精度求解大约需要 0.15 秒。

我现在已经对Eigens稀疏求解器进行了一些试验。我试过了：

SimplicialLLT
SimplicialLDLT
SparseQR
ConjugateGradient

和一些不同的排序。迄今为止最好的是

SimplicialLDLT

这需要0.35 - 0.5使用AMDOrdering.

例如，当我使用ConjugateGradient它时，大约需要6 s，0用作初始化。

解决问题的代码是：

    A_tot.makeCompressed();
     // Create solver
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>, Eigen::Lower, Eigen::AMDOrdering<int> > solver;
    // Eigen::ConjugateGradient<Eigen::SparseMatrix<float>, Eigen::Lower> cg;
    solver.analyzePattern(A_tot);
    t1 = omp_get_wtime();
    solver.compute(A_tot);
    if (solver.info() != Eigen::Success)
     {
         std::cerr << "Decomposition Failed" << std::endl;
         getchar();
     }
    Eigen::VectorXf opt = solver.solve(b_tot);
    t2 = omp_get_wtime();
    std::cout << "Time for normal equations: " << t2 - t1 << std::endl;

这是我第一次在 C++ 及其求解器中使用稀疏矩阵。对于此项目，速度至关重要，低于此速度是0.1 s最低要求。

我想得到一些反馈，哪些是最好的策略。例如，应该能够使用SuiteSparse和OpenMPin Eigen。您对这些类型的问题有何经验？例如，有没有办法提取结构？而应该conjugateGradient真的是慢？

编辑：

感谢您提出宝贵意见！今晚我在 Nvidia 上阅读了一些关于 cuSparse 的文章。它似乎能够进行分解一个偶数求解系统。特别是似乎可以重用模式等等。问题是这能有多快，可能的开销是多少？

鉴于我的矩阵 A 中的数据量与图像中的数据量相同，内存复制不应该是这样的问题。几年前我做了实时 3D 重建软件，然后你复制每一帧的数据，一个慢速版本仍然以超过 50 Hz 的频率运行。因此，如果分解速度要快得多，那么可能会加速吗？特别是项目的其余部分将在 GPU 上，所以如果可以直接在那里解决它并保留解决方案，我想这没有什么缺点。

Cuda 领域发生了很多事情，我并不是真正了解最新情况。

这是我找到的两个链接：基准？, MATLAB GPU

Answer 1

您的矩阵非常稀疏，对应于 2D 域上的离散化，因此预计SimplicialLDLT这里是最快的。由于稀疏模式是固定的，调用analyzePattern一次，然后factorize代替compute。这应该会节省一些毫秒。此外，由于您在常规网格上工作，您还可以尝试使用NaturalOrdering（不是 100% 确定，您必须进行工作台）绕过重新排序步骤。如果这仍然不够快，您可以搜索为天际线矩阵量身定制的 Cholesky 求解器（在这种情况下，Cholesky 分解更简单，因此速度更快）。