Ортогональная диагонализация двухдиагональных и якобиевых матриц с гарантированной оценкой точности тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

ВВЕДЕНИЕ.

ГЛАВА I. Вывод формул алгоритмов исчерпывания трехдиаго-нальных симметрических и двухдиагональных матриц. Анализ погрешностей, связанных с неточностью решения уравнений относительно параметров, определяющих ортогональные преобразования вращения.

§ I. Алгоритм исчерпывания симметрической трехдиагональной матрицы. •

§ 2. Алгоритм исчерпывания двухдиагональной матрицы.

ГЛАВА 2. Решение уравнений для параметров, задающих ортогональные преобразования вращения, и их машинное вычисление.

§ 3. Решение уравнений для параметров, определяющих цреобразования вращения.

§ 4. Машинное вычисление параметров, задающих преобразования вращения.

ГЛАВА 3. Машинное вычисление цреобразованных матриц.

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

§ 6. Машинное вычисление преобразованной двухдиагональной матрицы

ГЛАВА 4. Общие схемы алгоритмов исчерпывания трехдиагональных симметрических и двухдиагональных матриц.

§ 7. Общая схема исчерпывания симметрической трехдиагональной матрицы

§ 8. Общая схема исчерпывания двухдиагональной матрицы.ИЗ

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

УН- \ , а вторая клетка имеет порядок I и просто совпадает с } . Циклическое повторение процесса исчерпывания позволяет цривести исходную матрицу к диагональной форме. Этот алгоритм хорошо известен и широко освещен в литературе, см., например, [I] • Вариант для сингулярного разложения двухдиагональной матрицы описан в [2] .

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

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

4= с1{ к

1г к

I & г &1+0

Чти ^пги не имеющей нулевых элементов на побочных диагоналях,

Предположим, что нам известны собственное значение 1 этой матрицы и отвечающий ему собственный вектор пи) Т » так что имеют место равенства б^Н ^Сн> (I) V

Традиционно алгоритм исчерпывания состоит в следующем. Полу ложим и, - СС и подберем ортогональные матрицы вращения 2 к и к № )

Сг 1 1 о -л

Зс Сс так, чтобы вектор Сс имел вид:

Ц,СС)=(0,0,^,0 , и?, LLc+i,., Um)T •

Другими словами, параметры $с,Сс (+ t) находятся как решение уравнения tcU^-StUfO.

С)

При этом компонента Ц^с оказывается равной u£j?= Si uf^-t Сс U,i ■

Сс-1)

Одно из решений уравнения Cc^c-i 'Sc^'O имеет вид: s U£f LLC

Тогда ¡/[¡^rW и выражения для , Сс могут быть переписаны в виде: и (С-1) tt-C-i h= пШ ' ¿С= пш ' Hi и- С

Отметим, что знаменатель в этих формулах не равен нулю, так как, очевидно, Uf^O и, следовательно, и}^ -= |ju% + ui + — + uy 76 0- Поделим обе части I -го ( L « /, из равенств (I) на и последнего на ¿¿^ • в результате простых преобразований, использующих формулы (2), эти равенства принимают вид: hldrD + C-J^O, hhk+SbCzii^D+Cbtb-O, (3)

SmCnt-i * C-^Utn."}) - 0 .

С помощью этих равенств мы покажем, что матрица $ = = Ст,Спи-{"'С3С?ппС^ имеет вид а к. Ь I

4 ¿4 иПг1 с1щ--1 ¿тп с1/ О о Л. т.е. матрица А в результате применения к ней (т- {) -го подобных преобразований вращения приводится к клеточно-диа-гональной форме, в которой одномерная клетка совпадает с собственным значением X •

Обозначив С1'~ Сс^ С* жем сначала, что матрица Ас имеет ввд покак

1 4 к ¿г

4* и>с с Щ

4 т с1>т где для элементов , Ыс имеют место формулы:

5)

Для этого предположим, что матрица А ¿-I имеет такой же вид, если только значение индекса I заменить на 1*1 , и вычислим матрицу Ас = С с А С-{ С/ • Нетрудно проверить, что вычисленная по этой формуле матрица А^ имеет вид

4 4 4 4» 4 к

4 4

4 (¿с

Не т ¿Ли

4 ¿4/ я. где использованы обозначения: — н = Сс (Сс т-1 -*с(СеЬн&с ,

6)

Щ = Зс ($с иг^+ЬЬ-ЛЬ Ы^с-А^сси).

Подставляя выражения для ¿¿¿--^ } №¿-1 в формулы для % С-г , Не , иГс , находим:

Используя теперь равенства (3), видим, что 1 а Ъ&с имеют требуемые выражения:

Таким образом, предположение о том, что матрица АI имеет вид (5), обосновано.

Рассмотрим теперь матрицу А С ИРИ (--Кь , т.е. матрицу Ат- ~ А • Чтобы убедиться в том, что она имеет вид (4), выпишем формулы для элементов Ипи , 'ЬОЪь. : и-М1, пи " С т. [ЗщС-тг{

Выражение, заключенное в квадратные скобки, представляет собой в точности левую часть последнего из равенств (3). Следовательно, Ипь~ 0 , 10}^- 2 и матрица Д действительно имеет вид (4). Переход от матрицы А к подобной ей матрице /( вида (4) и составляет сущность алгоритма исчерпывания. Процесс приведения симметрической трехдиагональ-ной матрицы 4 к диагональному ввду, т.е. потение разложения где 2 ~ .диагональная матрица, содержащая собственные значения /4 , а ^ - ортогональная, состоит, таким образом, в выполнения (т-1) -го шагов, I -й из которых представляет собой применение алгоритма исчерпывания к симметричесI кой трехдиагональной матрице порядка т-С + { ( й ,

Особо подчеркнем, что при получении вида матрицы Д мы существенно пользовались всеми равенствами (3). В то же время последнее из них является в некотором смысле лишним. Конкретнее, для определения (¡71- {) -й пар неизвестных , С с достаточно первых (т- {) -го из этих равенств, а последнее должно выполняться автоматически, если ) точное собственное значение. Если же } мало отличается от точной величины собственного значения, то, казалось бы, это последнее уравнение должно быть почти выполненным. В действительности же бывают случаи, когда это не так, т.е. последнее уравнение имеет большую (по отношению к норме решения) невязку, хотя все остальные уравнения решены точно. По-видимому, первым обратил на это внимание Уилкинсон ( [з] , стр. 286). Один из наиболее показательных в этом отношении примеров рассмотрен в [4] , стр. 4.

На практике параметры $с , Сс подбираются именно из первых -го равенств (3). При этом из-за того, что последнее уравнение может иметь большую невязку, значения элементов ¿¿/^ , Ь/гп могут сильно отличаться от их точных значений, в частности, значение элемента и,пи может оказаться недостаточно малым для того, чтобы им можно было пренебречь. Таким образом, первая причина неустойчивости алгоритма исчерпывания связана с недооценкой роли последнего уравнения, которое обычно считается выполненным.

На первый взгляд эту причину можно устранить, вычислив собственный вектор так, чтобы хорошо удовлетворялись все уравнения (I) (а значит, и уравнения (3), так как они являются следствием (I)), включая последнее, и определял параметры в с » Ос по формулам (2), т.е. непосредственно через компоненты собственного вектора. Однако простое преобразование этих формул к виду показывает, что параметры $£ , Сс определяются не столько самими компонентами собственного вектора, сколько их отношениями Не/. В таком случае ясно, что даже чрезвычайно точное определение собственного вектора в метрике евклидова пространства не гарантирует еще правильности выбора параметров ¿с , С-с , поскольку оно все равно может приводить к неправильным отношениям компонент этого вектора. Это происходит, например, в том случае, когда среди этих компонент наряду с большими есть очень малые. Таким образом, вторая причина неустойчивости связана с недостаточно точным определением отношений компонент собственного вектора.

Задача вычисления отношений компонент собственного вектора с высокой точностью решена в [4] . Мы существенно опираемся на результаты этой работы. В главе 2 диссертации показывается, что использование этих отношений позволяет добиться достаточно точного удовлетворения соотношений (3). При этом мы находим некоторые "идеальные" параметры ,

Сс ( $с + Сс~1 ), которые удовлетворяют соотношениям (3) приближенно, т.е. с некоторыми возмущениями коэффициентов ¿4* , 4с • Приводится оценка этих возмущений. Мы называем параметры , Сс идеальными по той причине, что они в точности удовлетворяют соотношению $1 + Сс - { л используются только при выводе формул алгоритма, т.е. только в теоретических рассмотрениях. Реально же в памяти машины хранятся лишь некоторые их "машинные" приближения $1 , , мало (в смысле относительной погрешности) отличающиеся от , Сс . Чтобы различие меящу параметрами £с »

Сс и , Сс стало понятным, скажем, что $с , С с определяются в результате точных вычислений по формулам (7), в которых вместо Сеч используется С'с-ч , а ^'• ,

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

Во-первых, от величины погрешности параметров С'с зависит величина возмущений коэффициентов (¿с * ¿с » с которыми выполняются равенства (3) для параметров ¿с , С с . Во-вторых, параметры ¿с , Сс используются при построении явных формул, задающих элементы , преобразованной матрицы А (см. (4)). Понятно, что при реальных вычислениях в этих формулах можно использовать лишь приближения ¿с » С'с . Поэтому от точности этих величин зависит точность определения преобразованной матрицы А . Наконец, важно обеспечить, чтобы приближенная матрица о о

С'с -£

4' а ч была близка к ортогональной, ведь результатами работы алгоритма, кроме преобразованной матрицы, являются и матрицы, задающие требуемые преобразования вращения. Алгоритм вычисления величин йс » » обеспечивающий их высокую относительную точность, описан в § 4 (гл. 2). § 3 той же главы посвящен непосредственно решению уравнений (3) относительно параметров Сс •

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

Цри вычислении элементов преобразованной матрицы по полученным формулам неизбежно возникновение погрешностей. Эти погрешности имеют двоякую природу. Во-первых, как уже отмечалось, вместо участвующих в этих формулах параметров ^с » Сс мы можем использовать только их машинные приближения , С'с • Во-вторых, погрешности возникают и при непосредственном выполнении арифметических операций, предписывавмых этими формулами. Оба эти источника погрешностей учитываются в § 5 (гл. 3),

Проблемы, аналогичные описанным, встают и в алгоритме исчерпывания двухдиагональной матрицы А а* о о

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

О А

А* О

Алгоритм состоит в определении ортогональных матриц вращения I

Сг1 О

Сс с Сс I

С, о к к 1 так, чтобы матрица ^С^-^Щ^-С^С, имела вид 1 1 4я* к ч "г

N-1

8)

4-/ о ±б

Циклическое повторение алгоритма исчерпывания позволяет получить сингулярное разложение матрицы А

-fZQ*, где 2 - диагональная матрица, содержащая сингулярные числа матрицы /4 , а ф , @ - ортогональные матрицы,

1.

Аналогично случаю трехдиагональнои матрицы нетрудно показать, что параметры Р Сс , $с , Сс (Яс+С1 - { , ¿¡¿г+ У') должны удовлетворять равенствам 1 =

9) $£+4Сс^ ? ^с+Ас-и ~О Сыв* С" Оч

Первое из этих равенств можно рассматривать как уравнение для определения параметров £2 » С г » второе - как уравнение для , Сг и т.д. При этом последнее равенство остается неиспользованным. Если б известно с некоторой, пусть даже достаточно малой, погрешностью, то это неиспользованное уравнение может иметь большую невязку. Это приводит к тому, что элемент матрицы /) , стоящий в позиции теоретически равный нулю, в практических расчетах оказывается недостаточно малым для того, чтобы им можно было пренебречь. Поэтому встает задача научиться достаточно точно решать уравнения (9). При решении этой задачи мы снова опираемся на возможность предварительного определения отношений компонент собственного вектора с высокой точностью (см. [4] ). В данном случае речь идет о собственном векторе трехдиагональной симметрической матрицы о

1 0 4 4 о аг аг о ^

• •

V/ о 4 6« о а» ан о отвечающем ее почти собственному значению б • Отметим, что матрица Т получается из матрицы £ с помощью одноименных перестановок строк и столбцов. Задачам определения параметров ^ , , » » удовлетворяющих (приближенным) равенствам (9), и вычисления их достаточно точных машинных приближений > С'с , // , посвящены

§§ 3, 4 (гл. 2), где они решаются параллельно с аналогичными вопросами для случая трехдиагональной матрицы. Отметим, что мы смогли с достаточной точностью решить уравнения (9) только в случае, когда б является наибольшим сингулярынм числом матрицы /4 . Понятно, что это ограничение диктует такой порядок приведения матрицы к диагональному виду, согласно которому исчерпывание каждой очередной матрицы выполняется с наибольшим для этой матрицы сингулярным числом.

Поскольку уравнения (9) могут быть решены только приближенно, т.е. они выполняются с некоторыми возмущенными коэффициентами , , то матрица уже не может иметь вида (8), а оказывается заполненной. Однако, как показано в диссертации, эта заполненная матрица может быть представлена в виде суммы двух матриц, одна из которых имеет требуемую форму (8), а вторая - матрица погрешностей - имеет малую норму. Для элементов ¿^ , первой матрицы получены простые явные формулы. Вывод этих формул и оценка нормы матрицы погрешностей составляют содержание § 2 (гл. I). § 6 (гл. 3) посвящен анализу погрешностей, возникающих при непосредственном вычислении элементов ас , 1с по получениям формулам. В результате этого анализа находится оценка, показывающая, насколько вычисленная (т.е. рассматриваемая уже как таблица машинных чисел) преобразованная матрица отличается от некоторой матрицы, ортогонально эквивалентной исходной матрице А .

В §§ 7 и 8 (гл. 4) приведены общие схемы алгоритмов исчерпывания соответственно: трехдиагональных симметрических и двухдиагональных матриц.

Повторим 1фатко основные результаты диссертации.

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

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

Основные результаты диссертации полностью опубликованы в работах [б - 8] .

Автор выражает глубокую благодарность своему научному руководителю члену-корреспонденту АН СССР, профессору С.К.Годунову за постановку задачи и постоянное внимание к работе автора.

1. Golub G., Kaban W# Calculating the singular values and pseudoinverse of a matrix. J.SIAM Numer. Anal., Ser B, V.2, No 2, 1965, pp.205-224.

2. Уилкинсон Дж.Х. Алгебраическая проблема собственных значений. Москва, Наука, 1970, 564 с.

3. Годунов К., Костин В.И., Митченко А.Д. Вычисление собственного вектора симметрической трехдиагональной матрицы. Новосибирск, Б.и., 1983, 58 с. (Препринт/ИМ СО АН СССР, 44).

4. Годунов К. Решение систем линейных уравнений. Новосибирск, Наука, 1980, 177 с.

5. Митченко А.Д, Алгоритмы исчерпывания трездщагональных симметрических и двухдиагональных матриц. Новосибирск, Б.И., 1984, 44 с. (Препринт/ИМ СО АН СССР, J 59).

6. Митченко А.Д. Учет вычислительных погрешностей в алгоритме исчерпывания симметрической трех,диагональной матрицы. Новосибирск, Б.и., 1984, 44 с. (Препринт/ИМ СО АН СССР, 60).

7. Митченко А.Д. Учет вычислительных погрешностей в алгоритме исчерпывания двухдиагональной матрицы. Новосибирск, Б.И., 1984, 46 с. (Прецринт/ИМ СО АН СССР, 61).

8. Воеводин Б.В. Вычислительные основы линейной алгебры. Москва, Наука, 1977, 304 с.