Дж. Деммель - Вычислительная линейная алгебра (1156793), страница 51
Текст из файла (страница 51)
«Разделяй-и-властвуй». В настоящее время это самый быстрый метод вычисления всех собственных значений и собственных векторов симметричной трехдиагональной матрицы порядка п, большего 25. (ЬАРАСК-реализация метода, называемая ивьечб, переключается на ЯВ-итерацию, если п < 25.) В наихудшем случае алгоритм «разделяй-и-властвуй» требует 0(пз) флопов, но константа, упрятанная в этом символе, на практике оказывается весьма малой. На большой выборке случайных тестовых матриц в среднем затрачивается лишь 0(пз з) флопов, а для некоторых специальных распределений собственных значений даже 0(пз) флопов.
Теоретически алгоритм может быть реализован за 0(п . 1ояг и) флопов, где р — неболыпое целое число [131]. Эта сверхбыстрая реализация использует быстрый многополюсный метод (РММ) [124], первоначально придуманный для совершенно другой задачи (а именно вычисления электрических сил взаимодействия в ансамбле из п заряженных частиц).
Вследствие алгоритмической сложности этого сверхбыстрого метода, в настоящее время следует выбирать ь)В-итерацию, если вычисляются все собственные значения, и алгоритм «разделяй-и-властвуй» (без использования РММ), если вычисляются все собственные значения и собственные векторы. 4. Бисекция и обратная итерация. Висекцию можно использовать для вычисления только некоторого подмножества собственных значений симметричной трехдиагональной матрицы, скажем, тех, что расположены на отрезке [а, Ь] или [а;, гт, ]. Она обходится лишь в 0(ий) флопов, где й — число требуемых собственных значений. Поскольку для ЯВ-итерации нужны 0(пз) флопов, бисекция может быть гораздо быстрее при й « и.
Соответствующие собственные векторы могут быть найдены с помощью обратной итерации (алгоритм 4.2). В наилучшем случае, когда собственные значения «хорошо разделены» (смысл этого термина мы подробно объясним позже), обратная итерация также стоит лишь 0(пй) флопов. Это значительно меньше, чем стоимость ЯВ;итерации или алгоритма «разделяй-и-властвуй» (без ЕММ), даже при необходимости вычислять все собственные значения и собственные векторы (т. е.
при й = п). Однако в наихудшем случае, когда имеются обширные кластеры собственных значений, трудоемкость обратной итерации возрастает до О(пйз) флопов; при этом не гарантируется качество вычисленных собственных векторов (хотя на практике они почти всегда вычисляются с хорошей точностью). Поэтому при необходимости вычисления всех (или большей части) собственных значений и собственных векторов Майаь проверяет, является ли аргумент команды езй симметричной матрицей и, если да, использует симметричный алгоритм. 224 Глава 5. Симметричная проблема собственных значений следует выбирать алгоритм «разделяй-и-властвуй» или чьей-алгоритм, особенно прн наличии кластеров собственных значений.
Бисекция и обратная итерация реализованы как опции в ЬАРАСК-программе аауечх. В настоящее время ведутся активные исследования по проблеме близких собственных значений в обратной итерации. Они могут превратить ее в самый быстрый метод вычисления всех собственных векторов (если опустить теоретическое сравнение с алгоритмом «разделяй-и-властвуй» на основе РММ) ~105, 203, 83, 201, 176, 173, 175, 269). Однако программ, реализующих этот улучшенный вариант обратной итерации, пока нет. 5.
Метод Якоби. Этот метод, восходящий к 1846 г., исторически является старейшим для проблемы собственных значений. Обычно он много медленнее любого из перечисленных выше методов, имея трудоемкость 0(пз) флопов с большой константой внутри символа О. В то же время интерес к методу сохраняется, потому что подчас он дает гораздо более точные результаты, чем другие методы. Причина состоит в том, что метод способен иногда достигать уровня относительной точности, описываемого в разд. 5.2.1; при этом малые собственные значения вычисляются много точнее, чем в конкурирующих алгоритмах ~82).
Мы обсудим это замечательное свойство метода Якоби в разд. 5.4.3, где показано, как следует вычислять сингулярное разложение. В последующих разделах мы более подробно рассмотрим названные алгоритмы. В разд. 5.3.6 представлены результаты экспериментального сравнения их производительности. 5.3.1. Трехдиагональпая С)В;итерации Вспомним, что ч4В.-алгоритм для несимметричных матриц состоит из следую- щих двух этапов: 1. Применить к заданной матрице А алгоритм 4.6, вычисляющий ортогональную матрицу Я, такую, что цА(,Р = Н имеет верхнюю форму Хессенберга. 2. Применить ь)Гх-итерацию к матрице Н (как это описано в разд.
4.4.8). В результате будет получена последовательность верхних хессенберговых матриц Н = Не, Ны Нз, ..., сходящаяся к вещественной форме Шура. Наш первый алгоритм для симметричной проблемы собственных значений имеет точно такую же структуру: 1. Пользуясь модификацией алгоритма 4.6, описанной в разд.
4.4.7, вычислить для заданной матрицы А = Ат ортогональную матрицу ч, такую, что 1ЗАЧг = Т имеет трехдиагональпую форму. 2. Применить ()Е-итерацию к матрице Т. В результате будет получена последовательность трехдиагональных матриц Т = Те, Т„Тз, ..., сходящаяся к диагональному виду. То обстоятельство, что матрицы Т, в ч4В-итерации остаются трехдиагональными, можно обосновать так: поскольку матрица ЯАсэт одновременно симметричная и верхняя хессенбергова, она должна быть и нижней хессенберговой, т.е.
трехдиагональной. Именно это обстоятельство делает каждый шаг ь)В- итерации очень дешевым. Приведем подсчет числа операций: 5.3. Алгоритмы для симметричной проблемы собственных значений 225 1. Приведение А к симметричной трехдиагональной матрице Т стоит «пг + 0(пг) флопов или г пг + 0(пг) флопов, если нужны и собственные векторы.
2. Один шаг трехдиагональной ОВ-итерации с одинарным сдвигом («вытеснение выступа») стоит бп флопов. 3. При вычислении всех собственных значений матрицы Т требуется в среднем только 2 ЯВ-шага на одно собственное значение, что в совокупности дает бпг флопов. 4. Вычисление всех собственных значений и собственных векторов матрицы Т стоит 6пг + 0(пг) флопов. 5. Общая стоимость вычисления всех собственных значений матрицы А (без вычисления собственных векторов) составляет «пг + 0(пг) флопов. 6.
Общая стоимость вычисления всех собственных значений и собственных векторов матрицы А составляет 88~не + 0(п ) флопов. Нужно еще описать, как выбираются сдвиги при реализации ь1 К-алгоритма. Пусть на»чм шаге имеем матрицу а| Ь| Ь, Ь„, Ьь-1 аь Простейшим выбором было бы правило о, = а„. Это ЯВ-итерация с одинарным сдвигом, обсуждавшаяся в равд.
4.4.8. Оказывается, что почти для всех матриц она сходится кубически (это показано в следующем разделе). К сожалению, существуют примеры несходимости этого метода [197, с. 76]. Поэтому, чтобы обеспечить глобальную сходимость, используется несколько более сложная стратегия: в качестве сдвига о, выбирается то из собственных значений подматрицы Ь , которое ближе к а„.
Это правило выбора ( а„1 Ь„ ь — 1 аь сдвигов называется стратегией Уилкинсона. Теорема 5.8 (Уилкинсон). ЧЯ-итерация со стратегией Уилкинсона сходится глобально, причем, по меньшей мере, с линейной скоростью. Длл почти всех матриц метод асимптотически сходится с кубической скоростью. Доказательство этой теоремы можно найти в (197]. Метод реализован ВАРАСК-программой ввуеч. Если нужны только собственные значения, то внутренний цикл алгоритма можно организовать более эффективно (ввсехХ; см. также (104, 200]) по сравнению со случаем, когда вычисляются и собственные векторы (ввьеср.). Пример 5.7. Приведем иллюстрацию сходимости трехдиагональной ЯВ- итерации, примененной к следующей трехдиагональной матрице (ее диагонали 226 Глава 5.
Симметричная проблема собственных значений схематически изображены столбцами): .24929 1. 263 1.263 .96880 †.82812 .48539 -3.1883 †.91563 —.82812 -3.1883 Те — — Сг!д!а8 В таблице показано поведение последнего внедиагонального элемента матриц То последнего диагонального элемента и разности межлу последним диагональным элементом и его предельным значением (т.е. собственным значением о — — 3.54627). Из последнего столбца таблицы очевидна кубическая сходимость ошибки к нулю. В этот момент (з = 3) мы имеем матрицу 1.9871 .77513 .77513 1.7049 — 1.7207 -1.7207 .64214 6 1. 10-зз 6 1, 10-зз -3.5463 Тз —— Спей ай Ее очень малые элементы (4,3) и (3,4) заменяются нулями.
Этот прием, называемый дефляцией (понижением порядка), численно устойчив, так как Тз претерпевает возмущение с нормой всего лишь 6.1 10 зз. Теперь цВ-итерация применяется снова, но уже только к ведущей главной 3 х 3-подматрице матрицы Тз. Этот процесс повторяется, пока не будут найдены все собственные значения (см. НОМЕРАОЕ/МаС!аЬ/Сг!д!ЯВ..ш). 5.3.2. В!4-итерация Вспомним, что, согласно нашему анализу в рззд. 4.4, на каждом шаге ЯВ; алгоритма неявно выполняется обратная итерация. Исследуем эту связь между двумя алгоритмами более подробно в том случае, когда в качестве сдвига в обратной итерации выбирается отношение Рэлея. Такой метод будем называть ВХг-итерацией (от Вау!е!8)з ЯпоС!епС 1СегаСюп).
Алгоритм 5.1. ВО„-итерацияз для заданного вектнора хо (!!хо~)з = 1) и задаваемого пользователем критерия астапова Со! вьтолнять итерации ро = р(хо,А) = .а.,4хвв з=0 5.3. Алгоритмы для симметричной проблемы собственных значений 227 гереа1 Уг=1А — Р; г1) 'х; з * =у*ЛЫЬ р, = р(х;,А) 1=4+1 пока метод не сойдется (т.
е. пока не будет вьгполнено условие ((Ахг — р,х,Йз < ьо1) Если условие останова удовлетворено, то по теореме 5.5 число р; удалено от некоторого собственного значения матрицы А не более чем на 1о!. Пусть в ЯВ;итерации используются сдвиги в, = а„„, а КЯ-итерация применяется к начальному вектору хо = [О,..., О, 11т. Обсуждавшуюся в разд. 4.4 связь между ЯВ;алгоритмом и обратной итерацией можно использовать, чтобы показать, что последовательности чисел ог и р„генерируемые этими двумя алгоритмами, совпадают 1см.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
