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

Цель исследования: выявить зависимость между сходимостью рассматриваемых численных методов и размерностью СЛАУ.

Материалы и методы исследования

Вычисления проводились на программной платформе .NET Framework 4.5.1, в качестве IDE использовалась Microsoft Visual Studio 2017. Для реализации метода GMRES за основу был взят github проект участника Dan Lavrushko [1]. LU декомпозиция выполнена с помощью библиотеки Math.NET Numerics [2].

Постановка задачи

Рассмотрим задачу решения СЛАУ A•x = b итерационным методом GMRES [3] и LU разложением (декомпозицией) [4], где: A – невырожденная, плотная, несимметричная матрица размером m×m; x – вектор неизвестных размером m; b – свободные члены размером m (1).

(1)

Метод GMRES

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

Метод подразумевает использование алгоритма итераций Арнольди [5] с целью нахождения приближенных собственных векторов матрицы посредством построения ортонормированного базиса подпространства Крылова [6] для осуществления аппроксимации решения с минимальной невязкой.

Сам алгоритм метода GMRES был предложен в 1986 г. Юсефом Саадом и Мартином Х. Шульцем. GMRES стал обобщением метода MINRES Криса Пейджда и Майкла Сондерса, также является частным случаем метода DIIS Питера Пулайе, применимым к нелинейным системам.

Из сложности процесса Арнольди, минимизации невязки и матрично-векторного произведения складывается сложность метода GMRES.

Увеличение объема требуемой памяти для покрытия стоимости вычисления собственных векторов на каждой итерации составляет O(n2), где n – номер итерации. Очевидно, что многие вычислительные системы не в состоянии удовлетворить такие потребности алгоритма, поэтому для устранения проблемы производится перезапуск вычислений после заданного значения k итераций, начиная с xk в качестве начального предположения. Такой подход признан модификацией исходного GMRES алгоритма и обозначается с припиской Restarted. При использовании Restarted GMRES возникает существенный недостаток: перезапущенное подпространство часто близко к более раннему подпространству, поэтому метод страдает от стагнации конвергенции. Решением проблемы служит повторное использование подпространства Крылова методами GCROT и GCRODR, которые относятся к типу GCRO.

Сложность решения задачи минимизации невязки (минимальных квадратов) составляет O(m2), где m – размер матрицы.

Сложность получения решения матрично-векторным произведением составляет O(nm).

Математическое описание алгоритма GMRES. На каждом k-м шаге решение xk рассматривается в виде , где столбцы υ_k матрицы V находятся с помощью процесса Арнольди для построении ортонормированного базиса подпространства Крылова , вектор y_k представляет собой решение системы , где Hk – верхняя матрица Хессенберга размера k×k, элементы которой формируются в процессе ортогонализации Арнольди, , , и, наконец, e1 – канонический вектор размерности k.

Ортонормированный базис {υ1, υ2,... υk} подпространства вычисляется с помощью ортогонализации Арнольди по формулам

(2)

После завершения вычислений процессом Арнольди получаем ортонормированный базис подпространства Крылова и матрицу, расширенную последним вычисленным вектором : .

Hk – верхняя Хессенбергова матрица, элементы которой равны коэффициентам ортогонализации (3).

(3)

Выполнив несложные преобразования, можно получить выражение для невязки на k-й итерации: . Очевидно, что для задачи минимум достигается при .

Для решения последней системы матрица приводится к верхнему треугольному виду с помощью k вращений Гивенса (4).

(4)

где , .

Метод LU разложения

LU разложение основывается на методе последовательного исключения Гаусса. Метод LU разложения подразумевает представление исходной матрицы A в виде произведения нижнетреугольной матрицы L на верхнетреугольную матрицу U (5). Традиционно, либо матрица L, либо U нормализуются до получения единичной диагонали. Несмотря на то, что метод LU разложения предназначен для матриц любой (прямоугольной) формы, чаще всего он применяется к квадратным матрицам A•x = b.

(5)

Хорошим свойством LU декомпозиции является то, что для хранения опорных точек достаточно O(m) требуемой памяти, так как элементы матрицы L выше главной диагонали равны нулю (сама матрица L имеет единичную диагональ, в соответствии с классической схемой нормализации), поэтому должны быть сохранены только субдиагональные элементы. Точно так же с матрицей U – элементы ниже главной диагонали равны нулю, поэтому должны быть сохранены только главная диагональ и супердиагональные элементы (6).

(6)

Таким образом, LU представление матрицы A занимает столько же места, сколько исходная матрица. Исходную матрицу не требуется хранить, так как вся информация о ней расположена в полученных треугольных L и U матрицах.

Различают две реализации алгоритма: наивную и блочно-матричную факторизацию.

Наивная реализация LU разложения плохо масштабируется в зависимости от размера СЛАУ. Когда размерность матрицы A превышает 64×64, она не вписывается в типичную кэш-память ЦП L1, вследствие чего алгоритм затрачивает большую часть времени исполнения на ожидание поступления данных из кэш-памяти L2/L3 или ОЗУ.

Блочно-матричный подход является лучшим решением для решения СЛАУ больших размерностей. При его применении большая матрица разбивается на мелкие блоки и выполняется пошаговая факторизация. Алгоритм LU разложения при таком подходе будет состоять из нескольких шагов.

На первом шаге от входящей матрицы m×m производится отделение ведущей подматрицы m×k (7).

(7)

где , , .

На следующем шаге производится LU факторизация наименьшей матрицы m×k, в то время как остальные части остаются неизменными (8).

(8)

Сравнительный анализ итерационного метода GMRES с прямым LU разложением

Следующий шаг – обновление части матрицы результатами исходной LU. Верхний блок A21 умножается на inv(L11), а конечная матрица A22 получает более сложное обновление (9).

(9)

В результате имеется k столбцов L-фактора в L1, первые k строк U-фактора в матрицах U1 и U21, а оставшаяся часть, подлежащая обработке, A, сохраняется в конечных mk-строках/столбцах. На этом этапе остается только применить алгоритм к конечным строкам/столбцам A, пока не получится факторизация полной матрицы.

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

Частными случаями LU разложения являются: метод прогонки (используется для трехдиагональных матриц), факторизация Дулитла (если у матрицы L на главной диагонали все единицы), факторизация Крита (если у матрицы U на главной диагонали все единицы). Когда U = LT (или L = UT), это называется разложением Холецкого.

Результаты исследования и их обсуждение

Программно реализовав описанные выше методы, был проведен численный эксперимент, результаты которого представлены в таблице (в соответствии с таблицей).

Для сравнения поведения методов на СЛАУ больших размерностей на каждом шаге сгенерирована произвольная плотная несимметричная матрица квадратной формы A и массив свободных членов b (1).

Решение итерационным методом GMRES находилось с заранее заданной точностью 1E-3.

Погрешность решения рассмотренными численными методами [7] для каждой размерности СЛАУ находилась по формуле

(10)

Для проведения сравнительного анализа была найдена максимальная невязка между решениями, полученными представленными методами:

(11)

Полученные результаты могут быть использованы при решении задач, поставленных в работах [8–10].

Выводы

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