Богачев - Практикум на ЭВМ. Методы решения линейных систем и нахождения собственных значений (947493), страница 22
Текст из файла (страница 22)
По доказанному выше это означает, что все матрицы Аь почти треугольные (трехдиагональные). Методы нахождения собственных значений К.Ю.Богачев Замечание 3. Описанные ниже приемы не только ускоряют сходимость алгоритмов нахождения собственных значений, но и расширяют множество матриц, для которых они сходятся. Другими словами, при использовании этих приемов алгоритмы нахождения собственных значений часто работают для матриц, для которых не выполнены условия приведенных теорем о сходимости алгоритмов.
~7. БЯ АЛГОРИТМ 105 Замечание 4. Без использования описанных ниже приемов скорость сходи- мости алгоритмов нахождения собственных значений оказывается весьма низкой, а множество матриц, для которых они применимы, достаточно узким. Поэтому эти приемы всегда применяются при вычислениях по этим алгоритмам. ~ 7.3.1. Исчерпываиие матрицы На й-ом шаге матрица Ай имеет вид (й) а„. (1) а2; (й) ат,п-2 (й) (й) ат и, (й) ае„, (й) а,п (й) а2п (й) ат; (й) а1-~-1,1 (й) а;и 2 (й) а;и, (й) а;и а ап2пе а( ) ап — 1,п — 2 (й) 2 (й) ап — 1,п (й)' апп Если для некоторого г, 1 = 1,2,..., и — 1 выполнено условие ~ат „,;~ ( е~)Ай~( где е точность, с которой требуется найти собственные значения, то полагаем а;ттт = О.
Поскольку собственные значения являются непрерывными функциями элементов матрицы, то собственные значения матрицы Ай изменятся при этой замене на величину порядка е. Новая матрица Ай имеет блочную структуру: (й) (й) а,; а„, (й) (й) а2 1 а2 пй1 (й) ат — 2 (й) а2, — 2 (й) (й) а,„, ап (й) (й) а2, — 1 (й) а — 2 (й) ат-)-1, — 2 (й) атп, (й) а;„,„, (й) а;и (й) а;,и а( ) и — 2,п — 2 а( ап — 1,п — 2 а( ) ап — 2,п — 1 а( ап — 1,п — 1 (й)' айп, и ее множество собственных значений равно объединению множеств собственных значений подматриц (а 1 ),1 — 12; Е Мт и (а 1 ),,1 — 1+1,1.1.2 „Е Мп 1.
Следо(й) (й) вательно, можно применить алгоритм к каждой из этих подматриц по отдельности. Методы нахождения собственных значений К.Ю.Богачев (й) (й) а„а12 (й) (й) а2, а22 (й) а22 (й) (й) ап агй (й) (й) аи а22 (й) а32 (й) (й) а;; а1 111 (й) О а;„,; , (й) а1-~-2,1-1-1 а( ап — 2,п — 1 а( ) ап — т,п — 1 (й)' ап,п — 1 (й) ап 2„ (й) ' ап, „ (тс) апп ~7. ХВ АЛГОРИТМ 106 1асто для уменыпения вычислительных затрат вместо условия (а;~,;( < (ь) е//Аь!/„используют условие !а) ),;/ < я/!Ад//, ~ 7.3.2. Сдвиги Для матрицы А — ЛХ собственные значения равны Л; — Л, Л; — собственные значения матрицы А.
Скорость сходимости поддиагональных элементов а; „,; к ну(й) (Л... — Л''1 лю для этой матрицы по теореме о сходимости алгоритма есть О Если Л близко к Л;+~, то скорость сходимости элемента а; „- к нулю будет очень (ь) высокой. Модифицированный Х Л-алгоритм, основанный на этой идее, выглядит следующим образом. Будем строить для матрицы А Е М„последовательность (Аь) матриц Аь Е М„по следующим правилам: 1) А~ — — А; 2) для всех й = 1,2,... матрица А~+~ получается из матрицы Аь следующим образом: а) определяем требуемый сдвиг зь (его оптимальный выбор — отдельная задача), б) строим ХЛ-разложение матрицы Аь — зь1: Аь — зь1 = ХьХХь, в) вычисляем матрицу Аьч.~ как произведение матриц В~ и Хч, плюс зь1: А~+~ — — ЯьХь + зь1.
Здесь мы предполагаем, что для каждого й = 1,2,... ХВ-разложение матрицы Аь — зь1 существует, т.е. для нее выполнены условия теоремы 1. Если это не так, то надо изменить значение сдвига зь. На практике условия теоремы 1 не проверяют, а выполняют алгоритм построения ХВ-разложения. Если в алгоритме требуется осуществить деление на О, то немного изменяют значение сдвига зь и заново выполняют алгоритм построения ХВ-разложения.
Нетрудно проверить, что матрица Аь+~ подобна Аь. Аь+~ — — В~Х,ь+ зь1 = (Х~'Х~)(Л~Х~+ зь1) = Х~'(1~В~)Х~+ зьХь| = Х~'(1~Я+ зь1)Хь — — Х~'А~1~ и, следовательно, все матрицы Аь, й = 1, 2,... имеют те же собственные значения, что и матрица А. 3 7.3.3. Практическая организация вычислений в 1Гг алгоритме Пусть требуется определить все собственные значения матрицы А е М„с точностью е. Вначале приводим матрицу к почти треугольному виду А~ унитарным подобием одним из алгоритмов, описанных в 3 1.14 и 3 1.15. К.Ю.Богачев Методы нахождения собственных значений '38.
МЕТОД ХОЛЕЦКОГО 107 Затем к матрице А1 применяем ХВ-алгоритм со сдвигами. На шаге й в качестве сдвига эь возьмем а„„, т.е. эь — — а„„. Поскольку а„„-+ Л„, то гь (ь) , (ь) (ь) является приближением к Л„и скорость сходимости к нулю элемента а„„, будет очень высокой. Как только на некотором шаге й будет выполнено условие (а„„„) ( е)(А(), в качестве Л„берем а„„и применяем алгоритм к подматрице (ь) (ь) (а; ),д — цг „,,„1 Е М„| на 1 меньшей размерности.
Так поступаем до тех пор, пока размерность матрицы не станет равной 2. Для этой матрицы собственные значения определяются как решения соответствующего квадратного уравнения. 3 8. МЕТОД ХОЛЕЦКОГО Метод Холецкого позволяет находить все собственные значения самосопряженной положительно определенной матрицы А Е М„за вдвое меныпее количество арифметических операций, чем Х)т-алгоритм.
3 8.1. Разложение Холецкого, используемое в методе Холецкого Теорема 1. (О разложении Холецкого) Если матрица А = А* > 0 — самосопряженнол положительно определенная, то она допускает представление А = ХХ*, где матрица Х вЂ” нижняя треугольнол с вещественными положительными элементами на главной диагонали.
3 8.1.1. Алгоритм построения разложения Холецкого для произвольной самосопряженной матрицы Алгоритм построения разложения Холецкого самосопряженной матрицы А Е М„очень похож на алгоритм построения разложения Холецкого, используемого в методе Холецкого решения линейных систем с самосопряженной матрицей (см. стр. 36). Пусть для самосопряженной матрицы А = (а; ) (А* = А, т.е. а; = а,;) требуется найти нижнюю треугольную матрицу Х = (1; ) с вещественными положительными элементами на главной диагонали (1а > 0 для всех 1 = 1,..., п), такую, что А = Х Х*, т.е.
~1;ь1;ь = а;, г=1 1,) =1,...,п. Методы нахождения собственных значений К.Ю.Богачев Доказательство. Согласно теореме 1.10.1 и замечаниям 1.10.2, 1.10.3 матрица А допускает представление в виде А = В*В, где Л вЂ” верхняя треугольная матрица с вещественными положительными элементами на главной диагонали. Положим Х = Л* е ЬТ(п). Тогда А = ХХ." и матрица Х --- нижняя треугольная с вещественными положительными элементами на главной диагонали. Теорема доказана. Р. МКТОД ХОЛКЦКОГО Поскольку матрицы А и АА* — самосопряженные, то уравнение с номером О,г) получается из уравнения с номером (1, Я путем комплексного сопряжения и не дает ничего нового. Поэтому система (1) эквивалентна системе 1;й11,=а;, й=! (2) Таким образом, (2) представляет собой систему из п(п+ 1)/2 уравнений с п(п+ 1)/2 неизвестными 111, г' ) 1 (напомним, Л е 1Т(п) и 11 = 0 при 1' < 1, при этом 1йй > О).
Получим формулы для решения системы (2), которые и составляют алгоритм нахождения разложения Холецкого. Перепишем (2) в виде 11й1,й+ ~~> 11й11й — — а;,, (3) Поскольку матрица Х -- нижняя треугольная, то 1;й — — 0 при 4 < й и вторая из сумм в (3) равна нулю. Следовательно, система (2) эквивалентна следующей ~ 1;й1 й + 1нД; = а;,, й=1 1<1, 1,1=1,...,П, 1й1 й+1,Д; = а,;, й=! Выделим в здесь отдельно случай 1 = 1: 1 — 1 аи — Е ~11й~ й=! 1 — ! а; — 2.' 1,й1,й й=1 Отсюда получаем расчетные формулы: 1 — 1 а.! ~ ~1. ~г й=1 (4) 1 — 1 111 = (агч — х,„, 11й1,й)/1н й=1 Процесс вычислений по этим формулам строится следующим образом: вначале вычисляются неизвестные элементы первого столбца матрицы Л: 1п =,,/а11,, 1, ! — — а, ! /111,,1 = 2,..., и; Методы нахождения собственных значений К.Ю.Богачев (здесь считается, что сумма вида 2 й:, равна нулю, если верхний предел сум- мирования меньше нижнего; это позволяет не рассматривать отдельно случай 1 = 1).
Применим к этим уравнениям комплексное сопряжение и учтем, что а; = ач и Ц! - вещественный элемент: ()8. МЕТОД ХОЛЕЦКОГО 109 4,=г:;,:~~„~, 12 — — (а,г — 12111)/122, 1 = 3,...,п; затем по формулам (4) при г = 3 вычисляются неизвестные элементы третьего столбца матрицы Т и так далее. Замечание 1. Организация хранения матриц А и Ь в памяти.
Поскольку для самосопряженной матрицы а; = а;, то можно вместо всей матрицы А хранить только ее нижний треугольник: а„, 1 > г', 2,1 = 1,...,и. Формулы (4) таковы, что при вычислении элемента 1; используются значения элемента а; и вычисленных Ранее элементов 1 ь, Й ( 1.
Это позволЯет хРанить нижнюю треугольную матрицу Б на месте нижнего треугольника матрицы А: 1,,= а;, г > г, 2,1=1,...,п. Оценка количества арифметических операций в алгоритме построения разложения Холецкого Вычисление элемента 1а по формулам (4) требует одной операции извлечения корня и 1 — 1 мультипликативных и столько же адцитивных операций. При фиксированном г' = 1,...,п вычисление элементов 1зч для всех у = 1+ 1,...,п по формулам (4) требует 1+~,",ч1(2 — 1) = (п — 1)(1 — 1)+1 мультипликативных и ~ " 1,1(2 — 1) = (и — 1)(1 — 1) адцитивных операций. Следовательно, вычисление всех элементов матрицы Л требует и операций извлечения корня, 2," 1((и— 1)(1 — 1) + (1 — 1) + 1) = пз/6+ О(пг) (п — + оо) мультипликативных и ~," ((ив 2)(г — 1) + (1 — 1)) = пз/6+ О(пг) (и -+ со) аддитивных операций (подробное вычисление см.