Участник:EasyBreezy/Стабилизированный метод бисопряженных градиентов (BiCGSTAB)

Основные авторы описания:

Свойства и структура алгоритма

Общее описание алгоритма

Стабилизированный метод бисопряженных градиентов, так же известный как BiCGSTAB (Biconjugate gradient stabilized method), был предложен голландским математиком Хенком ван дер Ворстом в 1992 году. Статья, в которой описывался метод, стала самым цитируемым научным трудом в области математики в 1990х годах.

Метод появился, как ответ на задачу о построении способа решения систем линейных алгебраических уравнений произвольного вида. На момент появления BiCGSTAB уже существовали методы решения СЛАУ произвольного типа, в том числе и в комплексном пространстве, например метод бисопряженных градиентов. Однако существующие на тот момент методы либо обладали недостаточной скоростью сходимости, либо использовали в вычислениях тяжелые операции, такие как транспонирования матрицы, которые усложняли распараллеливание на многопроцессорных машинах.

BiCGSTAB является итерационным методом решения систем линейных алгебраических уравнений. Мощь метода заключается в том, что его можно применять к матрицам общего вида, в том числе с комплексными коэффициентами. В процессе выполнения BiCGSTAB не использует для вычисления плохо распараллеливаемых операций, не допускает накопления погрешностей округления и нестабильного поведения невязки. Эти факторы выгодно отличают его от предшественников: метода бисопряженных градиентов и квадратичного метода сопряженных градиентов, из которого BiCGSTAB фактически и вытекает. К преимуществам так же можно добавить более высокую в общем случае скорость сходимости относительно ближайших конкурентов. Тем не менее, для некоторых видов СЛАУ BiCGSTAB может уступать более простым методам.

Благодаря свойствам алгоритма область применения метода весьма широка. Фактически, он используется во всех задачах, где требуется решать СЛАУ, в том числе СЛАУ большой размерности, к которым сводится большое количество численных методов, применяемых при решении задач математического моделирования.

BiCGSTAB относится к классу так называемых проекционных методов. Основная идея методов данного класса: искать итерационным способом наилучшее приближение решения в некотором подпространстве пространства R^n, и чем шире подпространство, в котором ищется приближение решения, тем ближе будет приближенное решение к точному на каждой итерации. В качестве такого подпространства используются пространства Крылова, порожденные матрицей системы и невязкой первого приближения: [math]\displaystyle{ K^m (A, r^0) = span(r^0, Ar^0, A^2r^0, ..., A^{m-1} r^0). }[/math]

Изложим основную идею любого проекционного метода, заложенную в том числе и в BiCGSTAB. Для простоты выкладок и рассуждений здесь и далее будем считать, что матрица системы и ее правая часть вещественны. Без ограничения общности все дальнейшие рассуждения могут быть перенесены и на комплексный случай. Итак, пусть решается система [math]\displaystyle{ Ax=b }[/math] с невырожденной матрицей [math]\displaystyle{ A \in R^{n \times n} }[/math] и ненулевым вектором [math]\displaystyle{ b \in R^n }[/math]. Имеются два подпространства [math]\displaystyle{ K_n }[/math] и [math]\displaystyle{ L_n }[/math]. Требуется найти вектор [math]\displaystyle{ x \in K_n }[/math], который обеспечивает решение системы, оптимальное относительно подпространства L, то есть требуется выполнение условия Петрова-Галеркина: [math]\displaystyle{ (Ax,z)=(b,z), \forall z \in L_n }[/math]. В данном случае вектор x называется обобщенным решением системы.

Условие Петрова-Галеркина можно записать и через определение невязки. Пусть [math]\displaystyle{ r = b - Ax }[/math], тогда условие перепишется в виде: [math]\displaystyle{ (b - Ax, z)=(r, z),\forall z \in L_n }[/math]. Это равносильно тому, что вектор невязки ортогонален пространству L.

Пусть для СЛАУ известно некоторое приближение [math]\displaystyle{ x^0 }[/math] точного решения [math]\displaystyle{ x^* }[/math]. Данное приближение требуется уточнить поправкой [math]\displaystyle{ \delta: b - Ax^0 + \delta \perp L }[/math]. Невязка уточненного решения имеет вид: [math]\displaystyle{ r_{new} = b - Ax^0 + \delta = r^0- A\delta }[/math]. Условие Петрова-Галеркина тогда можно переписать так: [math]\displaystyle{ (r_{new}, z)= r^0 - (A\delta, z) = 0, \forall z \in L }[/math]

Пусть [math]\displaystyle{ dim K = dim L = m, V = (v_1, ..., v_n) \in R^{n \times m}, W = (w_1, ... w_n) \in R^{n \times m} }[/math], где V и W – матрицы, столбцы которых образуют базисы пространств K и L соответственно. Пусть так же [math]\displaystyle{ \delta = Vy, z = Wq }[/math], где y – коэффициенты разложения вектора поправки по базису пространства K, а q – коэффициенты разложения вектора z в пространстве L. Тогда новое приближение решения можно записать в виде: [math]\displaystyle{ x^1 = x^0 + \delta= x^0 + Vy }[/math]. Из условия ортогональности Петрова-Галеркина получаем, что: [math]\displaystyle{ (r - AVy, Wq) = (W^Tr^0 - AVy, q) = 0 }[/math]. Поскольку [math]\displaystyle{ q \neq 0 }[/math], то получаем СЛАУ для вектора y: [math]\displaystyle{ W^TAVy = W^Tr^0 }[/math]. Если предположить, что матрица [math]\displaystyle{ W^TAVy }[/math] невырождена, то вектор y однозначно определяется и равен: [math]\displaystyle{ y = W^TAV^{-1}W^Tr^0 }[/math]. Подставляя найденные коэффициенты разложения в выражение для первого приближения получаем: [math]\displaystyle{ x^1 = x^0 + \delta = x^0 + W^TAV^{-1}W^Tr^0y }[/math]

Подведем итог. Общий порядок шагов любого проекционного метода следующий:

1. Выбрать пару подпространств K и L
2. Найти базис V пространства K и базис W пространства L;
3. Вычислить невязку
4. Найти [math]\displaystyle{ y = W^T AV^{-1} W^T r^0 }[/math]
5. Положить [math]\displaystyle{ x = x + Vy }[/math]
6. Если не достигнут критерий остановки итерационного процесса, то повторить алгоритм, начиная с первого шага.

Как будет показано далее, в BiCGSTAB не требуется вычислять матрицу W^T AV^{-1} в явном виде.