Журнал вычислительной математики и математической физики, 2019, T. 59, № 8, стр. 1331-1339

ВВЕДЕНИЕ

В настоящее время итерационные методы решения систем линейных алгебраических уравнений (СЛАУ) приобретают все большую популярность. Это связано с тем, что во многих случаях необходимо решать СЛАУ огромной размерности N и с плотными матрицами. Подобные задачи возникают, например, при численном решении многомерных интегральных уравнений. Ясно, что в этих случаях практически невозможно решать СЛАУ прямыми методами, типа метода Гаусса, поскольку тогда для получения решения необходимо выполнить $T\sim {{N}^{3}}$ арифметических операций и хранить в памяти компьютера порядка $M\sim {{N}^{2}}$ чисел. Для итерационных методов указанные характеристики алгоритмов T и M оцениваются формулами $T\sim L{{T}_{A}}$ и $M\sim {{M}_{A}}$, где ${{T}_{A}}$ – число арифметических операций для умножения матрицы на вектор, L – количество итераций для получения решения с заданной точностью, а ${{M}_{A}}$ – число различных элементов матрицы.

Используя равномерную сетку для дискретизации многомерных интегральных уравнений и алгоритмы быстрого дискретного преобразования Фурье, для решения многих задач математической физики получаются СЛАУ для которых ${{T}_{A}}\sim N$ и ${{M}_{A}}\sim N$ [1]. Тогда главной проблемой численного решения становится уменьшение числа итераций L, требуемое для получения решения с заданной точностью.

Существует большое число итерационных методов. Их можно разделить на два класса: стационарные и нестационарные методы. К первому классу относятся методы, в которых итерационные параметры определяются до выполнения итераций и в дальнейшем не изменяются. Ко второму классу относятся методы, в которых итерационные параметры определяются в процессе выполнения итераций и зависят от номера итерации. К стационарным относятся методы простой итерации, чебышёвский итерационный метод и другие [2], [3]. К нестационарным относятся методы сопряженных градиентов, минимальных невязок и другие [3]–[5]. Условия применимости стационарных методов для решения уравнений, как правило, более жесткие по сравнению с нестационарными методами. Тем не менее практически у каждого итерационного метода есть класс задач, где его применение наиболее эффективно.

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

1. СИСТЕМЫ УРАВНЕНИЙ С ДЕЙСТВИТЕЛЬНЫМИ МАТРИЦАМИ

Рассмотрим сначала систему N линейных алгебраических уравнений с действительными коэффициентами матрицы и правой части вида

Систему линейных уравнений (1.1) можно записать в матричном виде: где A – матрица коэффициентов системы уравнений, а x и f суть N-мерные векторы. Ниже будем полагать, что определитель матрицы не равен нулю и поэтому система уравнений имеет единственное решение при любой правой части. Никаких других ограничений на матрицу A налагать не будем.

Определим функции

Очевидно, что $\Phi (x)$ является квадратичной функцией переменных и имеет минимум, который определяется из условий Из (1.3), (1.4) имеем Из условий экстремума (1.4) и (1.5) получим Поскольку определитель матрицы системы уравнений не равен нулю, то равенства (1.6) имеют место, если все ${{r}_{n}}$ равны нулю. Принимая во внимание обозначения в (1.5), это значит, что минимум функции $\Phi (x)$ будет при значении вектора x, который является решением исходной системы уравнений (1.2).

Далее, для решения системы уравнений (1.1) будем использовать алгоритм нахождения минимума функции $\Phi (x)$. Будем применять метод градиентного спуска. Из (1.5) следует, что градиент функции $\Phi (x)$ в любой точке x определяется формулой

В матричном виде (1.7) можно записать как где $g$ – вектор в направлении градиента функции $\Phi (x)$ в точке x, ${{A}^{{\text{т }}}}$ – матрица, транспонированная к исходной матрице $A$, $r$ – вектор невязки в точке x. Теперь в направлении вектора $g$ будем осуществлять спуск до минимального значения функции, т.е. искать минимум функции $\Phi (x + hg)$ по отношению к действительной величине h. Функцию $\Phi (x)$ можно представить в виде скалярного произведения Тогда Приравнивая производную функции $\Phi (x + hg)$ к нулю, получаем для значения h следующую формулу:

Таким образом, процедура итерационного метода градиентного спуска состоит в следующем. Выбирается какое-то начальное значение вектора ${{x}_{0}}$. Затем с использованием формулы (1.11) находится значение вектора следующего приближения ${{x}_{1}}$, затем ${{x}_{2}}$ и т.д. до тех пор, пока не получим требуемую точность. Используя (1.8), (1.11), окончательно запишем итерационную процедуру в векторном виде:

Поскольку $(A{{A}^{{\text{т }}}}{{r}_{k}},{{r}_{k}}) = {{\left\| {{{A}^{{\text{т }}}}{{r}_{k}}} \right\|}^{2}}$, то итерационную процедуру (1.12) можно записать в следующем эквивалентном виде:

Из (1.13) следует, что последовательные невязки итераций связаны соотношением

Далее имеем Таким образом, получаем Поскольку ${{\left\| {{{A}^{{\text{т }}}}{{r}_{k}}} \right\|}^{4}} = {{\left| {({{A}^{{\text{т }}}}{{r}_{k}},{{A}^{{\text{т }}}}{{r}_{k}})} \right|}^{2}} = {{\left| {({{r}_{k}},A{{A}^{{\text{т }}}}{{r}_{k}})} \right|}^{2}}$, то из (1.15) имеем

Отметим важное обстоятельство. Поскольку $\left| {\det A} \right| = \left| {\det {{A}^{{\text{т }}}}} \right| > 0$, то, так как ${{r}_{k}} \ne 0$ (а равенство нулю может быть только, если на k-м шаге итерационной процедуры мы получили точное решение), имеем из (1.15), (1.16), что $0 \leqslant \left\| {{{r}_{{k + 1}}}} \right\|{\text{/}}\left\| {{{r}_{k}}} \right\| < 1$. Значит, невязка решения уменьшается на каждой итерации и в пределе получим точное решение. Заметим, что если векторы ${{r}_{k}}$ и $A{{A}^{{\text{т }}}}{{r}_{k}}$ коллинеарные, то на (k + 1)-м шаге итерационной процедуры мы получим точное решение. Можно дать геометрическую интерпретацию этому условию. Если векторы u и $\text{v}$ коллинеарные, то $u = \lambda \text{v}$, где $\lambda $ – константа и поэтому векторы ${{A}^{{ - 1}}}u$ и ${{A}^{{ - 1}}}\text{v}$ тоже коллинеарные. Значит, если векторы ${{r}_{k}}$ и $A{{A}^{{\text{т }}}}{{r}_{k}}$ коллинеарные, то и векторы ${{A}^{{ - 1}}}{{r}_{k}}$ и ${{A}^{{ - 1}}}(A{{A}^{{\text{т }}}}{{r}_{k}}) = {{A}^{{\text{т }}}}{{r}_{k}}$ также коллинеарные. Тогда и векторы $({{x}_{k}} - x)$ и ${{g}_{k}}$ коллинеарные, где x – решение уравнения (1.2), а ${{g}_{k}}$ – вектор в направлении градиента функции $\Phi (x)$в точке ${{x}_{k}}$. Значит, в этом случае луч, проведенный из точки ${{x}_{k}}$ в направлении (–${{g}_{k}}$), проходит через точку x и поэтому на (k + 1)-м шаге итерационной процедуры получим точное решение. Если же луч проходит вблизи точки x, то на этом шаге мы сильно приблизимся к решению.

2. СИСТЕМЫ УРАВНЕНИЙ С КОМПЛЕКСНЫМИ МАТРИЦАМИ

Итерационная процедура (1.13) была получена для случая, когда матрица и правая часть действительны. Однако во многих задачах матрица системы уравнений и правая часть комплексные. Рассмотрим этот общий случай. Пусть задана система уравнений вида

где $H$ и $f$ – комплексные матрица и вектор соответственно. Будем также предполагать, что определитель матрицы не равен нулю. Никаких других ограничений на матрицу H не налагается.

Представим комплексную матрицу $H$ и векторы z и f в виде

где $A$, $B$ – действительные матрицы, а $x$, $y$ и ${{f}^{{(1)}}}$, ${{f}^{{(2)}}}$ – действительные векторы. Тогда система (2.1) может быть записана в виде

Разделяя действительную и мнимую части в (2.3), получаем, что система уравнений (2.1) эквивалентна следующей системе уравнений с действительными матрицами и векторами вида:

Размерность системы уравнений (2.4) в два раза больше размерности комплексной системы (2.1). Для системы (2.4) будем использовать итерационную процедуру (1.13).

Невязка k-й итерации для системы (2.4) описывается формулой

Используя (2.2), получаем, что невязка ${{r}_{k}} = r_{k}^{{(1)}} + ir_{k}^{{(2)}}$ может быть представлена в виде ${{r}_{k}} = H{{z}_{k}} - f$, где ${{z}_{k}} = {{x}_{k}} + i{{y}_{k}}$. Далее, вектор в направлении градиента функции $\Phi (x,y)$ в точке $({{x}_{k}},{{y}_{k}})$определяется, согласно (1.8), формулой Несложно видеть, что Поскольку $H* = {{A}^{{\text{т }}}} - i{{B}^{{\text{т }}}}$, где $H{\text{*}}$ – матрица, сопряженная к H, из (2.6), (2.7) следует, что вектор ${{g}_{k}} = g_{k}^{{(1)}} + i{\kern 1pt} g_{k}^{{(2)}}$ может быть представлен в виде ${{g}_{k}} = H{\text{*}}{{r}_{k}}$. Норма любого вектора где u, $\text{v}$ суть N-мерные действительные векторы, представима в виде ${{\left\| \omega \right\|}^{2}} = {{\left\| u \right\|}^{2}} + {{\left\| \text{v} \right\|}^{2}}$. Если $\omega $ представить в виде комплексного вектора $\omega = u + i\text{v}$, то получим такое же значение для нормы.

Учитывая вышеизложенное, итерационную процедуру для системы (2.4) можно записать для комплексной системы (2.1) в следующем виде:

Невязки итераций (2.8) будут уменьшаться в соответствии с формулой а итерации будут сходиться к решению (2.1).

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

Рассмотрим обобщение описанного выше метода.

Сначала будем рассматривать систему N линейных алгебраических уравнений с действительными коэффициентами матрицы и правой части вида (1.2). В описанном выше методе новое значение итерации ${{x}_{{k + 1}}}$ ищется в направлении градиента ${{g}_{k}}$ в точке ${{x}_{k}}$. Однако при этом не используется информация о предыдущих итерациях.

Определим следующую итерационную процедуру:

где ${{t}_{k}}$ и ${{h}_{k}}$ – действительные итерационные параметры. В (3.1) ${{x}_{0}}$ – начальное приближение; ${{t}_{0}} = 0$, поскольку ${{x}_{{ - 1}}}$ не существует; ${{h}_{0}}$ определяется, согласно (1.12), формулой Мы полагаем, что ${{x}_{0}}$ не является точным решением уравнения (1.2) и поэтому ${{r}_{0}} \ne 0$. В противном случае процесс нахождения решения заканчивается.

Итерационные параметры при $k \geqslant 1$ будем искать из условия минимума квадрата нормы невязки ${{r}_{{k + 1}}}$. Из (3.1), учитывая, что ${{r}_{m}} = A{{x}_{m}} - f$, получаем

Определим функцию $F({{t}_{k}},{{h}_{k}}) = ({{r}_{{k + 1}}},{{r}_{{k + 1}}})$. Тогда итерационные параметры определяются из условий минимума функции $F({{t}_{k}},{{h}_{k}})$

Из (3.3), (3.4) получим, что итерационные параметры определяются из решения системы двух линейных уравнений с двумя неизвестными:

где $\Delta {{r}_{k}} = {{x}_{k}} - {{x}_{{k - 1}}}$. Из (3.5) следует, что определитель системы определяется формулой

Докажем, что определитель (3.6) не равен нулю при всех значениях $k \geqslant 1$. Из (3.6) ясно, что равенство нулю определителя будет иметь место, если векторы $\Delta {{r}_{k}}$ и $A{{g}_{k}}$ коллинеарные, т.е. $\Delta {{r}_{k}} = \lambda A{{g}_{k}}$, где $\lambda $ – любое число. Поскольку матрица A невырожденная, то тогда ${{x}_{k}} - {{x}_{{k - 1}}} = \lambda {{g}_{k}}$. Кроме этого, определитель будет обращаться в ноль, если $A{{g}_{k}} = 0$, а значит, поскольку матрицы A и ${{A}^{{\text{т }}}}$ имеют обратные, то и ${{r}_{k}} = 0$, это означает, что на предыдущем шаге итерации мы получили точное решение (1.2). Определитель (3.6) также обращается в ноль, если $\Delta {{r}_{k}} = 0$, т.е. ${{x}_{k}} = {{x}_{{k - 1}}}$. Однако это также может иметь место, если на предыдущем шаге итерации получено точное решение (1.2).

Итерационный параметр ${{h}_{0}}$ определяется из условия минимума функции

Тогда имеем

Из (3.1) при k = 0 получим, что $ - {{h}_{0}}{{g}_{0}} = {{x}_{1}} - {{x}_{0}}$. Тогда из (3.7) следует, что векторы $({{x}_{1}} - {{x}_{0}})$ и ${{g}_{1}}$ ортогональны и поэтому не могут быть коллинеарными. Значит, согласно вышеизложенному, определитель (3.6) не равен нулю для k = 1.Таким образом, система уравнений (3.5) имеет единственное решение при k = 1.

Докажем, что система (3.5) имеет единственное решение при $k > 1$, если на предыдущем шаге мы не получили точное решение уравнения (1.2), т.е. ${{r}_{k}} \ne 0$. Доказательство будем проводить методом математической индукции. Пусть при значении k система (3.5) имеет единственное решение. Докажем, что тогда и при значении $(k + 1)$ система (3.5) будет также иметь единственное решение.

Из условий минимума (3.4) при значении k следует

Учитывая (3.3), из первого равенства (3.8) получаем Из второго равенства (3.8) имеем Из (3.1), (3.9), (3.10) получим

Таким образом, из (3.11)) следует, что векторы ${{x}_{{k + 1}}} - {{x}_{k}}$ и ${{g}_{{k + 1}}}$ ортогональны и поэтому не могут быть коллинеарными. Значит, согласно вышеизложенному, при значении $(k + 1)$ система (3.5) будет иметь единственное решение. Поскольку при $k = 1$ рассматриваемая система уравнений имеет единственное решение, то итерационная процедура (3.1), (3.5) является корректной, так как определитель (3.6) не равен нулю при всех значениях $k \geqslant 1$.

Рассмотренная итерационная процедура является обобщением (1.12), для которой, поскольку $\left\| {{{r}_{{k + 1}}}} \right\|{\text{/}}\left\| {{{r}_{k}}} \right\| < 1$, итерации сходятся к решению (1.2). Значит, и итерации (3.1), (3.2), (3.5) также сходятся к решению (1.2).

Отметим важное обстоятельство. Итерации (3.1), (3.5) изменяются в направлении суперпозиции векторов разности последовательных итераций $({{x}_{k}} - {{x}_{{k - 1}}})$ и вектора градиента ${{g}_{k}}$, которые ортогональны. Это позволяет избегать “мертвых” зон, которые могут иметь место в итерационной процедуре, основанной только на использовании вектора градиента ${{g}_{k}}$. Также следует отметить, что, в силу ортогональности векторов $({{x}_{k}} - {{x}_{{k - 1}}})$ и ${{g}_{k}}$, сходимость итерационной процедуры (3.1), (3.5) будет всегда лучше итерационной процедуры (1.13).

Теперь рассмотрим комплексную систему уравнений (2.1).

Скалярное произведение двух действительных векторов

где ${{u}_{1}}$, ${{u}_{2}}$, ${{\text{v}}_{1}}$, ${{\text{v}}_{2}}$ суть N-мерные действительные векторы, представимо в виде $({{\omega }_{1}},{{\omega }_{2}}) = ({{u}_{1}},{{u}_{2}}) + ({{\text{v}}_{1}},{{\text{v}}_{2}})$. Если ${{\omega }_{1}}$, ${{\omega }_{2}}$ представить в виде комплексных векторов $\omega _{1}^{c} = {{u}_{1}} + i{{\text{v}}_{1}}$, $\omega _{2}^{c} = {{u}_{2}} + i{{\text{v}}_{2}}$, то имеем где справа в (3.12) стоит скалярное произведение комплексных векторов.

Учитывая изложенное в разд. 2 и (3.12), итерационную процедуру (3.1), (3.2), (3.5) для действительной системы (2.4) можно записать для комплексной системы (2.1) в следующем виде:

Обоснование сходимости итерационных процедур (2.8) и (3.13), (3.14) к решению уравнения (2.1) не требовало того, чтобы оператор H был матрицей. Поэтому очевидно, что итерации также сходятся к решению, если рассматривается операторное уравнение (2.1), например, интегральное уравнение, в гильбертовом пространстве. Единственное требование – оператор H должен иметь ограниченный обратный оператор.

4. ЧИСЛЕННЫЕ РЕЗУЛЬТАТЫ

Рассмотрим применение описанных методов для численного решения СЛАУ. Для итерационных методов сходимость итераций зависит от расположения собственных значений матрицы на комплексной плоскости. Чем больше отношение Q максимального модуля собственного значения к минимальному модулю, тем хуже сходимость.

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

В табл. 1 приведены зависимости количества итераций от размерности задачи N = 1000 и диапазона собственных значений Q для чистого градиентного (ЧГ) и модифицированного градиентного методов (МГ). Итерации заканчивались по достижении относительной невязки решения $\delta = \left\| r \right\|{\text{/}}\left\| f \right\|\sim {{10}^{{ - 5}}}$. Вычисления проводились для различных правых частей и для заданного Q разных значений собственных чисел на комплексной плоскости. Количество требуемых итераций мало отличалось друг от друга и в таблицах приводятся усредненные значения.

Как видно из табл. 1 модифицированный градиентный метод значительно более эффективен, по сравнению с чистым градиентным методом.

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

В табл. 2 приведены зависимость количества итераций от размерности СЛАУ N и диапазона собственных значений Q по достижении относительной невязки решения $\delta = \left\| r \right\|{\text{/}}\left\| f \right\|\sim {{10}^{{ - 5}}}$.

Как видно из табл. 2, сходимость итераций при заданном Q практически не зависит от размерности матрицы N.

В табл. 3 приведена зависимость количества итераций от размерности задачи и диапазона собственных значений в зависимости от точности решения.

Из приведенных численных результатов видно, что, даже для очень большого разброса диапазона собственных значений, сходимость итераций модифицированного метода градиентного спуска оказывается достаточно хорошей. Полученные данные могут являться ориентиром при численном решении конкретных задач.

ЗАКЛЮЧЕНИЕ

Рассмотрены итерационные методы градиентного спуска для решения линейных уравнений. Проведено математически строгое обоснование сходимости итераций к решению. Отметим, что, в отличие от большинства итерационных методов, рассмотренные методы не требуют никаких дополнительных условий, кроме того, что оператор уравнения имеет ограниченный обратный оператор. Проведенные численные эксперименты показали эффективность предложенного модифицированного итерационного метода градиентного спуска.