Прямые методы для разреженных матриц. О. Эстербю, З.Златев (984134), страница 15
Текст из файла (страница 15)
Обычно этого можно добиться, присоединяя итерационное уточнение к общей к-шаговой схеме, как это было сде- линО в гл, 4 для гауссова исключения. Поэтому к двум эта'- пам, описаииым а й 5.2, мы побааим еще и следующий: " Оеиииииыиии ее ираиеаии И. Прим, парии, г,=с — В,у„1=1(1) ц — 1„ д., = Нг,.„1 = 1 (1) с1 — 1; у =у +6~, 1=1(1)ц — 1. Для окончания итерационного процесса нужно принять какой-нибудь критерий остановки (см. 15, 75) ); полученный вектор у„будет рассматриваться как приближение к у. Наконец пользунсь свнзью между х и У, по в~ктоРУ Уч вычисляем х.
Предположим на время, что формулы (2.10) и (4.1)— (4,3) могут быть реализованы без ошибок округлений. Введем вектор длЯ любого 1 ДОкйзйтВльстВО, Из формул (4.1) — (4.4) и (2.13) Быте" кает, что У! У=у~ 1+ п1 1 — У =У +НА«у — У; 1)+в) — у (4.7) = 1' «У; ~ У) + Н8 Теорема 5.25. Если ٠— ЙОВАВдОВитВльцость ВВкторОВ, ОПРВдВЯЯВЯпЯ ФОРЯЦлой (4.2), тО д1=Р 'д~ (4,8) Доказатсльство.
Утверждение теоремы сразу следует из равенства А = А-(+ ((1( — А-)) =А( )+ Ж(у; ( — у() =Р(1! ( (4,9) Теорема 5,26. Если р(Р) - 1, то у = у, + ~ б, — (Н~,)-' Нз откуда следует (4.1О). И Следствие 5.27. При р(Р) ~ 1 итерационный процесс (4,1) — (4.3) сходится к точному решению задачи (2.1), если выполнено одно из следующих трех условий: а=О, К=В)~, 1-1 =Йв'„ где Н вЂ” произвольная ') матрица из Р" )б Р. Следствие 5.28. ИтграЦионный арнесс (4.1) — (4.3) схо- дитсн, если !!т!!( !, еде !! !! обозначает матричную норм~ порождсннЦк) Выбранной в8кторной нормой.
Замечание 5.29. Условие (4.16) следует охарактеризовать как чисто теоретическое. В (4.14) (4.15) (4.16) для любоио фиксироданноио натуральноио Доказательство, Символ р (Р) обозначает спектральный радиус Р, и из р (Р) ~ 1 вытекает 1ип Р' = О„„~, Р' = (1 — Р) '. (4.11) ) -+со )=О Поэтому, согласно (4.6) и (2.13), !!т у,=у+ ~ р)Не=у+(НВ,) 'Не. (4.!2) ! -+ОО )неО Из (4.3) выводим Практическая ценность границы (4.21) невелика, поскольку теоретические значения для 1(гп, и) обычно очень грубы и приводят к сильной переоценке ~~Р~~.
Но (4.21) и (4.17) выявляют важную зависимость между числом обусловленности ((4,22) или (4.23)) и барьером, показывая еще раз, что при большом х нужно выбирать Т поменьше. Для матриц класса Р2(гп, п,с, г,а) можно управлять числом обусловленности, меняя а. Взаимовлияние Т и а для зтих матриц демонстрирует табл. 5,1.
Через таха обознаТаблнЧа Ы. Максимальное значение а для матрицы Р2 (22, 22, 11, 2сб), допускакрщее успешное решение задачи при данном значении Т Ортотональные нреобрааоаанна Гауссово исключение сомэ 1х — х!1 1х — х1пп чена наибольшая степень 2, допускающая успешное решение задач~. Гауссово исключение выполняетс~ усовершенство' ванным вариантом программы 51КЬМ 193,94), а ортогональ- ные преобразования — программой 1Л.ЫО1 195~. Параметр . СОИ1.р есть оценка числа обусловленности, вычисляемая фортранной подпрограммой„содержащейся в 11О9~. о,б. Ортогоиальиые преобразования В ~ 5.7 мы подробней изучим детали реализации двух- шаговой схемы, основанной на ортогональных преобразованиях.
Пока же мы обсудим ортогональное Й08 разложение п1 )( и матрицы А. Как в учебниках, так и на практике очень популярны два м~~ода: ~снов~н~~й на плос~~х вращениях метод Гивенса 142,431 и метод Хаусхолдера 15Ц, использующий элементарные отражения. Вычислительная стоимость методов Гивенса и Хаусхолдера, измеряемая числом операций умножения и извлечения квадратного корня, для случая плотных матриц приведена в табл. 5,2. Из нее ясно, почему начиная с 1959 г.'~ метод Хаускпллера был более популярен. ') Статья Хауояоанера 'ьбт1 а которой аоераеа быто предложено копользовать элементарные отражения для ортогональной триангуляризацин матрицы, была опубликована в 1958 г, — Прим. перев.
Умножения Метод Кеядретные корни ХиусхОлдер В последние годы ситуация изменилась: благодаря результатам работ 136, 37, 38, 50~ вычислительная стоимость метода Гивенса снизилась приблизительно до уровня стоимости метода Хаусхолдера. Поскольку для разреженных матриц ме- тоД Гивенса подходит больше, мы обсудим его более подробно. Ортогональное приведение состоит из и больших шагов, на каждом из которых все поддиагональные элементы некоторого столбца преобразуются в нули. Каждый большой шаг составляют несколько малых шагов — плоских вращений," каждое из них аннулирует Один элемент. Если имеется $1~ таких элементов„то К-Й большой шаг состоит из з1, малых шагов.
Б обычном методе Гивенса малый шаг можно описать ПРОИЗВЕДЕНИЕМ 1, Е 1, К+1''' 1, и а1. Р1, к+1 ° ° где один из элементов а11, а11 будет аннулирован, а 4= =01 —— 1, В д~йствительности первым двум сомножителям соответствуют п1,','п1 матрицы с единицами в диагональных позициях всех прочих строк, поэтому мы показываем только элементы, участвующие в вычислениях.
Ясно, что для каждого элемента а..., затронутого преобразованием (5.1), нужны два умножения. Чтобы уменьшить работу, Джентльмен предложил иную факторизацию произведения двух первых матриц в (5.1): ~т1 При последующих умножениях мы используем только последнюю матрицу из (5.2). Легко видеть, что в этой схеме Для ~аждог~ а... ~р~бу~~~я ~ол~ко Одно ум~Ожени~ ~и Одно сложение) . Мы еще должны решить, какой из элементов — ап, или ап,— аннулировать, и дать формулы для коэффициентов сс, р и '~.
Если ~! аф ~~ ~1~а~~~' то а~к .'= О, причем а~к/а к с' = М~Я' Б обоих случаях Приведенные формулы показывают, что, помимо экономии и умнОжениях, мы можем еще и избежать извлечения квадратных корней — которые были необходимы при определении у и а в обычном методе Гивенса, — если будем хранить велич~и~ Й'. в~~сто д, Вот ~оч~му этот метод называю~ методом Гивенса без квадратных корней, Кроме того, скоро мы увидим, что в наших вычислениях используется матрица В~, а не В (см. (3.15)). Для величин д~ устанавливаются начальные значения: 1,' = 1, 1 = 1 (1) гп. (5,7) Если задача предусматривает веса, то начальными значениями будут квадраты весов. На каждом малом шаге пересчитываются два коэффицеинта д~. д',:=д,'"~?', д';:=д» ~' (5.8) (см. (5.2), где у2 выражается формулой (5.6) ).
1 Из алгоритма (5,3) — (5.6) следует, что — <: у' < 1, поэтому'~ элементы матрицы В' убывают, но не слишком быстро. Для очень большой задачи имеет смысл во избежание : машинных нулей контролировать величины ф и при необходимОсти изменять масштабы. Рассмотрим еще раз данные табл. 5.2. При указанных выше модификация~ метод Джентльмена — Гивенса требует цримерно такой же вычислительной работы, как и метод Хаусхолдера; ОднакО В случае плотных матриц по-прежнему нет причин для того, чтобы предпочесть первый метод второму (вспомним, например, проблему машинных нулей для коэффициентов 4).
Можно упомянуть, что трапецеидально- 1 треугольное Ш-разложение обходится всего лишь в — 1пп 3 арифметических операций. Тем не менее обычно предпочитают ортогональные методы, поскольку они считаются более устойчивыми '1. При работе с разреженными матрицами сохранение разреженности становится важным вопросом, и в этом отношении метод Хаусхолдера менее привлекателен. Если в 1-м большом шаге метода Хаусхолдера участвуют за+ 1 строк, то каждая преобразованная строка есть линейная комбинация всех за+ 1 прежних строк; следовательно, ее портрет соответствует объединению этих за + 1 строк (если пренебречь возможностью точных взаимных сокращений).
В случае метода Гивенса портрет двух строк, комбинируемых в ~~ло~ шаге, после преобразо~~~ия буде~ соответствовать об~единению этих строк'1. Если одна из строк не принимает участия в последующих малых шагах (внутри рас- х о и х к х м о в х х Рис. 5.4. Метод Хаусхолдера. После первого большого шага — 9 новых ненулевых элементов. Рис. 5.5.
Метод Гивенса. После первого большого шага — 7 новых нену- левых элементов. Рис. 5.6. Гауссово исключение. После первого большого шага — 3 новых ненулевых элемента. м Смысл последних двух фаз не вполне ясен. Само по себе 1ЛЗ-разложение не помогает решить задачу наименьших квадратов. Если же оно проводится как первый шаг метода Питерса — Уилкинсона 1см. пример 5,21), то к его вычислительной стоимости нужно прибавить стоимость формирования и последующей факторизации матрицы 1.т1. системы нормальных уравнений.
О~но~ительно сравнительных дост~~н~~~ ортогональных методов и метода Питерса — Уилкинсона в случае разреженных задач нет полной ясности — см„например, упомянутый в предисловии обзор переводчика. — Прим. перев. '~ До преобразования. — Прин. перев, сматриваемого большого шага), то никаких новых ненулевых элемеитов В НЕЙ НЕ ПО~~И~СИ. Для полноты картины заметим, что в трапецеидальнотреугольном разложении в ведущей строке не возникает заполнения. Таким образом, это разложение является наилучшим методом в отношении сохранения разреженности. Проиллюстрируем простым примером появление новых элементов.
На рис. 5.3 — 5.6 показаны исходная (квадратная) матрица и матрицы, полученные из нее после первого большого шага соответственно методов Хаусхолдеря, Гивенса и Гаусса. о.б. Стратегия выбора главного элемента Как мы видели в ф 5.5, гауссово исключение обычно дает меньшее заполнение, чем ортогональные методы, а среди этих пОслеДних метОД Гивенса слеДует преДпочесть метОДу Хаус холдера (см. также теОретические результаты Даффа и Рида 121] и Зльфинга 1311). Обсуди~ тепер~ стратегию выбора главного эл~м~нт~ В метОДе Гивенса, направленную на тО, чтОбы, по ВОзмОжности, уменьшить размер заполнения.