Доклады Российской академии наук. Математика, информатика, процессы управления, 2020, T. 495, № 1, стр. 38-43

О МЕТОДАХ МОМЕНТОВ В ПОДПРОСТРАНСТВАХ КРЫЛОВА

В. П. Ильин^1, 2, *

¹ Институт вычислительной математики и математической геофизики Сибирского отделения Российской академии наук
Новосибирск, Россия

² Новосибирский государственный университет
Новосибирск, Россия

^* E-mail: ilin@sscc.ru

Поступила в редакцию 02.06.2020
После доработки 02.06.2020
Принята к публикации 11.11.2020

DOI: 10.31857/S2686954320060223

Полный текст (PDF)

Аннотация

Рассматриваются методы моментов в подпространствах Крылова для решения симметричных систем линейных алгебраических уравнений. Предложено семейство итерационных алгоритмов, основанное на обобщенной ортогонализации Ланцоша с выбором исходного вектора ${{v}^{0}}$ независимо от начальной невязки. Данный подход позволяет на одном наборе базисных векторов экономично решать серии систем линейных алгебраических уравнений с одинаковой матрицей, но разными правыми частями, а также реализовывать обобщенные методы моментов, сводящиеся к блочным крыловским алгоритмам с использованием совокупности линейно независимых исходных векторов $v_{1}^{0},...,v_{m}^{0}.$ Повышение производительности реализаций алгоритмов достигается за счет сокращения числа матричных умножений и эффективного распараллеливания векторных операций. Показана возможность расширения применимости методов моментов с использованием предобусловливания на различные классы алгебраических систем: знакопеременных, несовместных, несимметричных и комплексных, в том числе неэрмитовых.

Ключевые слова: методы моментов, подпространства Крылова, параметрическая ортогонализация Ланцоша, алгоритмы сопряженных направлений

В 1958 г. ленинградский математик Ю.В. Воробьёв опубликовал книгу о методах моментов в пространствах Галёркина [1], связанных со многими проблемами теории операторов и прикладной математики [2]. Следствием изложенных им результатов явилось обобщение метода сопряженных градиентов для решения систем линейных алгебраических уравнений (СЛАУ), предложенного в работах К. Ланцоша, М. Хестенеса и Е. Штифеля, см. монографии [3–6]. В частности, Ю.В. Воробьёвым построены экономичные алгоритмы для решения СЛАУ со знакопеременным спектром и разными правыми частями, но с одинаковой матрицей. В 2013 г. авторы книги [3] исследовали и применили методы моментов для решения комплексных СЛАУ, а также для редукции моделей многомасштабных динамических систем.

Целью данной работы является развитие метода моментов в подпространствах Крылова с обобщением на алгоритмы сопряженных направлений [5], а также на предобусловленные итерационные процессы для решения вещественных СЛАУ. Рассмотренные алгоритмы несложно переносятся на комплексные эрмитовы системы.

1. ПОЛИНОМИАЛЬНОЕ ПРЕДСТАВЛЕНИЕ ПРИБЛИЖЕННОГО РЕШЕНИЯ СЛАУ

Рассмотрим симметричную положительно определенную (с.п.о.) СЛАУ в виде

(1)

$Au = f,\quad A = {{A}^{ \top }} = {{\mathcal{R}}^{{N,N}}},\quad u,f \in {{\mathcal{R}}^{N}},$

где для векторов введено евклидово параметризованное скалярное произведение

${{(u,v)}_{\gamma }} \equiv ({{A}^{\gamma }}u,v)$ = ${{u}^{ \top }} \cdot {{A}^{\gamma }}v$, ${{(u,u)}_{\gamma }} = {\text{||}}u{\text{||}}_{{^{\gamma }}}^{2},$

причем γ равно 0, или 1, или 2.

Пусть для произвольного вектора ${{v}^{0}}$ совокупность линейно независимых векторов ${{v}^{{k + 1}}} = A{{v}^{k}},$ k = 0, 1, …, n – 1, образует n-мерное подпространство Крылова

(2)

${{\mathcal{K}}_{n}}({{v}^{0}},A) = \operatorname{Span} \{ {{v}^{0}},{{v}^{1}},...,{{v}^{{n - 1}}}\} ,$

так что элемент ${{E}_{n}}{{v}^{n}} = A{{v}^{{n - 1}}}$ есть проекция ${{v}^{n}}$ в подпространство ${{\mathcal{K}}_{n}}$, а E_n – соответствующий оператор проектирования в смысле введенного скалярного произведения (⋅, ⋅)_γ.

Сформулируем следующую проблему моментов: по заданному набору n + 1 линейно независимых элементов ${{v}^{0}},{{v}^{1}},...,{{v}^{n}}$ из гильбертова пространства ${{\mathcal{R}}^{N}}$ построить линейный оператор A_n, определенный в n-мерном подпространстве ${{\mathcal{K}}_{n}}$, удовлетворяющий условиям:

(3)

$\begin{gathered} {{v}^{1}} = {{A}_{n}}{{v}^{0}},...,{{v}^{{n - 1}}} = {{A}_{n}}{{v}^{{n - 2}}}, \\ {{E}_{n}}{{v}^{n}} = {{A}_{n}}{{v}^{{n - 1}}} = A_{n}^{n}{{v}^{0}}. \\ \end{gathered} $

Поскольку ${{E}_{n}}{{v}^{n}}$ есть проекция элемента ${{v}^{n}} \in {{\mathcal{R}}^{N}}$ в подпространство ${{\mathcal{K}}_{n}}$, то разность ${{v}^{n}} - {{E}_{n}}{{v}^{n}}$ ортогональна любому элементу ${{v}^{k}} \in {{\mathcal{K}}_{n}}$. Отсюда, используя представление

${{E}_{n}}{{v}^{n}} = - {{a}_{0}}{{v}^{0}} - {{a}_{1}}{{v}^{1}} - ... - {{a}_{{n - 1}}}{{v}^{{n - 1}}},$

после скалярного умножения обеих его частей на ${{A}^{\gamma }}{{v}^{k}}$, для вектора коэффициентов $a = {{({{a}_{0}},...,{{a}_{{n - 1}}})}^{ \top }}$ получаем систему

(4)

$Ga = w \equiv {{({{({{v}^{n}},{{v}^{0}})}_{\gamma }},...,{{({{v}^{n}},{{v}^{{n - 1}}})}_{\gamma }})}^{ \top }},$

с невырожденной матрицей Грама G = $\{ ({{v}^{k}}$, ${{v}^{l}}{{)}_{\gamma }}\} \in \mathcal{R}{{}^{{n,n}}}$. В силу этого СЛАУ (4) единственным образом определяет матричный многочлен n-й степени, аннулирующий элемент ${{v}^{0}}$:

(5)

${{P}_{n}}({{A}_{n}}){{v}^{0}} = (A_{n}^{n} + {{a}_{{n - 1}}}A_{n}^{{n - 1}} + ... + {{a}_{0}}I){{v}^{0}} = 0,$

а корни соответствующего скалярного полинома P_n(λ) являются собственными числами оператора A_n. Здесь и далее индекс γ у матриц, векторов и многочленов для краткости опускаем.

Рассмотрим теперь представление решения уравнения

(6)

${{A}_{n}}{{u}^{n}} = {{f}^{n}};\quad {{u}^{n}},{{f}^{n}} \in {{\mathcal{K}}_{n}},$

правая часть которого имеет вид

(7)

$\begin{gathered} {{f}^{n}} = \sum\limits_{k = 0}^{n - 1} {{{b}_{k}}{{v}^{k}}} = \sum\limits_{k = 0}^{n - 1} {{{b}_{k}}A_{n}^{k}{{v}^{0}}} = {{F}_{{n - 1}}}({{A}_{n}}){{v}^{0}}, \\ {{F}_{{n - 1}}}(\lambda ) = {{b}_{{n - 1}}}{{\lambda }^{{n - 1}}} + {{b}_{{n - 2}}}{{\lambda }^{{n - 2}}} + ... + {{b}_{0}}. \\ \end{gathered} $

Умножая (7) скалярно на ${{v}^{{_{0}}}},...,{{v}^{{_{{n - 1}}}}},$ формируем СЛАУ для вектора коэффициентов:

(8)

$\begin{gathered} Gb = \overline f \equiv {{({{({{v}^{0}}f)}_{\gamma }},{{({{v}^{1}}f)}_{\gamma }},...,{{({{v}^{{n - 1}}},f)}_{\gamma }})}^{ \top }}, \\ b = {{({{b}_{0}},{{b}_{1}},...,{{b}_{{n - 1}}})}^{ \top }}. \\ \end{gathered} $

Решение системы (6) может быть выражено с помощью многочлена

(9)

${{S}_{{n - 1}}}({{A}_{n}}) = A_{n}^{{ - 1}}\left[ {{{F}_{{n - 1}}}({{A}_{n}}) - {{{{P}_{n}}({{A}_{n}}){{F}_{{n - 1}}}(0)} \mathord{\left/ {\vphantom {{{{P}_{n}}({{A}_{n}}){{F}_{{n - 1}}}(0)} {{{P}_{n}}(0)}}} \right. \kern-0em} {{{P}_{n}}(0)}}} \right],$

в результате чего при условии P_n(0) = a₀ ≠ 0 приближенное решение проблемы моментов выражается с помощью коэффициентов многочленов a_k и b_k следующим образом:

(10)

$\begin{gathered} {{u}^{n}} = {{S}_{{n - 1}}}({{A}_{n}}){{v}^{0}} = \left( {{{b}_{1}} - {{{\bar {a}}}_{1}}} \right){{v}^{0}} + ... + \\ \, + \left( {{{b}_{{n - 1}}} - {{{\bar {a}}}_{{n - 1}}}} \right){{v}^{{n - 2}}} - {{{\bar {b}}}_{0}}{{v}^{{n - 1}}}, \\ {{{\bar {a}}}_{k}} = {{{{a}_{k}}{{b}_{0}}} \mathord{\left/ {\vphantom {{{{a}_{k}}{{b}_{0}}} {{{a}_{0}},}}} \right. \kern-0em} {{{a}_{0}},}}\quad {{{\bar {b}}}_{0}} = {{{{b}_{0}}} \mathord{\left/ {\vphantom {{{{b}_{0}}} {{{a}_{0}}.}}} \right. \kern-0em} {{{a}_{0}}.}} \\ \end{gathered} $

2. РЕКУРСИВНОЕ РЕШЕНИЕ ПРОБЛЕМЫ МОМЕНТОВ

Повышение экономичности алгоритма, представленного формулой (10), достигается путем ортогонализации векторов ${{v}^{0}},{{v}^{1}},\,\,...,{{v}^{n}}$, предложенной К. Ланцошем и лежащей в основе метода сопряженных градиентов. Рассмотрим этот подход в обобщенной форме, приводящей к методам сопряженных направлений [5]. Более конкретно, применяем параметрическую A^γ-ортогонализацию, т.е. строим совокупность векторов p^k, обладающих свойствами ${{({{p}^{k}},{{p}^{l}})}_{\gamma }}$ = ${{({{p}^{k}},{{p}^{k}})}_{\gamma }} \cdot {{\delta }_{{k,l}}},$ где γ = 0, 1, 2, а δ_k,l – символ Кронекера. Сами формулы имеют вид

(11)

$\begin{gathered} {{p}^{0}} = {{v}^{0}},\quad {{p}^{1}} = A{{v}^{0}} - \bar {\alpha }_{0}^{{(\gamma )}}{{v}^{0}};\quad k = 1,2,...: \\ {{p}^{{k + 1}}} = (A - \bar {\alpha }_{k}^{{(\gamma )}}I){{p}^{k}} - \bar {\beta }_{k}^{{(\gamma )}}{{p}^{{k - 1}}}, \\ \bar {\alpha }_{k}^{{(\gamma )}} = {{{{{({{p}^{k}},{{p}^{k}})}}_{{\gamma + 1}}}} \mathord{\left/ {\vphantom {{{{{({{p}^{k}},{{p}^{k}})}}_{{\gamma + 1}}}} {{\text{||}}{{p}^{k}}{\text{||}}}}} \right. \kern-0em} {{\text{||}}{{p}^{k}}{\text{||}}}}_{\gamma }^{2},\quad \bar {\beta }_{k}^{{(\gamma )}} = {{{\text{||}}{{p}^{k}}{\text{||}}_{\gamma }^{2}} \mathord{\left/ {\vphantom {{{\text{||}}{{p}^{k}}{\text{||}}_{\gamma }^{2}} {{\text{||}}{{p}^{{k - 1}}}{\text{||}}}}} \right. \kern-0em} {{\text{||}}{{p}^{{k - 1}}}{\text{||}}}}_{\gamma }^{2}. \\ \end{gathered} $

При этом векторы p^k выражаются через матричные многочлены, связанные между собой трехчленными соотношениями:

(12)

$\begin{gathered} {{p}^{k}} = {{P}_{k}}(A){{v}^{0}},\quad \beta _{{ - 1}}^{{(\gamma )}} = 0,\quad {{P}_{0}}(A) = 1, \\ {{P}_{{k + 1}}}(A) = (A - \bar {\alpha }_{k}^{{(\gamma )}}I){{P}_{k}}(A) - \bar {\beta }_{k}^{{(\gamma )}}{{P}_{{k - 1}}}(A). \\ \end{gathered} $

Поскольку векторы p⁰, …, pⁿ ^{– 1} образуют A^γ-ортогональный базис в ${{\mathcal{K}}_{n}}$, в силу свойств оператора проектирования E_n и рекурсий (12) имеем

${{P}_{n}}({{A}_{n}}){{v}^{0}} = {{E}_{n}}{{P}_{n}}(A){{v}^{0}} = {{E}_{n}}{{p}^{n}} = 0.$

Отсюда следует, что полученные полиномы P_n совпадают с многочленами из (5), корни уравнения P_n(λ) = 0 являются собственными числами оператора A_n, а в ортогональном базисе p⁰, …, pⁿ ^{– 1} оператор проектирования из ${{\mathcal{R}}^{N}}$ в ${{\mathcal{K}}_{n}}$ является ортогональным проектором, т.е.

${{E}_{n}} = {{V}_{n}}{{V}_{n}}{{\,}^{ \top }},\quad E_{n}^{2} = {{E}_{n}},\quad {{V}_{n}} \in {{\mathcal{R}}^{{N,n}}},$

где V_n есть матрица, столбцы которой составлены из векторов p^k.

Замечание 1. В работе [7] было предложено другое обобщение ортогонализации Ланцоша:

${{p}^{0}} = {{v}^{0}},\quad {{c}_{1}}{{p}^{1}} = B{{p}^{0}} - {{\alpha }_{0}}{{p}^{0}},\quad k = 2,3,...:$

(13)

$\begin{gathered} {{c}_{{k + 1}}}{{p}^{{k + 1}}} = B{{p}^{k}} - {{\alpha }_{k}}{{p}^{k}} - {{\beta }_{k}}{{p}^{{k - 1}}}, \\ {{\alpha }_{k}} = {{({{p}^{k}},CB{{p}^{k}})} \mathord{\left/ {\vphantom {{({{p}^{k}},CB{{p}^{k}})} {({{p}^{k}},C{{p}^{k}})}}} \right. \kern-0em} {({{p}^{k}},C{{p}^{k}})}}, \\ \end{gathered} $

${{\beta }_{k}} = {{({{p}^{{k - 1}}},CB{{p}^{{k - 1}}})} \mathord{\left/ {\vphantom {{({{p}^{{k - 1}}},CB{{p}^{{k - 1}}})} {({{p}^{{k - 1}}},C{{p}^{{k - 1}}})}}} \right. \kern-0em} {({{p}^{{k - 1}}},C{{p}^{{k - 1}}})}},$

где B и C – некоторые перестановочные симметричные матрицы, а c_k – произвольные ненулевые константы. Данные соотношения обеспечивают C-ортогональность векторов, т.е. $({{p}^{k}},C{{p}^{n}}) = 0$ при k ≠ n. Формулы (12) следуют из (13) при B = A и С = A^γ, а в [7] рассматривался более конкретно случай B = AA^⊤ = C для несимметричных СЛАУ.

Используя рекуррентные свойства построенных многочленов, проведем преобразования приближенного решения uⁿ из (10). Вводя обозначение ${{R}_{n}}(\lambda ) = {{{{P}_{n}}(\lambda )} \mathord{\left/ {\vphantom {{{{P}_{n}}(\lambda )} {{{P}_{n}}(0)}}} \right. \kern-0em} {{{P}_{n}}(0)}},$ для данных полиномов вследствие (12) получаем следующую рекурсию:

(14)

$\begin{gathered} {{R}_{{k + 1}}}(\lambda ) = (\lambda - \bar {\alpha }_{k}^{{(\gamma )}}){{R}_{k}}(\lambda ){{{{P}_{k}}(0)} \mathord{\left/ {\vphantom {{{{P}_{k}}(0)} {{{P}_{{k + 1}}}(0)}}} \right. \kern-0em} {{{P}_{{k + 1}}}(0)}} - \\ \, - \bar {\beta }_{{k - 1}}^{{(\gamma )}}{{{{R}_{{k - 1}}}(\lambda ){{P}_{{k - 1}}}(0)} \mathord{\left/ {\vphantom {{{{R}_{{k - 1}}}(\lambda ){{P}_{{k - 1}}}(0)} {{{P}_{{k + 1}}}(0)}}} \right. \kern-0em} {{{P}_{{k + 1}}}(0)}}. \\ \end{gathered} $

Отсюда с помощью обозначений

$\begin{gathered} {{\alpha }_{k}} = - {{{{P}_{k}}(0)} \mathord{\left/ {\vphantom {{{{P}_{k}}(0)} {{{P}_{{k + 1}}}(0)}}} \right. \kern-0em} {{{P}_{{k + 1}}}(0)}}, \\ {{\beta }_{k}} = \beta _{k}^{{(\gamma )}}{{R}_{{k - 1}}}(\lambda ){{P_{{k - 1}}^{2}(0)} \mathord{\left/ {\vphantom {{P_{{k - 1}}^{2}(0)} {P_{k}^{2}(0)}}} \right. \kern-0em} {P_{k}^{2}(0)}}, \\ {{Q}_{k}}(\lambda ) = {{ - [{{R}_{{k + 1}}}(\lambda ) - {{R}_{k}}(\lambda )]} \mathord{\left/ {\vphantom {{ - [{{R}_{{k + 1}}}(\lambda ) - {{R}_{k}}(\lambda )]} {{{\alpha }_{k}}\lambda }}} \right. \kern-0em} {{{\alpha }_{k}}\lambda }}, \\ \end{gathered} $

приходим к двучленным рекурсиям

(15)

$\begin{gathered} {{R}_{{k + 1}}}(\lambda ) = {{R}_{k}}(\lambda ) - \lambda {{\alpha }_{k}}{{Q}_{k}}(\lambda ), \\ {{Q}_{{k + 1}}}(\lambda ) = {{R}_{{k + 1}}}(\lambda ) + {{\beta }_{k}}{{Q}_{k}}(\lambda ). \\ \end{gathered} $

Вводя наряду с этими многочленами соответствующие элементы гильбертова пространства ${{r}^{k}} = {{R}_{k}}(A){{v}^{0}},$ ${{q}^{k}} = {{Q}_{k}}(A){{v}^{0}},$ получаем аналогичные векторные соотношения

(16)

$\begin{gathered} {{r}^{{k + 1}}} = {{r}^{k}} - {{\alpha }_{k}}A{{q}^{k}},\quad {{q}^{{k + 1}}} = {{r}^{{k + 1}}} + {{\beta }_{k}}{{q}^{k}}{\text{,}} \\ {{q}^{0}} = {{r}^{0}} = {{v}^{0}}. \\ \end{gathered} $

На основе свойств введенных выше полиномов приходим к следующим ортогональным свойствам полученных векторных последовательностей:

(17)

$\begin{gathered} {{({{r}^{k}},{{r}^{{k{\kern 1pt} '}}})}_{{\gamma - 1}}} = {{\sigma }_{k}}{{\delta }_{{k,k{\kern 1pt} '}}},\quad {{\sigma }_{k}} = {\text{||}}{{r}^{k}}{\text{||}}_{{\gamma - 1}}^{2}, \\ {{({{q}^{k}},{{q}^{{k{\kern 1pt} '}}})}_{\gamma }} = {{\rho }_{k}}{{\delta }_{{k,k{\kern 1pt} '}}},\quad {{\rho }_{k}} = {{({{q}^{k}},{{q}^{k}})}_{\gamma }}. \\ \end{gathered} $

Отсюда и из (16) получаем формулы для итерационных параметров:

(18)

${{\alpha }_{k}} = {{{{\sigma }_{k}}} \mathord{\left/ {\vphantom {{{{\sigma }_{k}}} {{{\rho }_{k}}}}} \right. \kern-0em} {{{\rho }_{k}}}},\quad {{\beta }_{k}} = {{{{\sigma }_{{k + 1}}}} \mathord{\left/ {\vphantom {{{{\sigma }_{{k + 1}}}} {{{\sigma }_{k}}}}} \right. \kern-0em} {{{\sigma }_{k}}}}.$

Определим теперь выражение для приближенного решения (10). Так как элементы r⁰,…, rⁿ ^{– 1} образуют в ${{\mathcal{K}}_{n}}$ ортогональный базис, мы можем записать:

${{f}^{n}} \equiv {{E}_{n}}f = {{F}_{{n - 1}}}(A){{v}^{0}} = \sum\limits_{k = 0}^{n - 1} {{{{{r}^{k}}{{{(f,{{r}^{k}})}}_{\gamma }}} \mathord{\left/ {\vphantom {{{{r}^{k}}{{{(f,{{r}^{k}})}}_{\gamma }}} {{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}}}} \right. \kern-0em} {{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}}}.} $

Поскольку отсюда следует

(19)

$\begin{gathered} {{F}_{n}}(0) = \sum\limits_{k = 0}^n {{{{{{(f,{{r}^{k}})}}_{\gamma }}{{R}_{k}}(0)} \mathord{\left/ {\vphantom {{{{{(f,{{r}^{k}})}}_{\gamma }}{{R}_{k}}(0)} {{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}}}} \right. \kern-0em} {{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}}} = {{F}_{{n - 1}}}(0) + f_{n}^{n},} \\ {\text{ }}f_{n}^{n} = {{{{{(f,{{r}^{n}})}}_{\gamma }}} \mathord{\left/ {\vphantom {{{{{(f,{{r}^{n}})}}_{\gamma }}} {{\text{||}}{{r}^{n}}{\text{||}}_{\gamma }^{2}}}} \right. \kern-0em} {{\text{||}}{{r}^{n}}{\text{||}}_{\gamma }^{2}}}, \\ \end{gathered} $

то из выражений (10), (16) получаем

$\begin{gathered} {{u}^{{k + 1}}} = {{A}^{{ - 1}}}[{{f}^{{k + 1}}} - {{F}_{k}}(0){{r}^{{k + 1}}}] = \\ \, = {{A}^{{ - 1}}}[{{f}^{k}} - {{F}_{{k - 1}}}(0){{r}^{k}} + {{\alpha }_{k}}({{F}_{{k - 1}}}(0) + f_{k}^{r})A{{q}^{k}}] = \\ \, = {{u}^{k}} + {{\alpha }_{k}}{{\gamma }_{k}}{{q}^{k}}, \\ \end{gathered} $

где введено обозначение γ_n = F_n(0). Таким образом, решение проблемы моментов определяется двучленными соотношениями

${{u}^{{k + 1}}} = {{u}^{k}} + {{\alpha }_{k}}{{\gamma }_{k}}{{q}^{k}},\quad {{\alpha }_{k}} = {{{{{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}} \mathord{\left/ {\vphantom {{{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}} {(A{{q}^{k}},{{q}^{k}})}}} \right. \kern-0em} {(A{{q}^{k}},{{q}^{k}})}}}_{\gamma }},$

(20)

$\begin{gathered} {{r}^{{k + 1}}} = {{r}^{k}} - {{\alpha }_{k}}A{{q}^{k}},\quad {{\gamma }_{k}} = {{\gamma }_{{k - 1}}} + f_{k}^{k}, \\ {{q}^{{k + 1}}} = {{r}^{{k + 1}}} + {{\beta }_{k}}{{q}^{k}},\quad {{\beta }_{k}} = {{{\text{||}}{{r}^{{k + 1}}}{\text{||}}_{\gamma }^{2}} \mathord{\left/ {\vphantom {{{\text{||}}{{r}^{{k + 1}}}{\text{||}}_{\gamma }^{2}} {{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}}}} \right. \kern-0em} {{\text{||}}{{r}^{k}}{\text{||}}_{\gamma }^{2}}},{\text{ }} \\ \end{gathered} $

${{r}^{0}} = {{q}^{0}} = {{v}^{0}},\quad {{\gamma }_{0}} = {{(f,{{v}^{0}})_{\gamma }^{2}} \mathord{\left/ {\vphantom {{(f,{{v}^{0}})_{\gamma }^{2}} {{\text{||}}{{v}^{0}}{\text{||}}_{\gamma }^{2}}}} \right. \kern-0em} {{\text{||}}{{v}^{0}}{\text{||}}_{\gamma }^{2}}}.$

Эти формулы упрощаются, если положить ${{v}^{0}}$ = f, поскольку при этом γ₀ = … = γ_n = 1. Полагая ${{v}^{0}} = f - {\text{ }}А{{u}^{0}}$ для произвольного начального вектора u⁰, на любой n-й итерации имеем вектор невязки rⁿ = f – Аuⁿ, т.е. приходим к методам сопряженных направлений [4, 5].

3. НЕКОТОРЫЕ СВОЙСТВА ИТЕРАЦИОННОГО МЕТОДА МОМЕНТОВ

Алгоритм (20) позволяет экономично решать такие актуальные задачи, как серии СЛАУ с одинаковой матрицей A, но разными правыми частями f^(m), m = 1, …, M, которые определяются последовательно, т.е. новый m-й вектор f^(m) может быть вычислен только после решения предыдущей (m – 1)-й системы. При этом первая СЛАУ решается по формулам (20), а последующие искомые векторы u^(m) – с помощью уже найденных коэффициентов α_n и векторов qⁿ. Перевычислять надо только величины γ_n, что значительно сокращает объем арифметических операций, особенно ввиду отсутствия матрично-векторных умножений. Для приближенных решений uⁿ вместо (20) можно использовать их представление через векторы qⁿ, которые в силу свойств ортогональности (17) образуют базис в ${{\mathcal{K}}_{n}}$ и удовлетворяют соотношениям

(21)

$\begin{gathered} {{A}_{n}}{{q}^{k}} = \alpha _{k}^{{ - 1}}[ - {{\beta }_{{k - 1}}}{{q}^{{k - 1}}} + (1 + {{\beta }_{k}}){{q}^{k}} - {{q}^{{k + 1}}}], \\ {{A}_{n}}{{q}^{0}} = \alpha _{0}^{{ - 1}}[(1 + {{\beta }_{0}}){{q}^{0}} - {{q}^{1}}], \\ {{A}_{n}}{{q}^{{n - 1}}} = \alpha _{{n - 1}}^{{ - 1}}[ - {{\beta }_{{n - 1}}}{{q}^{{n - 2}}} + (1 + {{\beta }_{{n - 1}}}){{q}^{{n - 1}}}]. \\ \end{gathered} $

Подставляя искомое представление решения

(22)

${{u}^{n}} = \sum\limits_{k = 0}^{n - 1} {{{\nu }_{k}}{{q}^{k}}} = {{u}^{{n - 1}}} + {{\nu }_{{n - 1}}}{{q}^{{n - 1}}},\quad {{u}^{{ - 1}}} = 0,$

в уравнение (6), с помощью соотношений (17) и (20) получаем следующее выражение:

$\begin{gathered} \sum\limits_{k = 0}^{n - 1} {{{\nu }_{k}}\alpha _{k}^{{ - 1}}[ - {{\beta }_{{k - 1}}}{{q}^{{k - 1}}} + (1 + {{\beta }_{k}}){{q}^{k}} - {{q}^{{k + 1}}}]} = \sum\limits_{k = 0}^{n - 1} {f_{k}^{{}}{{q}^{k}}} , \\ {{f}_{k}} = {{{{{(f,{{q}_{k}})}}_{\gamma }}} \mathord{\left/ {\vphantom {{{{{(f,{{q}_{k}})}}_{\gamma }}} {\left\| {{{q}_{k}}} \right\|}}} \right. \kern-0em} {\left\| {{{q}_{k}}} \right\|}}_{\gamma }^{2}. \\ \end{gathered} $

Приравнивая коэффициенты при одинаковых векторах, получаем трехдиагональную систему уравнений для нахождения ν_k из (22):

(23)

$\begin{gathered} - \,\,\alpha _{{k - 1}}^{{ - 1}}{{\nu }_{{k - 1}}} + \alpha _{k}^{{ - 1}}(1 + {{\beta }_{k}}){{\nu }_{k}} - \alpha _{{k + 1}}^{{ - 1}}{{\beta }_{k}}{{\nu }_{{k + 1}}} = f_{k}^{{}}, \\ {{\nu }_{{ - 1}}} = {{\nu }_{n}} = 0,\quad k = 0,1,...,n - 1. \\ \end{gathered} $

Подчеркнем, что в представлении (22) для решения uⁿ коэффициенты ν_k зависят от n, и для их нахождения требуется решать последовательность систем возрастающих порядков вида (23). Это можно сделать экономично с помощью метода прогонки [5], “стандартные” формулы которого в данном случае имеют следующий вид:

$\begin{gathered} {{\operatorname{ж} }_{1}} = {{{{f_{1}^{q}} \mathord{\left/ {\vphantom {{f_{1}^{q}} t}} \right. \kern-0em} t}}_{1}}{\text{,}}\quad {{\operatorname{ж} }_{k}} = (f_{k}^{q} - {{x}_{k}}{{\operatorname{ж} }_{{k - 1}}}{\text{)}}{{s}_{k}}{\text{,}} \\ {{s}_{k}} = {{({{t}_{k}} - {{x}_{k}}{{\operatorname{ж} }_{{k - 1}}}{\text{)}}}^{{ - 1}}}, \\ \end{gathered} $

(24)

${{\theta }_{1}} = {{{{{{y}_{1}}} \mathord{\left/ {\vphantom {{{{y}_{1}}} t}} \right. \kern-0em} t}}_{1}},\quad {{\theta }_{k}} = {{y}_{k}}{{s}_{k}},\quad k = 1,2,...,n - 1,$

$\begin{gathered} {{\nu }_{{n - 1}}} = {{\operatorname{ж} }_{{n - 1}}},\quad {{\nu }_{k}} = {{\theta }_{k}}{{\nu }_{{k + 1}}} + {{\operatorname{ж} }_{k}}, \\ k = n - 2,...,1. \\ \end{gathered} $

Здесь ради краткости введены обозначения для коэффициентов СЛАУ (23):

$\begin{gathered} {{t}_{k}} = \alpha _{k}^{{ - 1}}(1 + {{\beta }_{k}}),\quad {{x}_{k}} = \alpha _{{k - 1}}^{{ - 1}},\quad {{y}_{k}} = \alpha _{{k + 1}}^{{ - 1}}{{\beta }_{k}}, \\ k = 0,1,...,n - 1;\quad {{x}_{1}} = {{y}_{{n - 1}}} = 0. \\ \end{gathered} $

Величины æ_k и θ_k при разных n вычисляются одинаковыми рекурсиями, отличающимися только длиной. Реализация (22) требует одну компоненту решения из (24): ${{\nu }_{{n--}}}_{1} = {{{\unicode{230} }}_{{n--}}}_{1}.$

Вопросы устойчивости и точности метода прогонки (24) для решения трехдиагональной СЛАУ (23) играют ключевую роль в эффективности всего алгоритма. Отметим, что рассматриваемый итерационный процесс допускает обобщения и на более широкие типы симметричных алгебраических систем: знаконеопределенных, сингулярных и несовместных, – для которых в классических крыловских подпространствах развиты методы MINRES, MINRES-QLP и другие, использующие для решения трехдиагональных СЛАУ типа (23) алгоритмы QR- и QLP-разложений (см. [8] и цитируемую там литературу).

Рассматриваемый итерационный процесс является обобщением известных алгоритмов сопряженных направлений, поскольку исходный вектор ν⁰ подпространства Крылова ${{\mathcal{K}}_{n}}({{v}^{0}},A)$ никак не связан с решаемой СЛАУ. Однако методы моментов наследуют вариационные свойства алгоритмов минимальных итераций, или ошибок [9], а также сопряженных градиентов и сопряженных невязок [4, 5] в силу следующего из (20) представления для векторов r^k:

(25)

${{r}^{k}} = {{r}^{0}} - {{\alpha }_{0}}A{{q}^{0}} - {{\alpha }_{1}}A{{q}^{1}} - ... - {{\alpha }_{{k - 1}}}A{{q}^{{n - 1}}},$

где векторы обладают свойствами ортогональности (17). Умножая обе части (25) скалярно на вектор A^{γ – 2}r^k, при значениях α_k из (20) получаем выражения

(26)

$\begin{gathered} {{\varphi }_{\gamma }}({{r}^{k}},{{r}^{k}}) = {{({{r}^{k}},{{r}^{k}})}_{{\gamma - 2}}} = \\ \, = {{({{r}^{0}},{{r}^{0}})}_{{\gamma - 2}}} - \sum\limits_{l = 0}^{k - 1} {{{{{{({{r}^{0}},{{q}^{l}})_{{\gamma - 1}}^{2}} \mathord{\left/ {\vphantom {{({{r}^{0}},{{q}^{l}})_{{\gamma - 1}}^{2}} {({{q}^{l}},{{q}^{l}})}}} \right. \kern-0em} {({{q}^{l}},{{q}^{l}})}}}}_{\gamma }},} \\ \end{gathered} $

причем данный функционал достигает максимума в подпространстве ${{\mathcal{K}}_{k}}({{v}^{0}},A)$.

В силу оптимальных свойств функционала (26) можно также с помощью чебышевской методологии исследовать скорость сходимости и получить оценку числа итераций n(ε), необходимых для выполнения условия

(27)

${{\varphi }_{\gamma }}({{r}^{n}},{{r}^{n}}) \leqslant {{\varepsilon }^{2}}{{\varphi }_{\gamma }}({{r}^{0}},{{r}^{0}}),\quad \varepsilon \ll 1.$

Очевидно, что соответствующее неравенство имеет вид

(28)

$n(\varepsilon ) \leqslant 1 + 0.5\unicode{230} \ln \left| {{\varepsilon \mathord{\left/ {\vphantom {\varepsilon 2}} \right. \kern-0em} 2}} \right|,$

где æ – число обусловленности матрицы A (см. [4, 5]).

Остановимся на вопросах предобусловливания СЛАУ при использовании методов моментов для уменьшения числа æ в оценке (28). Для сохранения симметричности естественно использовать двустороннее предобусловливание, приводящее уравнение (1) к виду

(29)

$\bar {A}\bar {u} = \bar {f},\quad \bar {A} = {{C}^{ \top }}AC,\quad \bar {u} = {{C}^{{ - 1}}}u,\quad \bar {f} = {{C}^{ \top }}f,$

где C – некоторая легко обращаемая невырожденная матрица. В этом случае все предыдущие рассуждения повторяются для предобусловленной СЛАУ с точностью до замены обозначений $\bar {u},\bar {f},\bar {A}$. Примером могут служить методы неполной факторизации [5] для системы с матрицей A = = D + L + U, где D – блочно-диагональная, а L и U – нижняя и верхняя треугольные части матрицы. При этом предобусловленная СЛАУ определяется как

$\begin{gathered} \bar {A}\bar {u} = {{(I + \bar {L})}^{{ - 1}}}(\bar {D} + \bar {L} + \bar {U}){{(I + \bar {U})}^{{ - 1}}}\bar {u} = \\ \, = \bar {f} = {{(I + \bar {L})}^{{ - 1}}}G_{L}^{{ - 1}}f, \\ \end{gathered} $

(30)

$\bar {L} = G_{L}^{{ - 1}}LG_{U}^{{ - 1}},\quad \bar {D} = G_{L}^{{ - 1}}DG_{U}^{{ - 1}},$

$\begin{gathered} \bar {U} = G_{L}^{{ - 1}}UG_{U}^{{ - 1}},\quad \bar {u} = (I - \bar {U}){{G}_{U}}u, \\ {{G}_{L}}{{G}_{U}} = G = D - \overline {L{{G}^{{ - 1}}}U} , \\ \end{gathered} $

где $\overline {L{{G}^{{ - 1}}}U} $ – некоторая аппроксимация матрицы $L{{G}^{{ - 1}}}U$. Рассматриваемые в (30) предобусловливатели являются специальным случаем широко распространенных ILU-разложений (см. [4, 5]).

При заданном с.п.о. предобусловливателе B построение предобусловленного процесса Ланцоша представляется как метод ортогонализации, примененный к предобусловленной СЛАУ

$\bar {A}\bar {u} = \bar {f} = {{B}^{{ - 1/2}}}f,\quad \bar {A} = {{B}^{{ - 1/2}}}A{{B}^{{ - 1/2}}},\quad \bar {u} = {{B}^{{1/2}}}u.$

При этом в формулах (11) для новых векторов ${{w}^{k}} = {{B}^{{ - 1/2}}}{{p}^{k}}$ получим рекуррентное соотношение

${{w}^{{k + 1}}} = A{{B}^{{ - 1}}}{{w}^{k}} - \bar {\alpha }_{k}^{{(\gamma )}}{{w}^{k}} - \bar {\beta }_{k}^{{(\gamma )}}{{w}^{{k - 1}}},$

в котором фактически B^1/2 не требуется (см. подробнее [8]).

Замечание 2. Отметим также полезность таких простых предобусловливаний, как масштабирование или бинормализация СЛАУ [10], которые соответствуют в (27) диагональным матрицам C и могут употребляться в сочетании с другими предобусловливателями.

Замечание 3. Ю.В. Воробьёвым в [1] рассмотрен обобщенный метод моментов, заключающийся в выборе вместо одного исходного вектора ${{v}^{0}}$ заданного набора линейно независимых m векторов $v_{1}^{0},...,v_{m}^{0}$. В методах решения СЛАУ такие подходы, перспективные с точки зрения распараллеливания, называются блочными крыловскими алгоритмами, см. [11]. Описанные методы без труда обобщаются на их блочные варианты. Получаемые при этом алгоритмы формально допускают мультипредобусловливание, т.е. использование на одной итерации по несколько предобусловливающих матриц и направляющих векторов $p_{1}^{k},...,p_{m}^{k}$, см. [12]. Однако обоснование этого подхода для симметричных СЛАУ требует дополнительных ислледований.

Список литературы

Воробьёв Ю.В. Метод моментов в прикладной математике. М.: Физматлит, 1958.
Ахиезер Н.И. Классическая проблема моментов. М.: Физматлит, 1961.
Liesen J., Strakos Z. Krylov subspace methods. Principles and analysis. Oxford: Oxford University Press, 2013.
Saad Y. Iterative methods for sparse linear systems. 2nd. ed. SIAM, 2003.
Ильин В.П. Методы и технологии конечных элементов. Новосибирск: ИВМиМГ СО РАН, 2007.
Olshanskii M.A., Tyrtyshnikov E.E. Iterative methods for linear systems theory and applications. SIAM, Philadelphia, 2014.
Craig E.G. The N-step iteration procedure // J. Math. Phis. 1955. V. 8. P. 64–73.
Choi S.-C.T., Paige C.C., Saunders M.A. MINRES-QLP: Krylov subspace method for indefinite or singular symmetric systems. SIAM, arXiv:1003.4042[math.NA]. 27 Mar 2015.
Фридман В.М. Метод минимальных итераций с минимальными ошибками для систем линейных алгебраических уравнений с симметричной матрицей // ЖВМиМФ. 1962. Т. 2. № 2. С. 341–342.
Livne O.E., Golub J.H. Scaling by binormalazation // Numer. Algorithms. 2004. V. 35(1). P. 97–120.
Nikishin A.A., Eremin A.Y. Variable block CG algorithms for solving large sparse symmetric positive definite linear systems on parallel computers // SIAM J. Matrix Anal. Appl. 1995. V. 16. P. 1135–1153.
Il’in V. P. Multi-preconditioned domain decomposition methods in the Krylov subspaces // LNCS, 10187, Springer, 2017. P. 95–107.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Доклады Российской академии наук. Математика, информатика, процессы управления

Архивы выпусков Информация о журнале Отправить рукопись в журнал