Лучший способ решения разреженных линейных систем в графическом процессоре Возможно?

В настоящее время я работаю над проектом, где нам нужно решить

|Ax - b|^2,

В этом случае, A это очень разреженная матрица и A'A имеет не более 5 ненулевых элементов в каждом ряду.

Мы работаем с изображениями и размером A'A является NxN где N — количество пикселей. В этом случае N = 76800, Мы планируем пойти в RGB и тогда измерение будет 3Nx3N,

В матлаб решения (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 а также OpenMP в Eigen, Что вы думаете об этих типах проблем? Есть ли способ извлечения структуры, например? И должен conjugateGradient действительно так медленно?

Редактировать:

Спасибо за ценные комментарии! Сегодня вечером я немного читал о cuSparse на Nvidia. Похоже, что можно делать факторизацию и даже решать системы. В частности, кажется, можно использовать шаблон и так далее. Вопрос в том, насколько быстро это может быть и каковы возможные накладные расходы?

Учитывая, что объем данных в моей матрице A такой же, как в изображении, копирование памяти не должно быть такой проблемой. Я сделал несколько лет назад программное обеспечение для 3D-реконструкции в реальном времени, а затем вы копировали данные для каждого кадра, и медленная версия все еще работает с частотой более 50 Гц. Так что, если факторизация намного быстрее, возможно ли это ускорение? В частности, остальная часть проекта будет на графическом процессоре, поэтому, если кто-то сможет решить его там напрямую и сохранить решение, я полагаю, это не является недостатком.

Много произошло в области Cuda, и я не совсем в курсе.

Вот две ссылки, которые я нашел: Benchmark?, MatlabGPU