§ 46. Ускорение QR-алгорифма

Одним из важнейших факторов ускорения -алгорифма является инвариантность процесса (45.1) к правой почти треугольной форме матрицы. В самом деле, пусть А правая почти треугольная. Тогда можно считать, что матрица есть произведение матриц вращения. Но в этом случае матрица будет снова правой почти треугольной. Конечно, такой же вид будут иметь все матрицы

Если матрица А полная и не обладает какой-либо спецификой, то реализация одного шага процесса (45.1) требует выполнения около арифметических операций. Если же А правая почти треугольная, то на одном шаге нужно выполнить только операций. Поэтому каждый шаг процесса (45.1) будет осуществляться примерно в раз быстрее, если предварительно матрицу А привести подобным преобразованием к правой почти треугольной форме. Такое преобразование было описано и исследовано в § 32.

Предположим, что элемент правой почти треугольной матрицы А равен нулю. Разобьем А на такие четыре клетки, чтобы диагональные клетки были квадратными и клетка в левом верхнем углу имела порядок Но тогда клетка в нижнем левом углу будет нулевой,

Процесс (45.1) инвариантен в этой форме, причем на всех его шагах клетка в нижнем правом углу

преобразуется независимо от остальных клеток. Следовательно, нет особого смысла в применении -алгорифма непосредственно к таким матрицам.

Если мы интересуемся только собственными значениями матрицы то достаточно определить собственные значения клеток Даже при определении собственных векторов матрицы (46.1) целесообразно сначала определить собственные значения клеток Поэтому каждый раз, когда какие-нибудь из поддиагональных элементов правой почти треугольной матрицы по тем или иным причинам обратятся в нули, мы будем продолжать применение -алгорифма только к соответствующим матрицам меньших размеров.

Это, конечно, затрудняет использование -алгорифма для определения собственных векторов. Однако мы покажем в дальнейшем, что, вычислив собственные значения, можно весьма эффективно получить собственные векторы.

Всюду в дальнейшем мы будем считать, что -алгорифм применяется к правым почти треугольным матрицам с ненулевыми поддиагональными элементами.

Если какое-нибудь собственное значение йгакой матрицы является кратным, то оно должно входить только в один канонический ящик Жордана. В самом деле, пусть собственному значению К соответствует более одного ящика Жордана матрицы Тогда ранг матрицы должен быть не больше Но он заведомо не меньше так как поддиагональные элементы отличны от нуля и, следовательно, отличен от нуля минор, расположенный в первых столбцах и последних строках. Полученное противоречие подтверждает справедливость высказанного свойства.

Из этого свойства, в частности, вытекает, что если правая почти треугольная матрица с ненулевыми поддиагональными элементами имеет простую структуру, то все ее собственные значения различны. Как показывает пример матрицы (44.20), мы не можем надеяться на то, что наличие очень близких собственных значений обязательно приведет к появлению очень малых поддиагональных элементов.

Рассмотрим теперь любую последовательность чисел Снова построим ортогональные матрицы и

правые треугольные матрицы по рекуррентным формулам

Если А правая почти треугольная матрица, то аналогичный вид будут иметь все матрицы Как и раньше, легко показать, что матрицы ортогонально подобны исходной матрице при этом

Предположив, что при некоторой, вообще говоря, зависящей от нумерации собственных значений выполняются неравенства

получим, что поддиагональные элементы матрицы суть величины такого порядка:

Числа в (46.2) называются сдвигами. Выбирая их подходящим образом, можно в значительной мере ускорить сходимость -алгорифма.

По существу, все известные стратегии выбора сдвигов состоят из двух различных этапов. На первом этапе каким-либо способом обеспечивается ваметное уменьшение одного из поддиагональных элементов матрицы Теоретическое обоснование этого этапа для матриц общей структуры, как правило, отсутствует, но проводится экспериментальное подтверждение его эффективности. Малость поддйагонального элемента позволяет начать на втором этапе целенаправленный выбор сдвигов. При этом обеспечивается дальнейшее уменьшение поддйагонального элемента с существенно большей скоростью.

Пусть элемент не превосходит по модулю малого числа Покажем, как выбрать теперь сдвиги, чтобы скорость убывания элементов стала не менее, чем квадратичной.

Разобьем матрицу на клетки таких же размеров, как в матрице (46.1). Если обозначить это разбиение через

то диагональные клетки являются правыми почти треугольными, а клетка имеет лишь один ненулевой элемент порядка в верхнем правом углу.

Вычислим каким-либо способом, например, по формуле (26.5), характеристический многочлен матрицы Принимая во внимание неравенства (46.4) при заключаем, что согласно лемме 11.2 величины будут иметь порядок по крайней мере при и не будут малыми при Обозначим через корни многочлена и выполним дополнительно шагов процесса (46.2), взяв в качестве сдвигов. Рассмотрим соответствующие неравенства (46.4) при Для малых совокупности собственных значений при будут одинаковыми.

Но в этом случае

Циклически повторяя описанный выбор сдвигов через каждые шагов процесса (46.2), мы обеспечим не менее, чем квадратичную скорость убывания поддиагональных элементов Как только элемент в позиции станет достаточно малым, его можно заменить нулем, не потеряв при этом существенно в точности, и продолжать применение -алгорифма к матрицам меньших размеров.

Практическая реализация процесса ускорения в соответствии с изложенной схемой не всегда оказывается эффективной. Предположим, что вещественная матрица А имеет комплексные собственные значения. В этом случае среди корней вещественного многочлена могут оказаться комплексно сопряженные пары. Но тогда некоторые из промежуточных матриц будут комплексными, несмотря на то, что матрица является вещественной. Появление комплексных матриц весьма нежелательно как с точки зрения времени счета, так и с точки зрения использования памяти ЭВМ. Поэтому покажем сейчас, как можно вычислять матрицу минуя получение промежуточных матриц Для этого нам потребуется Лемма 46.1. Пусть для матрицы А выполнены унитарно подобные преобразования

причем матрицы правые почти треугольные с ненулевыми поддиагональными элементами. Если первые столбцы матриц совпадают, то существует такая диагональная матрица с элементами, равными по модулю единице, что

Доказательство. По существу, нам достаточно показать, что если унитарная матрица и правая почти треугольная матрица С с ненулевыми поддиагональными элементами связаны соотношением

то при заданном первом столбце матрицы 7, с определяются в основном однозначно.

Обозначим через столбцы матрицы через элементы матрицы С. Приравнивая первые столбцы левой и правой частей в (46.6), имеем

Так как унитарна, то

где — произвольное комплексное число, по модулю равное единице. Приравнивая вторые столбцы в (46.6), получим

откуда вытекает, что

Здесь снова произвольное комплексное число, по модулю равное единице.

Продолжая процесс, устанавливаем, что все столбцы матрицы начиная со второго, определяются однозначно с точностью до умножения на комплексные числа, по модулю равные единице. и подтверждает справедливость соотношений (46.5).

Теперь перейдем к обоснованию процесса вычисления матрицы Ясно, что

для некоторой ортогональной матрицы и правой треугольной матрицы

Предположим, что мы каким-либо способом нашли такую ортогональную матрицу которой первый столбец совпадает с первым столбцом и при этом матрица является правой почти треугольной. Возможны две ситуации. Если все поддиагональные элементы матрицы С отличны от нуля, то согласно лемме 46.1

где -диагональная матрица с элементами, равными по модулю единице. Модули соответствующих элементов матриц и с одинаковы, поэтому безразлично, с какой из них продолжать -алгорифм. Мы будем продолжать его с матрицей С. Если же какие-то поддиагональные элементы матрицы с равны нулю, то это более благоприятный случай, так как можно продолжать -алгорифм с матрицами меньших размеров.

Теоретически матрицу можно получить как произведение матриц вращения при исключении поддиагональных элементов например, по столбцам сверху вниз. Будем считать, что

В силу строения матриц вращения первая строка произведения совпадает с первой строкой произведения следовательно, с первой строкой матрицы Но матрицы определяются лишь первым столбцом которого могут быть отличны от нуля только первые элементов. Поэтому в действительности первый столбец матрицы совпадает с первым столбцом произведения

Вычислим первый столбец матрицы и определим по нему соответствующие матрицы вращения Получим, далее, матрицу

и приведем ее с помощью алгорифма, описанного в § 32, к ортогонально подобной правой почти треугольной матрице. Это и будет нужная нам матрица С.

Действительно, новая почти треугольная матрица получена как произведение вида Согласно построению ортогональная матрица имеет первые столбец и строку, совпадающие с первыми столбцом и строкой единичной матрицы. Следовательно, первый столбец матрицы будет совпадать с первым столбцом матрицы Последняя матрица не только ортогональная, но и является матрицей подобного преобразования, приводящего к правой почти треугольной, т. е. обладает свойствами матрицы

Матрица в в (46.7) имеет много нулевых элементов и отличается от правой почти треугольной тем, что почти

все элементы главного минора порядка могут быть отличны от нуля. Специальный вид матрицы В легко учесть при ее приведении к почти треугольной форме. На всех этапах приведения промежуточные матрицы будут отличаться от правой почти треугольной тем, что почти все элементы некоторого минора порядка опирающегося на главную диагональ, могут быть отличны от нуля. По мере выполнения процесса приведения этог минор будет перемещаться по диагонали вниз.

В целом прямое получение матрицы из матрицы требует выполнения примерно такого же объема вычислений, как и последовательное ее получение за шагов процесса (46.2) с вещественными сдвигами. При этом не нужно находить корни многочленов. Однако объем дополнительной памяти ЭВМ, необходимой для хранения результатов промежуточных вычислений, быстро растет с уменьшением

Описанное ускорение сходимости -алгорифма особенно эффективно, когда . В этом случае очередной сдвиг совпадает с последним диагональным элементом матрицы и заведомо будет вещественным. Поэтому прямое получение матрицы не является обязательным, а квадратичное убывание поддиагонального элемента будет происходить на каждом шаге процесса (46.2). Если же вещественная матрица А имеет комплексные собственные значения, прямое получение матрицы может оказаться необходимым при условии, конечно, что мы хотим выполнять лишь вещественные вычисления. В этом случае прямое получение матрицы наиболее эффективно при Обязательное появление малых поддиагональных элементов с меньшим значением связано в основном с наличием у матрицы А канонических ящиков Жордана больших размеров.

Подводя итоги, перечислим еще раз рассмотренные приемы ускорения времени счета при реализации -алгорифма.

1. Предварительное приведение матрицы А к правой почти треугольной форме. Без этого приведения -алгорифм обычно не применяется.

2. Использование сдвигов для повышения скорости убывания поддиагональных элементов. При наличии комплексных собственных значений наиболее эффективно прямое получение матриц из (46.2).

3. Замена малых поддиагональных элементов нулями. позволяет продолжать применение -алгорифма к матрицам меньших размеров.

В приемах ускорения остается не ясным только один вопрос. Именно, как начинать выбор сдвигов, чтобы достаточно быстро получить малый поддиагональный элемент в позиции с наибольшим по возможности значением За исключением некоторых специальных классов матриц, пока на него не получено убедительного ответа. Стратегии выбора сдвигов, гарантирующие убывание поддиагональных элементов, обычно оказываются слишком медленными.

Мы опишем сейчас одну из практических процедур выбора сдвигов. Она эффективна, хотя ее применение связано с некоторым риском.

Процесс начинается с вычисления матрицы при Предположим, что уже получены матрицы Проверяем выполнение неравенства

Если оно справедливо, то находим матрицу беря Если же это неравенство не имеет места, то. проверяем выполнение другого неравенства:

где — клетки второго порядка матриц находящиеся в нижнем правом углу. Если это неравенство справедливо, то вычисляем характеристический многочлен матрицы и находим матрицу используя прямой способ ее получения. Если же и последнее неравенство не имеет места, то находий матрйцу беря Конечно, при первой же возможности заменяем малые поддиагональные элементы нулями и продолжаем применение -алгорифма к матрицам меньших размеров.

Применение этой процедуры показало, что среднее число итераций на каждое собственное значение, как правило, не превосходит 5. Процесс исключительно устойчив. Анализ ошибок округления в нем полностью охватывается тем анализом, который был выполнен ранее. Если

тать, что на каждое собственное значение матрицы А требуется не более пяти итераций, то вычисленные собственные значения будут точными для некоторой возмущенной матрицы причем

УПРАЖНЕНИЯ

(см. скан)