5

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

|Ax - b|^2.

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

我们正在处理图像，其中 NA'A是NxN像素数。在这种情况下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是最低要求。

我想得到一些反馈，这将是最好的策略。例如，一个应该能够使用SuiteSparseand OpenMPin Eigen。你对这些类型的问题有什么经验？例如，有没有一种提取结构的方法？真的应该conjugateGradient那么慢吗？

编辑：

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

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

Cuda 领域发生了很多事情，而我并不是最新的。

这是我找到的两个链接：基​​准？, MatlabGPU

4

1 回答 1

5

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

于 2017-02-08T15:38:38.903 回答