Метод стабилизации бисопряженного градиента
 | Эта статья может быть слишком технической для понимания большинства читателей . ( Май 2015 г. ) |
В числовой линейной алгебре , метод стабилизации бисопряженного градиента часто сокращаемый как BiCGSTAB , представляет собой итерационный метод, разработанный Х.А. ван дер Ворстом для численного решения несимметричных линейных систем . Это вариант метода бисопряженных градиентов (BiCG), который имеет более быструю и плавную сходимость, чем исходный BiCG, а также другие варианты, такие как метод квадратов сопряженных градиентов (CGS). Это метод подпространств Крылова . В отличие от исходного метода BiCG, он не требует умножения на транспонирование системной матрицы.
Алгоритмические шаги
[ редактировать ]Непредварительный BiCGSTAB
[ редактировать ]В следующих разделах ( x , y ) = x Т y обозначает скалярное произведение векторов. Чтобы решить линейную систему Ax = b , BiCGSTAB начинает с начального предположения x 0 и действует следующим образом:
- р 0 = б - Ах 0
- Выберите произвольный вектор r̂ 0 такой, что ( r̂ 0 , r 0 ) ≠ 0 , например, r̂ 0 = r 0
- ρ 0 = ( r̂ 0 , r 0 )
- р 0 = р 0
- Для i = 1, 2, 3, …
- v = Ар я -1
- α = ρ i −1 /( r̂ 0 , v )
- час знак равно Икс я -1 + α п я -1
- s знак равно р я -1 - α v
- Если h достаточно точен, т. е. если s достаточно мало, установите x i = h и выйдите из программы.
- т = Как
- ω знак равно ( т , s )/( т , т )
- х я = час + ω с
- р я знак равно s - ω т
- Если x i достаточно точен, т. е. если r i достаточно мало, то закройте
- ρ я знак равно ( р̂ 0 , р я )
- β знак равно ( ρ я / ρ я -1 )( а / ω )
- п я знак равно р я + β ( п я -1 - ω v
В некоторых случаях выбор вектора r̂ 0 случайным образом улучшает численную стабильность. [1]
Предварительно подготовленный BiCGSTAB
[ редактировать ]Предварительные условия обычно используются для ускорения сходимости итерационных методов. Чтобы решить линейную систему Ax = b с предобусловливателем K = K 1 K 2 ≈ A , предобуславливание BiCGSTAB начинается с начального предположения x 0 и продолжается следующим образом:
- р 0 = б - Ах 0
- Выберите произвольный вектор r̂ 0 такой, что ( r̂ 0 , r 0 ) ≠ 0 , например, r̂ 0 = r 0
- ρ 0 = ( r̂ 0 , r 0 )
- р 0 = р 0
- Для i = 1, 2, 3, …
- у = К -1
2К -1
1 п я -1 - v = Есть
- α = ρ i −1 /( r̂ 0 , v )
- час знак равно х я -1 + α у
- s знак равно р я -1 - α v
- Если h достаточно точен, то x i = h и выходим.
- г = К -1
2К -1
1 с - т = Это
- ω = ( К −1
1 т , К −1
1 с )/( К −1
1 т , К −1
1 т ) - х я = час + ω z
- р я знак равно s - ω т
- Если x i достаточно точен, выйдите
- ρ я знак равно ( р̂ 0 , р я )
- β знак равно ( ρ я / ρ я -1 )( а / ω )
- п я знак равно р я + β ( п я -1 - ω v
- у = К -1
Эта формулировка эквивалентна применению непредобусловленного BiCGSTAB к явно предварительно обусловленной системе.
- Ãx̃ = b̃
с Ã = K −1
1 А К -1
2 , x̃ = K 2 x и b̃ = K −1
1 б . Другими словами, в этой формулировке возможна как левая, так и правая предобуславливание.
Вывод
[ редактировать ]BiCG в полиномиальной форме
[ редактировать ]BiCG направления поиска pi и и p̂i следующих остатки ri r̂i В : и обновляются с использованием отношений рекуррентных
- п я знак равно р я -1 + β я п я -1 ,
- p̂ я знак равно r̂ я −1 + б я p̂ я −1 ,
- р я знак равно р я -1 - α я Ap я ,
- р̂ я = р̂ я −1 − а я А Т π̂ я .
Константы α i и β i выбраны равными
- α я знак равно ρ я /( p̂ я , Ap я ) ,
- β я знак равно ρ я / ρ я -1
где ρ i = ( r̂ i −1 , r i −1 ), так что остатки и направления поиска удовлетворяют биортогональности и двусопряжённости соответственно, т. е. для i ≠ j ,
- ( р̂ я , р j ) знак равно 0 ,
- ( п̂ я , Ap j ) знак равно 0 .
Несложно показать, что
- р я знак равно п я ( А ) р 0 ,
- r̂ i = P i ( A Т ) р̂ 0 ,
- п я +1 знак равно Т я ( А ) р 0 ,
- p̂ i +1 = T i ( A Т ) р̂ 0
где P i ( A ) и Ti A ( ) . — i -й степени от A полиномы Эти полиномы удовлетворяют следующим рекуррентным соотношениям:
- п я ( А ) знак равно п я -1 ( А ) - α я А Т я -1 ( А ) ,
- Т я ( А ) знак равно п я ( А ) + β я +1 Т я -1 ( А ) .
Вывод BiCGSTAB из BiCG
[ редактировать ]Нет необходимости явно отслеживать остатки и направления поиска BiCG. Другими словами, итерации BiCG могут выполняться неявно. В BiCGSTAB желательно иметь рекуррентные отношения для
- р̃ я знак равно Q я ( А ) п я ( А ) р 0
где Q i ( A ) = ( I − ω 1 A )( I − ω 2 A )⋯ ( I − ω i A ) с подходящими константами ω j вместо r i = P i ( A ) r 0 в надежде, что Q i ( A ) обеспечит более быструю и плавную сходимость в r̃ i , чем в r i .
Из рекуррентных соотношений для Pi A ( ) что и Ti , ( A ) и определения Q i ( A ) следует
- Q я ( А ) п я ( А ) р 0 знак равно ( я - ω я А )( Q я -1 ( А ) п я -1 ( А ) р 0 - α я А Q я -1 ( А ) Т я -1 ( А ) р 0 ) ,
что влечет за собой необходимость рекуррентного соотношения для Q i ( A ) T i ( A ) r 0 . Это также можно вывести из соотношений BiCG:
- Q я ( А ) Т я ( А ) р 0 знак равно Q я ( А ) п я ( А ) р 0 + β я +1 ( я - ω я А ) Q я -1 ( А ) п я -1 ( А ) р 0 .
Аналогично определению r̃ i , BiCGSTAB определяет
- п̃ я +1 знак равно Q я ( А ) Т я ( А ) р 0 .
Записанные в векторной форме рекуррентные соотношения для p̃ i и r̃ i имеют вид
- п̃ я знак равно р̃ я -1 + β я ( я - ω я -1 А ) п̃ я -1 ,
- р̃ я знак равно ( я - ω я А )( р̃ я -1 - α я А п̃ я ) .
Чтобы получить рекуррентное соотношение для x i , определите
- s я знак равно р̃ я -1 - α я А п̃ я .
Тогда рекуррентное соотношение для r̃ i можно записать как
- r̃ я знак равно r̃ я -1 - α я А p̃ я - ω я Как я ,
что соответствует
- Икс я знак равно Икс я -1 + α я п̃ я + ω я s я .
Определение констант BiCGSTAB
[ редактировать ]Теперь осталось определить константы BiCG α i и β i и выбрать подходящую ω i .
В BiCG β i = ρ i / ρ i −1 с
- ρ я знак равно ( р̂ я -1 , р я -1 ) знак равно ( п я -1 ( А Т ) р̂ 0 , п я -1 ( А ) р 0 ) .
не отслеживает явно r̂i или . ri формуле , ρi Поскольку BiCGSTAB не может быть немедленно вычислено по этой Однако это может быть связано со скаляром
- ρρi = ( Q i −1 ( А Т ) р̂ 0 , п я -1 ( А ) р 0 ) знак равно ( р̂ 0 , Q я -1 ( А ) п я -1 ( А ) р 0 ) знак равно ( р̂ 0 , р я -1 ) .
В силу биортогональности r i −1 = P i −1 ( A ) r 0 ортогонально U i −2 ( A Т ) r̂ 0 , где U i −2 ( A Т ) — любой многочлен степени i − 2 из A Т . Следовательно, только члены высшего порядка P i −1 ( A Т ) и Q i −1 ( A Т ) материя в скалярном произведении ( P i −1 ( A Т ) r̂ 0 , п я -1 ( А ) р 0 ) и ( Q я -1 ( А Т ) р̂ 0 , п я -1 ( А ) р 0 ) . Старшие коэффициенты P i −1 ( A Т ) и Q i −1 ( A Т ) равны (−1) я -1 α 1 α 2 ⋯ α i −1 и (−1) я -1 ω 1 ω 2 ⋯ ω i −1 соответственно. Отсюда следует, что
- ρ я знак равно ( α 1 / ω 1 )( α 2 / ω 2 )⋯( α я -1 / ω я -1 ) ρ̃ я ,
и таким образом
- β я знак равно ρ я / ρ я -1 знак равно ( ρ̃ я / ρ̃ я -1 )( α я -1 / ω я -1 ) .
Простую формулу для αi . можно вывести аналогичным образом В БиКГ,
- α я знак равно ρ я /( p̂ i , Ap я ) знак равно ( п я -1 ( А Т ) r̂ 0 , п я -1 ( А ) р 0 )/( Т я -1 ( А Т ) р̂ 0 , А Т я -1 ( А ) р 0 ) .
Как и в случае выше, только члены высшего порядка P i −1 ( A Т ) и T i −1 ( A Т ) имеют значение в скалярных произведениях благодаря биортогональности и двусопряжению. Бывает, что P i −1 ( A Т ) и T i −1 ( A Т ) имеют одинаковый ведущий коэффициент. Таким образом, их можно заменить одновременно с Q i −1 ( A Т ) в формуле, что приводит к
- α i знак равно ( Q i −1 ( А Т ) r̂ 0 , п я -1 ( А ) р 0 )/( Q я -1 ( А Т ) r̂ 0 , А Т я -1 ( А ) р 0 ) знак равно ρ̃ я /( r̂ 0 , А Q я −1 ( А ) Т я −1 ( А ) р 0 ) знак равно ρ̃ я /( r̂ 0 , Ap̃ я )
Наконец, BiCGSTAB выбирает ω i для минимизации r̃ i = ( I − ω i A ) s i в 2 -норме как функции ω i . Это достигается, когда
- (( я - ω я А ) s я , As я ) знак равно 0 ,
дающий оптимальное значение
- ω я знак равно ( As я , s я )/( As я , As я ) .
Обобщение
[ редактировать ]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 .
См. также
[ редактировать ]Ссылки
[ редактировать ]- ^ Шутроп, Крис; Бунккамп, Ян тен Тидже; Дейк, Ян ван (июль 2022 г.). «Исследование надежности решателей BiCGStab и IDR для уравнения реакции адвекции-диффузии» . Коммуникации в вычислительной физике . 32 (1): 156–188. doi : 10.4208/cicp.oa-2021-0182 . ISSN 1815-2406 .
- Ван дер Ворст, ХА (1992). «Bi-CGSTAB: быстрый и плавно сходящийся вариант Bi-CG для решения несимметричных линейных систем». СИАМ J. Sci. Стат. Вычислить. 13 (2): 631–644. дои : 10.1137/0913035 . hdl : 10338.dmlcz/104566 .
- Саад, Ю. (2003). «§7.4.2 БИКГСТАБ» . Итерационные методы для разреженных линейных систем (2-е изд.). СИАМ. стр. 231–234 . ISBN 978-0-89871-534-7 .
- ^ Гуткнехт, МХ (1993). «Варианты BICGSTAB для матриц с комплексным спектром». СИАМ J. Sci. Вычислить. 14 (5): 1020–1033. дои : 10.1137/0914062 .
- ^ Слейпен, GLG; Фоккема, ДР (ноябрь 1993 г.). «BiCGstab( l ) для линейных уравнений, включающих несимметричные матрицы с комплексным спектром» (PDF) . Электронные труды по численному анализу . 1 . Кент, Огайо: Кентский государственный университет: 11–32. ISSN 1068-9613 .