Прямые методы для разреженных матриц. О. Эстербю, З.Златев (984134), страница 14
Текст из файла (страница 14)
Из (1.7) и (1.8) следует, что (1.7) (1.8) (1.9) 1 Ю А $ гп п, х е=- й"Х'. (1.1О) И Замечание 5.9. Условие (1.8) существенно для методов, обсуждение которых мы начинаем. Если гапЕА <. п, то следует Обра~и~ься к друГим методам, например к сингулярному разложению (см. 144, 45, 751).
И Лемма 5.1О. Если выполигны ославил (1.7), (1.8), то А+ =(АТА) 'Ат. (1.11) И :Замечание 5.11. Если выполняются (1.7) и (1.8), то линейную задачу наи~еньш~х квадратов (1.5) ~ожно перефор- МУлиРОватЬ так: Решить систему Ах ~ Ь вЂ” Г при условии Атг = О. Таким образом„задача (1.5) эквивалентна линейной си- СТЕМЕ ПОРЯДКа П1+ П Б.2. Обилий 1~-ииыовый нр~®ой ®етод Учитывая неизбежные ошибки вычислений, мы заменим линейную задачу наименьших квадратов более слабой. Задача 5.12.
Бай~и ~риб~иже~ие х ~~ И" ' к псевдорешению х =А+Ь. В этом параГрафе будет Определена Общая вычислитель" ная ~~е~а, включаюгцая в себя м~о~и~ так ~азы~ае~ы~ прямые м~тоды для ~и~~йн~х задач н~име~ьших ~~~дра~~в, $,2ь Обилий й-аагоеый арямМ метод Предположим, что задачу 5.12 можно заменить сле.-- дующей, ЗадаЧа 5.13. НаЙТИ ВЕКТОР у Е:- =Р ь ТаКОЙ, Ч7О' у=В» с, (2.1)» где р Ъ~ яь р Е- =Мь я Е="- ""Мь В, =й'~', гагйВ» =ц, В, и с могут быть вычислены по А и Ь, Приближение у» к вектору у можно вычислить посред;- ством следующего двухзтапного процесса: Зтап 1 — Обобщенное разложение.
В~~~с~~~~ В» = л»߻+ Еь 1= 1(1) ~„1~ а= К„(2.7)~ где Риф — матрицы перестановОк, Е; — матрицы вОзмуще: ний (ошибок). Предположим, что В» разложена в произведение В; = С;С;1:»ь 1= 1 (1) К. (2.8) Если 1~."> 1, положим В»+» = С»ТС»Сь 1=1(1) К вЂ” 1. (2.9) Потребуем, чтобы для матриц 0; произведения В~»2;, 1 = = 1(1) К легко вычислялись при любом векторе т„Кроме того, разложение матрицы В», должно быть таким, чтобы легко вычислялись произведения В»+к. Никаких других ограничений на матрицы В», С», Сь 1»» мы не накладываем, за исключением того, что размерности должны быть согласованы так, чтобы все указанные произведения могли быть вычислены. Зтап 2 — Обобщенная подстановка. Вычислить Ы-» К-» Т у = П ЯЭ" Я В'Р' П Р С, с=Нс, (2.10) » ьюь » ьь~ Замечание 5.14.
В формуле (2;1О) и последующих выражениях мы придерживаемся соГлашения ЦА;=А~А~+1 ... Аь (2.1 1) ЗВАЛ ~ ЦА;=1, Ы<~. Н~ М;Р,Е10~ П оЭ~) ~" 1 М-1 Т м~-(П о~0+)оЛ~Р~(П Р,'с~) . ДокйзйтВльстВО. Если 1 = 1, то Р=1 — НВ =1 — Я В+Р В =6~,(1 — В+,Р,ВД,)~, =Я~В~ ЕА1 =МА РР1 Е,Я, =М.~Р,Е~Я» — — Н~. + т -+ т т т т и Поскольку В имеет полный столбцоный ранг, то В+В =1.-Прим. иир86. (2.12) И Если у1 — приближение к вектору у (см. (2.1)), то с помощью у~ и свойства (2.6) можно получить приближение х к псевдорешению х.
Следовательно, нужно доказать, что у1 буде~ Хорошим приближением к у, если матрицы возмущений Еь 1= 1(1)1, ~малыэ в том или ином смысле, Дли матрицы К~БР' риз формулы (2.10) положим Г =1 — НВ,, Г, 1 = К"". (2.13) Если 1~ - 2, то„используя (2.13), (2.1О), (2,16), (2.11), (2.7), (2.8), (2,9), (2,15) и (2.12), получим!1 Е = 1 — М1В! — — 1 — М2С! Р1В! = 1 — М~С! (В! — Е1) Я~ = 1 — М2С! С1СЛ1Ф + М2С! Р!Р! Е1Я! т — т т т т =! — марв~ (П яэ!) + н~.
Это было началом доказательства по индукции, а индуктивный переход записывается следующей выкладкой (2 =1:. :: 1~ — 1): м~В; (.П яю~') = м)+~с~"Р|В~ (П яэ~ ) т = М!+1С» (В1 — Е») Ф П Я10 1 1 =м»~Вн-~ ПяВ ) — н~. 1 ° ! Для окончания доказательства воспользуемся (2.12), (2.7), (2.16), (2.15) и предположением, что В1, и В! имеют полный стОл бцо вы й р а н Г: Р=! — И~В~(П яэ~) + К н~ = ! — (П ЯВ~ ) ЯкВ»РьВ~(П ЯЭ~~) + К н! =1 — П ЯР! ЯЖ(В!! — Ек) Я~к Ц ©В! +,Г, Н! ! аа! 1 ° 1 1=1 К-1 = н + ~„' н, = ~„" н,. (2,1 Следствие 5.16. Если разлОжение ВыпОлняется без Ошибок, т. е.
Е! — — О, 1=1(1)1, то Н! =О, 1=1(1)К и Н=в!+. Если к тому же подстановка выполняется без ощибок ОкруГлениЙ, то у1 = у. В Определение 5.17. Вычислительная схема, состоящая из этапа 1 и этапа 2, называется общим М-шаговым прямым ме*одом, нли М-шаговой вычисли~е~ьн~Й схе~~й для решения задачи 5.13. И Теперь мы пр~~едем ~ес~ь примеров х~рошо известных и часто используемых прямых методов, которые могут рассматриваться как специальные случаи общей к-шаговой вычислительной схемы. Большинство из них являются одношаговыми методами, а для к = 1 общая схема сводится к соотношениям: Следовательно, нужно указать С~, С1, В~ и проверить, что В1 У леГко вычисляется для любОГО к.
Пример 5.18. Если 1п = и, тО классическое ГауссовО исключение получается из общей схемы, если пОложить М =1 и С~ =1, 01 =14 При этом (3.4а) есть преобразование'" задачи 5,12 в задачу 5.13; (3.46) указывает метод, а (3,4в) — связь между у~ и х. Матрицы 1.а и 11а суть треугольные множители А, вычислен- ные посредством гауссова исключения. Пример 5.19. Пусть Гп =. и. Предположим, что система нормальных уравнений ~> решается каким-либо симметричным вариантом гауссова исключения. Эта схема получается, если положить Е = 1 и 5.8.
Сйе~ийльйа~е ~~фчйи 06~1еао ~негода В, = Р,ВЯ~+ Е„ В, =С1С~Э„ у=НС =Я~В~ Р~с. + В1=А А, с=А Ь, у=х, С1 = 1.с, С~ = Вс, 01 = 1.с, т х= у~. (3.5а) (3.56) (3.5в) Обозначим через О( ~~ (у2 ~~~~ ~ ° е ~ ~(уп ~~>' О сингулярные числа матрицы А (т. е. квадратные корни из собственных значениЙ АТА). Спектральное число ОбусловленнОсти матрицы А равно ~((А) =11А ~И~А+ $= о~/(т,. (3.7) Легко видеть, что х (В1) = (х (А))~.
Это одна из причин, почему в общем случае нормальные уравнения нельзя рекомендовать как метод решения линейных задач наименьших квадратов. Кроме того, если матрица А — большая, разреженная и не слишком хорошо обусловленная, то (3.8) может сильно ограничить выбор барьера„ а АтА может уже не быть очень разреженной 16„71. Пример 5.2О. Пусть т > и, Предположим, что расширенная линейная система решается посредством гауссова исключения. Этот метод получается при ~( = 1 и (сравни (1.12) ) С1 = 1 „С~ = 1, В» — — 13,; х = последние и КООрдинат уь (3.9в) В1 — Это так называемая расширенная матрица, Я 1 а и 13,— ее треугольные множители, вычисленные в гауссовом исключении, Бьорк 13~ показал, что ~(В)= Д (А) ~= 4~Я; х(В,) ~(х(А)Р, аР-.О,/ /2.
(3.10) (3.11) В 122) сообщается, что при и=1 получены хорошие результаты. Однако, как следует из (3,11), знЯчение с(, следует выбирать осмотрительно, Заметим, что при а < о„/~2 число обусловленности х(В~) снова возрастает (см. также 161). Й Пример 5.21. Метод Питерса — Уилкинсона 1641 можно пОлучить из об~цей схемы при 1 = 2 и С1 =1, С2 — — В„ 1~2 С1=О, 0~=8; В~ А, С~= й„ х у! ° Здесь Й а= Р ~ " — ортогональная матрица, т. е.
Йтй = 1, Π— диагональная и Х и матрица; 8 — верхняя треугольная и Х и матрица и К08 — ортогональное разложение А. Если оно Вычисляется методами Хаусхолдера или Гивенса, то 0 1, И Пример 5.23. Другой вариант ортогонального разложения получается, если взять 1 = 2, Р1 = Р2 = 91 = Я2 = 1 и исользовать разложение РАЯ+ Е = КВ8 из примера 5.22. Олагаем и АтА Ат1 С1 — — 1, С1=АтРтК, 01=38Ят; С,=Я„С,=8т, О,=В; Х= у1. (3.14а) (3.14б) (3.14в) (3.14г) Здесь 1.р — гпмнп нижняя трапецеидальная матрица с единичной диагональю, Цр — верхняя треугольная и,", и матрица и 1.р11р А+Е1. 1ФЭ,Л., есть Вычисленное разложение сим метРичной матРиЦы В2 = 1 р1 р.
Можно ожидать, что точность вычислений на втором шаге будет приблизительно такой же, как на первом, если м р.а) О~н(А). Ма| не можем нонааать поденного ооотно. шения, Одна~~ звристические с~ображе~ия у~а~~ва~~, что плОхая обуслОВленнОс~~ А будет скОрей всего отражат~ся в 1.1р, а 1., будет, как правило„хорошо обусловленной 113, 64).
Численная практика подтверждает, что обычно метод устойчив (см., например, 1221). В этом случае у~ = Нс = Я ~ О~+ЯяВг Р2С~ Р1с =(Ю~') '~ЯЗ'В) 'С=Як 'В '~Я') 'Я'АтЬ. (3.15) Заметим, что все матрицы, необходимые для построения у„вычисляются на первом шаге, а матрица К не участвует в фактических вычислениях и, следовательно, не нуждается в хранении.
В случае метода, работающего с плотными матрицами, мы могли бы разместить под диагональю А информацию, необходимую для восстановления К, т. е. можно было бы хранить Й (см. 1761). Если, однако, А — большая и разреженная и применяется метод для разреженных матриц, то мы можем получить значительную экономию памяти '>, так как Р зачастую гораздо менее разрежена, чем А. Этот факт подчеркивается в 14,6) и используется в программе Златева и Нильсена 1961 (см.
~ 5.7). Отметим, правда, что приходится хранить и А, и Ь. Для матриц возмущений Е1 и Е2 справедливы следующие выражения: Е~ — — В| — Р1ВД~ = (АТРТ1~) (ОБЛЕТ) — АТА = А" Рт (808 — РАЯ) Ят = АТРтЕЯт, (3.16) Е В Р В,О ~Цт1.) АтртР О (~т1~трт ~тАтрт) Р ~Етр (3.17) Теперь мы сосредоточим свое внимание на линейных задачах наименьших квадратов с большими и разреженными матрицами коэффициентов А.
Чтобы эффективно решать такие задачи, нужно, по возможности, минимизировать память и время, Для достижения этой цели мы применим технику разреженных матриц„возьмем подходящую стратегию выбора главного элемента, зададим надлежащий коэффициент устойчивости и большой барьер (см., например, 111,66,831). Эти попытки использовать и сохранить разреженность— в особенности последнее — часто будут сопровождаться потерей точности, которую нужно попытаться как-то восстановить.