Лучший способ решения разреженных линейных систем на C++ - GPU Возможно?

Question

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

|Ax - b|^2.

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

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

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

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

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

Многое произошло в области Куды, и я не в курсе последних событий.

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

Answer 1

Ваша матрица является чрезвычайно скудны и соответствует дискретности на 2D области, поэтому ожидается, что SimplicialLDLT является самым быстрым здесь. Поскольку шаблон разреженности фиксирован, вызовите analyzePattern один раз, а затем factorize вместо compute. Это должно сэкономить несколько миллисекунд. Более того, поскольку вы работаете с обычной сеткой, вы также можете попытаться обойти шаг повторного упорядочения, используя NaturalOrdering (не на 100% уверен, вам нужно скамейку). Если это еще не достаточно быстро, вы можете найти решение Cholesky, адаптированное к строкам горизонта (факторизация Cholesky намного проще и, следовательно, быстрее в этом случае).