Как эффективно решить линейную систему с лапласианом + диагональной матрицей?

В моей реализации алгоритма обработки изображений мне нужно решить большую линейную систему вида A*x=b, где:

  • Матрица A=L+D представляет собой сумму матрицы Лапласа L и диагональной матрицы D
  • Матрица Лапласа L разреженная, около 25 отличных от нуля в строке.
  • Система большая, в ней столько неизвестных, сколько пикселей во входном изображении (обычно > 1 миллиона).

Матрица Лапласа L не меняется между последовательными запусками алгоритма; Я могу построить эту матрицу в предварительной обработке и, возможно, вычислить ее факторизацию. Диагональная матрица D и правосторонний вектор b изменяются при каждом запуске алгоритма.

Я пытаюсь выяснить, какой самый быстрый способ решения системы во время выполнения; Я не против потратить время на предварительную обработку (например, для вычисления факторизации L).

Моя первоначальная идея заключалась в том, чтобы предварительно вычислить факторизацию Холецкого для L, затем обновить факторизацию во время выполнения значениями из D (обновление ранга 1 с помощью cholupdate) и быстро решить проблему с обратной подстановкой. К сожалению, факторизация Холецкого не такая разреженная, как исходная L-матрица, и только загрузка ее с диска уже занимает 5,48 с; для сравнения, для прямого решения системы с обратной косой чертой требуется 8,30 с.

Учитывая форму моих матриц, есть ли какой-либо другой метод, который вы бы порекомендовали для ускорения решения во время выполнения, независимо от того, сколько времени это занимает во время предварительной обработки?


person Asmo    schedule 29.01.2013    source источник
comment
Попробуйте SVD, хотя я не знаю, насколько он быстр для такого большого объема данных...   -  person Ander Biguri    schedule 29.01.2013
comment
Я не понимаю, почему применение SVD здесь было бы хорошей идеей, учитывая форму моих матриц.   -  person Asmo    schedule 30.01.2013
comment
да, возможно вы правы, я не эксперт. На всякий случай, если вы не учли это :)   -  person Ander Biguri    schedule 31.01.2013
comment
Взгляните на это: заголовок scicomp.stackexchange.com/questions/1001/   -  person johnish    schedule 31.01.2013


Ответы (1)


Предполагая, что вы работаете с сеткой (поскольку вы упоминаете изображения — хотя это не гарантируется), что вас больше интересует скорость, чем точность (поскольку 5s кажется уже слишком медленным для 1 миллиона неизвестных), я вижу несколько вариантов.

Во-первых, забудьте о точных методах, таких как Холецкий (+переупорядочивание). Даже если они позволяют хранить факторизацию и повторно использовать ее для нескольких правых, вам, вероятно, придется хранить гигантские матрицы, которые кажутся неразрешимыми в вашем случае (надеюсь, вы переупорядочиваете строки/столбцы с обратным Катхиллом Макки или чем-то еще). в противном случае - это сильно разрежает факторизацию).

В зависимости от ваших граничных условий я бы сначала попробовал Matlab poisolv, который решает задачу Пуассона с использованием БПФ и возможных перепроекций, если вы хотите, чтобы граничные условия Дирихле были вместо периодических. Это очень быстро, но может не подходить для вашей проблемы (вы упоминаете наличие 25 nnz для матрицы Лапласа + идентичность: почему? - это матрица Лапласа высокого порядка, и в этом случае вас может больше интересовать точность, чем то, что я предполагаю или это на самом деле другая проблема, чем та, которую вы описываете?).

Затем вы можете попробовать многосеточные решатели, которые очень быстро работают с изображениями и гладкими задачами. Вы можете использовать простой метод релаксации для каждой итерации и каждого уровня мультисетки или использовать более сложные методы (например, предварительно обусловленный парный уровень сопряженного градиента). В качестве альтернативы вы можете выполнить более простой предварительно подготовленный сопряженный градиент (или даже SSOR) без многосеточной сетки, и если вас интересует только приблизительное решение, вы можете остановить итерации до полной сходимости.

Мои аргументы в пользу итеративных решателей:

  • вы можете остановиться до сходимости, если вам нужна приблизительная задача
  • вы по-прежнему можете повторно использовать другие результаты для инициализации вашего решения (например, если ваши разные прогоны соответствуют разным кадрам видео, то использование решения предыдущего кадра в качестве инициализации следующего имеет смысл).

Конечно, прямой решатель, для которого вы можете предварительно вычислить, сохранить и сохранить факторизацию, также имеет смысл (хотя я не понимаю вашего аргумента в пользу обновления ранга 1, если ваша матрица постоянна), поскольку остается сделать только обратную подстановку в время выполнения. Но, учитывая, что это игнорирует структуру проблемы (регулярная сетка, возможный интерес к результатам с ограниченной точностью и т. д.), я бы выбрал методы, разработанные для этих случаев, такие как методы, подобные Фурье, или многосеточные методы. Оба метода могут быть реализованы на графическом процессоре для получения более быстрых результатов (напомним, что графические процессоры скорее приспособлены для работы с изображениями/текстурами!).

Наконец, вы можете получить интересные ответы от scicomp.stackexchange, который больше ориентирован на численный анализ.

person nbonneel    schedule 01.02.2013
comment
+1 за упоминание об изменении порядка! Я бы просто добавил к вашему списку неполные факторизации Холецкого, в остальном совершенные! - person Dr_Sam; 04.02.2013
comment
спасибо :) У меня сложилось впечатление, что неполные факторизации Холецкого для предобусловливания были намного медленнее, чем простой предобусловливатель Якоби, лишь для незначительного улучшения скорости сходимости; не так ли? - person nbonneel; 06.02.2013