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

Question

我目前正在进行一个项目，我们需要解决

CCD_ 1。

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

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

在matlab中，使用doubles求解(A'A)(A'b)大约需要0.15s。

我现在已经用|Ax - b|^20稀疏解算器做了一些实验。我试过：

SimplicialLLT
SimplicialLDLT
SparseQR
ConjugateGradient

以及一些不同的排序。到目前为止最好的是

SimplicialLDLT

其使用CCD_ 12占用大约CCD_。

例如，当我使用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是最低要求。

我想得到一些关于哪种策略是最好的反馈。例如，假定能够在CCD_ 19中使用CCD_ 17和CCD_ 18。你对这些类型的问题有什么经验？例如，有没有提取结构的方法？conjugateGradient真的应该那么慢吗？

编辑：

感谢您的宝贵意见！今晚我在英伟达上读了一些关于cuSparse的文章。它似乎可以做因式分解，甚至可以求解系统。特别是，似乎可以重用模式等等。问题是这可能有多快，可能的开销是多少？

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

在库达领域发生了很多事情，而我并不是真正的最新。

以下是我发现的两个链接：Benchmark？，MatlabGPU

Answer 1

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