CC BY
28
DOI: https://doi.org/10.23670/IRJ.2018.78.12.002
О ЧИСЛЕННОМ РЕШЕНИИ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ С ПЛОХО
ОБУСЛОВЛЕННЫМИ МАТРИЦАМИ
Научная статья
Рябов В.М.1, Бурова И.Г.2' *, Кальницкая М.А.3, Малевич А.В.4, Лебедева А.В.5, Борзых А.Н.6
1 ORCID: 0000-0002-1364-5428;
2 ORCID: 0000-0002-8743-1377;
3 ORCID: 0000-0001-5671-5070;
4 ORCID: 0000-0003-0753-4621;
5 ORCID: 0000-0001-9026-0292;
6 ORCID: 0000-0001-5489-911X;
i, 2, з, 4, 5, 6 Санкт-Петербургский государственный университет, Санкт-Петербург, Россия
* Корреспондирующий автор (i.g.burova[at]spbu.ru)
Аннотация
Представлены результаты численного решения систем линейных алгебраических уравнений (СЛАУ) с симметричными и несимметричными плохо обусловленными матрицами методом регуляризации. Рассматриваются положительно определенные, а также осцилляционные матрицы. В статье показано, что для регуляризации вычислительного процесса по методу Тихонова достаточно заменить матрицу Ап системы матрицей An + аЕ п где Е п - единичная матрица, а а - некоторое положительное число (параметр регуляризации), которое стремится к нулю.
Ключевые слова: плохо обусловленные системы линейных алгебраических уравнений, матрицы Гильберта, параметр регуляризации.
ON NUMERICAL SOLUTION OF SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS WITH ILL-
CONDITIONED MATRICES
Research article
Ryabov V.M.1, Burova I.G.2' *, Kalnitskaya M.A.3, Malevich A.V.4, Lebedeva A.V.5, Borzykh A.N.6
1 ORCID: 0000-0002-1364-5428;
2 ORCID: 0000-0002-8743-1377;
3 ORCID: 0000-0001-5671-5070;
4 ORCID: 0000-0003-0753-4621;
5 ORCID: 0000-0001-9026-0292;
6 ORCID: 0000-0001-5489-911X;
1 2, 3 4 5, 6 St. Petersburg State University, St. Petersburg, Russia
* Corresponding author (i.g.burova[at]spbu.ru)
Abstract
The results of the numerical solution of systems of linear algebraic equations (SLAE) with symmetric and asymmetrical ill-conditioned matrices by the regularization method are presented in the paper. Positive-definite and oscillatory matrices are considered. The article shows that in order to regularize the computational process according to the Tikhonov method, it is enough to replace the system matrix A_n with the matrix A_n + aE_n where E_n is the identity matrix, and а is some positive number (regularization parameter) that tends to zero.
Keywords: ill-conditioned systems of linear algebraic equations, Hilbert matrices, regularization parameter.
Введение
При решении различных задач возникает необходимость решать плохо обусловленные СЛАУ с положительно определенными симметричными матрицами. Такие системы линейных алгебраических уравнений возникают, например, когда функцию f (x) аппроксимируют алгебраическими многочленами в метрике пространства L 2( 0 , 1 ): если мы возьмем полином в виде и определим погрешность аппроксимации как
то из условия минимизации придем к СЛАУ с матрицей Гильберта. Такие системы линейных алгебраических уравнений также возникают при решении обыкновенных дифференциальных уравнений методом Ритца, что приводит к матрицам Грама. Эти матрицы размерности П являются симметричными и положительно определенными, но при неограниченном увеличении n наименьшее собственное значение может стремиться к нулю, что приводит к неустойчивости решения.
Обычно для получения надежного решения используют методы регуляризации. Общей стратегией является использование стабилизатора Тихонова [1] или его модификаций [2], [3], [4], [5], [6], [7]. В этой статье рассматриваются особенности численного решения СЛАУ с положительно определенной симметричной матрицей с использованием регуляризации. В следующих разделах показано, что для регуляризации вычислительного процесса по методу Тихонова достаточно заменить матрицу системы матрицей , где - единичная матрица, а -
положительное число (параметр регуляризации), стремящееся к нулю. Таким образом, мы уменьшаем число обусловленности системы линейных алгебраических уравнений, что увеличивает устойчивость.
Постановка задачи
Пусть A - невырожденная вещественная квадратная матрица порядка п . В этом случае решение СЛАУ A z = / существует и единственно. Известны различные модификации метода Гаусса для решения СЛАУ, например, метод Гаусса с выбором ведущего элемента (в столбце или в матрице) и другие. Предположим, что число обусловленности с оп d (A ) = | | A | | ■ | | A " 1 | | матрицы A велико, т.е. матрица системы уравнений плохо обусловлена. Решение плохо обусловленных СЛАУ методом Гаусса не всегда дает удовлетворительное решение. Например, пусть
/0.0000001 333 555\ / 888 \ А = 33333333 1 70 ,/= 33333404 V55555555 70 32/ V55555657/
Число обусловленности матрицы A приближенно равно 107. Решая эту систему методом Гаусса без перестановок (с помощью программы, написанной на C++ с вещественными числами типа double), получим
, что существенно отличается от точного решения . Подобные
примеры были рассмотрены в [6]. Эти примеры показывают, что в процессе решения необходимо избегать деления на малые по абсолютной величине элементы. Избежать эту ситуацию помогает использование модифицированного метода Гаусса, заключающегося в выборе ведущего элемента, являющегося наибольшим по абсолютному значению элементом в столбце (стратегия Уилкинсона) или по всей матрице (стратегия полного упорядочения Жордана) [2]. Применение метода Гаусса с выбором ведущего элемента по столбцам дает решение z = ( 1 . 0, 1 .0, 1 .0) т. Если система уравнений является плохо обусловленной, например, в случае СЛАУ с матрицей Гильберта порядка П с элементами h = 1 / (, + j — 1), то практически невозможно получить приемлемое решение СЛАУ с
помощью известных методов (прямыми методами, такими как метод Гаусса, метод квадратного корня, итерационными методами и др.).
В таблице 1 приведены числа обусловленности матриц Гильберта порядков от 3 до 20, полученные с помощью пакета Maple (Digits=50). В таблице 2 представлены решения СЛАУ A z = / с матрицами Гильберта #„, полученные методом Гаусса (стратегия Уилкинсона). Здесь представлены решения, полученные при решении систем уравнений для , вычисленные в C++ с числами типа double. Решения, представленные в таблице 2, далеки от
истинных решений. В следующих разделах представлен метод регуляризации, с помощью которого получены решения исходной системы с матрицами Гильберта при использовании нормы
П
| \HJ | = max,. X \hij\.
j=i
Таблица 1 - Числа обусловленности матриц Гильберта Нп
п cond (Нп) п cond (Нп)
3 748 12 4.11541 016
4 28375 13 1.3244 1 018
5 9.4366 1 05 14 4.5378 1 019
6 2.9070 1 07 15 1.5392 1 02 1
7 9.85191 08 16 5.0628 1 022
8 3.3873-1 010 17 1.6808 1 024
9 1.0996 1 012 18 5.7661 1 025
10 3.5357-1 013 19 1.9258 1 02 7
11 1.2337 1 015 20 6.2836 1 028
Таблица 2 - Решения СЛАУ с матрицами Нп.
п=10 п=12 п=14 п=20
1.000000 1.000000 1.000000 1.000000
1.000000 1.000003 0.999990 1.000075
0.999998 0.999916 1.000310 0.997232
1.000020 1.001125 0.996220 1.043007
0.999903 0.991859 1.019948 0.657461
1.000267 1.035385 0.981738 2.485524
0.999563 1.000421 0.999779 0.902315 1.175419 0.795768 0.675197 2.922896 -4.498100 -2.218021 2.012694 11.914026 -20.728838 5.732215 34.233277 -42.153623 18.727738 2.721547 -16.803379 49.405743 -57.516277 33.288344 -5.798746
1.0000480 1.148663 10.515965
0.938525 1.011022 -9.428208 8.094763 -1.740964 1.460245
Численное решение плохо обусловленной системы уравнений
Различные подходы к решению систем уравнений с плохо обусловленными матрицами известны (см. [1-7]). В данной работе для получения приемлемого решения СЛАУ рассматривается применение модификации метода регуляризации. Известный стандартный метод регуляризации Тихонова позволяет найти нормальное решение системы Аг =/. Этот метод основан на поиске элемента, на котором функционал Ма (г, А,/) = | | Аг — / | | 2 + а | | г | | 2 достигает наименьшего значения для фиксированного положительного а. Для этого необходимо решить уравнение Эйлера , где - сопряженная матрица. Решение уравнения Эйлера зависит от числа
обусловленности матрицы А * А. Это число может быть очень большим. Если матрица А симметрична и положительно определена, мы предлагаем найти нормальное решение другим способом. В данной работе предлагается найти нормальное решение путем решения системы уравнений, для которой число обусловленности значительно меньше. Пусть матрица СЛАУ
Az=f (1)
является симметричной и положительно определенной (например, матрица Гильберта). В этом случае существует единственный положительно определенный корень из матрицы В = л/А, т.е. положительно определенная матрица В такая, что В2 = А. Установим связь между собственными значениями и собственными векторами матриц А и В: пусть ц и х — собственное значение и собственный вектор матрицы В:
Б=цх (2)
Умножив (2) на Б, получаем
Ах=цБх (3)
Используя (2), перепишем (3) как Ах = ц2х. Это равенство означает, что собственные векторы матриц А и В одинаковы, а собственные значения матрицы А равны квадратам собственных значений матрицы Б. Умножив (1) на , получаем
Вг=В"1/ (4)
Запишем для (4) уравнение Эйлера для минимизации сглаживающего функционала Ма (г, В ,В " 1 /) = | | Вг — В—1/2 + аг2: оно имеет вид
(Б*Б+аБК=Б*(Б-14 а>0 (5)
В случае симметричной матрицы матрица является самосопряженной, поэтому мы получаем уравнение
(А+аЕ) Ъ а=f (6)
Формально в исходной системе осуществляется сдвиг, но фактически это метод регуляризации для уравнения (1). В памяти компьютера числа обычно хранятся с некоторыми погрешностями. Будем считать, что матрица и вектор
даны приближенно, т.е. вместо матрицы А и правой части / известны Л8 и такие, что | | А — А § | | < 5, | | / — /§ | | < 5. Говорить о сходимости можно лишь при неограниченном увеличении точности исходной информации, т.е. при 5^0. Однако практически мы не можем неограниченно уменьшать 5 , из-за чего результаты могут быть далеки от желаемых решений. Решая задачу (6), получаем приближенное решение системы (1).
Замечание. Методы вычисления квадратного корня матрицы описаны в [8-11]. Если матрица А симметрична, нам не нужно вычислять матрицу . Применение регуляризации непосредственно к уравнению (1) просто увеличит число обусловленности результирующей системы, что невыгодно. В случае несимметричной положительно определенной матрицы уравнение Эйлера для системы (1) имеет вид (5). Сходимость метода регуляризации изложена в [1]. В отличие от стандартного подхода процедуры регуляризации представление (5) имеет меньшее число обусловленности, что очень важно. Но, в отличие от симметричных матриц, необходимо знать матрицу В в случае несимметричной матрицы. Последнее требование усложняет применение данного метода на практике.
Существуют классы несимметричных матриц, все собственные значения которых положительны. Покажем, как в этом случае можно осуществить процесс регуляризации в форме (5), а не в традиционной форме (Л Л + аЕ)г = Л /, а> 0, что, как уже упоминалось выше, существенно уменьшает число обусловленности решаемой системы. Рассматривается класс осцилляционных матриц (см. [12]), все собственные значения которых положительны и попарно различны, а соответствующие собственные векторы обладают особыми свойствами колебаний.
Пусть матрица А „ — осцилляционная, а V есть матрица собственных векторов А п, Л = [Я 1 ,. . ,Д„] — диагональная матрица собственных чисел Ап . Тогда Ап V = VАn или Ап = V ЛV " 1. Пусть Л 1 = [лД!,. . .,,/Я^] . Положим В = VЛ1V" 1. Понятно, что . Таким образом, процесс регуляризации может быть
осуществлен в форме (5). Обобщенная матрица Вандермонда, то есть матрица вида ,
, относится к классу осцилляционных матриц. Для регуляризации СЛАУ с такими матрицами применима описанная выше схема.
Результаты применения метода регуляризации
Проведена серия вычислительных экспериментов по применению метода регуляризации для решения СЛАУ с матрицами Гильберта порядков п=2,3,...,20. В таблицах 3, 4 показаны результаты применения метода регуляризации для различных параметров для решения возмущенной СЛАУ , где матрица
исходной системы Яп = ( /гу) 1 есть матрица Гильберта порядка п . Решение возмущенной системы получено методом Гаусса с выбором ведущего элемента по столбцам. Точное решение исходной невозмущенной системы есть п -мерный вектор из единиц: . Вычисляя решение возмущенной СЛАУ при
различных значениях , находим значение параметра, при котором погрешность решения имеет минимальное значение. Параметр а, соответствующий решению с наименьшей нормой, назовем оптимальным. Таким образом, для решения системы уравнений (1) необходимо решить несколько систем уравнений для различных а. В таблицах 3, 4
полужирным шрифтом выделено наименьшее значение нормы погрешности для заданного п, соответствующее оптимальному значению параметра возмущения а. В последней строке этих таблиц приведены нормы погрешностей решений, полученных без регуляризации. Погрешность решения вычислялась с помощью евклидовой нормы.
Таблица 3 - Нормы погрешностей решений СЛАУ с матрицами Гильберта порядка п = 1 0 , 1 2 методом (6) при
различных значениях параметра а
а 71 = 10 тг = 12
ю-12 0.13 1 0"5 0.22 ■ 1 0" 5
10"11 0.14 • 1 0 " 6 0.69 • 1 0 " 7
10-ю 0.16 0.19
10"8 0.16 0.20
10"7 0.58 ■ 1 0 " 3 0.58
10"6 0.17 ■ 1 0 " 2 0.19 ■ 1 0 " 2
10"5 0.56 0.63
10-4 0.18 0.19
10"3 0.56 0.61
10"2 0.1799 0.1985
10"1 0.5788 0.6355
Нормы пог] решностей решений без применения метода регуляризации
0.71 ■ 1 0"3 0.3306
Точность арифметики с плавающей запятой можно охарактеризовать машинным е, т.е. наименьшим положительным числом с плавающей запятой е, для которого 1 + е >1 (см. [6]). Мы можем вычислить машинное е. Например, в 64-битном С++ переменные типа double дают е = 2.2 ■ 1 0 " 1 6.
Таблица 4 - Нормы погрешностей решения СЛАУ с матрицами Гильберта порядков п = 1 4,2 0 методом (6) при _различных значениях параметра а_
а 71 = 14 тг = 20
10"12 0.19 ■ 1 0"4 0.17 ■ 1 0"4
10"11 0.27 0.23 • 1 0 " 6
10-ю 0.20 • 1 0 " 6 0.25
10"8 0.20 0.25
10"7 0.65 0.78
10"6 0.21 ■ 1 0" 2 0.25 ■ 1 0" 2
10"5 0.66 0.80
10-4 0.21 0.25
10"3 0.67 0.81
10"2 0.2153 0.2581
10"1 0.6872 0.8220
Погрешности решении без регуляризации
17.0703 105.2819
Теорема Тихонова утверждает, что теоретически, по мере уменьшения а, регуляризованное решение улучшается, но в практических расчетах для достаточно малых а (в пределах точности машины в C++) погрешности округления и число обусловленности матрицы имеют значительное влияние. Это можно увидеть, изучив результаты, представленные в начале таблиц 3, 4.
Заключение
В данной работе представлены результаты численного решения СЛАУ с положительно определенными симметричными (или несимметричными, но почти симметричными) плохо обусловленными матрицами модифицированным методом регуляризации. Показано, что решение СЛАУ с матрицами Гильберта с использованием метода регуляризации может быть существенно улучшено.
Конфликт интересов Conflict of Interest
Не указан. None declared.
Список литературы / References
1. Тихонов А. Н. Методы решения некорректных задач / А. Н. Тихонов, В. Я. Арсенин. - М.: Наука. Главная редакция физико-математической литературы, 1979. Изд. 2-е. - 284 с.
2. Фаддеев Д. К. Вычислительные методы линейной алгебры / Д. К. Фаддеев, В. Н. Фаддеева. - М.: Физматгиз, 1960. - 656 с.
3. Фаддеев Д.К. О плохо-обусловленных системах линейных уравнений / Д. К. Фаддеев, В. Н. Фаддеева // Журнал вычислительной математики и математической физики. - 1961. - Т. 1. - №3. - С. 412-417.
4. Гавурин М.К. О плохо-обусловленных системах линейных алгебраических уравнений / М. К. Гавурин // Журнал вычислительной математики и математической физики. - 1962. - Т. 2. - №3. - С. 387-397.
5. Гавурин М. К. Применение полиномов Чебышева при регуляризации некорректных и плохо обусловленных уравнений в гильбертовом пространстве / М. К. Гавурин, В. М. Рябов // Журнал вычислительной математики и математической физики. - 1973. - Т. 13. - №6. - С. 1599-1601.
6. Форсайт Дж. Численное решение систем линейных алгебраических уравнений / Дж. Форсайт, К. Молер. Перевод с английского В.П. Ильина и Ю.И. Кузнецова. Под редакцией Г.И. Марчука. - М.: Мир, 1969. - 167 с.
7. Кабанихин С. И. Обратные и некорректные задачи / С.И. Кабанихин. - Новосибирск: Сибирское научное издательство, 2009. - 457 с.
8. Bjorck A. A Schur method for the square root of a matrix / A. Bjorck, S. Hammarling // Linear Algebra Appl. - 1983. -V. 52/53, P. 127-140.
9. Higham Nicholas J. Functions of matrices: Theory and computation / J. Higham Nicholas - Philadelphia: Society for Industrial and Applied Mathematics, 2008. - 425 p.
10. Deadman E. Blocked Schur algorithms for computing the matrix square root / E. Deadman, N.J. Higham, R. Ralha // Lecture Notes in Computer Science, V. 7782, Springer-Verlag. - 2012. P. 171-182.
11. Higham Nicholas J. Newton's Method for the Matrix Square Root / J. Higham Nicholas // Mathematics of computation. - 1986. - V. 46. - №174. - P. 537-549.
12. Гантмахер Ф. Р. Осцилляционные матрицы и ядра и малые колебания механических систем / Ф.Р. Гантмахер, М.Г. Крейн // Москва-Ленинград: Государственное издательство технико-теоретической литературы [ГИТТЛ], 1950. 359 с.
Список литературы на английском языке/ References in English
1. Tikhonov A. N. Solutions of Ill-Posed Problems / A. N. Tikhonov, V. Ya.Arsenin - Halsted, New York, 1977
2. Faddeev D. K. Vychislitelnye metody linejnoj algebry [Computational methods of linear algebra] / D. K. Faddeev, V. N. Faddeeva. - Moskow: Fizmatgiz, 1960. - 656 P. [in Russian]
3. Faddeev D. K. O plokho-obuslovlennykh sistemakh linejnykh uravnenij [On ill-posed systems of linear equations] / D. K. Faddeev, V. N. Faddeeva // Zhurnal vychislitelnoj matematiki i matematicheskoj fiziki [Journal of Computational Mathematics and Mathematical Physics]. - 1961. - V. 1. - №3. - P. 412-417. [in Russian]
4. Gavurin M. K. O plokho-obuslovlennykh sistemax linejnykh algebraicheskikh uravnenij [On ill-posed systems of linear algebraic equations] / M. K. Gavurin // Zhurnal vychislitelnoj matematiki i matematicheskoj fiziki [Journal of Computational Mathematics and Mathematical Physics]. - 1962. - V. 2. - №3. - P. 387-397. [in Russian
5. Gavurin M. K. Primenenie polinomov Chebysheva pri regulyarizacii nekorrektnykh i plokho obuslovlennykh uravnenij v gilbertovom prostranstve [Application of Chebyshev polynomials in the regularization of ill-posed and ill-conditioned equations in a Hilbert space] / M. K. Gavurin, V. M. Ryabov // Zhurnal vychislitelnoj matematiki i matematicheskoj fiziki [Journal of Computational Mathematics and Mathematical Physics]. - 1973. - V. 13. - №6. - P. 1599-1601. [in Russian
6. Forsythe George E. Computer Solution of Linear Algebraic Systems / George E.Forsythe & Cleve B.Moler - Prentice-Hall, Englewood Cliffs, N. J. 1967
7. Kabanikhin S. I. Obratnye i nekorrektnye zadachi [Inverse and incorrect tasks] / S. I. Kabanikhin. - Novosibirsk: Sibirskoe nauchnoe izdatelstvo [Siberian Scientific Publishing], 2009. - 457 P. [in Russian
8. Bjorck A. A Schur method for the square root of a matrix / A. Bjorck, S. Hammarling // Linear Algebra Appl. - 1983. -V. 52/53, P. 127-140
9. Higham Nicholas J. Functions of matrices: Theory and computation / J. Higham Nicholas - Philadelphia: Society for Industrial and Applied Mathematics, 2008. - 425 P.
10. Deadman E. Blocked Schur algorithms for computing the matrix square root [Text] / E. Deadman, N. J. Higham, R. Ralha // Lecture Notes in Computer Science, V. 7782, Springer-Verlag. - 2012. P. 171-182.
11. Higham Nicholas J. Newton's Method for the Matrix Square Root / J. Higham Nicholas // Mathematics of computation. - 1986. - V. 46. - №174. - P. 537-549.
12. Gantmakher F. R. Oscillation matrices and kernels and small vibrations of mechanical systems. / F. R.Gantmakher, M.G.Krein - Dept. Commerce USA. Joint Publ. Service - 1961 (Translated from Russian)
