CC BY
101
УДК 519.61 Вестник СПбГУ. Математика. Механика. Астрономия. 2019. Т. 6 (64). Вып. 4
МБС 65Р22
О численном решении систем линейных алгебраических уравнений с плохо обусловленными матрицами
А. В. Лебедева, В. М. Рябов
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Лебедева А. В., Рябов В. М. О численном решении систем линейных алгебраических уравнений с плохо обусловленными матрицами // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2019. Т. 6 (64). Вып. 4. С. 619626. https://doi.org/10.21638/11701/spbu01.2019.407
Рассматриваются системы линейных алгебраических уравнений (СЛАУ). Если матрица системы невырождена, то существует единственное решение системы. В вырожденном случае система может не иметь решения или иметь бесконечно много решений. В этом случае вводится понятие нормального решения. Случай невырожденной квадратной матрицы теоретически можно считать хорошим в смысле существования и единственности решения. Однако в теории вычислительных методов невырожденные матрицы подразделяют на две категории: «плохо обусловленные» и «хорошо обусловленные». Плохо обусловленными называют матрицы, для которых решение системы уравнений практически является неустойчивым. Одной из важных характеристик практической устойчивости решения системы линейных уравнений является число обусловленности. Обычно для получения надежного решения используют методы регуляризации. Общей стратегией является использование стабилизатора Тихонова или его модификаций, либо представление искомого решения в виде ортогональной суммы двух векторов, один из которых определяется устойчиво, а для поиска второго необходима некая процедура стабилизации. В настоящей статье рассматриваются методы численного решения СЛАУ с положительно определенной симметричной матрицей или с матрицей осцилляционного типа с использованием регуляризации, приводящие к СЛАУ с уменьшенным числом обусловленности.
Ключевые слова: система линейных алгебраических уравнений, некорректные задачи, плохо обусловленные задачи, число обусловленности, метод регуляризации.
1. Общая теория решения некорректных и плохо обусловленных уравнений изложена в фундаментальных работах [1—3], вопросы реализации общей теории применительно к конкретным прикладным задачам — в книге [4] и многочисленных статьях.
Рассмотрим конечномерный случай решения системы линейных алгебраических уравнений (СЛАУ). Дана система
Аг = и, (1)
где А — вещественная квадратная матрица, и — заданный вещественный вектор и 2 — искомый вектор размерности п. Матрица и правая часть в уравнении (1) могут
(¡5 Санкт-Петербургский государственный университет, 2019
быть заданы приближениями к ним As и us такими, что
||ДЛ|| = ||A - As|| < S, ||Ди|| = ||u - us|| < S.
Эти погрешности повлекут ошибку Дх определения вектора х. В качестве векторной нормы используем евклидову норму, а в качестве матричной — подчиненную ей норму [5]
\\А\\ = ^(Ш- (2)
Введем относительные погрешности
_ ||ДЛ|| Л -UM Л -ИМ
\\а\\ ' INI ' INI •
Заметим, что эти погрешности присутствуют практически всегда, поскольку вычисления, как правило, проводятся с конечным числом значащих цифр.
Если матрица A невырождена, то существует единственное решение системы (1) и можно ввести в рассмотрение число обусловленности ^(A) = cond(A) = ||A|| •
||A-11|.
Важность этой характеристики задачи видна из следующего результата.
Теорема 1 (см. [5, с. 111]). Для произвольных Su и достаточно малых Sa справедлива оценка погрешности
Следовательно, независимо от метода решения СЛАУ погрешность результата может оказаться примерно в ^(A) раз больше погрешности исходной информации.
Приведем пример классической задачи подобного типа: пусть непрерывная функция f аппроксимируется алгебраическим многочленом в метрике пространства L2[0,1]; будем искать многочлен в виде
pk xk
k=0
и определим погрешность аппроксимации как
1
Pn (x) = ^ Pk xk (3)
Условия ее минимизации
E = f (f (x) - Pn(x))2 dx.
0
dE
— =0, j = 0,1,..., n, dPj
приводят к СЛАУ с матрицей Гильберта
1
H
n+1
k + j + Ч kj=0
Эта матрица симметрична и положительно определена, но при возрастании п наименьшее собственное число матрицы Ип+\ стремится к нулю, что приводит к
n
быстрому росту числа обусловленности и неустойчивости решения СЛАУ (так, cond(Hю) « 3.5 • 1013). Следовательно, на этом пути вряд ли удастся получить приближение с малой погрешностью.
Р. Беллман говорил: «Теорию матриц можно считать арифметикой высшей математики. Решение линейной системы является самой фундаментальной в анализе проблемой. К этой проблеме пытаются сводить многие другие задачи и иногда, к своему ужасу, достигают успеха».
2. Напомним, что задача решения уравнения (1) называется корректно поставленной, если ее решение существует для любой правой части и, оно единственно и непрерывно зависит от и. Задача называется плохо обусловленной, если ее решение неустойчиво или велико ее число обусловленности.
Такие системы возникают, например, в теории и практике приближения функций, при решении задач математической физики вариационными методами [6].
Будем считать, что уравнение (1) задано в гильбертовом пространстве H и его приближенное решение разыскивается в виде линейной комбинации нескольких первых элементов линейно независимой и полной в энергетическом пространстве Ha координатной системы {wk}CC=i. Матрица СЛАУ в таком случае есть матрица Грама [7]
Гп = ([wk,Wj Wj .
При этом «качество» СЛАУ определяется разумным выбором координатных систем [7].
Определение 1. Множество элементов гильбертова пространства называется минимальной в этом пространстве системой, если вычеркивание любого элемента этого множества сужает натянутое на него пространство.
Приведем важный пример неминимальной системы. По теореме Г. Мюнца последовательность {xpk} (k = 1, 2,...), 0 ^ pi < Р2 < рз ..., полна в ¿2[0,1], если ряд
СО ..
Е1
расходится.
Очевидно, вычеркивание любого элемента этой системы не влияет на сходимость ряда и не сокращает натянутого на остальные элементы подпространства, потому эта система элементов не минимальна.
Матрица Грама симметрична и положительно определена и ее собственные числа Аук положительны; запишем их в порядке возрастания
0 < а1п) < А2п) < ••• < АПп).
Определение 2. Множество элементов гильбертова пространства называется сильно минимальной в этом пространстве системой, если
inf а(п) = lim а(п) = А0 > 0.
Из определения (2) следуют равенства
ЦГ„|| = А("), НГ-1!^-¿у.
А1
Следовательно, нормы обратных матриц Грама, порожденных сильно минимальной системой, ограничены.
Определение 3. Множество элементов гильбертова пространства называется почти ортонормированной в этом пространстве системой, если существуют такие постоянные Ао и Ло, что для любых п и т ^ п выполняются неравенства
Ас < Ат < Ло.
Очевидно, для почти ортонормированной системы нормы матриц Гп, Г—1 ограничены, а вместе с ними ограничены и числа еопй(Гп) = || Гп || • || Г—11|. В частности, любая ортонормированная система является и почти ортонормированной при Ао = Ло.
При решении конкретной задачи во избежание плохой обусловленности получаемой СЛАУ необходимо использовать координатную систему с наилучшими свойствами.
В рассмотренном выше примере построения аппроксимирующего многочлена в виде (3) в пространстве ^[0,1] была использована неминимальная система {хк}^=о, что привело к СЛАУ с плохо обусловленной матрицей Гильберта. Если же взять в качестве координатной систему многочленов Лежандра {Ь^(2х — 1)}£=о и искать многочлен в виде Рп(х) = ^к=о СкЕи(2х — 1), то придем к СЛАУ с диагональной матрицей, решение которой не составляет никакого труда.
3. В общем случае СЛАУ (1) может не иметь решения. Тогда целесообразно ввести в рассмотрение псевдорешения и нормальные решения уравнения (1).
Определение 4. Псевдорешением уравнения (1) называется любой элемент г, минимизирующий величину ЦАг — и||.
Определение 5. Нормальным 'решением уравнения (1) называется псевдорешение с наименьшей нормой.
Очевидно, нормальное решение существует и единственно; если уравнение (1) имеет единственное решение, то нормальное решение совпадает с ним.
Для нахождения нормального решения необходимы процедуры регуляризации, в которых параметры регуляризации выбираются в зависимости от точности исходной информации и обусловленности исходной задачи.
Далее описаны некоторые алгоритмы, которые формально укладываются в общую схему регуляризации А. Н. Тихонова [1], но имеют меньшую обусловленность.
В методе Тихонова вводится функционал
Ма(А,г,и) = ||Аг — и ||2 + а||г||2, а> 0, (4)
и ставится задача нахождения элемента га, минимизирующего этот функционал. Такой элемент однозначно определяется как решение уравнения Эйлера
(А* А + аЕ)га = А* и. (5)
Теорема 2 (см. [1, с. 119]). Пусть г° есть нормальное решение системы Аг = и и вместо вектора и мы имеем вектор и$ такой, что ||и —ид || ^ 5. Пусть далее в1(5) и в2(5) — какие-либо непрерывные и положительные на [0, 52] функции, монотонно стремящиеся к нулю при 5 ^ 0 и такие, что
с
Тогда для любой положительной на [0, ö2] функции а = a(ö), удовлетворяющей условию
с
(6)
векторы za(Ä), являющиеся 'решением уравнения Эйлера
(A*A + a(ö)E)za(<5) = A*u, (7)
сходятся к нормальному решению системы Az = и при ö ^ 0, т. е.
lim ||za(<5) - zo|| =0.
Замечание. Если система (1) имеет единственное решение, то неравенство (6) можно заменить на такое:
^«(¿к/ад, (8)
откуда, в частности, при ßi(S) = (S) = S из (8) находим a(S) = 6. Аналогично из (6) при ßi(S) = /?2(S) = VS находим а(6) = у/5.
Переход от уравнения (1) к уравнению Эйлера (7) можно представить в виде двух шагов: сперва переходим от уравнения Az = и к уравнению
A* Az = A*u, (9)
а затем делаем сдвиг
(A*A + aE)z = A*u.
Уравнение (9) в отличие от исходного всегда разрешимо, но его число обусловленности практически возводится в квадрат по сравнению с исходным. Эта же черта характерна и для уравнения Эйлера при малых а. Наша задача заключается в построении алгоритмов регуляризации уравнения (1), лишенных этого недостатка для некоторого класса матриц.
4. Регуляризация СЛАУ с симметричной положительно определенной матрицей. Для симметричной положительно определенной матрицы A существует единственный положительно определенный корень из нее, т. е. симметричная матрица B такая, что B2 = A.
Пусть ^ и x — собственное число и собственный вектор матрицы B, т.е.
Bx = y«x. (10)
После умножения этого равенства слева на B получим Ax = ^Bx. Используя (10), получим Ax = y«2x. Это равенство означает, что собственные векторы матриц A и B совпадают, а собственные числа матрицы B суть корни из собственных чисел матрицы A.
Умножая исходное уравнение Az = и слева на матрицу B-1, получим
Bz = B-1u.
Запишем функционал (4) для этого уравнения:
Ma(B,z,B-1u) = ||Bz - B-1u||2 + а||z||2, а > 0,
уравнение Эйлера (5) для него примет вид
(В*В + аЕ)га = В*(В—1и) (11)
или в силу самосопряженности В
(А + аЕ)га = и. (12)
Итак, метод регуляризации А.Н.Тихонова для уравнения (1) с симметричной положительно определенной матрицей вместо общего уравнения (5) может быть сведен к более простой системе (12) (фактически это метод регуляризации М. М. Лаврентьева). Основное преимущество такого подхода состоит в том, что число обусловленности при этом не возрастает.
5. Регуляризация СЛАУ в случае осцилляционных матриц. Существуют классы несимметричных матриц, все собственные значения которых положительны. Покажем, как в этом случае можно осуществить процесс регуляризации в форме (11), а не в традиционной форме, что, как уже упоминалось выше, существенно уменьшает число обусловленности решаемой системы. Рассмотрим класс осцилля-ционных матриц, появляющихся при исследовании малых колебаний механических систем, все собственные значения которых положительны и попарно различны, а соответствующие собственные векторы обладают особыми свойствами колебаний (см. [8, 9]). Справедлива следующая
Теорема 3 (см. [9, с. 400]).
1. Осцилляционная матрица всегда имеет п различных положительных собственных чисел А1 > А2 > . . . , Ап > 0.
2. У первого собственного вектора матрицы, отвечающего наибольшему собственному числу, все координаты отличны от нуля и одного знака. У собственного вектора, отвечающего второму по величине собственному числу, в ряду координат имеется точно одна перемена знака, и вообще в ряду координат к-го собственного вектора, соответствующего к-му собственному числу, имеется точно к — 1 перемена знака (к = 1, 2,..., п).
Пусть А — осцилляционная матрица, а V есть матрица ее собственных векторов, Л = [А1,..., Ап] —диагональная матрица собственных чисел матрицы А. Тогда AV = УА или А = УАУ-1. Пусть Л1 = [\ДГ,... ,л/К]. Положим В = УЛ^-1. Понятно, что В2 = VЛ1V —1 = VЛ2V—1 = А.
Следовательно, В есть корень из А, и собственные векторы этих матриц совпадают.
Таким образом, процесс регуляризации может быть осуществлен в форме (11).
Число обусловленности задачи (11) существенно меньше, чем у стандартной формы метода регуляризации (5). Отметим, что матрица В*В-1 в правой части уравнения (11) невелика (в случае симметрии А она равна единичной матрице) и не сильно увеличивает погрешность вектора и. Методы вычисления необходимого нам квадратного корня матрицы описаны в работе [10]. Заметим, что для осцилляцион-ных матриц с успехом можно использовать метод Ньютона в форме
Вк+1 = ^(Вк+В^А), к = 0,1,..., В0=А.
Обобщенная матрица Вандермонда, то есть матрица вида (а^ )пк=1, 0 < а1 < • • • < ап, 61 < • • • < Ьп, относится к классу осцилляционных матриц.
Эти матрицы плохо обусловлены. Так, для обычной матрицы Вандермонда = k — 1, ai = i/n) числа обусловленности растут как 8n [11]. В работе [12] представлена программа для решения системы линейных алгебраических уравнений с положительно определенной матрицей методом регуляризации.
Литература
1. Тихонов А.Н., Арсенин В. Я. Методы решения некорректных задач. М.: Наука, 1979.
2. Лисковец О. А. Вариационные методы решения неустойчивых задач. Минск: Наука и техника, 1981.
3. Иванов В. К., Васин В. В., Танана В. П. Теория линейных некорректных задач и ее приложения. М.: Наука, 1978.
4. Кабанихин С. И. Обратные и некорректные задачи. Новосибирск: Сибирское научное издательство, 2009.
5. Воеводин В. В., Кузнецов Ю. А. Матрицы и вычисления. М.: Наука, 1984.
6. Михлин С. Г. Вариационные методы в математической физике. М.: Наука, 1970.
7. Михлин С. Г. Численная реализация вариационных методов. М.: Наука, 1966.
8. Гантмахер Ф. Р., Крейн М. Г. Осцилляционные матрицы и ядра и малые колебания механических систем. М.; Л.: ГИТТЛ, 1950.
9. Гантмахер Ф. Р. Теория матриц. М.: Наука, 1967.
10. Higham N. J. Functions of matrices: Theory and computation. Philadelphia: Society for Industrial and Applied Mathematics, 2008.
11. Gautschi W. The Condition of a Matrix Arising in the Numerical Inversion of the Laplace Transform // Mathematics of Computation. 1969. Vol. 23, No. 105. P. 109-118.
12. Бурова И. Г., Рябов В.М., Кальницкая М.А., Малевич А. В., Демьянович Ю. К. Программа для решения системы линейных алгебраических уравнений с положительно определенной матрицей методом регуляризации. Патент № 2018661356, 6.09.2018.
Статья поступила в редакцию 20 апреля 2019 г.;
после доработки 2 июня 2019 г.; рекомендована в печать 13 июня 2019 г.
Контактная информация:
Лебедева Анастасия Владимировна — канд. физ.-мат. наук, доц.; a.v.lebedeva@spbu.ru Рябов Виктор Михайлович — д-р физ.-мат. наук, проф.; v.ryabov@spbu.ru
On the numerical solution of system of linear algebraic equations with ill-conditioned matrices
A. V. Lebedeva, V. M. Ryabov
St. Petersburg State University, Universitetskaya nab., 7-9, St. Petersburg, 199034, Russian Federation
For citation: Lebedeva A. V., Ryabov V. M. On the numerical solution of system of linear algebraic equations with ill-conditioned matrices. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2019, vol. 6(64), issue 4, pp. 619-626. https://doi.org/10.21638/11701/spbu01.2019.407 (In Russian)
The system of linear algebraic equations (SLAE) is considered. If the matrix of the system is non-degenerate, then there is a unique solution to the system. In a degenerate case, the system may not have a solution or have infinitely many solutions. In this case, the concept of a normal solution is introduced. The case of a non-degenerate square matrix can theoretically be considered good in terms of existence and uniqueness of the solution, but in the theory of computational methods, nondegenerate matrices are divided into two categories: "ill-conditioned" and "well-conditioned". Badly called matrices for which the solution
of the system of equations is practically unstable. One of the important characteristics of practical solution stability A system of linear equations is a condition number. Usually, regularization methods are used to obtain a reliable solution. A common strategy is to use Tikhonov's stabilizer or his modifications, or the representation of the desired solution in the form of orthogonal sums of two vectors, one of which is determined stably, and for searching the second requires some stabilization procedure. In this article the methods of numerical solution of SLAE are considered with a positive defined symmetric matrix or oscillating matrix type using regularization, leading to SLAEs with a reduced conditionality number.
Keywords : a system of linear algebraic equations, ill-posed problems, ill-conditioned problems, condition number, regularization method.
References
1. Tikhonov A.N., Arsenin V. Ya., Solutions of Ill-Posed Problems (Winston, 1977)
2. Liskovets O.A., Variational methods for solving unstable problems (Nauka i Tehnika Publ., Minsk, 1981). (In Russian)
3. Ivanov V. K., Vasin V. V., Tanana V. P., Theory of linear ill-posed problems and its applications (Nauka Publ., Moscow, 1978). (In Russian)
4. Kabanikhin S.I., Obratnye i nekorrektnye zadachi (Sibirskoe nauchnoe izdatelstvo, Novosibirsk, 2009). (In Russian)
5. Voevodin V.V., Kuznetsov Yu.A., Matrices and computations (Nauka Publ., Moscow, 1984). (In Russian)
6. Mikhlin S. G., Variational methods in mathematical physics (Pergamon Press, 1982).
7. Mikhlin S. G., The numerical performance of variational methods (Wolters-Noordhoff Publishing, Groningen, 1971).
8. Gantmakher F. R., Krein M.G., Oscillation matrices and kernels and .small vibrations of mechanical systems (AMS Chelsea Publ., Providence, R. I., 2002).
9. Gantmakher F. R., The Theory of Matrices (Chelsea Publishing Company, 1989).
10. Higham N. J., Functions of matrices: Theory and computation (Society for Industrial and Applied Mathematics, Philadelphia, 2008).
11. Gautschi W., "The Condition of a Matrix Arising in the Numerical Inversion of the Laplace Transform", Mathematics of Computation 23(105), 109-118 (1969).
12. Burova I.G., Ryabov V. M., Kalnitskaya M.A., Malevich A. V., Demjanovich Yu.K., Program for solving a system of linear algebraic equations with a positively defined matrix by the regularization method (Patent No. 2018661356, September 6, 2018). (In Russian)
Received: April 20, 2019 Revised: June 2, 2019 Accepted: June 13, 2019
Authors' information:
Anastasia V. Lebedeva — a.v.lebedeva@spbu.ru Victor M. Ryabov — v.ryabov@spbu.ru
