Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 16
Текст из файла (страница 16)
Чтобы матрица P была положительно определенной, необходимо и достаточно, чтобы все ее главные миноры былиположительны.Достаточное условие. Диагональное преобладание, т. е. свойство∀i = 1, . . . , n :p(i, i) >nXj=1,j6=i|p(i, j)| ,влечет положительную определенность матрицы P = [p(i, j)].6.2Квадратные корни из P и алгоритмы ХолесскогоОпределение 6.2.Если матрица P может быть представлена какP = SS T(6.1)с квадратной матрицей S, то S называют квадратным корнем из P [97].Квадратные корни матриц, когда они существуют, определяются равенством(6.1) неединственным образом, так как, если S удовлетворяет равенству (6.1)и T есть любая ортогональная (T T T = I и T T T = I) матрица, то ST такжеудовлетворяет (6.1).916 Разложения ХолесскогоЗамечание 6.1. Хотя определение 6.2, т.
е. формула (6.1), частоиспользуется, оно не является универсальным. Обобщения заключаются вотказе от квадратной формы матрицы S или (если она остается квадратной) в допущении комплексных значений ее элементов. В последнем случаепишут P = SS H , где S H обозначает комплексно сопряженную при транспонировании матрицу. Ограничения в определении квадратного корня изматрицы могут заключаться в различных требованиях к ее форме. Например, можно требовать симметричности S = S T или эрмитовости S = S H . Вслучае симметричности, т. е. при S = S T , квадратный корень из матрицыP обозначают S 1/2, т. е. пишут P = SS 1/2.
Если же требовать, чтобы S в(6.1) имела треугольную форму, то в (6.1) приходим к четырем вариантамразложения Холесского [13, 15].Определение 6.3. Разложением Холесского принято называть любоеиз следующих представлений положительно определенной матрицы P :P = LLT ,P = L̄DL̄T ,P = UUT,P = Ū DŪ T ,(6.2)где L — нижняя треугольная матрица с положительными элементами надиагонали, U — верхняя треугольная матрица с положительными элементами на диагонали, L̄ — нижняя треугольная матрица с единичными элементами на диагонали, D — диагональная матрица с положительными элементами на диагонали, Ū — верхняя треугольная матрица с единичнымиэлементами на диагонали.Теорема 6.1 (Нижнее треугольное разложение Холесского).
Симметрическая матрица P > 0 имеет разложение P = LLT , где L — нижняятреугольная матрица. Это разложение на сомножители с положительнымидиагональными элементами в L дается следующим алгоритмом.Для j = 1, . . . , n − 1 рекуррентно выполнять цикл, образованный следующим упорядоченным набором выражений:L(j, j) = P (j, j)1/2,L(k, j) = P (k, j)/L(j, j),k = j + 1, . . . , n,(L(n, n) = P (n, n)1/2.k = j + 1, . . . , n P (i, k) := P (i, k) − L(i, j)L(k, j)i = k, .
. . , nЗамечание 6.2. Приведенная формулировка алгоритма используетобозначения элементов матриц P и L для наглядности. Это не должно создавать впечатления, что данные матрицы хранятся в памяти по отдельности; наоборот, в лабораторном проекте все приведенные вычисления должны926.2 Квадратные корни из P и алгоритмы Холесскогопроводиться в одном и том же массиве. Это значит, что элементы матрицы Pдолжны замещаться элементами матрицы L по мере вычисления последних.Замечание 6.3. Легко видеть, что если L1 и L1 суть два различныхварианта факторизации матрицы P > 0, тоL1 = L2 diag [±1, .
. . , ±1],т. е. в приведенном алгоритме можно брать L(j, j) = ±P (j, j)1/2. Обычнорекомендуют выбирать единственое решение, отвечающее положительнымдиагональным элементам матрицы L.Внимательное прочтение теоремы 6.1 обнаруживает, что данный алгоритм требует вычисления n квадратных корней, что замедляет вычисления.Этот недостаток устраняется, если перейти ко второму варианту разложения Холесского из общего списка (6.2).Следствие 6.1 (Нижнее треугольное разложение Холесского без операции квадратного корня). Симметрическая матрица P > 0 имеет разложение P = L̄DL̄T , где L̄ — нижняя треугольная матрица с единичной диагональю и D = diag (d1, . .
. , dn). Элементы матриц L̄ и D даются следующималгоритмом.Для j = 1, . . . , n − 1 рекуррентно выполнять цикл, образованный следующим упорядоченным набором выражений:dj = P (j, j),L̄(j, j) = 1,(k = j + 1, . . . , n d(n) = P (n, n)),P (i, k) := P (i, k) − L̄(i, j)P (k, j)i = k, . .
. , nL̄(n, n) = 1.L̄(k, j) = P (k, j)/d(j),k = j + 1, . . . , nДанный результат легко получается из теоремы 6.1 с использованием обозначений dj = L(j, j)2 и L̄(i, j) = L(i, j)/L(j, j).Следующий аналог теоремы 6.1 формулирует алгоритм третьей версииразложения Холесского из списка (6.2).Теорема 6.2 (Верхнее треугольное разложение Холесского).
Симметрическая матрица P > 0 имеет разложение P = U U T , где U — верхняятреугольная матрица. Это разложение на сомножители с положительнымидиагональными элементами в U дается следующим алгоритмом.936 Разложения ХолесскогоДля j = n, n − 1, . . . , 2 рекуррентно выполнять цикл, образованный следующим упорядоченным набором выражений:1/2U (j, j) = P (j, j) ,U (k, j) = P (k, j)/U (j, j),k = 1, . . . , j − 1,(U (1, 1) = P (1, 1)1/2.k = 1, .
. . , j − 1 P (i, k) := P (i, k) − U (i, j)U (k, j)i = 1, . . . , kАналогично следствию 6.1, зафиксируем последний из четырех вариантовразложения Холесского (6.2).Следствие 6.2 (Верхнее треугольное разложение Холесского без операции квадратного корня).Симметрическая матрица P > 0 имеет разTложение P = Ū DŪ , где Ū — верхняя треугольная матрица с единичнойдиагональю и D = diag (d1, . . . , dn ).
Элементы матриц Ū и D даются следующим алгоритмом.Для j = n, n − 1, . . . , 2 рекуррентно выполнять цикл, образованный следующим упорядоченным набором выражений:dj = P (j, j),Ū (j, j) = 1, d(1) = P (1, 1),Ū (k, j) = P (k, j)/d(j),k = 1, . . . , j − 1,(Ū (1, 1) = 1.k = 1, . . . , j − 1 P (i, k) := P (i, k) − Ū (i, j)Ū(k, j)dji = 1, . . . , k6.3Программная реализация алгоритмов ХолесскогоВ справочных целях включаем реализацию трех из четырех приведенныхвыше алгоритмов на языке FORTRAN [15]. Эти примеры реализации могутпомочь студентам написать свои собственные программы на других языкахвысокого уровня при выполнении лабораторного проекта № 5, описание которого дано ниже в подразд.
6.7 [97, 15].946.3 Программная реализация алгоритмов ХолесскогоВерхнетреугольное разложение Холесского:P (N, N ), P = U U T , P > 0, U — верхнетреугольная матрица.c j = n, n − 1, . . . , 2DO 5 J = N, 2, −1U (J, J) = SQRT(P (J, J))α = 1./U (J, J)DO 5 K = 1, J − 1U (K, J) = α ∗ P (K, J)β = U (K, J)DO 5 I = 1, K5 P (I, K) = P (I, K) − β ∗ U (I, J)U (1, 1) = SQRT(P (1, 1))Замечание 6.4. Матрица U должна замещать P в компьютерной памяти. Нижние части матриц U и P не должны использоваться вовсе (памятьдля них не выделяется).
В любом случае верхнетреугольная часть матрицыP теряется, так как на ее месте появляется верхнетреугольная часть матрицы U , т. е. все вычисления ведутся в одном и том же массиве P .Верхнетреугольное без√· разложение Холесского:P (N, N ), P = Ū DŪ T , P > 0, Ū — верхнетреугольная матрицас единичной диагональю, D — диагональная матрица.c j = n, n − 1, .
. . , 2DO 5 J = N, 2, −1D(J) = P (J, J)α = 1./D(J)DO 5 K = 1, J − 1β = P (K, J)Ū (K, J) = α ∗ βDO 5 I = 1, K5 P (I, K) = P (I, K) − β ∗ Ū (I, J)D(1) = P (1, 1)Замечание 6.5. В любом случае верхнетреугольная часть матрицыP теряется, так как на ее месте появляется верхнетреугольная часть матрицы Ū , при этом единицы, соответствующие диагонали матрицы Ū , толькоподразумеваются, а на их месте пишутся диагональные элементы матрицыD.
Для поддиагональной части массива P память не выделяется.956 Разложения ХолесскогоПредыдущие замечания 6.4 и 6.5 свидетельствуют, что фактически массив, выделяемый для исходной матрицы P и одновременно для получаемыхна ее месте результатов разложений Холесского, должен быть оформленкак одномерный массив. Размер этого массива, очевидно, равен N (N + 1)/2элементов. Напишем для предыдущего алгоритма его эквивалентную «одномерную» версию.Верхнетреугольное без√· «одномерное» разложение P = Ū DŪ T :Одномерный массив P (N (N + 1)/2) соответствует P = Ū DŪ T .JJ = N (N + 1)/2JJN = JJc j = n, n − 1, . .
. , 2DO 5 J = N, 2, −1α = 1./P (JJ)KK = 0c JJN = следующийJJN = JJ − Jдиагональный элементDO 4 K = 1, J − 1c JJN + K = (K, J)β = P (JJN + K)P (JJN + K) = α ∗ βDO 3 I = 1, Kc KK + I = (I, K)3 P (KK + I) = P (KK + I) − β ∗ P (JJN + I) c KK = K(K − 1)/24 KK = KK + Kc JJ = J(J − 1)/25 JJ = JJN6.4Разложение Холесского: ijk-формыРазложение Холесского симметричной положительно определенной матрицы P может быть получено в результате незначительных измененийбазовых LU -разложений квадратной матрицы A, изложенных в подразд. 3.5.При этом симметрия матрицы P используется для сокращения числа действий примерно вдвое.
Способ хранения матрицы P должен быть компактным, т. е. в одномерном массиве хранится по строкам (или по столбцам)только нижняя (или верхняя) треугольная часть матрицы P вместе с диагональю.В той же последовательности, как выше изложены (см. подразд. 3.5)ijk-формы L̄U -разложения матрицы A, приведем ijk-формы L̄DL̄T и LLT разложений Холесского матрицы P > 0. Из них видно, сколь незначительны966.4 Разложение Холесского: ijk-формытребуемые изменения. В каждом алгоритме объединены оба разложения,при этом те изменения, что относятся к LLT -разложению, заключены вскобки.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
