Содержание к диссертации
ВВЕДЕНИЕ 3
ГЛАВА I. Вывод формул алгоритмов исчерпывания трехдиаго-нальных симметрических и двухдиагональных матриц. Анализ погрешностей, связанных с неточностью решения уравнений относительно параметров, определяющих ортогональные преобразования вращения.
I. Алгоритм исчерпывания симметрической трехдиаго-
нальной матрицы 18
2. Алгоритм исчерпывания двухдиагональной матрицы.. 32 ГЛАВА 2. Решение уравнений для параметров, задающих ортогональные преобразования вращения, и их машинное вычисление.
3. Решение уравнений для параметров, определяющих
преобразования вращения 55
4. Машинное вычисление параметров, задающих преоб
разования вращения 66
ГЛАВА 3. Машинное вычисление преобразованных матриц.
5. Анализ погрешностей, возникающих при вычислении
элементов преобразованной трехдиагональной
матрицы 77
6. Машинное вычисление преобразованной двухдиаго
нальной матрицы 93
ГЛАВА 4. Общие схемы алгоритмов исчерпывания трехдиаго-
нальных симметрических и двухдиагональных матриц.
7. Общая схема исчерпывания симметрической трехдиа
гональной матрицы 103
8. Общая схема исчерпывания двухдиагональной
матрицы ИЗ
ПРИЛОЖЕНИЕ. Численный пример 126
ЛИТЕРАТУРА 135
Введение к работе
Работа посвящена выводу нового варианта формул и анализу погрешностей при реализации так называемых алгоритмов ортогонального исчерпывания трехдиагональннх симметрических и двухдиагональных матриц, с помощью которых достигается приведение этих матриц к диагональному виду. Если мы знаем собственное значение 2 симметрической трехдиагональной матрицы порядка Ш и соответствующий ему собственный вектор, то алгоритм исчерпывания позволяет с помощью подобных ортогональных преобразований вращения привести эту матрицу к клеточно-диагональной форме, в которой одна клетка является снова симметрической трехдиагональной матрицей порядка
У1Ь- і , а вторая клетка имеет порядок I и просто совпадает с } . Циклическое повторение процесса исчерпывания позволяет привести исходную матрицу к диагональной форме. Этот алгоритм хорошо известен и широко освещен в литературе, см., например, [і] Вариант для сингулярного разложения двухдиагональной матрицы описан в [2] .
При практическом применении известных реализаций алгоритмов исчерпывания иногда обнаруживается численная неустойчивость этих алгоритмов, состоящая в том, что окончательная форма преобразованной матрицы отличается от желаемой. В дис-вертации выясняются причины этой неустойчивости и приводятся способы их устранения. В результате разработаны численно устойчивые варианты алгоритмов исчерпывания трехдиагональннх симметрических и двухдиагональных матриц с гарантированной оценкой точности.
Продемонстрируем проблемы, возникающие при реализации алгоритмов исчерпывания, на примере симметрической трехдиа-тональной матрицы &; Si So &. о* А-
Vftu-i &№-{ От не имеющей нулевых элементов на побочных диагоналях, щФО
Предположим, что нам известны собственное значение 1 этой матрицы и отвечающий ему собственный вектор Ьі~(М'і, и*г}"ч LLnu) » так что имеют место равенства V/n, lltrv-i4- (cLtn.-})tlnb ~0 .
Традиционно алгоритм исчерпьшания состоит в следующем. По- а) лежим и, - U. и подберем ортогональные матрицы вращения (2 4= am) і Sc Cc 'ЧІ так, чтобы вектор Ца)^Си^(С~і} имел вид: (lcc>- '=(Qt0,~,0,1lftU*+i,~.,Utn,)T.
Другими словами, параметры SdCc (Sc+Cc^i) находятся как решение уравнения
При этом компонента Ц, оказывается равной
Сс-1) Одно из решений уравнения Cc&c-i "Si'U^O имеет вид: о № г- Ui
Тогда uf= l/JU^T^W и выражения для S^ , С і могут быть переписаны в виде: к* п(С) > ссs ш ' (2)
Отметим, что знаменатель в этих формулах не равен нулю, так как, очевидно, ІС^Ф 0 и, следовательно, и,і ~ = ]1и% + и1 + — + иУ 0- Поделим обе части і -го ( я /, lr.>,Wl-i ) из равенств (I) на Uc+i и последнего -на Ц(^ в результате простых преобразований, использующих формулы (2), эти равенства принимают вид:
С помощью этих равенств мы покажем, что матрица /f = * tn. ^CtruCnu-rC^CzACiC\ "CtnriCtn, имеет вид V*. &i Ох Vtnrl С^т^і Ь(п,-і unt-i dms-i 0 о X\ т.е. матрица А в результате применения к ней (tlt-{) -го подобных преобразований вращения приводится к клеточно-диа-гональной форме, в которой одномерная клетка совпадает с собственным значением /I
Обозначив &і~СіС[гі"С'ьСі&СгСі>'"Сс-іС , пока- жем сначала, что матрица fit і имеет вид
к h h і 4 ^ 44 d- і
От dm где для элементов & , ІІЇ имеют место формулы:
Для этого предположим, что матрица пі-і имеет такой же вид, если только значение индекса I заменить на 1-і , и вычислим матрицу 4с =Сс^С-{Сс Нетрудно проверить, что вычисленная по этой формуле матрица л имеет вид к 4
S. L L /г w% t/j
4 4 4 Щн 4 4 ^ ч-4
4vl*i ^W ^+і ирг
m (Л'гп где использованы обозначения: (6) k &Ъ = Л & ^-y+ьс* Ah Ct(№c-A+Ctdt).
Подставляя выражения для U>~t f Ьґс-і в формулы для %і'-і » ^^ » иГс- , находим:
Используя теперь равенства (3), видим, что %с-і~0 * а & , l&t имеют требуемые выражения:
Таким образом, предположение о том, что матрица ft і имеет вид (5), обосновано.
Рассмотрим теперь матрицу А С при L-Ль , т.е. матрицу /4^ =* /} . Чтобы убедиться в том, что она имеет вид (4), выпишем формулы для элементов tlnu , 1^7п* : &т> "'SttvlSm Cnu-i Vnt + Zftv (cLtru" A /J >
Выражение, заключенное в квадратные скобки, представляет собой в точности левую часть последнего из равенств (3). Следовательно, Ыщ.-О , 1ВДУ, «.. потение ^з-ложения где 2 ~ .диагональная матрица, содержащая собственные значения /4 , a tf* - ортогональная, состоит, таким образом, в выполнении (tn-{) -го шагов, і -й из которых представляет собой применение алгоритма исчерпывания к симметричес-кой трехдиагональной матрице порядка im-L + i ( L~ і ,
2,-.,m-l ).
Особо подчеркнем, что при получении вида матрицы Д мы существенно пользовались всеми равенствами (3). В то же время последнее из них является в некотором смысле лишним. Конкретнее, для определения (т- і) -й пар неизвестных ^ , С с достаточно первых (fti- {) -го из этих равенств, а последнее должно выполняться автоматически, если } точное собственное значение. Если же } мало отличается от точной величины собственного значения, то, казалось бы, это последнее уравнение должно быть почти выполненным. В действительности же бывают случаи, когда это не так, т.е. последнее уравнение имеет большую (по отношению к норме решения) невязку, хотя все остальные уравнения решены точно. По-видимому, первым обратил на это внимание Уилкинсон ( [з] , стр. 286). Один из наиболее показательных в этом отношении примеров рассмотрен в [4] , стр. 4.
На практике параметры $с , Сс подбираются именно из первых (tit-і) -го равенств (3). При этом из-за того, что последнее уравнение может иметь большую невязку, значения элементов ilm, , Ь/щ могут сильно отличаться от их точных значений, в частности, значение элемента U^m-может оказаться недостаточно малым для того, чтобы им можно было пренебречь. Таким образом, первая причина неустойчивости алгоритма исчерпывания связана с недооценкой роли последнего уравнения, которое обычно считается выполненным.
На первый взгляд эту причину можно устранить, вычислив собственный вектор так, чтобы хорошо удовлетворялись все уравнения (I) (а значит, и уравнения (3), так как они являются следствием (I)), включая последнее, и определяя параметры $с , Сс по формулам (2), т.е. непосредственно через компоненты собственного вектора. Однако простое преобразование этих формул к виду , ж !ЩШ' 'Ї^Ш показывает, что параметры $с » Сі определяются не столько самими компонентами собственного вектора, сколько их отношениями Не/^Н . В таком случае ясно, что даже чрезвычайно точное определение собственного вектора в метрике евклидова пространства не гарантирует еще правильности выбора параметров Sc » Сс » поскольку оно все равно может приводить к неправильным отношениям компонент этого вектора. Это происходит, например, в том случае, когда среди этих компонент наряду с большими есть очень малые. Таким образом, вторая причина неустойчивости связана с недостаточно точным определением отношений компонент собственного вектора.
Задача вычисления отношений компонент собственного вектора с высокой точностью решена в [4] . Мы существенно опираемся на результаты этой работы. В главе 2 диссертации показывается, что использование этих отношений позволяет добиться достаточно точного удовлетворения соотношений (3). При этом мы находим некоторые "идеальные" параметры Sc » Сс ( St + Cc~i ), которые удовлетворяют соотношениям (3) приближенно, т.е. с некоторыми возмущениями коэффициентов d-i , щ Приводится оценка этих возмущений. Мы называем параметры Sc » Сс идеальными по той причине, что они в точности удовлетворяют соотношению Sc -+Ci - { и используются только при выводе формул алгоритма, т.е. только в теоретических рассмотрениях. Реально же в памяти машины хранятся лишь некоторые их "машинные" приближения $'с ,
СІ , мало (в смысле относительной погрешности) отличающиеся от Si . Сс . Чтобы различие между параметрами Sc »
Сс и Sc » &' стало понятным, скажем, что Sc » С с определяются в результате точных вычислений по формулам (7), в которых вместо Сс'і используется с'с-ч » a si »
С'с представляют собой результаты машинной реализации этих же формул. Обеспечение высокой относительной точности при вычислении S'c » Сс важно по нескольким причинам. Во-первых, от величины погрешности параметров С'с зависит величина возмущений коэффициентов dc # щ » с которыми выполняются равенства (3) для параметров Sc , Сс . Во-вторых, параметры Sc » Сс используются при построении явных формул, задающих элементы (Ц , . преобразованной матрицы А (см. (4)). Понятно, что при реальных вычислениях в этих формулах можно использовать лишь приближения Sc » Сі . Поэтому от точности этих величин зависит точность определения преобразованной матрицы А . Наконец, важно обеспечить, чтобы приближенная матрица а,- а - si с' была близка к ортогональной, ведь результатами работы алго- ритма, кроме преобразованной матрицы, являются и матрицы, задающие требуемые преобразования вращения. Алгоритм вычисления величин S[ » С'с » обеспечивающий их высокую относительную точность, описан в 4 (гл. 2). 3 той же главы посвящен непосредственно решению уравнений (3) относительно параметров Sc * Сс '
В главе I на основе соотношений (3) выводятся явные формулы для элементов преобразованной матрицы fi . В традиционно используемом алгоритме исчерпывания формулы для этих элементов имеют, как нетрудно видеть из (6), рекуррентный характер, что, конечно, затрудняет анализ погрешностей, возникающих при их реализации на вычислительной машине. Мы же приводим окончательные формулы для этих элементов, так что сами преобразования вращения выполняются не в процессе вычислений, а аналитически. Поскольку уравнения (3) могут быть решены только приближенно, то преобразованная матрица уже не будет иметь клеточно-диагональной формы (4), а оказывается заполненной. Мы показываем, что эта заполненная матрица может быть представлена в виде суммы двух матриц, одна из которых имеет требуемый вид (4), а другая, называемая матрицей погрешностей, имеет малую норму. Выводу формул для элементов преобразованной матрицы и оценке нормы матрицы погрешностей посвящен I.
При вычислении элементов преобразованной матрицы по полученным формулам неизбежно возникновение погрешностей. Эти погрешности имеют двоякую природу. Во-первых, как уже отмечалось, вместо участвующих в этих формулах параметров Sc » Сс мы можем использовать только их машинные приближения Sj $ С'с Во-вторых, погрешности возникают и при непосредственном выполнении арифметических операций, предписывав- мых этими формулами. Оба эти источника погрешностей учитываются в 5 (гл. 3),
Проблемы, аналогичные описанным, встают и в алгоритме исчерпывания двухдиагональной матрицы А- &j Ьг О Q>N-{ ON і-
Коротко опишем этот алгоритм, предполагая элементы Сії , отличными от нуля. Пусть нам известно сингулярное число этой матрицы. Напомним, что сингулярными числами матрицы А называются N наибольших собственных значений составной матрицы S*
О А А* О
Алгоритм состоит в определении ортогональных матриц вращения і Сг- О і 4 Сс . d- о і так, чтобы матрица Ь*СыСы-ГСъЪ№№"СЛ имела вид >4- a,. 4-/
Циклическое повторение алгоритма исчерпывания позволяет получить сингулярное разложение матрицы А где 2 - диагональная матрица, содержащая сингулярные числа матрицы /t , а ф , Q - ортогональные матрицы,
Аналогично случаю трехдиагональноі матрицы нетрудно показать, что параметры Sc>Сс ,Sc , CL- (Sc+СІ-і > &г+ С*= і) должны удовлетворять равенствам ^C+iH^C^i Sc+iSrSz 6% $rh a* +CcJc+rO, > і = г,зг..,Мч, 0+ґ"Ьг - Sj-Sz X*-^.
Первое из этих равенств можно рассматривать как уравнение для определения параметров *>2 , и% , второе - как уравнение для Sz , Сг «Jo » W и т.д. При этом последнее равенство остается неиспользованным. Если б известно с некоторой, пусть даже достаточно малой, погрешностью, то это неиспользованное уравнение может иметь большую невязку. Это приводит к тому, что элемент матрицы Д , стоящий в позиции {Н~І7Н)) теоретически равный нулю, в практических расчетах оказывается недостаточно малым для того, чтобы им можно было пренебречь. Поэтому встает задача научиться достаточно точно решать уравнения (9). При решении этой задачи мы снова опираемся на возможность предварительного определения отношений компонент собственного вектора с высокой точностью (см. [4] ). В данном случае речь идет о собственном векторе трехдиагональной симметрической матрицы
S, т-
CL. N-i О і і отвечающем ее почти собственному значению б Отметим, что матрица Т получается из матрицы S с помощью одноименных перестановок строк и столбцов. Задачам определения параметров ^. , С^ , j , ^- , удовлетворяющих (приближенным) равенствам (9), и вычисления их достаточно точных машинных приближений s' , с'с » S/ » ?/ посвящены 3, 4 (гл. 2), где они решаются параллельно с аналогичными вопросами для случая трехдиагональной матрицы. Отметим, что мы смогли с достаточной точностью решить уравнения (9) только в случае, когда б является наибольшим сингулярним числом матрицы /f . Понятно, что это ограничение диктует такой порядок приведения матрицы А к диагональному виду, согласно которому исчерпывание каждой очередной матрицы выполняется с наибольшим для этой матрицы сингулярным числом.
Поскольку уравнения (9) могут быть решены только прибли женно, т.е. они выполняются с некоторыми возмущенными коэффи циентами & , -^- , то матрица fi уже не может иметь вида (8), а оказывается заполненной. Однако, как показано в диссертации, эта заполненная матрица может быть представлена в виде суммы двух матриц, одна из которых имеет требуемую форму (8), а вторая - матрица погрешностей - имеет малую нор му. Для элементов el; , А. первой матрицы получены прос тые явные формулы. Вывод этих формул и оценка нормы матрицы погрешностей составляют содержание 2 (гл. I). 6 (гл. 3) посвящен анализу погрешностей, возникающих при непосредствен- вой вычислении элементов cLc , Je по полученным форму- лам. В результате этого анализа находится оценка, показываю щая, насколько вычисленная (т.е. рассматриваемая уже как таб лица машинных чисел) преобразованная матрица отличается от некоторой матрицы, ортогонально эквивалентной исходной мат рице А
В 7 и 8 (гл. 4) приведены общие схемы алгоритмов исчерпывания соответственно: трехдиагональных симметрических и двухдиагоналъннх матриц.
Повторим кратко основные результаты диссертации.
Разработан новый вариант алгоритма исчерпывания симметрической трехдиагональной матрицы. Это исчерпывание может выполняться с произвольным собственным значением. Получена оценка погрешности, допускаемой при реализации этого алгоритма на вычислительной машине.
Разработан алгоритм исчерпывания двухдиагональнои матрицы. В случае, когда это исчерпывание выполняется с наибольшим сингулярным числом, получена гарантированная оценка точности.
Основные результаты диссертации полностью опубликованы в работах [б - 8] .
Автор выражает глубокую благодарность своему научному руководителю члену-корреспонденту АН СССР, профессору С.К.Го-дунову за постановку задачи и постоянное внимание к работе автора.
