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

Дефляция До сих пор предполагалось, что все д, различны и все и, отличны от нуля. Если это не так, вековое уравнение у(Л) = О имеет (с вертикальных асимптот, где (с ( и, а потому й корней. Однако оказывается, что остальные и — )с собственных значений могут быть определены без каких-либо усилий: если д, = дс 1 или и, = О, то легко показать, что 4 является собственным значением и для матрицы В+ ринг (см.

вопрос 5.16). В такой ситуации мы говорим о дефляции. На практике выбирается некоторое пороговое значение и дефляция для числа 4 регистрируется, если в смысле этого порога 4 достаточно близкб к 4ч.1 либо и, достаточно мало. Оказывается, что в практических вычислениях дефляция происходит весьма часто. В экспериментах со случайными плотными матрицами, собственные значения которых равномерно распределены, дефляция имела место для более чем 15% собственных значений наибольшей матрицы Ю + риит.

Когда же собственные значения случайных плотных матриц сходились к нулю со скоростью геометрической прогрессии, дефляция наблюдалась для более чем 85% собственных значений! Важно научиться извлекать выгоду из этого обстоятельства, чтобы сделать алгоритм по-настоящему быстрым (59, 210). Основной выигрыш от использования дефляции состоит не в том, что убыстряется решение векового уравнения — этот этап в любом случае стоит лишь 0(п ) операций. Выигрыш заключается в ускорении матричного умножения на последнем шаге алгоритма. Действительно, если и; = О, то соответствующий собственный вектор есты-й столбец е; единичной матрицы (см. вопрос 5.16).

Это означает, что е; является г-м столбцом в матрице сг', поэтому при формировании матрицы Я посредством левого умножения на Я1 и Ят вычисление 1-го столбца не требует никаких затрат. Аналогичное упрощение имеет место в случае 4 = дсы. При дефляции многих собственных значений устраняется большая часть работы, связанной с матричным умножением.

Это отчетливо проявляется в численньгх экспериментах, обсуждаемых в равд. 5.3.6. Решение векового уравнения Предположим, что некоторое число и;, хотя и мало, все же недостаточно мало для того, чтобы была зарегистрирована дефляция. В этом случае применение метода Ньютона к решению векового уравнения встречается с затруднениями. Вспомним, что пересчет приближенного решения Л уравнения у'(Л) = О в методе Ньютона основан на следующих положениях: 1. Вблизи точки Л = Л функция у(Л) аппроксимируется линейной функцией 1(Л); график последней есть прямая линия, касающаяся графика функ- цииДЛ) приЛ=Л.

2. В качестве Л сы берется нуль этого линейного приближения, т. е. ((Л ~1) =О. Функция, показанная на рис. 5.2, не доставляет видимых трудностей методу Ньютона, поскольку вблизи каждого своего нуля у(Л) достаточно хорошо ап- 234 Глава 5. Симметричная проблема собственных значений 9 $99 1.996 2 2.009 2.01 0 2 4 Рис. 5.3. График функции 11Л) = 1+ е „+ зе + зе + ~~ проксимируется линейными функциями.

Однако рассмотрим график функции на рис. 5.3. Она получена из функции на рис. 5.2 заменой значения .5 для и~ на .001. Это новое значение недостаточно мало для того, чтобы вызвать дефляцию. График функции в левой части рис. 5.3 визуально не отличим от ее вертикальных и горизонтальных асимптот, поэтому в правой части укрупненно воспроизведен фрагмент графика, прилегающий к вертикальной асимптоте Л = 2. Видно, что график слишком быстро «выполняет поворотв и для большей части значений Л почти горизонтален. Поэтому, применяя метод Ньютона почти к любому начальному приближению Ло, мы получаем линейное приближение 1(Л) с почти горизонтальным графиком и малым положительным угловым коэффициентом.

В результате Лз является отрицательным числом, огромным по абсолютной величине, которое совершенно бесполезно в качестве приближения к истинному корню. Чтобы найти выход из этого положения, можно модифицировать метод Ньютона следующим образом: раз 1(Л) нельзя хорошо приблизить линейной функцией 11Л), попробуем взять в качестве приближения какую-нибудь другую простую функцию Ь(Л).

Нет ничего особого именно в прямых линиях: для метода Ньютона вместо ЦЛ) можно взять любое приближение Ь(Л), значения и нули которого легко вычисляются. Функция 1(Л) имеет полюсы в точках д; и 4~~, которые определяют ее поведение в соответствующих окрестностях. Поэтому при поиске корня в интервале (4+~,4) естественно выбрать функцию Ь1Л), также имеющую эти полюсы, т.

е. функцию вида Ач. — Л Константы сы сз и сз, обеспечивающие, что Ь(Л) есть приближение к ДЛ), можно выбрать несколькими способами. Мы изложим несколько упрощенный вариант способа, используемого ВАРАСК-программой в1ае04 ~172, 45], но вначале отметим, что если сз, сз и сз уже известны, то уравнение Ь(Л) = О легко решается относительно Л, поскольку сводится к эквивалентному квадратному уравнению с~ (сЦ-~ — Л) + сз(д~ — Л) + сз (А — Л) (4+~ — Л) = О.

5.3. Алгоритмы для симметричной проблемы собственных значений 235 Пусть Л вЂ” приближенное значение корня. Определим сы сз и сз так, чтобы П + + сз — — Ь(Л) 7'(Л) = 1+ р~~» 2-~-1 йь — л для Л в окрестности Л, Заметим, что 2 изь у(л) = 1 + р ~~' " + р '5 ' " = — 1 + ~,(л) + р,(л). а=1 " 2=2Ч-1 Если Л Е (4~.ы б2), то 2рг(Л) есть сумма положительных слагаемых, а 2рз(Л)— сумма отрицательных.

Поэтому и 2р1(Л), и 2р2(Л) могут быть вычислены с высокой точностью; однако при их сложении вполне вероятно взаимное уничтожение верных разрядов и потеря относительной точности в сумме. Возьмем числа с1 и сы такие, что функция С2 Ь1(Л) = с2 + удовлетворяет условиям Ь1(лу) = ф~(лб) и Ьг(Л1) = 2р~~(лб). (5.15) Это означает, что гипербола, являющаяся графиком функции Ь1(Л), касается графика функции 2р1(Л) при Л = Л . Два условия в (5.15) — это обычные условия метода Ньютона, за исключением того, что вместо прямой в качестве приближения используется гипербола.

Легко проверить, что с2 — — ф',(Лз)(г(2— Лб)2 и сг = 2р1(лб) — 2р', (Лб)(г1, — Л ) (см. вопрос 5.17). Подобным же образом выбираем сз и с2 так, чтобы функция С2 Ь2(л) — = с2 + удовлетворяла условиям 2+1 Ь2(лз) = 2р2(лу) и Ь',(Л,) = 2р2(Л,). (5.16) Наконец, полагаем Ь(Л) = 1 + Ь,(Л) + Ь2(Л) С2 С2 (1 + с1 + с2) + + 2+1 СГ С2 и; — л а;. — л Пример 5.8. Если для функции на рис. 5.3 начать с Ло — — 2.5, то 1 1111, 10-2 1 1111, 10-2 Ь(Л) — + + 1. График этой функции визуально не отличим от графика функции ДЛ), показанного в правой части рис. 5.3.

Решая уравнение Ь(Л1) = О, получаем Л1 = 2.0011. Это значение имеет 4 верных десятичных разряда. Продолжая таким же образом, находим Л2 с 11 верными разрядами и Лз, все 16 разрядов которого верны. О Алгоритм, реализованный ЬАРАСК-программой в1ае44, является небольшой модификацией только что описанного метода (который назван в [172] 236 Глава 5. Симметричная проблема собственных значений алгоритмом среднего пути).

В очень обширных численных тестах ЕАРАСК- программе требовалось в среднем две-три итерации для того, чтобы найти собственное значение с полной машинной точностью, и ни разу не понадобилось более семи итераций. Устойчивое вычисление собственных векторов Как только из векового уравнения найдены собственные значения сп матрицы Р+ риит, можно вычислить собственные векторы, пользуясь простой формулой (Р— а11) 'и из леммы 5.2.

К сожалению, вычисления по этой формуле могут быть неустойчивы (59, 90, 234); так будет, в частности, когда собственные значения а1 н аг+1 очень близки. Интуитивно проблема состоит в том, что выражения (Р— о11) 'и и (Р— с»1«11) 'и «очень похожи», в то время как требуется получить ортогональные собственные векторы. Более точно, если пг н пг+1 очень близки, оба близки к находящемуся между ними числу д«. Поэтому имеет место заметная потеря верных знаков при вычислении д1 — сг, или с(1 — сг, »1, либо при решении векового уравнения методом Ньютона.

В любом из этих случаев числа д1-с»1 и д1 — а«+1 могут содержать большие относительные ошибки, вследствие чего вычисленные собственные векторы (Р— с»,1) 'и и (Р— сгсь11) 'и очень неточны и далеки от ортогональности. Первоначальные попытки разрешить эту проблему (90, 234) предусматривали использование арифметики двойной точности при решении векового уравнения (в предположении, что входные данные имеют обычную точность), чтобы разности д, — аг и д1 — аг+1 могли быть вычислены с хорошей точностью. Однако, если уже входные данные имеют двойную точность, это означает, что требуется переход к учетверенной точности, недоступной для многих компьютеров и языков или, по меньшей мере, доступной лишь дорогой ценой.

Как указано в равд. 1.5, учетверенную точность можно моделировать посредством двойной точности (234, 204). Это можно сделать переносимым и относительно эффективным образом, если округления в основной арифметике с плавающей точкой выполняются достаточно точно. В частности, требуется, чтобы, исключая переполнения и машинные нули, выполнялось соотношение й(а а 5) = (а а Ь)(1+б), где ~б) = 0(е); см.