Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 17
Текст из файла (страница 17)
В приводимых ниже алгоритмах явно не указано, когда элементыматриц D, L̄ и L должны замещать соответствующие элементы исходнойматрицы P . Такие замещения могут происходить сразу, а могут откладываться до того момента, когда элементы матрицы P станут ненужными длядальнейших вычислений. В этом отношении не все ijk-формы одинаковоэкономичны в реализации, и для каждой из них вопрос о скорейшем замещении исходных элементов матрицы P нужно решать отдельно.Два алгоритма Холесского: разложения LLT и L̄DL̄Tс немедленными модификациями1) kij-алгоритм1/2(l11 = p11 )Для k = 1 до n − 1Для i = k + 1 до nlik = pik /pkk(lik = pik /lkk )Для j = k + 1 до ipij = pij − lik pjk(pij = pij − lik ljk )1/2(lk+1,k+1 = pk+1,k+1)2) kji-алгоритм1/2(l11 = p11 )Для k = 1 до n − 1Для s = k + 1 до nlsk = psk /pkk (psk /lkk )Для j = k + 1 до nДля i = j до npij = pij − lik pjk(pij = pij − lik ljk )1/2(lk+1,k+1 = pk+1,k+1)Четыре алгоритма Холесского: разложения LLT и L̄DL̄Tс отложенными модификациями3) jki-алгоритм1/2(l11 = p11 )Для j = 2 до nДля s = j до nls,j−1 = ps,j−1/pj−1,j−1(ls,j−1 = ps,j−1/lj−1,j−1)Для k = 1 до j − 1Для i = j до npij = pij − lik pjk(pij = pij − lik ljk )1/2(lj,j = pj,j )4) jik-алгоритм1/2(l11 = p11 )Для j = 2 до nДля s = j до nls,j−1 = ps,j−1/pj−1,j−1(ls,j−1 = ps,j−1/lj−1,j−1)Для i = j до nДля k = 1 до j − 1pij = pij − lik pjk(pij = pij − lik ljk )1/2(lj,j = pj,j )976 Разложения Холесского5) ikj-алгоритм1/2(l11 = p11 )Для i = 2 до nДля k = 1 до i − 1li,k = pi,k /pk,k(li,k = pi,k /lk,k )Для j = k + 1 до ipij = pij − lik pjk(pij = pij − lik ljk )1/2(li,i = pi,i )6) ijk-алгоритм1/2(l11 = p11 )Для i = 2 до nДля j = 2 до ili,j−1 = pi,j−1/pj−1,j−1(li,j−1 = pi,j−1/lj−1,j−1)Для k = 1 до j − 1pij = pij − lik pjk(pij = pij − lik ljk )1/2(li,i = pi,i )Замечание 6.6.
Приведенные алгоритмы L̄DL̄T и LLT -разложенийХолесского матрицы получены из соответствующих ijk-алгоритмов L̄U разложения матрицы A (см. подразд. 3.5). Для получения Ū DŪ T и U U T разложений Холесского матрицы удобно исходить из Ū L-разложения матрицы A, если для него предварительно построить ijk-алгоритмы. Этопостроение нетрудно выполнить, если учесть, что Ū L-разложение соответствует измененному (инверсному) порядку исключения переменных.
В этомслучае модификация системы уравнений начинается с последней переменнойпоследнего уравнения.Суммируя вышеизложенное по ijk-формам алгоритмов Холесского, полученных из ijk-форм алгоритмов Гаусса, имеем 24 разновидности разложенийсимметричной положительно определенной матрицы P :6 ijk-форм для P = L̄DL̄T ,6 ijk-форм для P = LLT ,6 ijk-форм для P = Ū DŪ T ,6 ijk-форм для P = U U T .6.5Разложение Холесского: алгоритмы окаймленияКак и для LU -разложения, для разложения Холесского в любых еговариантах (6.2) существует еще один класс алгоритмов, — так называемыематрично-векторные алгоритмы, объединяемые идеей окаймления.
Получение этих алгоритмов базируется на блочном перемножения матриц, участвующих в разложении. Здесь полностью применимы принципы, изложенныев подразд. 6.2.986.5 Разложение Холесского: алгоритмы окаймленияПокажем, как выводятся такие матрично-векторные алгоритмы на примере одного из четырех вариантов разложения Холесского (6.2), а именно,варианта P = LLT .
Пользуясь этим справочным материалом, любой студентсможет самостоятельно построить родственные алгоритмы для других трехвариантов разложения. Для этого поделим все матрицы в данном вариантена блоки, выделяя в каждой матрице j-ю строку и j-й столбец. Тем самымразложение P = LLT будет представлено в блочной форме P11j ⇒ aTP31j⇓apjjbj⇓P13 L11T = TcljjbL31 dP33 L33T L11 j⇓cljjLT31 T ,dLT33(6.3)где фрагменты j-й строки и j-го столбца обозначены как векторы-столбцывыделенными символами a, b, c и d, а заглавные буквы обозначают матрицы. Нулевые элементы треугольных матриц не показаны.Перемножение матриц (6.3), выполняемое поблочно, дает девять соотношений относительно блок-элементов матриц P и L.
Пользуясь этим, рассмотрим два основных способа разложения матрицы P методом окаймления.Окаймление известной части разложения. Из указанных девятисоотношений возьмем только те, что окаймляют блок P11 = L11LT11, считая,что в этой части разложение уже сделано, т. е. что блок L11 уже вычислен.В силу симметрии P из трех окаймляющих произведений имеем только два:2a = L11c и pjj = cT c + ljj.(6.4)Отсюда сначала находим c как решение нижнетреугольной системы уравнений L11c = a; затем находим ljj = (pjj − cT c)1/2.Окаймление неизвестной части разложения. Из указанных соотношений возьмем те, что окаймляют блок P33 , считая, что до этого блокаразложение уже сделано, т.
е. что блоки L11, L31 и c уже найдены. В силусимметрии P из трех окаймляющих произведений имеем только два:2pjj = cT c + ljjи b = L31c + dljj .(6.5)Отсюда сначала находим ljj = (pjj − cT c)1/2; затем d = (b − L31c)/ljj .Существует два естественных способа реализации окаймления известнойчасти в LLT -разложении.996 Разложения ХолесскогоВ первом варианте треугольная система в (6.4) решается с помощьюстрочного алгоритма (аналог алгоритма на рис.
4.1 слева), во втором — спомощью алгоритма скалярных произведений (аналог алгоритма на рис. 4.1справа). Псевдокоды этих двух вариантов приведены на рис. 6.1.√√l11 = p11l11 = p11Для j = 2 до nДля j = 2 до nДля k = 1 до j − 1Для i = 2 до jljk = pjk /lkklj,i−1 = pj,i−1/li−1,i−1Для i = k + 1 до jДля k = 1 до i − 1pji = pji − ljk likpji = pji − ljk lik√√ljj = pjjljj = pjjРис. 6.1.
Алгоритмы окаймления известной части LLT -разложения: строчный (слева); алгоритм скалярных произведений (справа)Для окаймления неизвестной части в LLT -разложении также существуютдва естественных способа реализации выражений (6.5). Здесь основной операцией является умножение вектора на прямоугольную матрицу.Можно реализовать такие умножения посредством скалярных произведений или линейных комбинаций, что приводит к двум различным формам алгоритма, показанным на рис. 6.2, которые аналогичны алгоритмамДонгарры–Айзенштата на рис.
4.4.Для j = 1 до nДля k = 1 до j − 1pjj = pjj − ljk ljk√ljj = pjjДля k = 1 до j − 1Для i = j + 1 до npij = pij − lik ljkДля s = j + 1 до nlsj = psj /ljjДля j = 1 до nДля i = j + 1 до nДля k = 1 до j − 1pij = pij − lik ljkДля k = 1 до j − 1pjj = pjj − ljk ljk√ljj = pjjДля s = j + 1 до nlsj = psj /ljjРис. 6.2. Алгоритмы окаймления неизвестной части LLT -разложения: алгоритм линейных комбинаций (слева); алгоритм скалярных произведений (справа)Таким образом, выше показано, что алгоритмы окаймления в LU разложении (см. разд.
4) легко модифицируются для случая симметриче1006.6 Особенности хранения ПО-матрицы Pской положительно определенной матрицы P . Тогда мы имеем 4 вариантаразложения Холесского (6.2), 2 способа окаймления и 2 схемы вычисленийдля каждого алгоритма окаймления. Всего получается 16 вариантов алгоритмов окаймления для разложения Холесского симметрической положительно определенной матрицы. Добавляя к ним 24 разновидности ijk-форм,получаем 40 различных вычислительных схем разложений Холесского, которые и составляют весь набор вариантов (см.
подразд. 6.8) задания на лабораторный проект № 5 (см. подразд. 6.7).6.6Особенности хранения ПО-матрицы PКак уже отмечалось (см. стр. 96), особенностью данного проекта являетсяиспользование линейных (одномерных) массивов для хранения матрицы P .Так как матрица P симметрическая, то достаточно хранить только нижнюю (или верхнюю) треугольную часть этой матрицы вместе с диагональю.Причем для хранения заполненной матрицы используется один одномерныймассив, а для хранения разреженной — два.Хранение матрицы P может быть организовано по столбцам или по строкам в зависимости от используемого алгоритма разложения.Рассмотрим строчный вариант хранения нижней треугольной частизаполненной матрицы P .
В этом случае все элементы нижней треугольной матрицы записываются построчно в одномерный массив. Так как дляхранения первой строки матрицы требуется один элемент массива, для хранения второй строки — два элемента и т. д., то для хранения симметрической матрицы размера n требуется одномерный массив размера n(n + 1)/2.Положение (i, j)-го элемента матрицы P в массиве определяется по формулеk = (i − 1)i/2 + j .Аналогичным образом организуется хранение матрицы P по столбцам.Как уже отмечалось, для хранения разреженной матрицы P используются два одномерных массива. Так, хранение нижней треугольной частиматрицы P по строкам можно организовать следующим образом.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
