Дж. Деммель - Вычислительная линейная алгебра, страница 75
Описание файла
PDF-файл из архива "Дж. Деммель - Вычислительная линейная алгебра", который расположен в категории "". Всё это находится в предмете "квантовые вычисления" из 7 семестр, которые можно найти в файловом архиве МГУ им. Ломоносова. Не смотря на прямую связь этого архива с МГУ им. Ломоносова, его также можно найти и в других разделах. .
Просмотр PDF-файла онлайн
Текст 75 страницы из PDF
Мы можем видеть это по тому, что число и'; различных собственных значений в ряде случаев отличается от числа к; шагов, требуемых для сходимости, хотя, как было показано, теоретически эти числа должны совпадать. Все же 4 и к; были величинами одного порядка для всех систем. В действительности, если бы метод С6 был проведен в точной арифметике и полученные решения и невязки были сравнены с аналогичными величинами, вычисленными в машинной арифметике, то, скорей всего, результаты с ростом номера й итерации быстро стали бы совершенно различны.
Все же, если матрица А не слишком плохо обусловлена, результаты с плавающей точкой в конечном счете сойдутся к искомому решению системы, поэтому метод остается очень полезным. Тот факт, что точные и машинные (промежуточные) результаты могут быть совсем не похожи, интересен, но не препятствует практическому применению метода.
Вслед за открытием метода С6 было доказано, что в точной арифметике он дает точное решение системы за и шагов, поскольку невязка гк ы должна быть ортогональна к другим п ортогональным векторам гм..., г„, а потому равна О. Иначе говоря, метод С6 воспринимался как прямой, а не как итерационный. ззо Глава 6. Итерационные методы для линейных систем Поскольку на практике (точная) сходимость за и шагов не наблюдалась, метод стал рассматриваться как неустойчивый, а затем был на много лет заброшен. В конечном счете, было признано, что С6 — это очень хороший итерационный метод, часто дающий приближенные решения высокой точности уже после й « п шагов. Не так давно был выполнен тонкий обратный анализ ошибок, объясняющий наблюдаемое поведение метода С6 в арифметике с плавающей точкой и его возможное отличие от поведения в точной арифметике [123].
Картина сходимости метода может включать в себя протяженные участки «плоскогорья~, где 1)ть((т почти не уменьшается в течение многих итераций, перемежаемые периодами быстрой сходимости. Эту картину можно объяснить, показав, что метод С6, применяемый к системе Ах = Ь в арифметике с плавающей точкой, ведет себя в точностаи так, как тот же метод, применяемый в точной арифметике к системе Ах = Ь. Матрица А близка к А в следующем смысле: ее порядок намного больше порядка А, однако все собственные значения матрицы А расположены узкими кластерами возле собственных значений матрицы А.
Таким образом, плоскогорья в картине сходимости соответствуют тому, что много- член дю лежащий в основе метода С6, приобретает все больше и больше нулей возле кластеров собственньгх значений матрицы А. 6.6.5. Предобуславливание В предыдущем разделе мы видели, что скорость сходимости метода С6 зависит от числа обусловленности матрицы А или, более общо, от распределения ее собственных значений. Тем же свойством обладают и другие методы крыловского подпространства.
Предобуславливание означает, что система Ах = Ь заменяется системой М «Ах = М "Ь, где матрица М есть приближение к А со следующими свойствами: 1. М симметрична и положительно определена. 2. Матрица М 'А хорошо обусловлена или же имеет небольшое число собственных значений на периферии своего спектра. 3. Система Мх = Ь решается легко. Продуманный, зависящий от конкретной задачи выбор М может сделать число обусловленности матрицы М 'А много меньшим числа обусловленности матрицы А и тем самым в огромной степени ускорить сходимость.
В действительности, хороший предобуславливатель часто необходим для того, чтобы итерационный метод вообще сходился. Поэтому большая часть современных исследований по итерационным методам ориентирована на отыскание эффективных предобуславливателей (см. об этом также разд. 6.10). Мы не можем применить метод С6 непосредственно к системе М 1Ах = М 'Ь, поскольку матрица М "А в общем случае не является симметричной.
Мы построим предобусловленный метод сопряженных градиентов следующим образом. Пусть М = 1~ЛЯт — спектральное разложение матрицы М; положим М~~~: — ЯЛИ~Ц~. Заметим, что матрица Ми~а также симметрична и положительно определена и (М'~т)т = М. Теперь умножим систему М 'Ах = М Ч на матрицу М'~т, что дает новую симметричную положительно определен- 6.6. Методы крыловского подпросгравства ную систему (М з1гАМ 1~г)(Мз1гх) = М з~гЬ, или Ах = Ь. Отметим, что А и М 1А имеют одни и те же собственные значения, поскольку зти матрицы подобны (М 1А = М 1~гАМ1~г).
Теперь применим метод СС к системе Ах = Ь неявно, таким образом, чтобы не было необходимости в умножениях на М '~з. Это приводит к следующему алгоритму: Алгоритм 6.12. Предобусловленный алгоритм СС: й = О; хс = О; го = Ь; р1 = М 'Ь; ус = М 'го гереа1 й = й + 1 г = А.рь иь = (у„,гг 1)((р„г) хь = хь 1+ игра гг = гь 1 — иьг уь = М ~ге рьз.з — †(уг гг)((уь гг з) рг+1 = уь + Рг+зрг пока число 6гьбг не станет достаточно малым Теорема 6.9. Пусть А и М вЂ” симметричные положительно определенные матрицы.
Положим А = М з~зАМ1~з и 6 = М 11зЬ. Алгоритм СС в применении к системе Ах = Ь, т. е. алгоритм, описываемый формулами й=О;хо=О;гь=Ь;р1 =Ь; гереас й = й + 1 г=А рь ~„-т „-„) ~1д1тхг =ха з+йгрь гь = гь г — йьг аль+1 = (гьгь)~Я бг з) рг+1 = гг + рьетрг пока число Оггбг не станет достпапючно малым и алгоритм 6.1з связаны следующим образом: рг = Иг~ иг = иь, г=М хь = М~~~хг, гг = М ~ гь, Рг = М~' рг. Следовательно, векторы хг сходятся к произведению матрицы М 1~г и ре- шения системы Ах = 6, т. е. к вектору М 1~гА 16 = А 1Ь.
По поводу доказательства теоремы см. вопрос 6.14. Теперь мы опишем некоторые распространенные типы предобуславливате- лей. Заметим, что преследуемые нами две цели — минимизация числа обусло- вленности матрицы М 'А и сохранение простоты решения системы Мх = 6,— Глава б. Итерационные методы для линейных систем ° Если диагональные элементы матрицы А сильно различаются по величине, то можно использовать простой диагональный предобуславливатаель М = йая(аы,..., а„„).
Можно показать, что такой выбор дает число обусловленности произведения М 'А, не более чем в п раз превышающее минимальное значение по всем возможным диагональным предобуславливателям [244]. Этот выбор М называют еще предобуславливанием Якоби. ° Обобщая первый тип предобуславливания, станем рассматривать А как блочную матрицу с квадратными диагональными блоками Ап. .4ы .- .4гь 1 Аы Аьь ~ А= [ Тогда среди всех блочно-диагональных предобуславливателей с блоками Мн тех же размерностей, что и Ан, выбор Мн = А„с точностью до множителя й минимизирует число обусловленности матрицы М ИгАМ' '~г [69]. Этот выбор называют также блочным предобуславливанием Якоби.
° Подобно методу Якоби, метод ВБОК может быть использован в качестве (блочного)предобуславливателя. ° Под неполным разложением Холесского Ь.ьт матрицы А понимают приближение А — ЛЬг, в котором Ь ограничена требованием иметь конкретный тип разреженности, например тот, что имеет исходная матрица А. Другими словами, не допускается заполнения в ходе алгоритма Холесского. Произведение ЬЬ~ используется в качестве матрицы М. (Для несимметричных систем существует соответствующий неполный Ш предобуславливатель.) ° Если матрица А представляет уравнение (например, уравнение Пуассона) на физической области П, то используется декомпозиции области. До сих пор в случае уравнения Пуассона П была для нас единичным квадратом.
Более общо, область П может быть разбита на непересекающиеся (нли немного перекрывающиеся) подобласти П = О Йу, на каждой из которых уравнение может быть решено независимо от других подобластей. Пусть, например, решается уравнение Пуассона и подобластями являются квадраты или прямоугольники. Тогда подзадачи могут быть решены очень быстро вступают в противоречие друг с другом: выбор М = А минимизирует число обусловленности матрицы М гА, но оставляет систему Мх = о столь же трудной для решения, что и исходная система.
Выбор М = 1 делает решение системы Мх = о тривиальным, но не меняет числа обусловленности матрицы М 'А. Поскольку система с матрицей М решается во внутреннем цикле алгоритма, мы ограничимся обсуждением матриц М, для которых подобные системы легко решаются, и укажем те из них, что потенциально могут уменыпить число обусловленности произведения М зА. 333 6.6. Методы крыловского подорострглстве с помощью алгоритма РГТ [см. рэзд. 6.7). Эти подзадачи соответствуют использованию блочно-диагональной матрицы М (если подобласти не пересекаются) или произведения таких матриц (если подобласти перекрываются).
Более подробно эти вопросы обсуждаются в разд. 6.10. Многие предобуславливатели описанных типов реализованы в программных пакетах РЕТЯс [232] и РАВРВЕ (ИЕТ1 1В/зса1араск/рагргейаг.йв). 6.6.6. Другие алгоритмы крыловского подпространства для решения системы Ах = б До сих пор мы были сосредоточены на решении симметричных положительно определенных линейных систем и минимизации А «-нормы невязки. В этом разделе будут описаны методы для других типов линейных систем н даны советы о том, какой метод следует использовать, основанные на простых свойствах матрицы. Наши выводы суммируются на рис. 6.8 [15, 107, 136, 214].
Детали, более полные рекомендации по поводу выбора метода и программы можно найти в ХЕТЫВ/гешр1а1еэ. Всякую систему Ах = б можно преобразовать в симметричную положительно определенную систему АтАх = Атб [или систему ААту = б, х = Агу), называемую нормальными уравнениями. Это относится и к задаче наименьших квадратов ш1п, [[Ах — б[[т. Указанное преобразование позволяет применить метод СС, если есть возможность умножать векторы и на А, и на Ат. Число обусловленности каждой из матриц АтА и ААт равно квадрату числа обусловленности матрицы А, поэтому СС будет сходиться медленно, если А плохо обусловлена, и быстро, если А обусловлена хорошо (или же АтА имеет «хорошее» распределение собственных значений в смысле, обсуждавшемся в разд.
6.6.4). Если А — симметричная положительно определенная матрица, то вместо А "-нормы невязки можно минимизировать ее 2-норму. Эта установка определяет алгоритм минимальных невязок, илн М1ХВЕЯ [194]. Поскольку М1б1ВЕБ более трудоемок, чем СС, и часто менее точен вследствие численной неустойчивости, он не используется для положительно определенных систем. Однако, в отличие от СС, М1ХВ.ЕЯ может быть применен к симметричной незнакоопределенной матрице.