Scientific journal
Modern high technologies
ISSN 1812-7320
"Перечень" ВАК
ИФ РИНЦ = 0,899

COMPARATIVE ANALYSIS OF THE GMRES ITERATIVE METHOD WITH DIRECT LU DECOMPOSITION WHEN SOLVING THE SYSTEMS OF LINEAR EQUATIONS OF LARGE DIMENSIONS

Bogdanov A.E. 1 Torshina O.A. 1
1 Federal State Financed Educational Institution of Higher Education Magnitogorsk State Technical University named after Nosov
When choosing numerical methods for solving systems of linear algebraic equations, first of all, they are guided by such parameters as accuracy (does the method allow to obtain a solution with the required accuracy) and speed (depending on how quickly a solution is needed, the only relevant indicator is evaluated: time on a specific machine or the possibility of scalability of a method on a large computing cluster). Usually, it is necessary to choose a compromise between these parameters. The structure of the matrix (symmetric, dense or discharged, lamellar) and eigenvalues (whether they are all positive, their grouping, dispersion of values) are considered, equally with the above parameters, to select a specific method. Also, the important role is played by the implementation of the method (it is not always possible to find a good implementation for a suitable method). When it comes to systems of equations of large dimensions, the question on the convergence of methods for such systems appears. In this regard, to implement the optimal choice of a numerical method for solving systems of linear algebraic equations, the problem of evaluating various methods for a specific problem comes up. This article presents a comparative analysis of the iterative GMRES method with direct LU decomposition when solving systems of linear algebraic equations of large dimensions on dense asymmetric square matrices. For comparison of convergence, the error of the obtained solutions is found; the discrepancy between the solutions of the considered methods is calculated for all dimensions selected within the article.
numerical methods
systems of linear algebraic equations
GMRES
LU
iterative
direct

Рассмотрим численные методы 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).

bogdan01.wmf (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 рассматривается в виде bogdan02.wmf, где столбцы υk матрицы V находятся с помощью процесса Арнольди для построении ортонормированного базиса bogdan03.wmf подпространства Крылова bogdan04.wmf, вектор yk представляет собой решение системы bogdan05.wmf, где Hk – верхняя матрица Хессенберга размера k×k, элементы которой формируются в процессе ортогонализации Арнольди, bogdan06.wmf, bogdan07.wmf, и, наконец, e1 – канонический вектор bogdan08.wmf размерности k.

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

bogdan11.wmf (2)

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

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

bogdan15.wmf (3)

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

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

bogdan20.wmf (4)

где bogdan21.wmf, bogdan22.wmf.

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

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

bogdan23.wmf (5)

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

bogdan24.wmf (6)

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

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

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

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

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

bogdan25.wmf (7)

где bogdan26.wmf, bogdan27.wmf, bogdan28.wmf.

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

bogdan29.wmf (8)

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

Размерность

GMRES

LU

Невязка

Погрешность

Погрешность

100×100

2,0116088480075E-12

2,11060429475097E-12

2,0925483568135E-11

200×200

1,83254249558738E-13

5,04210124070489E-13

1,3944401189292E-13

300×300

1,95004482511071E-12

2,47861082279823E-12

3,47103679132488E-10

400×400

2,76235503824824E-13

1,1413797766996E-12

1,38888900380607E-13

500×500

4,48703350072223E-12

1,3246281079378E-11

1,64552815817842E-11

600×600

5,9363410406452E-13

2,72143384747093E-12

7,18314296932476E-14

700×700

1,42739638197509E-12

7,78304694683747E-12

2,00772731773213E-12

800×800

1,61541324545316E-12

8,03706747369854E-12

8,83848549904087E-12

900×900

5,80869227195123E-12

4,70100645431032E-11

3,67883501439792E-12

1000×1000

1,78459029722083E-12

1,50275970448393E-11

2,58693066967908E-12

2000×2000

8,9116957833455E-11

8,4298279402325E-10

4,32651248161164E-09

3000×3000

1,02300509221517E-10

2,25055203080853E-09

3,241424906264E-10

4000×4000

2,11655434842488E-10

5,10361178511771E-09

1,20999033015323E-09

5000×5000

2,45441089565206E-11

7,38593135178026E-10

9,10804764941986E-12

6000×6000

4,77423520551148E-11

1,13189903713514E-09

6,47268016962244E-10

7000×7000

9,69911363063168E-10

1,10518722842222E-08

1,12652671191427E-07

8000×8000

2,66089964273607E-11

1,36072345874094E-09

1,36598510280805E-11

9000×9000

4,99308489261675E-11

2,23105000993159E-09

1,92463156523104E-10

 

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

bogdan30.wmf (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] для каждой размерности СЛАУ находилась по формуле

bogdan31.wmf (10)

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

bogdan32.wmf (11)

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

Выводы

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