Метод двусопряженной градиентной стабилизации

Эта статья может быть слишком технической для понимания большинством читателей . Пожалуйста, помогите улучшить его, чтобы он был понятен неспециалистам , не удаляя технических деталей. ( Май 2015 г. ) ( Узнайте, как и когда удалить этот шаблон сообщения )

В численной линейной алгебре метод двусопряженной градиентной стабилизации , часто обозначаемый аббревиатурой BiCGSTAB , представляет собой итерационный метод, разработанный HA van der Vorst для численного решения несимметричных линейных систем . Это вариант метода двусопряженных градиентов (BiCG), который имеет более быструю и плавную сходимость, чем исходный BiCG, а также другие варианты, такие как метод квадрата сопряженных градиентов (CGS). Это метод подпространств Крылова .

Алгоритмические шаги [ править ]

Безусловный BiCGSTAB [ править ]

Чтобы решить линейную систему Ax = b , BiCGSTAB начинает с первоначального предположения x ₀ и действует следующим образом:

1. r ₀ = b - Ax ₀
2. Выберем произвольный вектор r̂ ₀ такой, что ( r̂ ₀ , r ₀ ) ≠ 0 , например, r̂ ₀ = r ₀ . ( x , y ) обозначает скалярное произведение векторов ( x , y ) = x ^T y
3. ρ ₀ = α = ω ₀ = 1
4. v ₀ = p ₀ = 0
5. Для i = 1, 2, 3,…
   1. ρ _i = ( r̂ ₀ , r _{i −1} )
   2. β = ( ρ _i / ρ _{i −1} ) ( α / ω _{i −1} )
   3. p _i = r _{i −1} + β ( p _{i −1} - ω _{i −1} v _{i −1} )
   4. v _i = Ap _i
   5. α = ρ _i / ( r̂ ₀ , v _i )
   6. h = x _{i −1} + α p _i
   7. Если h достаточно точное, установите x _i = h и выйдите
   8. s = r _{i −1} - α v _i
   9. t = Как
   10. ω _i = ( t , s ) / ( t , t )
   11. x _i = h + ω _i s
   12. Если x _i достаточно точен, то выйти
   13. r _i = s - ω _i t

Предварительно подготовленный BiCGSTAB [ править ]

Предобуславливатели обычно используются для ускорения сходимости итерационных методов. Чтобы решить линейную систему Ax = b с предварительным условием K = K ₁ K ₂ ≈ A , предварительное кондиционирование BiCGSTAB начинается с исходного предположения x ₀ и выполняется следующим образом:

1. r ₀ = b - Ax ₀
2. Выберем произвольный вектор r̂ ₀ такой, что ( r̂ ₀ , r ₀ ) ≠ 0 , например, r̂ ₀ = r ₀
3. ρ ₀ = α = ω ₀ = 1
4. v ₀ = p ₀ = 0
5. Для i = 1, 2, 3,…
   1. ρ _i = ( r̂ ₀ , r _{i −1} )
   2. β = ( ρ _i / ρ _{i −1} ) ( α / ω _{i −1} )
   3. p _i = r _{i −1} + β ( p _{i −1} - ω _{i −1} v _{i −1} )
   4. у = К −1
      2 п _я
   5. v _i = Ay
   6. α = ρ _i / ( r̂ ₀ , v _i )
   7. h = x _{i −1} + α y
   8. Если h достаточно точно, то x _i = h и выйти
   9. s = r _{i −1} - α v _i
   10. z = K −1
       2 s
   11. t = Az
   12. ω _i = ( K −1
       1 т , К −1
       1 с ) / ( К −1
       1 т , К −1
       1 т )
   13. х _я = h + ω _я z
   14. Если x _i достаточно точен, выйдите
   15. r _i = s - ω _i t

Эта формулировка эквивалентна применению безусловного BiCGSTAB к явно предварительно подготовленной системе.

Ãx̃ = b̃

с Ã = K −1
1 А К −1
2 , x̃ = K ₂ x и b̃ = K −1
1 б . Другими словами, с этой формулировкой возможны как левое, так и правое предварительное кондиционирование.

Вывод [ править ]

BiCG в полиномиальной форме [ править ]

В BiCG направления поиска p _i и p̂ _i и остатки r _i и r̂ _i обновляются с использованием следующих рекуррентных соотношений :

р _я = г _{я -1} + β _я п _{я -1} ,

p̂ _i = r̂ _{i −1} + β _i p̂ _{i −1} ,

r _i = r _{i −1} - α _i Ap _i ,

r̂ _i = r̂ _{i −1} - α _i A ^T p̂ _i .

Константы α _i и β _i выбраны равными

α _i = ρ _i / ( p̂ _i , Ap _i ) ,

β _i = ρ _i / ρ _{i −1}

где ρ _i = ( r̂ _{i −1} , r _{i −1} ), так что остатки и направления поиска удовлетворяют биортогональности и двусопряженности соответственно, т. е. для i ≠ j ,

( r̂ _i , r _j ) = 0 ,

( p̂ _i , Ap _j ) = 0 .

Несложно показать, что

r _i = P _i ( A ) r ₀ ,

r̂ _i = P _i ( A ^T ) r̂ ₀ ,

р _{я +1 знак} равно Т _я ( А ) г ₀ ,

p̂ _{i +1 знак} равно T _i ( A ^T ) r̂ ₀

где Р _я ( ) и Т _я ( ) являются I - й степени-многочлены A . Эти многочлены удовлетворяют следующим рекуррентным соотношениям:

P _i ( A ) = P _{i −1} ( A ) - α _i A T _{i −1} ( A ) ,

Т _я ( А ) = Р _я ( А ) + β _{я +1} Т _{я -1} ( А ) .

Получение BiCGSTAB от BiCG [ править ]

Нет необходимости явно отслеживать остатки и направления поиска BiCG. Другими словами, итерации BiCG могут выполняться неявно. В BiCGSTAB желательно иметь рекуррентные соотношения для

r̃ _i = Q _i ( A ) P _i ( A ) r ₀

где Q _i ( A ) = ( I - ω ₁ A ) ( I - ω ₂ A ) ⋯ ( I - ω _i A ) с подходящими константами ω _j вместо r _i = P _i ( A ) в надежде, что Q _i ( A ) обеспечит более быструю и плавную сходимость по r̃ _i, чем по r _i .

Из рекуррентных соотношений для P _i ( A ) и T _i ( A ) и определения Q _i ( A ) следует, что

Q _i ( A ) P _i ( A ) r ₀ = ( I - ω _i A ) ( Q _{i −1} ( A ) P _{i −1} ( A ) r ₀ - α _i A Q _{i −1} ( A ) T _{i −1} ( A ) r ₀ ) ,

что влечет за собой необходимость рекуррентного соотношения для Q _i ( A ) T _i ( A ) r ₀ . Это также может быть получено из отношений BiCG:

Q _i ( A ) T _i ( A ) r ₀ = Q _i ( A ) P _i ( A ) r ₀ + β _{i +1} ( I - ω _i A ) Q _{i −1} ( A ) P _{i −1} ( A ) г ₀ .

Аналогично определению r̃ _i , BiCGSTAB определяет

p̃ _{я +1 знак} равно Q _я ( А ) Т _я ( А ) г ₀ .

Записанные в векторной форме рекуррентные соотношения для p̃ _i и r̃ _i имеют вид

p̃ _i = r̃ _{i −1} + β _i ( I - ω _{i −1} A ) p̃ _{i −1} ,

r̃ _i = ( I - ω _i A ) ( r̃ _{i −1} - α _i A p̃ _i ).

Чтобы вывести рекуррентное соотношение для x _i , определим

s _i = r̃ _{i −1} - α _i A p̃ _i .

Тогда рекуррентное соотношение для r̃ _i можно записать как

r̃ _i = r̃ _{i −1} - α _i A p̃ _i - ω _i As _i ,

что соответствует

х _{я знак} равно х _{я -1} + α _я п _{̃ я} + ω _я с _я .

Определение констант BiCGSTAB [ править ]

Теперь осталось определить константы BiCG α _i и β _i и выбрать подходящий ω _i .

В BiCG β _i = ρ _i / ρ _{i −1} с

ρ _i = ( r̂ _{i −1} , r _{i −1} ) = ( P _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) .

Поскольку BiCGSTAB явно не отслеживает r̂ _i или r _i , ρ _i не вычисляется сразу по этой формуле. Однако это может быть связано со скалярным

ρ̃ _i = ( Q _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) = ( r̂ ₀ , Q _{i −1} ( A ) P _{i −1} ( A ) r ₀ ) = ( r̂ ₀ , r _{i −1} ) .

Благодаря биортогональности, г _{я -1} = Р _{я -1} ( ) г ₀ ортогонален U _{я -2} ( ^Т ) Г ₀ , где U _{я -2} ( ^Т ) является любой многочлен степени я - 2 в А ^T . Следовательно, только члены высшего порядка P _{i −1} ( A ^T ) и Q _{i −1} ( A^T ) имеют значение в скалярных произведениях ( P _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) и ( Q _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) . Старшие коэффициенты P _{i −1} ( A ^T ) и Q _i−1 ( A ^T ) равны (−1) ^{i −1} α ₁ α ₂ ⋯ α _{i −1} и (−1) ^{i −1} ω ₁ ω ₂ ⋯ ω _{i −1} соответственно. Следует, что

ρ _i = ( α ₁ / ω ₁ ) ( α ₂ / ω ₂ ) ⋯ ( α _{i −1} / ω _{i −1} ) ρ̃ _i ,

и поэтому

β _i = ρ _i / ρ _{i −1} = ( ρ̃ _i / ρ̃ _{i −1} ) ( α _{i −1} / ω _{i −1} ) .

Аналогичным образом может быть получена простая формула для α _i . В BiCG,

α _i = ρ _i / ( p̂ _i , Ap _i ) = ( P _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) / ( T _{i −1} ( A ^T ) r̂ ₀ , A T _{i −1} ( A ) r ₀ ) .

Как и в случае выше, только члены высшего порядка P _{i -1} ( A ^T ) и T _{i -1} ( A ^T ) имеют значение в скалярных произведениях благодаря биортогональности и двусопряженности. Бывает, что P _{i −1} ( A ^T ) и T _{i −1} ( A ^T ) имеют одинаковый старший коэффициент. Таким образом, их можно заменить одновременно на Q _{i −1} ( A ^T ) в формуле, что приводит к

α _i = ( Q _{i −1} ( A ^T ) r̂ ₀ , P _{i −1} ( A ) r ₀ ) / ( Q _{i −1} ( A ^T ) r̂ ₀ , A T _{i −1} ( A ) r ₀ ) = ρ̃ _i / ( r̂ ₀ , A Q _{i −1} ( A ) T _{i −1} (А ) r ₀ ) = ρ̃ _i / ( r̂ ₀ , Ap̃ _i ).

Наконец, BiCGSTAB выбирает ω _i, чтобы минимизировать r̃ _i = ( I - ω _i A ) s _i в 2- норме как функцию от ω _i . Это достигается, когда

(( I - ω _i A ) s _i , As _i ) = 0 ,

давая оптимальное значение

ω _i = ( As _i , s _i ) / ( As _i , As _i ) .

Обобщение [ править ]

BiCGSTAB можно рассматривать как комбинацию BiCG и GMRES, где за каждым шагом BiCG следует шаг GMRES ( 1 ) (т.е. GMRES перезапускается на каждом шаге) для исправления нерегулярного поведения конвергенции CGS, в качестве улучшения которого был разработан BiCGSTAB . Однако из-за использования полиномов минимальной остаточной степени один такое исправление может быть неэффективным, если матрица A имеет большие комплексные собственные пары. В таких случаях BiCGSTAB, вероятно, будет стагнировать, что подтверждается численными экспериментами.

Можно ожидать, что минимальные остаточные многочлены более высокой степени могут лучше справиться с этой ситуацией. Это дает начало алгоритмам, включая BiCGSTAB2 ^[1] и более общий BiCGSTAB ( l ) ^[2] . В BiCGSTAB ( l ) шаг GMRES ( l ) следует за каждыми l шагами BiCG. BiCGSTAB2 эквивалентен BiCGSTAB ( l ) с l = 2 .

См. Также [ править ]

• Метод двусопряженного градиента
• Метод сопряженного градиента в квадрате
• Метод сопряженных градиентов

Ссылки [ править ]

• Ван дер Ворст, HA (1992). «Bi-CGSTAB: быстрый и плавно сходящийся вариант Bi-CG для решения несимметричных линейных систем». SIAM J. Sci. Стат. Comput. 13 (2): 631–644. DOI : 10.1137 / 0913035 . hdl : 10338.dmlcz / 104566 .
• Саад, Ю. (2003). «§7.4.2 BICGSTAB» . Итерационные методы для разреженных линейных систем (2-е изд.). СИАМ. С.  231–234 . DOI : 10.2277 / 0898715342 . ISBN 978-0-89871-534-7.
• ^ Гуткнехт, MH (1993). «Варианты BICGSTAB для матриц с комплексным спектром». SIAM J. Sci. Comput. 14 (5): 1020–1033. DOI : 10.1137 / 0914062 .
• ^ Sleijpen, GLG; Фоккема, Д.Р. (ноябрь 1993 г.). «BiCGstab ( l ) для линейных уравнений, включающих несимметричные матрицы со сложным спектром» (PDF) . Электронные транзакции по численному анализу . Кент, Огайо: Государственный университет Кента. 1 : 11–32. ISSN 1068-9613 . Архивировано из оригинального (PDF) 06.06.2011 . Проверено 7 февраля 2010 .