Решение системы с обновлением диагонали малого ранга

9

Предположим, у меня есть оригинальная большая, разреженная линейная система: Ax0=b0 . Теперь у меня нет A1 как A слишком велика для разложения или любого вида разложения A , но предположим, что у меня есть решение x0 найденное с помощью итеративного решения.

Теперь я хочу применить небольшое ранговое обновление к диагонали A (измените несколько диагональных элементов): (A+D)x1=b0 где D - диагональная матрица с в основном 0 по диагонали и несколько ненулевых значений. Если бы у меня был A1 я бы смог воспользоваться формулой Вудбери, чтобы применить обновление к обратному. Тем не менее, у меня нет этого в наличии. Есть ли что-то, что я могу сделать, кроме как заново решить всю систему? Может быть, есть какой-то способ, которым я могу придумать предобусловливатель M который легко \ легче инвертировать, такой, что MA1A0 , так что все, что мне нужно сделать, если у меня есть x0 это применить M1 и итерационный метод будет сходиться в пару / несколько итераций?

КОСТИС
источник
Вы начинаете с хорошего предварительного кондиционера для и хотите знать, как его обновить? Какой рейтинг это обновление? ( Обновление ранга «мало» по сравнению с матрицей размером но не мало с точки зрения количества итераций.)A1000109
Джед Браун
AРазмер составляет от до , а обновление составляет <1000 (вероятно, <100) элементов. Я использую предкондиционер диагонального типа для A, который работает очень хорошо, поэтому обновление будет тривиальным, но мне было интересно, есть ли что-то лучшее, что я могу сделать, вместо разрешения новой системы с нуля. 106107
Костис
2
Решение одной системы мало о чем вам говорит. Если вы решаете одну и ту же систему несколько раз, обратная карта этих векторов (и / или связанных пространств Крылова) дает вам некоторую информацию, которую можно использовать для ускорения сходимости. Сколько систем вы решаете в каждом конкретном случае?
Джед Браун
В настоящее время я решение только для одного РИТ ( вектор) с каждой матрицы перед модификацией . bAA
Костис

Ответы:

4
  1. Сохраните в столбцах двух матриц и все векторы к которым вы применяли матрицу на предыдущих итерациях, и результаты .BCbjcj=Abj

  2. Для каждой новой системы (или , что является особым случаем ), приблизительно решить переопределенную линейную систему , например, выбрав подмножество строк (возможно, все) и используя плотный метод наименьших квадратов. Обратите внимание, что только выбранная часть должна быть собрана; так что это быстрая операция!(A+D)x=bAx=bD=0(C+DB)ybC+DB

  3. Положите . Это хорошее начальное приближение, с которого начинается итерация для решения . В случае, если дальнейшие системы должны быть обработаны, используйте матричные векторные произведения в этой новой итерации, чтобы расширить матрицы и в результирующей подсистеме.x0=By(A+D)x=bBC

Если матрицы и не помещаются в основную память, сохраните на диске и заранее выберите поднабор строк. Это позволяет вам хранить в ядре соответствующую часть и необходимую для формирования системы наименьших квадратов, и следующий может быть вычислен за один проход через с небольшим использованием памяти ядра.BCBBCx0B

Строки должны быть выбраны таким образом, чтобы они приблизительно соответствовали грубой дискретизации полной задачи. Достаточно взять в пять раз больше строк, чем общее число ожидаемых умножителей на матрицы.

Редактировать: Почему это работает? По построению матрицы и связаны соотношением . Если подпространство, натянутое на столбцы содержит точный вектор решения (редкая, но простая ситуация), то имеет вид для некоторого . Подставляя это в уравнение, определяющее получаем уравнение . Таким образом, в этом случае вышеуказанный процесс дает в качестве отправной точки , что является точным решением.BCC=ABBxxx=Byyx(C+DB)y=bx0=By=x

В общем, нельзя ожидать, что будет лежать в пространстве столбцов , но сгенерированная отправная точка будет точкой в ​​этом неясном пространстве, ближайшей к , в метрике, определяемой выбранными строками. Таким образом, это, вероятно, разумное приближение. По мере того, как обрабатывается больше систем, пространство столбцов увеличивается, и аппроксимация, вероятно, значительно улучшится, так что можно надеяться на сходимость за меньшее и меньшее количество итераций.xBx

Edit2: О сгенерированном подпространстве: если каждый решает каждую систему методом Крылова, векторы, используемые для получения начальной точки для второй системы, охватывают подпространство Крылова первой правой части. Таким образом, можно получить хорошее приближение всякий раз, когда это подпространство Крылова содержит вектор, близкий к решению вашей второй системы. В общем, векторы, используемые для получения начальной точки -й системы, охватывают пространство, содержащее подпространство Крылова первых правых частей.(k+1)k

Арнольд Ноймайер
источник
Спасибо, профессор Ноймайер. Я попробую это. Не могли бы вы дать мне краткое объяснение того, как это работает?
Костис
Кроме того, что, если я хочу решить одну и ту же систему для множества различных векторов RHS? т.е.Ax0=b0, Ax1=b1, Ax2=b2и т.д. Есть ли какая-либо информация, которую я могу использовать из предыдущих решений, чтобы ускорить последующие?
Костис
@Costis: разрешение с той же матрицей - это особый случай общей проблемы. По первому вопросу смотрите редактирование. D=0
Арнольд Ноймайер
@Costis: я добавил немного больше деталей к шагу 2. - Если вы пишете заявку, пожалуйста, пришлите мне препринт.
Арнольд Ноймайер
Спасибо за объяснение! Почему я не могу просто решить переопределенную систему(C+DB)ybиспользуя подход на основе факторизации QR и используя все строки, а не только подмножество? Я думаю, что по мере увеличения числа столбцов C и B, возможно, придется избавиться от некоторых строк, чтобы ускорить выполнение операции. Конечно, я напишу описание системы и отправлю вам по электронной почте. Я действительно думаю, что можно придумать схему для конкретного приложения, которая могла бы работать лучше, чем самый общий случай. Спасибо!
Костис