Об одной вычислительной реализации блочного метода Гаусса-Зейделя для нормальных систем уравнений



Цитировать

Полный текст

Аннотация

Статья посвящена модификации блочного варианта метода Гаусса-Зейделя для нормальных систем уравнений, который является одним из достаточно эффективных методов решения, в общем случае переопределенных, систем линейных алгебраических уравнений большой размерности. Основным недостатком методов, основанных на нормальных системах уравнений, является тот факт, что число обусловленности нормальной системы равно квадрату числа обусловленности исходной задачи. Этот факт отрицательно влияет на скорость сходимости итерационных методов, основанных на нормальных системах уравнений. Для повышения скорости сходимости итерационных методов, основанных на нормальных системах уравнений, при решении плохо обусловленных задач в настоящее время используются различные варианты предобуславливателей, позволяющие снизить число обусловленности исходной системы уравнений. Однако универсального предобуславливателя для всех задач не существует. Одним из эффективных подходов, позволяющих повысить скорость сходимости итерационного метода Гаусса-Зейделя для нормальных систем уравнений, является использование его блочного варианта. Недостатком блочного метода Гаусса-Зейделя для нормальных систем является тот факт, что на каждой итерации необходимо вычислять псевдообратную матрицу. Известно, что нахождение псевдообратной матрицы является достаточно сложной вычислительной процедурой. В настоящей работе предлагается процедуру псевдообращения матрицы заменить на задачу решения нормальных систем уравнений методом Холецкого. Нормальные уравнения, возникающие на каждой итерации метода Гаусса-Зейделя, имеют сравнительно невысокую размерность по сравнению с исходной системой. Приводятся результаты вычислительных экспериментов, демонстрирующие эффективность предлагаемого подхода.

Полный текст

Введение. В статье рассматривается блочный метод Гаусса-Зейделя для нормальных систем уравнений [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 Об одной вычислительной реализации блочного метода Гаусса-Зейделя. . . шения методом Холецкого нормальных систем уравнений для каждого отдельного блока. Приведены численные результаты, демонстрирующие эффективность данного подхода. Показано, что увеличение размерности блока, то есть увеличение числа столбцов в блоке, приводит, соответственно, к снижению числа требуемых итераций. Однако каждая итерация становится вычислительно более сложной, и в этом случае требуется определять оптимальный размер блока. Декларация о финансовых и других взаимоотношениях. Исследование не имело спонсорской поддержки. Все авторы принимали участие в разработке концепции статьи и в написании рукописи. Авторы несут полную ответственность за предоставление окончательной рукописи в печать. Окончательная версия рукописи была одобрена всеми авторами. Авторы не получали гонорар за статью.
×

Об авторах

Александр Иванович Жданов

Самарский государственный технический университет

Email: zhdanovaleksan@yandex.ru
(д.ф.-м.н., проф.; zhdanovaleksan@yandex.ru), декан, факультет дистанционного и дополнительного образования; заведующий кафедрой, каф. Высшей математики и прикладной информатики Россия, 443100, Самара, ул. Молодогвардейская, 244

Екатерина Юрьевна Богданова

Самарский государственный технический университет

Email: fwinter@yandex.ru
аспирант, каф. высшей математики и прикладной информатики Россия, 443100, Самара, ул. Молодогвардейская, 244

Список литературы

  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 с.

Дополнительные файлы

Доп. файлы
Действие
1. JATS XML

© Самарский государственный технический университет, 2016

Creative Commons License
Эта статья доступна по лицензии Creative Commons Attribution 4.0 International License.