Дж. Деммель - Вычислительная линейная алгебра, страница 76

Для этого случая может быть также использован алгоритм ЯУММЩ, предложенный Пэйджем и Сондерсом [194]; в нем на каждом шаге для невязки обеспечивается условие г».1» »[А, б). К сожалению, существует немного матриц, отличающихся от симметричных и допускающих алгоритмы типа СС, которые одновременно обладают такими свойствами: 1. минимизируется норма невязки [[ге[[э либо поддерживается условие ть1Хь[А, б); 2. во внутреннем цикле используется фиксированное, не зависящее от Й число скалярных произведений и операций типа эахру. По существу, алгоритмы с этими двумя свойствами возможны лишь для матриц вида ем(Т + о1), где Т = Тт [либо ТН = [НТ)т для некоторой симметричной положительно определенной матрицы Н), 6 — вещественное, а и— 334 Глава б. Итерационные методы для линейных систем комплексное число [102, 251].

Для таких симметричных или специальных несимметричных матриц А, как оказывается, можно найти короткие рекурсии типа рекурсий алгоритма Ланцоша, позволяющие вычислить ортогональный базис ]оы..., д„] подпространства Кь(А, Ь). То обстоятельство, что рекурсия, модифицирующая векторы оы содержит лишь несколько членов, означает высокую эффективность метода. Для несимметричных матриц А общего вида подобных коротких рекурсий не существует. Для таких матриц можно применить алгоритм Арнольди. Тогда вместо трехдиагональной матрицы Ть = Яь АЯн получим полноценную верхнюю хессенбергову матрицу Нь = ЬД АЯы Это разложение используется в алгоритме СМВЕЯ (Яепега11хеб ппшпппп геэ16па1) для выбора вектора хь = Яьуь 6 Кя(А, Ь), минимизирующего невязку ]]ть]]а = ]]Ь вЂ” Ахь]]т = ]]Ь - Аб),у,]], = ]]Ь вЂ” ЯН1„1~)Яьуь]]а согласно уравнению (6.30) = ]]Я~Ь вЂ” НЯ~Яьуя]]т, так как матрица О ортогональна, = "]]Ь]].- Й'„„"„' согласно уравнению (6.30), а также потому, что первым столбцом матрицы 0 = ]фы Я„] является вектор ЬДЬ]]т ег]]Ь]]т — Н Уь Н ьн В блоке Ньн отлична от нуля только первая строка, поэтому для элементов вектора уь получена (Й + 1) х Й-задача наименьших квадратов с верхней хессенберговой матрицей.

В силу этого свойства, 6)В;разложение, необходимое для решения задачи, может быть найдено посредством й плоских вращений за 0(йт) операций вместо 0(йе) в общем случае. Кроме того, для хранения Яь требуется 0(йп) слов памяти.

Один из способов ограничения трудоемкости и памяти состоит в том, чтобы обноелягаь метод. Это значит, что берется приближенное решение хь, вычисленное после й шагов, СМВЕЯ применяется к решению линейной системы Ад = гь = Ь вЂ” Ахю а затем находится новое приближенное решение хь + И. Описанный алгоритм называется СМВЕЯ(й). Однако даже этот алгоритм более трудоемок, чем СС, где стоимость внутреннего цикла вообще не зависит от Й. Другой подход к решению несимметричных линейных систем заключается в том, чтобы отказаться от задачи вычисления ортонормированного базиса в подпространстве Кь(А, Ь) и строить вместо этого неортогональный базис, в котором бы А снова приводилась к (несимметричной) трехдиагональной форме.

Этот подход называется несимметричным методом Ланцоша; он требует возможности формировать матрично-векторные произведения с участием как матрицы А, так и матрицы Ат. Это существенное требование, поскольку вычисление произведений вида Атя подчас затруднительно или вообще невозможно (см.

пример 6.13). Преимущество трехдиагональной формы — в том, что системы с трехдиагональными матрицами легче решать, чем хессенберговы. ззб б.?. Быстрое преобразование Фурье Симы нчналнА? Нет Доступна лн Ат? Определена лн А? Да Нет Нет Нухшо ян экономить память? Известны лн крайние собственные значения? Хорошо лн обусловлена А? Хорошо ли обусловлена А? Нет Да Нет Да Да Нет Нет Да Применять К нормальным МК уравнениям применитьСО Применить СО Применить СО с чебы- шеяекям ус ко нясм Применить ОМКЕБ Применить СОБ, иля ВЬСОБшь, ОМКЕБ к При ыешпь М1?ч КЕБ яяи какой-лябе Метод лая несимметричной А Рис. 8.8.

Дерево решений прн выборе итерационного алгоритма для системы Ах = Ь. Приняты аббревиатуры: ВЬССБгаь = стабилизированный метод бисопряженных градиентов, СгМН = метод квазиминнмальных невязок; ССБ = СС ейиагеб. 6.7. Быстрое преобразование Фурье В данном разделе символ 1 всегда обозначает ч/ — Т. Мы покажем вначале, как можно решить двумерное уравнение Пуассона, используя умножения на матрицу собственных векторов для Т~~. Прямолинейная реализация этого подхода стоила бы 0(А? ) = 0(п~?~) операций, что чересчур дорого. Затем мы объясним, как этот подход может быть реализован с помощью алгоритма ГГТ за 0(гч'г 1обтч') = 0(п1ояп) операций, что лишь множителем 1о8 и отличается от оптимального результата. Данный подход является дискретным аналогом метода рядов Фурье, применяемого к исходным дифференциальным уравнениям (6.1) или (6.6).

Позднее мы уточним характер этой аналогии. Недостаток же — в том, что базисные векторы могут быть очень плохо обусловленными. Иногда базисный вектор вообще не удается определить; это явление называется обрывом. Потенциальная эффективность данного подхода вызвала к жизни большое число исследований, посвященных тому, чтобы избежать этого типа неустойчивости или ослабить ее (методы Лаицоша с загллдывонием вперед), а также конкурирующие методы, среди которых отметим методы бисопрлзгсенных градисюпов и квазимин мааьных нсвлзок. Существуют и версии этих методов, не требующие произведений с участием матрицы Ат, в том числе методы СОЯ (соп)пйайе ягад!епФз гс?пагед) и стабилизированный метод бисопрлзшенных градиентов.

Ни один метод не является наилучшим для всех случаев. На рис. 6.8 показано дерево решений, дающее простые указания, какой метод следует попробовать первым, если более подробной информации о матрице А (например, что А возникает из уравнения Пуассона) не имеется. 336 Глава 6. Итерапяовные методы для линейных систем Пусть Тм = ЯЛЯт — спектральное разложение матрицы Тм, указанное в лемме 6.1. Будем работать с записью двумерного уравнения Пуассона в виде уравнения (6. 11): т )г+1гт =Ь Р.

Подставляя сюда Тм = ЯЛЯт и умножая слева на Ят, а справа на Я, получим Рт(РЛРт))гав+ Хт1г(ХЛИт)~ Ет(лгр)У или ЛГ + ГЛ = ЬгР', где у' = Я~уЯ и Р' = г ~ГЯ. Элемент (з,Ь) в последнем равенстве равен (Л)~ + Ъ"Л)уг = Л,';„+ оеьЛ, = ЬгУ,',. Определяя отсюда о'.„, находим Ь'1'.г Л.

+Л Это дает первый вариант нашего алгоритма. Алгоритм 6.13. Решение двумерного уравнения Пуассона с помощью спектрального разложения Тн = г ЛЯ~ г 1) Е" = ЯтРт, 2) для всех г' и Ь вычисяитпь и'г — — т.-+г~ —,' З)1 =г1 г' Стоимость шага 2 составляет ЗЖг = Зп операций; шаги 1 и 3 состоят из 4 матрично-матричньгх умножений с участием матриц Я и г"т = с; их стоимость при обычном алгоритме умножения равна 81г"г = 8пг1г операций.

В следующем разделе будет показано, что умножение на матрипу Я, по существу, есть то же самое, что вычисление дискретного преобразования Фурье. С помощью алгоритма ге Т это преобразование может быть проделано за 0(Лгг!ойдо) = 0(п1обп) операций. (Используя введенный в равд. 6.3.3 язык кронекеровых произведений и, в частности, спектральное разложение матрицы Тмхм, указанное в предложении 6.1, т. е.

т„„, = 1 З т„+ т З 1 = (г З г) (1 З Л+ Л З 1) (г З г) т, можно следующим образом переписать формулу, обосновывающую алгоритм 6.13: чес(У) = (тггхгг) ~ тес(Ь~г') =((ЯЗВ) (1ЗЛ+ЛЗ1) (ЕЗЯ)т) ' вес(ЬгР) =(ЯЗЯ) т (1ЗЛ+ЛЭ1) ' (г Зг) ' нес(ЬгР) = (Я З Я) ° (1ЗЛ+ Л З1) ~. (Е Э Я~) тес(Ь~г). (6.47) Мы утверждаем, что выполнение указанных матрично-векторных умножений справа налево математически зквивалентно алгоритму 6.13 (см. вопрос 6.9). Это заодно показывает, как распространить алгоритм на уравнения Пуассона более высокой размерности.) 337 6.7. Быстрое преобразование Фурье 6.7.1.

Дискретное преобразование Фурье В этом подразделе строки и столбцы матриц будут нумероваться числами от О до Х вЂ” 1, а не от 1 до Ю, как обычно. Определение 6.17. Дискретное преобразование Фурье (ПРТ) 1ч'-вектора х есть вектор у = Фх, где Ф вЂ” матрица порядка Ю, определяемая следующим г ~ образом. Пусть ы = е в = сов $ — зэ1пф есть примитивный корень Л-й степени из единицы. Тогда р ь = ыдь. Обратное дискретное преобразование Фурье (10РТ) вектора у есть вектор х = Ф ~у. Лемма 6.9. Матрица г Ф симметрична и унитарна, поэтому Ф чъ'и нФ*= нФ Итак, и ПРТ, и 10РТ являются всего лишь матрично-векторными умножениями и могут быть бесхитростно реализованы за 2Хз флопов.

Рассматриваемая операция называется дискретным преобразованием Фурье ввиду своей тесной математической связи с двумя другими типами анализа Фурье (см. таблицу) . Мы конкретизируем зту тесную связь двумя способами. Во-первых, мы покажем, как решить модельную задачу с помощью 1эРТ, а затем как решить уравнение Пуассона (6.1) посредством рядов Фурье. Этот пример послужит стимулом для отыскания быстрого метода умножения на матрицу Ф, потому что такой метод приведет к быстрому способу решения модельной задачи. Этот быстрый метод называется быстрым преобразованием Фурье, или РРТ.

Вместо 2Юз флопов в нем требуется гораздо меньшее число операций, приблизительно равное -Ю1ояэ Х. Алгоритм РРТ будет выведен с использованием з второго математического свойства, общего для различных видов анализа Фурье, а именно преобразования свертки в умножение. Доказательство. Ясно, что Ф что Ф Ф = Х 1. Имеем (ФФ)Б поскольку ы = ы '.