CC BY
39
ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION. NATURAL SCIENCE. 2021. No. 2
УДК 519.612
doi 10.18522/1026-2237-2021-2-11-15
РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ С ИСПОЛЬЗОВАНИЕМ ВЫБОРОЧНОЙ РЕГУЛЯРИЗАЦИИ
© 2021 г. В.Н. Лутай1
1Южный федеральный университет, Таганрог, Россия
SOLVING SYSTEMS OF LINEAR ALGEBRAIC EQUATIONS WITH THE USE OF AN SELECTIVE REGULARIZATION
V.N. Lutay1
Southern Federal University, Taganrog, Russia
Лутай Владимир Николаевич - кандидат технических наук, доцент, кафедра математического обеспечения и применения ЭВМ, Инженерно-технологическая академия, Южный федеральный университет, пер. Некрасовский, 44, г. Таганрог, 347922, Россия, e-mail: vnlu-tay@sfedu.ru
Vladimir N. Lutay - Candidate of Technical Sciences, Associate Professor, Department of Mathematical Support and Computer Applications, Academy for Engineering and Technologies, Southern Federal University, Nekrasovsky Lane, 44, Taganrog, 347922, Russia, e-mail: vnlutay@sfedu.ru
Рассматривается решение систем линейных алгебраических уравнений (СЛАУ), матрицы которых могут быть плохо обусловленными или неполного ранга. Обосновывается представление матрицы СЛАУ в виде произведения двух матриц, одна из которых отличается от исходной увеличенными значениями некоторых диагональных элементов. В качестве метода решения используется разложение исходной матрицы на треугольные составляющие по Гауссу или Холецкому с дополнительной операцией, которая заключается в увеличении минимальных (положительно определенные матрицы) или нулевых (сингулярные матрицы) диагональных членов треугольных матриц в процессе разложения. Для положительно определенного случая обсуждается стратегия выбора наименьшего диагонального члена в сочетании с общепринятым при решении СЛАУ выбором главного элемента. В случае сингулярных матриц добавление к 0 некоторого числа позволяет закончить процесс треугольного разложения и представить исходную матрицу в виде произведения двух матриц того же ранга. Показано, что если ранг исходной матрицы меньше размерности на единицу, то для её псевдоинверсии целесообразно использовать алгоритм Гревиля. Приводятся результаты вычислительных экспериментов.
Ключевые слова: плохо обусловленные матрицы, матрицы неполного ранга, треугольное разложение, увеличение диагональных членов, псевдообратная матрица.
The solution of systems of linear algebraic equations, which matrices can be poorly conditioned or singular is considered. As a solution method, the original matrix is decomposed into triangular components by Gauss or Cholesky with an additional operation, which consists in increasing the small or zero diagonal terms of triangular matrices during the decomposition process. In the first case, the scalar products calculated during decomposition are divided into two positive numbers such that the first is greater than the second, and their sum is equal to the original one. In further operations, the first number replaces the scalar product, as a result of which the value of the diagonal term increases, and the second number is stored and used after the decomposition process is completed to correct the result of calculations. This operation increases the diagonal elements of triangular matrices and prevents the appearance of very small numbers in the Gauss method and a negative root expression in the Cholesky method. If the matrix is singular, then the calculated diagonal element is zero, and an arbitrary positive number is added to it. This allows you to complete the decomposition process and calculate the pseudo-inverse matrix using the Greville method. The results of computational experiments are presented.
Keywords: ill-conditioned matrices, matrices of incomplete rank, triangular decomposition, increase of diagonal terms, pseudo-inverse matrix.
ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION. NATURAL SCIENCE. 2021. No. 2
Введение
Широко известным методом повышения устойчивости решения систем линейных алгебраических уравнений (СЛАУ) с плохо обусловленной и особенной матрицей является регуляризация по Тихонову [1], которая предлагает добавлять ко всем диагональным членам симметричной матрицы
А = 1а,Л в системе -'-'пхп
Ах = Ь (1)
положительное число Я, значение которого выбирается по некоторым дополнительным соображениям. Решение СЛАУ
(А + И)х = Ь (2)
с матрицей, которая становится неособенной или лучше обусловленной, чем исходная, является наилучшим псевдорешением.
В основе предлагаемого метода лежит модификация Ьи-разложения матрицы А или ее разложения по Холецкому, которая заключается в увеличении диагональных членов треугольной матрицы, вычисляемых в процессе разложения, на некоторую положительную величину. В [2] показано, что такая модификация треугольного разложения соответствует решению СЛАУ вида (А + Е)х = Ь, где Е - квадратная матрица и-го порядка, у которой ненулевыми являются к < п) диагональных членов, индексы которых совпадают с индексами увеличенных диагональных членов треугольных матриц, причем величина вц соответствует значению, на которое увеличивается а^. Количество вц определяется в процессе разложения. Таким образом, в отличие от (2) увеличиваются не все диагональные коэффициенты А, а значения вц для разных I могут быть различными.
Точное решение (1) может быть получено при решении двух систем
Мх = b, Zx = х,
где M = А + Е , Z = (I - М-1Е). Матрица Z имеет следующий вид:
-У1к У2к
(3)
г1 0
Z =
У11
0 1 -У21
0 0 1-У31
00
Уя1
Узк
1 - Уяк]
(4)
к раз использовать обратную подстановку в соответствующих алгоритмах.
Количество дополнительных операций, необходимых для решения системы с дополнительной операцией, зависит от к и оценивается в [2] при максимально возможном к=и-1, не большим, чем стандартное разложение.
Ниже приведены выражения из [3] для решения системы (1) с положительно определенной матрицей по Холецкому (6) и ЬЦ-разложением (7) (опущены формулы для вычисления недиагональных членов треугольных матриц) ННТ = А,
1
hu = (ßü - hik2)2, t G [2, n], Hz = b, LTx = z, LU = A,
(6)
Un dj,
Zfc =11 'ikUki, i G [2,n],
Столбцы у, могут быть получены при решении систем ЛАУ
Му, = е,1б[1,&], (5)
где е; - столбец матрицы Е, у которого есть ненулевой элемент.
Так как предполагается, что треугольные матрицы уже получены, то для решения (5) достаточно
Ьг = Ь, их = г, (7)
где Н, Ь - нижние треугольные матрицы, и - верхняя, последние строки в (6) и (7) носят название обратной подстановки соответствующего разложения.
Положительно определенные матрицы
В качестве средства увеличения диагональных коэффициентов в [2] предложена операция разделения скалярного произведения двух векторов на два числа: 81 и 82, первое из которых содержит старшие разряды скалярного произведения, а второе - младшие. Если количество разрядов мантиссы числа с плавающей точкой равно ^ а количество разрядов в 82 - х, то количество разрядов в 81 равно (1 - х) при одном и том же значении порядка обоих чисел. Число 81 продолжает участвовать в дальнейших операциях разложения, а 82 отсекается от вычислительного процесса и запоминается в матрице Е.
Ьи-разложение. В вычислительной практике общепринятым для Ьи-разложения является выбор максимального ведущего элемента, как правило, по столбцу. В результате элементы матрицы Ь близки к 1, а диагональные элементы матрицы и не возрастают [4]. Это позволяет выбрать отличную от предложенной в [2] стратегию инициализации операции разделения, а именно использование барьера - некоторого положительного числа е такого, что операция разделения выполняется в том случае, когда иц < е. Значение е можно выбрать из разных соображений, в частности, исходя из принятого представления - матрица хорошо обусловлена, если ее число обусловленности не превышает 103. Тогда е можно принять равным 10-3 или 10-4.
ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION. NATURAL SCIENCE. 2021. No. 2
В качестве примера выбрана верхняя матрица Хессенберга вида [5]
А =
п п — 1 п -2 3 2 1
п — 1 п — 1 п — 2 3 2 1
0 п — 2 п — 2 3 2 1
0 0 0 0 2 2 1
0 0 0 0 0 1 1-
для n=8 с числом обусловленности 2,75 • 105. Максимальный элемент матрицы L стандартного разложения без выбора ведущего элемента равен 6,85, с выбором - 1. Минимальное иц с выбором -2,710-5.
Результаты приведены в табл. 1 (все элементы матрицы были разделены на 8, значение т в операции разделения было принято равным 15 при t=17).
Таблица 1
Результаты применения операции разделения для LU-разложения / Results of applying the split operation for LU decomposition
8 k min ии cond (M) cond (Z)
- - 2,7-10-5 - -
0,001 1 0,005 1,2103 5,2-102
Как следует из табл. 1, решение плохо обусловленной СЛАУ сводится к получению матрицы U с существенно меньшим числом обусловленности, чем при стандартном разложении, вычислению обратными подстановками вектора х из (3) и вектора из (5) и решению системы c треугольной в данном случае матрицей (4).
Разложение Холецкого. Для проверки эффективности разложения Холецкого использовалась матрица Гильберта, известная своей плохой обусловленностью [6], размерностью 8, 9 и 10. Стандартное разложение Холецкого прекращается во всех трех случаях при вычислении h88. Результаты применения операции разделения приведены в табл. 2.
Таблица 2
Результаты применения операции разделения для разложения Холецкого / Results of applying the separation operation for the Cholesky decomposition
n k cond (M) cond (Z) min hn
8 1 1,1 ■ 109 1,9-106 2,6-10-2
9 2 1,1 ■ 109 1,7-107 2,9-10-2
10 3 1,2 109 9,2-07 1,210 -3
10 4 2,0-107 4,1 ■ 107 5,0-10-2
Данные первых трех строк табл. 2 получены при условии, что разделение чисел выполняется, когда подкоренное становится меньшим 0. Числа в четвертой строке получены при использовании барьера со значением 0,000001, что соответствует 0,001
для . Обусловленность матрицы М при этом уменьшается, но возрастает количество дополнительных операций.
Матрицы неполного ранга
Покажем, что (3) можно использовать для решения (1) и в том случае, когда квадратная матрица А размерностью п является сингулярной (ее ранг г < п). В этом случае наилучшее решение по методу наименьших квадратов выглядит следующим образом [7]: х = А+Ь, где А+ - псевдообратная матрица. Наиболее известным методом вычисления таких матриц является метод сингулярного разложения (SVD). Он является достаточно сложным в реализации: оценка количества операций - 12п3 [4]. Ниже рассматривается возможность вычисления псевдообратной матрицы предлагаемым в статье методом на примере матрицы с рангом п—1.
Так как А в нашем случае представляется в виде произведения двух матриц M и Z, то для получения А+ матрица М должна иметь г независимых столбцов, а Z — г независимых строк [8]. Положим, что матрица А имеет одну линейно зависимую строку. Тогда в процессе стандартного LU-разложения одна из строк треугольной матрицы U оказывается равной 0. Переставим ее вниз, прибавим к ее диагональному члену некоторое положительное число ипп и завершим треугольное разложение. При этом элементы последней строки матрицы L равны 0, за исключением 1пп, который равен 1. В А сделаем соответствующую перестановку строк и последнюю строку сделаем нулевой, так что М отличается от А только членом тпп = ипп.
При к = 1 в (5) еп = тпп, уп1 равно единице и последняя строка в матрице Z равна нулю. При этом все ненулевые строки Z линейно независимы, последний столбец является линейной комбинацией остальных и ранг Z по строкам и по столбцам равен п—1. Так как последние строки в А и Z нулевые, то без внесения искажений в произведение MZ можно сделать нулевым последний столбец, а значит, и последнюю строку матрицы М. Получившуюся таким образом матрицу обозначим W. Так как матрицы W и Z имеют ранг г, то можно записать следующее выражение: А = ШХ.
Тогда А+ = 2+Ш+.
При инвертировании положительно определенной матрицы M значение элемента шпп влияет только на тот столбец обратной матрицы, в котором он расположен. Если этот столбец сделан нулевым, то главная подматрица матрицы М-1 не изменится и будет равна [9]. Это значит, что вычисление Ш+сводится к обратной подстановке (7).
ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION. NATURAL SCIENCE. 2021. No. 2
Вычисление псевдообратной матрицы 7+ удобно выполнить алгоритмом Гревиля: к единичной квадратной матрице размерностью (п — 1) добавляется один столбец, который является линейной комбинацией остальных. Как следует из формул, приведенных в [8, с. 52, 53], вычисления в этом случае сводятся к умножению столбца на строку, определению нормы вектора и делению элементов вектора на положительное число.
Приведем пример работы алгоритма. В качестве матрицы неполного ранга выберем матрицу 4*4, последняя строка которой является линейно зависимой. После Ьи-разложения, обнаружения соответствующей строки в треугольной матрице, ее обнуления в А и добавления к ее диагональному члену единицы получим
А =
2 1 1 3 2 1 1 3
1 0 1 -1 M = 1 0 1 -1
0 1 2 3 0 1 2 3
0 0 0 0 0 0 0 1
1 0 0 -1/3
0 1 0 13/3
0 0 1 -2/3
0 0 0 0
Z =
Далее вычисляем W+ и Z+ :
W =
2110 1010 0120 0000
M
-1 _
1/3 1/3 -1/3 1/3 2/3 -4/3 1/3 -4/3 1/3 2/3 1/3 2/3
0
0
0
1
Z+ =
А+
Z+W+ =
1/3 1/3 -1/3 0 2/3 -4/3 1/3 0 -1/3 2/3 1/3 0 - 0 0 0 0 0,994 0,071 -0,011 0 0,071 0,076 0,142 0 -0,011 0,142 0,978 0 -0,016 0,213
0,3825 0,2295
0,0273 0,0164 -0,235 0,459
0,1475 -0,3115
0,033 0
0,3115 0
0,0492 0
0,377 0
0,0656 0
Матрица А+ совпадает с матрицей, полученной из исходной посредством алгоритма SVD (функция pinv из библиотеки linalg Python).
Для получения решения СЛАУ можно не вычислять полностью М-1, а воспользоваться решением следующих двух систем:
х = г+х, шх = ь.
Треугольное разложение Ш совпадает с разложением М, если в М сделать нулевым последний столбец матрицы и. Это позволяет воспользоваться обратной подстановкой, приняв значение хп равным произвольному числу, в том числе и 0: последний столбец нулевой и значение хп не влияет на остальные элементы этого вектора. При оценке трудоемкости операций можно ориентироваться на решение системы с помощью ЬЦ-2 3
разложения, а именно так как сложность использованных операций по Гревилю и умножения матрицы на вектор не превышает п2.
При дальнейшем уменьшении ранга матрицы можно использовать приведенные обоснования для представления сингулярной матрицы в виде произведения 2 матриц, но вычисление рекуррентным алгоритмом Гревиля становится слишком дорогим. Поэтому оценка сложности алгоритма при г < п — 1 требует дальнейших исследований.
Заключение
В работе получены следующие результаты:
1. Показано, что увеличение минимальных или нулевых диагональных членов треугольных матриц, получаемых при стандартных треугольных разложениях, приводит к представлению матрицы исходной системы в виде произведения двух матриц, одна из которых отличается от исходной некоторым количеством увеличенных диагональных членов.
2. Для плохо обусловленных матриц в качестве средства увеличения выбрана операция разделения скалярного произведения на два числа; обсуждается условие инициализации этой операции и ее эффективность для симметричных и несимметричных матриц.
3. Если исходная матрица сингулярна и ранг ее на единицу меньше размерности, то для ее представления в виде произведения двух матриц следует в процессе ЬЦ-разложения добавить некоторое положительное число к нулевому диагональному элементу матрицы и. В этом случае сложность решения СЛАУ можно сравнить со сложностью алгоритма Гаусса.
Литература
1. Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.: Наука, Главная редакция физ.-мат. лит-ры, 1979. 285 с.
2. Лутай В.Н. Повышение устойчивости треугольного разложения плохо обусловленных матриц //
ISSN 1026-2237 BULLETIN OF HIGHER EDUCATIONAL INSTITUTIONS. NORTH CAUCASUS REGION. NATURAL SCIENCE. 2021. No. 2
Сиб. журн. вычисл. математики. 2019. Т. 22, № 4. С. 473-481.
3. Воеводин В.В., Кузнецов Ю.А. Матрицы и вычисления. М.: Наука, Главная редакция физ.-мат. литры, 1984. 320 с.
4. Голуб Дж., Ван Лоун Ч. Матричные вычисления: пер. с англ. М.: Мир, 2010. 548 с.
5. Уилкинсон Дж. Алгебраическая проблема собственных значений. М.: Наука, Главная редакция физ.-мат. лит-ры, 1967. 564 с.
6. Форсайт Дж., Моллер К. Численное решение систем линейных алгебраических уравнений. М.: Мир, 1969. 168 с.
7. Гантмахер Ф.Р. Теория матриц. 5-е изд. М.: Физматлит, 2010. 560 с.
8. Альберт А. Регрессия, псевдоинверсия и рекуррентное оценивание : пер. с англ. М.: Наука, Главная редакция физ.-мат. лит-ры, 1977. 224 с.
9. Лоусон Ч., Хенсон Р. Численное решение задач метода наименьших квадратов : пер. с англ. М.: Наука, Главная редакция физ.-мат. лит-ры, 1986. 232 с.
References
1. Tikhonov A. N., Arsenin V. Ya. (1979). Methods for solving ill-posed problems. Moscow, Nauka Publ.,
Main Edition of the Physical and Mathematical Literature, 285 p. (in Russian).
2. Lutay V.N. (2019). Increasing the stability of triangular decomposition of ill-conditioned matrices. Siberian Journal of Numerical Mathematics, vol. 12, No. 4, pp. 388-394.
3. Voevodin V. V., Kuznetsov Yu.A. (1984). Matrixes and computations. Moscow, Nauka Publ., Main Edition of the Physical and Mathematical Literature, 320 p. (in Russian).
4. Golub G., Van Loan C. (2010). Matrix computations. Moscow, Mir Publ., 548 p. (in Russian).
5. Wilkinson J. (1967). The algebraic eigenvalue problem. Moscow, Nauka Publ., Main Edition of the Physical and Mathematical Literature, 564 p. (in Russian).
6. Forsythe G., Moler C. (1969). Computer solution of linear algebraic systems. Moscow, Mir Publ., 168 p. (in Russian).
7. Gantmaher F. R. (2010). Matrix theory. Moscow, Fizmatlit Publ., 560 p. (in Russian).
8. Albert A. (1977). Regression, pseudo-inversion and recurrent estimation. Moscow, Nauka Publ., Main Edition of the Physical and Mathematical Literature, 224 p. (in Russian).
9. Lawson C., Hanson R. (1986). Solving least squares problems. Moscow, Nauka Publ., Main Edition of the Physical and yMathematical Literature, 232 p. (in Russian).
Поступила в редакцию /Received
9 февраля 2021 г /February 9, 2021
