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

Abstract


This article focuses on the modification of the block option Gauss-Seidel method for normal systems of equations, which is a sufficiently effective method of solving generally overdetermined, systems of linear algebraic equations of high dimensionality. The main disadvantage of methods based on normal equations systems is the fact that the condition number of the normal system is equal to the square of the condition number of the original problem. This fact has a negative impact on the rate of convergence of iterative methods based on normal equations systems. To increase the speed of convergence of iterative methods based on normal equations systems, for solving ill-conditioned problems currently different preconditioners options are used that reduce the condition number of the original system of equations. However, universal preconditioner for all applications does not exist. One of the effective approaches that improve the speed of convergence of the iterative Gauss-Seidel method for normal systems of equations, is to use its version of the block. The disadvantage of the block Gauss-Seidel method for production systems is the fact that it is necessary to calculate the pseudoinverse matrix for each iteration. We know that finding the pseudoinverse is a difficult computational procedure. In this paper, we propose a procedure to replace the matrix pseudo-solutions to the problem of normal systems of equations by Cholesky. Normal equations arising at each iteration of Gauss-Seidel method, have a relatively low dimension compared to the original system. The results of numerical experimentation demonstrating the effectiveness of the proposed approach are given.

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
244, Molodogvardeyskaya st., Samara, 443100, Russian Federation
(Dr. Phys. & Math. Sci.), Dean, Faculty of the Distance and Additional Education; Head of Dept., Dept. of Higher Mathematics & Applied Computer Science

Ekaterina Yu Bogdanova

Samara State Technical University

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

References

  1. Saad Y. Basic Iterative Methods / Iterative Methods for Sparse Liner Systems. Philadelphia, PA, USA: SIAM, 2003. pp. 103-128. doi: 10.1137/1.9780898718003.ch4.
  2. Golub G. H., Van Loan C. F. Matrix Computations / Johns Hopkins Studies in Mathematical Sciences. Baltimore, London: Johns Hopkins University Press, 1996. xxvii+728 pp.
  3. Björck A. Linear Least Squares Problems / Numerical methods in matrix computations / Texts in Applied Mathematics, 59. Berlin: Springer, 2015. pp. 211-430. doi: 10.1007/ 978-3-319-05089-8_2.
  4. Young D., Rheinboldt W. Iterative Solutions of Large Linear Systems. New York: Academic Press, 1971. 572 pp. doi: 10.1016/c2013-0-11733-3.
  5. Ma A., Needell D., Ramdas A. Convergence Properties of the Randomized Extended Gauss-Seidel and Kaczmarz Methods // SIAM. J. Matrix Anal. Appl., 2015. vol. 36, no. 4. pp. 1590-1604. doi: 10.1137/15m1014425.
  6. Gill P. E., Murray W., Ponceleón D. B., Saunders M. A. Preconditioners for Indefinite Systems Arising in Optimization // SIAM. J. Matrix Anal. Appl., 1992. vol. 13, no. 1. pp. 292-311. doi: 10.1137/0613022.
  7. Benzi M. Preconditioning Techniques for Large Linear Systems: A Survey // Journal of Computational Physics, 2002. vol. 182, no. 2. pp. 418-477. doi: 10.1006/jcph.2002.7176.
  8. Benzi M., Tûma M. A comparative study of sparse approximate inverse preconditioners // Appl. Numer. Math., 1999. no. 2-3. pp. 305-340. doi: 10.1016/s0168-9274(98)00118-4.
  9. Bergamaschi L., Pini G., Sartoretto F. Approximate inverse preconditioning in the parallel solution of sparse eigenproblems // Numerical Linear Algebra with Applications, 2000. vol. 7, no. 3. pp. 99-116. doi: 10.1002/(sici)1099-1506(200004/05)7:3<99::aid-nla188>3.3.co;2-x.
  10. Benzi M., Joubert W. D., Mateescu G. Numerical experiments with parallel orderings for ILU preconditioners // Electronic Transactions on Numerical Analysis, 1999. vol. 8. pp. 88-114.
  11. Ильин В. П. Об итерационном методе Качмажа и его обобщениях // Сиб. журн. индустр. матем., 2006. Т. 9, № 3. С. 39-49.
  12. Gower R. M., Richtárik P. Randomized Iterative Methods for Linear Systems // SIAM. J. Matrix Anal. Appl., 2015. vol. 36, no. 4. pp. 1660-1690. doi: 10.1137/15m1025487.
  13. Strohmer T., Vershynin R. A Randomized Kaczmarz Algorithm with Exponential Convergence // J. Fourier Anal. Appl., 2009. vol. 15, no. 2. pp. 262-278. doi: 10.1007/s00041-008-9030-4.
  14. Жданов А. И., Сидоров Ю. В. Параллельная реализация рандомизированного регуляризованного алгоритма Качмажа // Комп. оптика, 2015. Т. 39, № 4. С. 536-541. doi: 10.18287/0134-2452-2015-39-4-536-541.
  15. Horn R. A., Johnson C. R. Matrix Analysis. Cambridge: Cambridge University Press, 1989. xviii+643 pp. doi: 10.1017/cbo9781139020411.
  16. Жданов А. И., Иванов А. А. Проекционный регуляризирующий алгоритм для решения некорректных линейных алгебраических систем большой размерности // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2010. № 21. С. 309-312. doi: 10.14498/vsgtu827.
  17. Малышев А. Н. Введение в вычислительную линейную алгебру. Новосибирск: Наука, 1991. 229 с.
  18. Гантмахер Ф. Р. Теория матриц. М.: Наука, 1966. 576 с.
  19. Беклемишев Д. В. Дополнительные главы линейной алгебры. М.: Наука, 1983. 336 с.

Statistics

Views

Abstract - 27

PDF (Russian) - 19

Cited-By


PlumX

Dimensions

Refbacks

  • There are currently no refbacks.

Copyright (c) 2016 Samara State Technical University

Creative Commons License
This work is licensed under a Creative Commons Attribution 4.0 International License.

This website uses cookies

You consent to our cookies if you continue to use our website.

About Cookies