On a computer implementation of the block Gauss-Seidel method for normal systems of equations

Full Text

Введение. В статье рассматривается блочный метод Гаусса-Зейделя для нормальных систем уравнений [1]. Данный метод является одним из достаточно эффективных итерационных методов решения, в общем случае переопределенных, систем линейных алгебраических уравнений (СЛАУ), которые возникают в многочисленных приложениях [2]. Пусть задана СЛАУ Ax = b, (1) где A ∈ Rm×n , m n, rank(A) = n. Необходимо вычислять в общем случае псевдорешение x∗ = A† b. Здесь A† - псевдообратная матрица Мура-Пенроуза матрицы A [3, p. 230]. Как хорошо известно [1, 3], вычисление x∗ эквивалентно решению нормальных систем уравнений: A Ax = A b. (2) Если m = n, то x∗ = A-1 b, то есть x∗ является решением исходной СЛАУ (1). Важной характеристикой итерационных методов является скорость сходимости, которая зависит от числа обусловленности исходной задачи. Плохая обусловленность приводит к значительному снижению скорости сходимости итерационных методов. Многие практические задачи имеют достаточно высокое число обусловленности, а у методов, основанных на решении нормальных систем уравнений, оно равно квадрату числа обусловленности исходной задачи: µ2 (A A) = µ22 (A), где µ2 (A) = A 2 · A† 2 - спектральное число обусловленности матрицы A, · 2 - спектральная матричная норма. Существуют различные подходы для увеличения скорости сходимости итерационных методов решения плохо обусловленных задач [4, 5]. Один из таких подходов основывается на предобуславливании исходной задачи [6-10]. Однако не для всех задач удается получить эффективную процедуру предобуславливания. При этом использование предобуславливателя приводит к сильному увеличению вычислительной сложности алгоритма, что крайне нежелательно при решении задач большой размерности. Для увеличения скорости сходимости итерационных методов Качмажа [11] и Гаусса-Зейделя для нормальных систем уравнений в последнее время также используется подход, основанный на идее рандомизации [12-14]. Однако применение такого подхода не всегда оказывается достаточно эффективным. В таких случаях наиболее перспективным способом увеличения скорости сходимости является блочный вариант этих методов. В данной статье рассматривается блочный вариант метода Гаусса-Зейделя для нормальных систем уравнений [3]. При использовании этого метода возникает необходимость вычисления псевдообратных матриц на каждой итерации, что, в свою очередь, резко повышает вычислительную сложность алгоритма. 731 Ж д а н о в А. И., Б о г д а н о в а Е. Ю. Целью данной работы является получение эффективной вычислительной реализации блочного метода Гаусса-Зейделя для нормальных систем уравнений. В работе предлагается заменить операцию псевдообращения матрицы на каждой итерации на задачу решения нормальных систем с помощью метода Холецкого [15]. Нормальные уравнения, возникающие на каждой итерации метода Гаусса-Зейделя, имеют сравнительно невысокую размерность по сравнению с исходной системой. Это позволяет снизить вычислительную сложность метода и тем самым повысить скорость сходимости алгоритма. 1. Описание метода. Рассмотрим блочный вариант метода Гаусса-Зейделя для нормальных систем уравнений (2): djk = A†j rkj , xjk+1 = xjk + ωdjk , rkj+1 = rkj - (3) ωAj djk , где j = 1, 2, . . . , s, k = 0, 1, . . ., A = (A1 A2 . . . As ); s - количество блоков матрицы A; Aj ∈ Rm×nj , sj=1 nj = n, 0 < ω < 2 - параметр релаксации. Алгоритм (3) можно записать с использованием «микроитераций» [16]: dk = A†j rk , xk+1 = xk + ωEj dk , rk+1 = rk - ωAj dk , (4) где j = j(k) = (k mod p) + 1; {j(k)}∞ k=0 - периодическая последовательность вида 1, 2, . . . , s, 1, 2, . . . , s, . . .; E = (E1 E2 . . . Es ); Ej ∈ Rn×nj ; E - единичная матрица порядка n. Если x0 - произвольный вектор, а r0 удовлетворяет условию согласования k→∞ r0 = b - Ax0 , то, как показано в [1], xk ---→ x∗ , где x∗ - псевдорешение системы уравнений (1). Для вычисления псевдообратной матрицы A†j можно использовать различные алгоритмы, например, сингулярное разложение [17] или метод Гревилля [18]. Данные методы вычисления псевдообратной матрицы имеют много недостатков, которые снижают скорость решения задачи. Пусть rank(Aj ) = nj , тогда псевдообратная матрица каждого блока определяется по формуле [19]: A†j = (Aj Aj )-1 Aj . (5) Однако данная процедура требует обращения матрицы Aj Aj , что, в свою очередь, сильно усложняет вычислительную процедуру. Из (4) и (5) следует, что вычисление dk можно свести к решению линейной системы уравнений Aj Aj dk = Aj rk . (6) Поскольку rank(Aj ) = nj , матрица Aj Aj - симметричная положительно определенная матрица и для решения линейной системы (6) можно использовать метод, основанный на разложении Холецкого: Aj Aj = Lj Lj , 732 Об одной вычислительной реализации блочного метода Гаусса-Зейделя. . . Lj y = Aj rjk , Lj dj = y. Здесь Lj - нижняя треугольная матрица (фактор Холецкого). Таким образом, предлагается вариант реализации блочного метода ГауссаЗейделя для нормальных уравнений, который не требует на каждой итерации вычисления псевдообратной матрицы, а dk вычисляется на основании решения системы (6) с использованием разложения Холецкого, что, в конечном счете, позволит снизить вычислительную сложность и повысить скорость итерационного метода. 2. Тестовые исследования. В данном пункте продемонстрируем преимущество предлагаемого вычислительного варианта блочного метода Гаусса- Зейделя для нормальных систем уравнений по сравнению с методом, требующим непосредственного вычисления псевдообратной матрицы. Численные тестовые исследования проводились с использованием программного математического пакета SciLab. Предлагаемый вариант блочного метода Гаусса-Зейделя проиллюстрируем на следующем примере: A ∈ Rm×n , A = (aij ), элементы aij имеют равномерный закон распределения на интервале (0; 10) и сформированы с помощью генератора случайных чисел, b ∈ Rm , элементы bi вектора b имеют равномерный закон распределения на интервале (0; 10). Проводилось численное исследование для различной размерности блоков матрицы A с целью решения вопроса о выборе оптимального размера блоков. В рассматриваемых примерах параметр релаксации ω = 1, m = 2200, n = 700, число обусловленности матрицы A равно µ2 (A) = 3.1256. В таблице сверху приведены результаты решения данной задачи блочным методом Гаусса-Зейделя с использованием разложения Холецкого. Здесь Ncol - размер блока (количество столбцов в блоке), k - количество всех выРезультаты численного решения нормальных систем уравнений методом Гаусса-Зейделя [Results for numerical solution of normal systems of equations by Gauss-Seidel method] Block size, Ncol Number of iterations, k Time (s) x - x∗ 2 Results for Cholesky decomposition 1 2 14 28 50 700 9800 8750 900 550 266 1 31.0156 26.9531 4.0312 3.8437 3.4375 9.1718 2.5 · 10-6 0.4 · 10-6 5.29 · 10-8 4.68 · 10-8 5.08 · 10-8 5.92 · 10-8 Results for Greville’s method 1 2 14 28 50 700 9800 8750 900 550 266 1 32.5937 26.8593 4.3281 4.4062 4.2187 31.5937 2.5 · 10-6 5.29 · 10-8 4 · 10-7 4.68 · 10-8 6.82 · 10-8 5.92 · 10-8 733 Ж д а н о в А. И., Б о г д а н о в а Е. Ю. полненных итераций, Time - время выполнения алгоритма в секундах, x∗ - точное псевдорешение, которое было найдено SVD-методом [3]. Условием остановки итерационного процесса является критерий: xk+1 - xk xk+1 2 2 10-5 . Проведен также аналогичный численный эксперимент при той же матрице A и том же векторе b, только поставленная задача была решена блочным методом Гаусса-Зейделя с использованием метода Гревилля для вычисления псевдообратной матрицы (результаты представлены в таблице снизу). Согласно данным таблицы, время решения задачи убывает до Ncol = 50, а потом начинает возрастать. Увеличение времени решения задачи методом Гаусса-Зейделя при использовании метода Гревилля для вычисления псевдообратной матрицы происходит более резко, чем в предложенном методе. При увеличении количества столбцов в блоке, то есть увеличении размера Ncol блоков, количество итераций, которые требуются для решения задачи, уменьшается. Более наглядно разницу между двумя численными алгоритмами можно увидеть на рисунке. Видно, что решение задачи блочным методом Гаусса- Зейделя с использованием разложения Холецкого выполняется быстрее, чем тот же алгоритм с использованием метода Гревилля для вычисления псевдообратной матрицы. При этом чем больше размерность блоков, тем больше разница во времени решения задачи между этими двумя алгоритмами. Время решения тестовых задач для нормальных систем уравнений методом Гаусса-Зейделя [Time of the test problems solution for normal systems of equations by Gauss-Seidel method] С увеличением размерности блоков возрастает сложность вычислений для каждого блока, и этот факт объясняет вид графиков на приведенном рисунке. Заключение. В данной статье предлагается оригинальная модификация блочного метода Гаусса-Зейделя для нормальных систем уравнений, которая заключается в замене операции псевдообращения матрицы на задачу ре734 Об одной вычислительной реализации блочного метода Гаусса-Зейделя. . . шения методом Холецкого нормальных систем уравнений для каждого отдельного блока. Приведены численные результаты, демонстрирующие эффективность данного подхода. Показано, что увеличение размерности блока, то есть увеличение числа столбцов в блоке, приводит, соответственно, к снижению числа требуемых итераций. Однако каждая итерация становится вычислительно более сложной, и в этом случае требуется определять оптимальный размер блока. Декларация о финансовых и других взаимоотношениях. Исследование не имело спонсорской поддержки. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами. Авторы не получали гонорар за статью.

About the authors

Alexander I Bogdanova

Samara State Technical University

Email: zhdanovaleksan@yandex.ru
(Dr. Phys. & Math. Sci.), Dean, Faculty of the Distance and Additional Education; Head of Dept., Dept. of Higher Mathematics & Applied Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

Ekaterina Yu Bogdanova

Samara State Technical University

Email: fwinter@yandex.ru
Postgraduate Student, Dept. of Higher Mathematics & Applied computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russian Federation

References

Supplementary files