Для заданных матрицы А, начального вектора x 0

x 0

1 , величины

2

сдвига s и точности ε положить k 0 и выполнить:

1.

Вычислить вектор y k 1 A sE 1 x k .

y k 1

2.

Вычислить приближенный собственный вектор x k 1

.

y k 1

2

~ k 1

x

3.

k 1 T

k 1

.

Определить приближенное собственное значение

Ax

4.

Если для векторов x k

и x k 1 выполняется условие

x k 1

x k

,

то сходимость

достигнута; в противном случае

положить

k k 1 и перейти к п.1.

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

5.8. Алгоритмы для симметричной проблемы собственных значений

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

5.8.1. Метод Якоби

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

чем другие алгоритмы, имея трудоемкость O Cn3 , где С – достаточно большая константа. Однако интерес к методу до сих пор сохраняется, т.к. он способен вычислять малые собственные значения и соответствующие им собственные векторы на много точнее, чем конкурирующие с ним алгоритмы. Кроме того, он легко распараллеливается.

Метод Якоби не предполагает предварительного приведения матрицы к трехдиагональной форме, а работает с исходной плотной матрицей.

По заданной симметричной матрице A A 0 в методе Якоби строится последовательность ортогонально подобных матриц A 1 , A 2 , которая

сходится к диагональной матрице, на диагонали которой стоят собственные
значения. Матрица A k 1	получается из A k по формуле A k 1 J k T A k J k ,
J k – ортогональная	матрица, называемая вращением Якоби. Таким
образом,

A k J k 1 T A k 1 J k 1 J k 1 T J k 2 T A k 2 J k 2 J k 1J k 1 T J 0 T A 0 J k 1 J 0 J T AJ

При подходящем выборе вращений матрица			A k для больших k будет
близкой к	диагональной матрице Λ.	Следовательно, можно записать,
J T A J	или J J T A. Поэтому	столбцы	матрицы J являются

приближенными собственными векторами для А.
	Матрица J k выбирается		так, чтобы		сделать нулевыми пару
внедиагональных элементов матрицы A k 1 J k T A k J k .
			i	j
	1
					0
			1
	J k R i, j,		cos	sin		i
			sin	cos		j
	0				1
					1
	Пусть θ – угол вращения. Он выбирает так, чтобы элементы i, j и j,i в
A k 1	стали нулевыми. Тогда, если обозначить c cos , s sin ,							t	s	tg ,

c

ctg 2 , можно получить следующий алгоритм для вращения Якоби в координатной плоскости i,j.

Если aij не слишком мал, т.е. больше погрешности вычисления:

1.	Вычислить		aii			a jj
					2aij
2.	Вычислить t				sign
					1 2
3.	Вычислить c					1
			1			t 2

4.Вычислить s c t

5.Применить вращение A R i, j, T A R i, j, .

6.Если нужны собственные векторы, J J R i, j, .

Стоимость применения вращения R i, j, к матрице А (или J) составляет O n операций, т.к. в А изменяются строки и столбцы с номерами i, j (а в J –

только столбцы i, j ). Тогда классический алгоритм Якоби можно записать:

Пока норма наддиагональной части матрицы А off A
						aij2	:
						1 i k n
1.	Выбрать i, j так, чтобы		aij		был наибольшим	по абсолютной
	величине внедиагональным элементом.
2.	Выполнить вращение Якоби в координатной плоскости i, j .

Если N	n n 1	– число наддиагональных элементов матрицы А, то
2

классический алгоритм Якоби сходиться, по меньшей мере, со скоростью

геометрической прогрессии, знаменатель которой не превосходит 1		1	.
		N

Асимптотически метод сходится квадратично.

Классический алгоритм Якоби не используется в практических

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

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

циклический метод Якоби – построчный циклический обход внедиагональных элементов матрицы А:

Пока норма наддиагональной части матрицы А off A
	aij2	:
	1 i k n
for i 1 to n 1
for j i 1 to n
Выполнить вращение для А в плоскости i, j

Если в течение всего внутреннего цикла процедура вращения дает только значения c 1, s 0 , то матрица А уже больше не изменяется. Циклический алгоритм Якоби также асимптотически сходится квадратично.

5.8.2. Трехдиагональная QR-итерация

Позволяет найти все собственные значения и, если необходимо, собственные векторы симметричной трехдиагональной матрицы. Нахождение всех собственных значений обходится за O n2 операций, что при достаточно

больших n составляет малую часть от 43 n3 операций для приведения исходной

плотной матрицы к трехдиагональному виду. Если необходимы и все собственные векторы, то QR-итерация в среднем требует больше, чем 6n3

операций, и остается самым быстрым алгоритмом только для матриц малого порядка ( n 25 ).

5.8.3. RQ-итерация

Этот алгоритм (Rayleigh Quotient Iteration) не предполагает первоначального приведения матрицы к трехдиагональному виду. Здесь в качестве сдвига в обратной итерации используется отношение Рэлея.

Алгорит м:

x 0

1 , величины

Для заданных матрицы А, начального вектора

x 0

точности ε положить k 0 и выполнить:

2

x 0 T Ax 0 .

1.

Вычислить отношение Рэлея 0

x 0 , A

x 0 T x 0

2.

Положить k 1.

3.

Вычислить вектор y k A k 1 E 1 x k 1 .

4.

Вычислить приближенный собственный вектор x k

y k

.

y k

2

5.

Определить отношение Рэлея k x k , A .

6.

Если выполняется условие

Ax k k x k

,

то сходимость

2

достигнута; в противном случае положить k k 1 и перейти к

п.3.

Если условие останова выполнено, то число k удалено от некоторого собственного значения матрицы А не более, чем на ε.

5.8.4. «Разделяй-и-властвуй»

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

O n3 операций; на большой выборке случайных тестовых матриц в среднем затрачивается только O n2.3 , а для некоторых специальных распределений собственных значений – O n2 операций.

5.8.5. Бисекция и обратная итерация

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

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

Рассмотрим трехдиагональную матрицу вида

				0	0	0	0
		1	1
	1		2	2	0	0	0
A								.	(5.12)
					n 2	n 1	n 1
0			0	0
	0		0	0	0	n 1	n
Обозначим через			1 , 2 , , n главные миноры						трехдиагональной

матрицы А. Для их вычисления можно получить рекуррентную формулу путем разложения по элементам j-го столбца

0

1;

1

1;

(5.13)

2

j

j

j 1

j 2

, j 2, n.

j 1

Выражение (5.13) называют последовательностью Штурма . Из нее

ивида матрицы А вытекают следствия:

1.Никакие два соседних минора матрицы (5.12) не могут одновременно равняться нулю;

2.Если минор j , 1 j n, равняется нулю, то соседние миноры отличны от нуля и имеют противоположные знаки;

3.Все собственные значения матрицы (5.14) простые.

Вычисления по формулам (5.13) устойчивые в случае арифметики с плавающей запятой, но даже для матриц малого порядка значения i при i,

близких к n, выходит за пределы диапазона допустимых для ЭВМ чисел.

Особенно это происходит, когда матрица имеет кластер близких собственных

значений.

Поэтому

последовательность

i

заменяют

такой

последовательностью:

i

q

i

,

i 1, n.

i 1

Полиномы qi

удовлетворяют соотношениям

2

q

,

q

i

, i 2, n.

i

i

1

1

qi 1

Пусть все собственные значения пронумерованы в порядке убывания, т.е.

1 2 n .

Определим

k-е

собственное

значение

независимо от

остальных.

n

Обозначим через

количество собственных значений матрицы А,

строго больших λ. Если известны такие числа a0 , b0 ,

что b0

a0 , n a0 k ,

n b0 k ,

то k

принадлежит

полуинтервалу

a0 , b0 .

Примем

теперь

операций, бисекция может быть гораздо

c0

a0

b0

и

вычислим n c0 . Если

n c0 k , то

k принадлежит

2

c0 , b0 . Если n c0 k , то

полуинтервалу

k принадлежит полуинтервалу

a

, c

. Поэтому всегда можно найти полуинтервал длины

b0 a0

, который

0

0

2

содержит k . Повторяя этот процесс, получим систему вложенных интервалов
as , bs , содержащих k , причем
b a	s	2 s b a		0	.
s		0
Это дает возможность локализировать			собственное значение k с

произвольной, наперед заданной точностью.

Бисекцию можно использовать для вычисления некоторого подмножества собственных значений, например, тех, которые расположены на отрезке a,b или с номерами i , i k . На ее реализацию не влияет наличие близких или

кратных собственных значений, даже большого их скопления. Она обходится только в O kn операций, где k – количество требуемых собственных значений.

Т.к. для QR-итерации нужно O n2

быстрее при k n . Соответствующие собственные векторы могут быть найдены с помощью обратной итерации. В наилучшем случае, когда нет близких собственных значений, обратная итерация может стоить также O kn операций. В наихудшем случае, когда имеются обширные кластеры собственных значений, трудоемкость обратной итерации возрастает до O k 2 n операций, при этом не гарантируется качество вычисленных собственных векторов (хотя на практике они почти всегда вычисляются с хорошей точностью).

Пример. Методом бисекций найти наибольшее собственное значение 1

	2	1	0
матрицы A	1	2	1	.
	0	1	2

Найдем сначала собственные значения классическим способом, решив характеристическое уравнение

	2	1	0
det A E	1	2	1		0 ,
	0	1	2

т.е.

2 3 2 2 0 .

Получим корни

1 2	2, 2 2,		1 2 2,

т.е.

1 2 3 .

Теперь применим метод бисекций для нахождения 1 . Оценим нормы

A

A

1 4 , т.е.

i

4 ,

i 1, 2, 3.

Примем a0 4 ,

b0 4 и построим

последовательность Штурма

0 1; 1 1 2 ;

2 2 1 12 0 2 2 1;

3 3 2 22 1 2 2 1 .

Вычислим значения последовательности на концах отрезка 4, 4 .

В точке -4: 0 1,	1 6 ,		2 35, 3 204 ; n 4 0 , n 4 3 ,
изменений знаков – 0, собственных значений, меньших -4 нет.
В точке 4: 0 1,	1 2, 2 3 , 3 5;							n 4 3,					n 4 0 , т.е.
имеются три собственных значения на полуинтервале 4, 4 . Вычислим
	c			a0 b0	0.
		0
			2
В точке 0: 0 1,	1 2 ,		2 3 ,		3 4 ; n 0 3,						т.е.		имеются три
собственных значения на полуинтервале 0, 4 . c							a1	b1		0 4			2 .
					1			2			2
В точке 2: 0 1,	1 0, 2 1, 3					0 ,		следовательно, 2 –

собственное значение, одно изменение знака, одно собственное значение,

больше 2, другое – больше 2; n 2 1.

n 2 1,

n 4 0 1, т.е.

1 –

наибольшее собственное

значение,

2, 4 . c

a2 b2

2 4

3.

1

2

2

2

В точке 3:

0 1, 1 1, 2 0 , 3

1; два

изменения знаков,

n 3 1, n 3 2 ,

n 3 1,

n 4 1,

1 3, 4 .

c3 3,5.

В точке 3,5:

0 1,

1 1,5,

2 1,25,

3 0,375; три изменения

знаков, n 3 1, n 3,5 0 1, n 4 1.

Следовательно, 1 3; 3,5 , примем 1 3,25 .