Безопасное применение итерационных методов на диагонально-доминантных матрицах

9

Предположим, что задана следующая линейная система

(1)Lx=c,
где представляет собой взвешенное лапласиан , как известно, положительно определенной с одномерным нуль - пространство , натянутое на , а перевод дисперсия , т. е. не меняет значение функции (производная которой ). Единственные положительные элементы в находятся на его диагонали, которая является суммой абсолютных значений отрицательных недиагональных элементов.Lsemi1n=(1,,1)RnxRnx+a1n(1)L

Я обнаружил в одной высоко цитируемой научной работе в этой области, что, хотя является диагонально доминирующим, такие методы, как Conjugate Gradient, Gauss-Seidl, Jacobi, все еще можно безопасно использовать для решения . Обоснование состоит в том, что из-за инвариантности перевода можно безопасно зафиксировать одну точку (например, удалить первую строку и столбец и первую запись из ), тем самым преобразовав в диагонально доминирующую матрицу. В любом случае, исходная система решается в полной форме с .Lnot strictly(1)LcLstrictly(1)LRn×n

Верно ли это предположение, и, если да, каково альтернативное обоснование? Я пытаюсь понять, как сближение методов все еще держится.

Если метод Якоби сходится с (1)Что можно сказать о спектральном радиусе? ρ итерационной матрицы D1(DL), где D диагональная матрица с записями Lпо диагонали? Являетсяρ(D1(DL)1Таким образом, отличается от общих гарантий сходимости для ρ(D1(DL))<1? Я спрашиваю это, так как собственные значения матрицы ЛапласаD1Lс теми, которые по диагонали должны быть в диапазоне[0,2],

Из оригинальной работы:

......................................

На каждой итерации мы вычисляем новый макет (x (t +1), y (t + 1)), решая следующую линейную систему:

(8)L·x(t+1)=L(x(t),y(t))·x(t)L·y(t+1)=L(x(t),y(t))·y(t)
Без потери общности мы можем зафиксировать местоположение одного из датчиков (используя степень перемещения локализованного напряжения) и получить строго диагонально доминирующую матрицу. Поэтому мы можем смело использовать итерацию Якоби для решения (8)

.......................................

Выше понятие «итерация» относится к базовой процедуре минимизации, и ее не следует путать с итерацией Якоби. Итак, система решается Якоби (итеративно), и затем решение покупается в правой части (8), но теперь для другой итерации базовой минимизации. Я надеюсь, что это проясняет вопрос.

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

Usero
источник
Не могли бы вы опубликовать ссылку или ссылку на цитируемую работу?
Джефф Оксберри
Его можно получить по адресу : citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.164.1421. Поскольку вы не должны читать всю работу, посмотрите на стр.7 (внизу). Я полагаю, что выбор итерационных решателей оправдан, но я считаю, что требуется лучшее (или, по крайней мере, другое) обоснование.
Usero
Интересно, являются ли эти ребята из того же сообщества, что и комбинаторные прекондиционеры?
shuhalo

Ответы:

5

Итерация Якоби может быть доказана сходящейся.

Первое, что вы должны убедиться, что cT1n=0, что является условием существования решения (полагаю L=LTиначе тебе нужно c(KerLT)) потому что ты сказал V0:=KerL=span{1n}, Мы будем использовать соглашение, котороеV0также является матрицей со столбцами, являющимися ее ортонормированной основой. В твоем случае,V0:=1n/n,

Тогда для ошибок итерации Якоби в исходной системе вы

e1=(ID1L)e0=(ID1L)(Pe0+V0a)=(ID1L)Pe0+V0a,
где P:=IV0V0 это ортогональная проекция на V1:=V0, Из приведенной выше итерации мы знаем, что
Pe1=P(ID1L)Pe0,

из которой у нас есть итерационная матрица S в V1,
S:=P(ID1L)P.
Не то S имеет одинаковые спектры (кроме нулей) со следующей матрицей
S~:=(ID1L)PP=(ID1L)P=(ID1L)(IV0V0)=ID1LV0V0.
Мы хотим, чтобы спектральный радиус S меньше, чем один, чтобы доказать сходимость.

Следующая цитата старая и хранится только для справки. Смотрите после для нового доказательства.

В твоем случае, V0V0=1n1n×n. И вы можете убедиться, что D1L+V0V0 строго по диагонали с использованием предположения, что записи Lположительны по диагонали и отрицательны в противном случае. Чтобы показать собственные значения D1L+V0V0 реальны, отметим, что матрица самосопряжена относительно внутреннего произведения <x,y>:=yTDx.

Если V0не в вашей конкретной форме, я не нашел ответа на вопрос о конвергенции. Может кто-то это прояснить?

Обратите внимание, что V0 собственный вектор, соответствующий собственному значению 1 из ID1L, Основываясь на наблюдении, мы называем теорему 2.1 из собственных значений обновленных матриц ранга 1 с некоторыми приложениями Джиу Дина и Ай-Хуэй Чжоу.

Теорема 2.1 Пустьu а также v быть двумя nвекторы столбцов такие, что u является собственным вектором A связано с собственным значением λ1, Тогда собственные значенияA+uvT находятся {λ1+uTv,λ2,,λn} считая алгебраическую кратность.

Тогда мы знаем, что спектры S~ такой же как ID1L кроме того, что собственное значение 1 в последнем сдвигается на 1в собственное значение ноль в первом. посколькуρ(ID1L)(1,1], у нас есть ρ(S~)(1,1),

Хуэй Чжан
источник
Спасибо за ответы. Нечто похожее было то, что я рассмотрел: а именно, с взвешенным лапласианом, структурированным какD1L выше, можно показать, что его собственные значения находятся в пределах [0,2)следовательно, со спектральным радиусом в пределах (0,2) (одно собственное значение больше 0и по крайней мере один 0). Следовательно, спектральный радиус итерационной матрицыID1L меньше чем 1следовательно, со сходящимся Якоби. Возможно, приведенное выше предположение о спектральном радиусеID1L (без учета 0) не безопасно?
Usero
Я думаю, что спектры D1L должен быть в [0,2]закрыто на 2, Я не знаю, как вы можете получить2Исключенный. С моей точки зрения, (теорема Гершгорина о круге) [ en.wikipedia.org/wiki/Gershgorin_circle_theorem] может дать только оценку, включающую2, Если это так, то оценка спектрального радиусаID1L является 1 с равенством, достижимым с векторами в ядре L, Я думаю, что сходимость, которую вы хотите, заключается в том, что в пространстве ортогонального дополненияV1как отмечено в приведенном выше ответе.
Хуэй Чжан
Вы можете взглянуть на лемму 1.7 (v) из math.ucsd.edu/~fan/research/cb/ch1.pdf Матрица D1L можно рассматривать как взвешенный лапласиан на полном графе, следовательно, с исключенным 2, Я предполагаю, что это достаточный аргумент для доказательства сходимости? ........... Требует ли ваш подход другой предварительной / последующей обработки итераций, кроме центрирования?c, Я спрашиваю, потому что вы представилиV0 И относительно спектров ID1LV0V0: учитывая, что спектральный радиус (sr) из ID1L является (0,1], Добавление 1n, даст sr<1, Разве это не достаточно хороший аргумент?
Usero
Привет, спасибо за указание на хорошую книгу. Но я обнаружил, что не могу быстро взглянуть. Что касается вашего последнего аргумента, он выглядит почти так же, как «ответ» выше. Просто будьте осторожны, вы не добавляете1n но 1n1n×nтак что это не простое дополнение к sr из ID1L, Как правило,srсуммы двух матриц не являются простой суммойsrс отдельных матриц.
Хуэй Чжан
Хорошо, что ты указал на это. Требует ли ваш подход другой предварительной / последующей обработки итераций за пределами центровки c. Я спрашиваю, потому что вы представилиV0и я подумал, что вы говорите о проецировании нулевого пространства. Если это так, является ли проекция пустым пространством действительно необходимой для сходимости?
Usero
5

Методы Крылова никогда явно не используют размерность пространства, в котором они итерируются, поэтому вы можете запускать их в особых системах, пока вы сохраняете итерации в ненулевом подпространстве. Обычно это делается путем проецирования пустого пространства на каждой итерации. Есть две вещи, которые могут пойти не так: первая встречается гораздо чаще, чем вторая.

  1. Предусловие неустойчиво применительно к сингулярному оператору. Прямые решатели и неполная факторизация могут обладать этим свойством. На практике мы просто выбираем разные предварительные кондиционеры, но есть более принципиальные способы создания предварительных кондиционеров для единичных систем, например, Zhang (2010). .
  2. На некоторой итерации, x находится в ненулевом подпространстве, но Axживет полностью в нулевом пространстве. Это возможно только с несимметричными матрицами. Немодифицированный GMRES ломается в этом сценарии, но см. Reichel и Ye (2005) для вариантов без поломок.

Для решения особых систем с использованием PETSc см. KSPSetNullSpace() . Большинство методов и предварительных кондиционеров могут решать особые системы. На практике малое нулевое пространство для PDE с граничными условиями Неймана почти никогда не является проблемой, если вы сообщаете решателю Крылова о нулевом пространстве и выбираете разумный предобусловливатель.

Судя по комментариям, вас особенно интересует Якоби. (Почему? Якоби полезен как многосеточный сглаживатель, но есть гораздо лучшие методы для использования в качестве решателей.) Якоби применил кAx=b не сходится, когда вектор b имеет компонент в нулевом пространстве Aоднако часть решения, ортогональная нулевому пространству, действительно сходится, поэтому, если вы проецируете нулевое пространство из каждой итерации, она сходится. В качестве альтернативы, если вы выбираете последовательныйbи начальное предположение, итерации (в точной арифметике) не накапливают компоненты в нулевом пространстве.

Джед браун
источник
Вы можете выполнить ортогональное изменение базиса так, чтобы на диагонали был ноль (найдите любую ортогональную матрицу Qв котором первый столбец является постоянным вектором). Под этим преобразованиемA1=QTAQ, матрица A1все еще остается симметричным положительным полуопределенным, но первая диагональная запись равна 0, поэтому прямое применение Якоби не удастся. посколькуA1плотный, вы не будете делать это на практике, но это показывает, что основа имеет значение. ЕслиZ является ортогональным базисом для нулевого пространства, спроектированный GMRES просто решает (IZ)P1Ax=(IZ)P1b.
Jed Brown
Hmm, it seems I replied to a comment that was deleted. I'll leave the comment here in case it's useful.
Jed Brown
Thanks for the answer, it's on much higher specialized level then I expected. Therefore, I'll need some guides on: 1) how to project out the null space at each iteration? 2) In my understanding, you stated that the Jacobi application to the system as stated primarily might not converge to the exact solution (i.e. the iterands are not getting better solution estimates). It is therefore suggested to choose different preconditioners? If so, does that practically imply a dynamic check on behaviour with diag(A), and change if problem occurs (with the above case of the linear system)?
usero
My 1) from above should be regarded as: given the Jacobi iteration with the system primarily posted, is it needed to project out the nullspace, and, if so, how could one incorporate it within the update Xk+1=D1(b(AD)Xk)? Postprocessing the iterate Xk+1, and considering the postprocessed version for Xk?
usero
1
In a reasonable basis, Jacobi should be stable. It can also use 1 on the diagonal if the diagonal matrix element is 0, the projection still removes the null space. Are you planning to use a Krylov method like CG or GMRES? If not, why not? If you are, you just need an orthogonal basis for the null space. You only have the constant mode in your null space, so an orthogonal projector into the null space is N=ZZT where Z is the column vector. The orthogonal projector that removes the null space is thus IN. (My first comment had a mistake, if Z is the basis, N=IZZT is the projector.)
Jed Brown