В настоящее время я работаю над проектом, где нам нужно решить
|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
Ваша матрица чрезвычайно разрежена и соответствует дискретизации в 2D-области, поэтому ожидается, что SimplicialLDLT
здесь самый быстрый Поскольку шаблон разреженности исправлен, звоните analyzePattern
один раз, а потом factorize
вместо compute
, Это должно сэкономить несколько миллисекунд. Более того, поскольку вы работаете с обычной сеткой, вы можете также попытаться обойти этап переупорядочения, используя NaturalOrdering
(не уверен на 100%, вы должны скамейку). Если это все еще не достаточно быстро, вы можете найти решатель Холески, специально разработанный для матриц горизонта (факторизация Холески намного проще и, следовательно, быстрее в этом случае).
Других решений пока нет …