Богачев - Практикум на ЭВМ. Методы решения линейных систем и нахождения собственных значений (947493), страница 8
Текст из файла (страница 8)
Алгоритм построения разложения Холецкого Построим алгоритм нахождения разложения из теоремы 1. Пусть для самосопряженной матрицы А (А* = А) требуется найти верхнюю треугольную матрицу К = (т,") с вещественными положительными элементами на главной диагонали (тн > О для всех г = 1,..., п) и диагональную матрицу Р с равными по модулю единице диагональными элементами (~4;~ = 1 для всех г = 1,...,и) такую, что А = К'РК. а Элемент (И,1) матрицы РК равен (РК)ь = ',„йтит; = гнить-, так как мат=1 трица Р— диагональная; элемент (г, Й) матрицы Л* равен (К*) ь = тм, элемент (1,Я матрицы К*РК равен (ЯРК);, = х,,(К*);ь(РК)ь, — — 1„г~,йььт~,.
Следоь=1 й=1 вательно, равенство А = К'РК дает нам уравнения тмпььть1 = ай, й=1 (4) г,у = 1,...,ть К.Ю.Богачев Точные методы ретения линейных систем Действительно, если матрица А положительно определена и для нее справедливо разложение А = К*РК, то для всякого х Е С", х ф О выражение (Ах, х) = (К*Р Кх, х) = (РКх, Кх) вещественно и положительно. Для всякого у е С" положим х = К 'у.
Тогда выражение (Ру, у) = (Ах, х) вещественно и положительно для всех у е С", у ф О, т.е. матрица Р положительно определена. Поскольку Р = йаи (с~п,..., д„„), то (Ру, у) = Ц', й," у у: = 2 ", й Д у, ~~. По доказанному, эта сумма вещественна и положительна. Выбирая здесь у = еь, й = 1,...,и, где еь — орты стандартного базиса, находим, что все Иьь, й = 1,..., и должны быть вещественными и положительными. Поскольку ~дьь~ = 1, то это означает, что дьд = 1, й = 1,..., п. 'ч1 О.
МЕТОД ХОЛЕЦКОГО (КВАДРАТНОГО КОРНЯ) 37 тй;Йййтй — — а;, й=1 (5) г < 1, г, 1 = 1,..., п. Таким образом, (5) представляет собой систему из п(п+1)~2 уравнений с п(п+ 1)/2 неизвестными т;, г' < 1 (напомним, Л Е ВТ(п) и т; = О при г > )) и п неизвестными Ййй, Й = 1,..., п (при этом тйй > О и Ййй Е 1 — 1, 1) ). Получим формулы для решения системы (5), которые и составляют алгоритм метода Холецкого. Перепишем (5) в виде тй;Йййгй, + ~ , 'гй;Йййтй, = а;,, г < 1', г,,1' = 1,...,п. (б) й=1 й=ьй1 Поскольку матрица  — верхняя треугольная, то тй; = О при Й > г и вторая из сумм в (б) равна нулю.
Следовательно, система (5) эквивалентна следующей ~ „гй1Йййтйэ+тцЙцт; = а;, г < 1, г,1 =1,...,п, (7) й=1 (здесь считается, что сумма вида З„й, равна нулю, если верхний предел сум- мирования меньше нижнего; это позволяет не рассматривать отдельно случай г = 1; также в (7) мы учли, что тц — вещественный элемент). Выделим в сум- ме (7) отдельно случай г' = 1' 1 — 1 т "Й" = а" — ) ~тй ~ Ййй й=1 1 — 1 гцЙцт;, = а;, — 2. гй1Йййтйэ, й=1 1=1,...,п, г<1, г,1=1,...,п, Отсюда получаем расчетные формулы: 1 — 1 Йц =з'К (ац — Е ~~йГЙйй), 1=1,...,п, (8) 1=1,...,п, т;, = (а;, — 2 11ц Йййтйу)((тайц), й=1 г<у, г,у=1,...,п.
Процесс вычислений по этим формулам строится следующим образом: вначале вычисляются неизвестные элементы первых строк матриц Р и В: Йп — — з18п ап„тп — — фа11~, ,г1 — — а11/(т11Йп), 1 = 2,...,п; К.Ю.Богачев Точные методы решения линейных систем Поскольку матрицы А и К*РК вЂ” самосопряженные, то уравнение с номером (~, г) получается из уравнения с номером (г, Я путем комплексного сопряжения и не дает ничего нового. Поэтому система (4) эквивалентна системе ~1 О. МЕТОД ХОЛЕЦКОГО (КВАДРАТНОГО КОРНЯ) потом по формулам (8) при г = 2 вычисляются неизвестные элементы вторых строк матриц Р и В: А. = 'в ( — 3 3'А ), " = ~М вЂ” ~ Ыч2, тй1 (а21 т12 г~11тц)I(т22~~22) р у 8~ р п~ затем по формулам (8) при г = 3 вычисляются неизвестные элементы третьих строк матриц В и В и так далее.
Замечание 4. Организация хранения матриц А, В и П в памяти. Замечание 5. Случай положительной матрицы А. В силу замечания 2 в этом случае разложение матрицы имеет более простой вид А = В*В. Формулы (8) тоже упрощаются: ю — 1 тн = (ан — Е ~тй'~'), й=й 1=1,...,п, (9) г — 1 тд = (а; — ~ ттитй,)~тйь й=1 г(1, з,~ =1,...,п.
Если матрица А еще и вещественная, то в формулах (9) можно убрать знак модуля и комплексного сопряжения: ~ — 1 тн — — (аи — ~ тй2;), й=1 1=1,...,п, 1 =( 1 — Е ттй1)/ й=1 1<у, 1,~=1,...,п. К.Ю.Богачев Точные методы ретения линейных систем Поскольку для самосопряженной матрицы а; = аО, то можно вместо всей матрицы А хранить только ее верхний треугольник: а;, г < ~, г, у = 1,..., и. Формулы (8) таковы, что при вычислении элемента т;, используются значения элемента а;, и вычисленных ранее элементов тй, Й < з. Это позволяет хранить вернюю треугольную матрицу Л на месте верхнего треугольника матрицы А: т„: — а;,, г ( у, 1,у = 1,..., и.
Матрица В обычно хранится в виде отдельного вектора с~ = (Ып,..., д„„)'. "811. МЕТОД ОРТОГОНАЛИЗАПИИ З 10.3. Оценка количества арифметических операций в алгоритме построения разложения Холецкого 8 11. МЕТОД ОРТОГОНАЛИЗАЦИИ Пусть требуется решить линейную систему Ах=Ь с матрицей А = (а; ) е М„и правой частью Ь = (Ьн...,Ь„)'. Пусть х = (хн..., х„)' точное решение этой системы. Введем следующие вектора из С" ~': У = (хн..., х„, 1) Е С™~, (2) К.10.Богачев Точные методы решения линейных систем 1.
Вычисление элемента ди, г' = 1,..., и по формулам (8) требует г — 1 мультипликативных и столько же адцитивных операций. Следовательно, вычисление всех элементов матрицы .0 требует 2„";,(г — 1) = п(п — 1)/2 = О(л~) мультипликативных и столько же аддитивных операций. 2. Вычисление элемента гн по формулам (8) требует одной операции извлечения корня и 1 — 1 мультипликативных и столько же аддитивных операций (умножение на 4; Е 1 — 1, 1) мы за операцию не считаем).
При фиксированном г' = 1,...,п вычисление элементов г;, для всех у = г+ 1,...,п по формулам (8) требует 1+ ~,".;~,(г — 1) = (и — г)(г — 1) + 1 мультипликативных и 2 ";„,(г — 1) = (и — г)(г — 1) аддитивных операций. Следовательно, вычисление всех элементов матрицы В требует и операций извлечения корня, ),",((н — г)(г — 1) + (г — 1) + 1) = пп/6 + 0(п~) (и -+ сс) мультипликативных и ~,",((и — г) (г — 1) + (г — 1)) = пз/6+ 0(п~) (п -+ сс) алдитивных операций (подробное вычисление см. при подсчете количества арифметических операций для алгоритма ИУ-разложения). На вычислительных машинах, имеющих аппаратную поддержку для вычисления трансцендентных функций вещественного аргумента, операция извлечения корня выполняется за время, по порядку равное времени деления двух вещественных чисел.
Поэтому на таких вычислительных системах операции извлечения корня в описанном выше алгоритме займут время, равное времени выполнения 0(п) мультипликативных операций. На вычислительных машинах, не имеющих аппаратной поддержки для вычисления трансцендентных функций вещественного аргумента, операции извлечения корня в описанном выше алгоритме займут время, равное времени выполнения 0(сопз1 и) = 0(п) мультипликативных операций. Таким образом, алгоритм построения разложения Холецкого требует для своего проведения выполнения ва/6+ 0(п2) (и — ~ сс) мультипликативных и столько же аддитивных операций, а в сумме — пз/3+ 0(п~) (и -+ сс) арифметических операций, т.е. асимптотически вдвое меныпе, чем в методе Гаусса или алгоритме построения ИУ-разложения.
~11. МЕТОД ОРТОГОНАЛИЗАПИИ 40 С использованием этих обозначений система (1) может быть записана в виде п условий ортогональности (УА) =О, 1=1,...,п, (4) где (, ) — обычное евклидово скалярное произведение в С"+'. Соотношения (4) дают новое определение понятию решения системы (1): решить систему (1) — это значит найти вектор У вида (2), ортогональный всем векторам (3). Будем строить последовательность подпространств С" и = Е(о) ~ Е(') ~ ~ Е(") (й — Ц~ (й 0 (Ай,е; ) (й 0 е, = е; —, ей, 1= 1+1,...,п.
(Ай, ей ) (5) Проверим, что Е(") состоит из векторов, ортогон альных векторам А„..., Ай. а) Базис пространства ЕОО получается как линейная комбинация базиса пространства Е(" '), все элементы которого ортогональны векторам Ай,..., Ай ~. Следовательно, вектора из Б ортогональны векторам Ай,..., Ай й.
(й) б) Проверим, что все элементы базиса пространства Е(й) ортогональны (й — й) (Ай,ей(" ')) (Ай, е(" ')) = О, 1 = й + 1,..., п. Лемма 1. Для всех й = О,..., п и всех 1 = 1+1,..., п+1 вектор е( имеет (й) не более 1+1 отличных от нуля компонент. Именно, могут быть отличными от нуля компоненты этого вектора с номерами 1,..., Й и 1. Доказательство. При й = 0 это следует из выбора е;, 1 = 1,..., и. Предпо(о) (й 1)' ложим, что утверждение леммы справедливо для е;, ю' = к,..., п.
Покажем, К.Ю.Богачев Точные методы решения линейных систем где Е(й) = (е(„)„,..., е(~+Д вЂ” (и — й+ 1)-мерное подпространство в Со~', состоя(й) (й) щее из векторов, ортогональных векторам Ай,..., Ай (в частности, все вектора одномерного подпространства Е(") ортогональны Ай,..., А„, т.е.
Е(") содержит решение, которое выделяется из всех других векторов этого подпространства тем, что его последняя координата равна 1). В качестве базиса начального пространства Е(о) возьмем стандартный базис С" "': е( = (1, О,..., 0),..., е = (О,..., О, 1).
Для всех й = 1,..., и базис ей „..., е„„каждого следующего подпространства Е строится по предыдущему базису е„,...,е„„, подпространства (й) (й й) (й 1) Е(" ') по формулам ~11. МЕТОД ОРТОГОНАЛИЗАПИИ 41 что оно верно для е;, г = 1+ 1,..., п. Действительно, по предположению у век(й) тора ей только компоненты 1,..., й могут быть отличны от нуля.