Методическое пособие для практических и лабораторных работ по курсу «Численные методы»

170
Вместо прямых итераций, требующих предварительного обращения исходной матрицы
A , обычно предпочитают строить ту же последовательность векторов
{ }
k
yr
, решая при
1, 2,...
k
=
линейные системы
1
k
k
Ay
y
−
=
r r
,
1, 2,...
k
=
Построение последовательности векторов, приближающих собственный вектор
n
xr по этой неявной формуле, и называют обратными
итерациями.
Важными достоинствами степенного метода являются его простота, отсутствие необходимости преобразования матрицы A.
Недостаток метода в применении ко многим прикладным задачам – зачастую медленная сходимость. Часто в приложениях
2
λ
оказывается близким к
1
λ
. Так как скорость сходимости степенного метода определяется величиной
1 1
2
≈
λ
λ
, то в этом случае сходимость будет очень медленной.
Чтобы преодолеть эту неприятность прибегают, в частности, к степенному
методу со сдвигом. Идея метода со сдвигом заключается в следующем: применяют степенной метод не к матрице A, а к матрице
E
A
A
σ
−
=
)
, где
σ
– некоторое число. Собственными значениями матрицы A
)
являются числа
σ
λ −
i
, получаемые сдвигом собственных значений
i
λ
на число
σ
. Если число
σ
λ −
1
по-прежнему остается максимальным по модулю, то следует попытаться подобрать
σ
так, чтобы сделать величину
σ
λ
σ
λ
−
−
≤
≤
1 2
max
i
n
i
минимальной. Если, например, все собственные значения положительны, то такой минимум достигается при
2 2
n
λ
λ
σ
+
=

9.3. Методы исчерпывания
После того как пара
(
)
1 1
, xr
λ
найдена, возникает задача о вычислении следующей пары
(
)
2 2
, xr
λ
, где
n
λ
λ
λ
>
>
>
3 2
,
2
xr соответствующий
2
λ
собственный вектор. Существует несколько способов избавления от уже вычисленных собственных чисел и соответствующих собственных векторов с целью избежать их повторного вычисления. Эти способы принято называть
исчерпыванием (методами исчерпывания). Приведем для примера один из методов исчерпывания, основанный на подобном преобразовании матриц.
Напомним, что подобные матрицы имеют одинаковые собственные числа и соответствующие им собственные вектора, которые всегда определяются с точностью до множителя, умножаются на соответствующие матрицы.

171
Например, пусть матрицы A и B подобны, т.е.
1
−
=
HAH
B
, где H – некоторая неособая матрица. Пусть справедливо равенство
x
x
A
r r
λ
=
, тогда верно равенство
x
H
x
HA
r r
λ
=
или
x
H
x
H
HAH
r r
λ
=
−
)
(
1
. Это выражение можно переписать в виде
y
y
B
r r
λ
=
, где
x
H
y
r r
=
. Последние равенства и устанавливают связь между собственными числами и собственными векторами подобных матриц.
Пусть задача
1 1
1 1
x
x
A
r r
λ
=
решена, т.е.
1
λ
и
)
,...,
(
)
(
1
)
1
(
1 1
n
x
x
x
=
r найдены, здесь
A
A
≡
1
Образуем матрицу




















−
−
=
1 0
0
...
...
0 0
1 0
0 0
1
)
1
(
1
)
(
1
)
1
(
1
)
2
(
1
)
1
(
1 1
x
x
x
x
x
H
n

Тогда
(
)
0
,...,
0
,
1
,
1 1
1
=
=
e
e
x
H
r r
r
. Равенство
1 1
1 1
x
x
A
r r
λ
=
преобразуем в эквивалентное
)
(
)
(
1 1
1 1
1 1
1 1
x
H
x
H
H
A
H
r r
λ
=
−
или
1 1
1 1
1 1
1
e
e
H
A
H
r r
λ
=
−
, откуда следует, что первый столбец матрицы
1 1
1 1
−
H
A
H
равен
1 1
er
λ
. Введем в рассмотрение матрицу




′
=
=
−
2 1
1 1
1 1
1 2
0
B
b
H
A
H
A
r v
λ
,
(9.9)
где
)
,...,
(
)
(
2
)
2
(
1 1
n
b
b
b
=
′
r и
2
B – матрица размерности
)
1
(
−
n
. Из (9.9) следует, что матрица
2
B имеет своими собственными значениями числа
n
λ
λ
,...,
2
, т.е. соответствующие собственные числа матрицы A . Для того, чтобы найти следующую пару
(
)
2 2
, xr
λ
нужно уже работать с матрицей
2
B . После того, как
2
λ
и
2
xr найдены, делаем соответствующие преобразования с матрицей
2
B и т.д., пока не «исчерпаем» все собственные числа исходной матрицы A.

172
Лекция 10. Метод Якоби решения полной проблемы собственных значений для симметричной матрицы. QR- метод.
Уточнение собственных чисел и векторов. Оценки собственных чисел. Теоремы Гершгорина.
10.1. Метод Якоби решения полной проблемы собственных
значений для симметричной матрицы
Рассмотрим задачу
x
x
A
r r
λ
=
(10.1)
Требуется решить полную проблему собственных значений: найти все собственные числа и соответствующие им все собственные векторы матрицы
A . Вспоминая свойства подобных матриц, поставим вопрос: нельзя ли подобными преобразованиями привести матрицу A к такой матрице, для которой собственные числа легко находятся?
Пусть A – симметричная вещественная матрица. Известно [6],
что такую матрицу всегда можно представить в виде

1
−
Λ
=
Q
Q
A
,
(10.2)

где Q – ортогональная матрица,
Λ
–
диагональная матрица у
которой на главной диагонали стоят собственные числа
матрицы A. Как мы знаем для ортогональной матрицы
T
Q
Q
=
−
1
,
поэтому равенство (10.2) равносильно следующему

Λ
=
AQ
Q
T
(10.3)

Перепишем (10.3) в виде

173

Q
AQ
Λ
=
(10.4)

Из равенства (10.4) следует, что столбцы матрицы Q являются
соответствующими собственными векторами матрицы A.
Формула (10.3) дает возможность строить алгоритмы
вычисления матрицы
Λ
,
отличающиеся между собой лишь
способами вычисления матрицы Q . В основе их лежит
следующий простой факт. Пусть каким-либо ортогональным
преобразованием Q
получили

Λ
=
AQ
Q
T
(10.5)

где матрица
)
(
ij
λ
=
Λ
мало отличается от
Λ
.
Собственные числа
матриц A и
Λ
совпадают. Если все
ij
λ
при
j
i
≠

равны нулю, то
собственные значения матрицы A равны диагональным
элементам матрицы
Λ
, т.е.
ii
λ
. Если же недиагональные
элементы
ij
λ
не равны нулю, но имеют малые значения, то
можно ожидать, что собственные значения
i
λ

близки к
i
λ
и
i
λ

можно взять за их приближенные значения. Поэтому на
основании
равенства
(10.3)
строят
последовательность
ортогональных приближений, позволяющих неограниченно
уменьшать модули недиагональных элементов матрицы A:

Λ

 →

∞
→
n
n
T
n
AW
W
,
n
n
Q
Q
Q
W
⋅⋅
⋅
=
2 1
.

174

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

2
)
(
∑
≠
=
j
i
ij
a
A
σ
(10.6)

Перейдем теперь непосредственно к изложению метода
вращений решения полной проблемы собственных значений
симметричной матрицы (
метода Якоби).
По заданной матрице A будем строить последовательность
)
( k
A
такую, чтобы матрица
)
1
(
+
k
A

получалась из матрицы
)
( k
A при
помощи преобразования подобия ортогональной матрицей
вращения
)
(
ϕ
ij
Q
:

столбец
й
j
столбец
й
i
строка
я
j
строка
я
i
Q
ij
−
−
−
−




























−
=
1 1
cos sin sin cos
1 1
)
(
O
M
M
L
L
L
M
O
M
L
L
L
M
M
O
ϕ
ϕ
ϕ
ϕ
ϕ
(10.7)

175
В матрице (10.7) выписаны только ненулевые элементы. Все
остальные элементы равны нулю.
Предположим, что уже построена матрица
( )
)
(
)
(
k
ij
k
a
A
=
.
Найдем в
ней максимальный по модулю недиагональный элемент. Пусть
это будет
)
( k
ij
a
.
Ввиду симметричности можно считать, что
j
i
<
.
Если таких элементов несколько, то берем любой из них. При
фиксированных индексах
j
i, строим матрицу вращения
)
(
)
(
)
(
)
(
k
k
ij
k
ij
Q
Q
ϕ
=
,
угол
)
( k
ϕ
определим позже. Образуем матрицу

)
(
)
(
)
(
)
1
(
k
ij
k
T
k
ij
k
Q
A
Q
A
=
+
(10.8)

Чтобы изучить вид элементов матрицы
)
1
(
+
k
A
,
предварительно
рассмотрим
вспомогательную
матрицу
)
(
)
(
)
(
)
(
)
(
k
ij
k
k
ij
k
Q
A
b
B
=
≡
.
Согласно определению матрицы
ij
Q (10.7), все столбцы матрицы
)
( k
B
кроме i -го и j -го будут такими же как и в
)
( k
A , элементы же
i -го и
j -го столбцов вычисляются по формулам

k
k
j
k
k
i
k
i
a
a
b
ϕ
ϕ
ν
ν
ν
sin cos
)
(
)
(
)
(
)
(
+
=
,
k
k
j
k
k
i
k
j
a
a
b
ϕ
ϕ
ν
ν
ν
cos sin
)
(
)
(
)
(
)
(
+
−
=
,
n
,...
2
,
1
=
ν
(10.9)

Аналогично, строки матрицы
)
(
)
(
)
1
(
k
T
k
ij
k
B
Q
A
=
+

кроме i -ой и j -ой
будут такими же как в
)
( k
B
,
а элементы i -ой и j -ой строк
вычисляются по формулам

k
k
j
k
k
i
k
i
b
b
a
ϕ
ϕ
ν
ν
ν
sin cos
)
(
)
(
)
(
)
1
(
+
=
+
,
(10.10)

176
k
k
j
k
k
i
k
j
b
b
a
ϕ
ϕ
ν
ν
ν
cos sin
)
(
)
(
)
(
)
(
+
−
=
,
n
,...
2
,
1
=
ν

Из (10.9) и (10.10), учитывая, что
)
(
)
(
k
ji
k
ij
a
a
=
получаем

(
)
(
)
(
)
2
sin
2 1
2
cos sin cos sin cos cos sin sin cos
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
(
)
1
(
k
k
ii
k
jj
k
k
ij
k
k
k
jj
k
k
ji
k
k
k
ij
k
k
ii
k
k
jj
k
k
ij
k
ij
a
a
a
a
a
a
a
b
b
a
ϕ
ϕ
ϕ
ϕ
ϕ
ϕ
ϕ
ϕ
ϕ
ϕ
−
+
=
+
−
+
+
+
−
=
+
=
+
(10.11)

Выбираем теперь угол
)
( k
ϕ
так, чтобы элемент
)
1
(
+
k
ij
a

обратился в
нуль. Находим, что

k
k
jj
k
ii
k
ij
k
P
a
a
a
tg
≡
−
=
)
(
)
(
)
(
)
(
2 2
ϕ
,

отсюда

k
k
arctgP
2 1
)
(
=
ϕ
.

Оценим теперь
)
(
)
1
(
+
k
A
σ
.
Учитывая симметричность матрицы A,
используя формулы (10.9) – (10.11), получаем

177

(
) ( )
( )
(
)
[
]
( )
( )
(
)
( )
( )
2 2
2 1
2 2
cos
2 2
sin
2 1
2 2
)
(
)
(
)
1
(
2
)
(
)
(
2
)
(
)
(
)
(
)
(
)
(
2
)
(
)
(
)
1
(
k
ij
k
k
ij
k
ij
k
k
k
ij
k
k
ii
k
jj
k
ij
k
k
a
A
a
a
A
a
a
a
a
A
A
−
=
+
−
=
=
+
−
+
−
=
+
+
σ
σ
ϕ
ϕ
σ
σ

Так как
)
( k
ij
a
согласно выбору максимальный по модулю элемент в
)
( k
A , то

( )
( )
2
)
(
)
(
)
1
(
k
ij
k
a
n
n
A
−
≤
σ
.

отсюда


( )
( )
)
1
(
)
(
2
)
(
−
≥
n
n
A
a
k
k
ij
σ
,

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

(
) ( )
( )
( )
( )
( )
)
(
)
(
)
(
2
)
(
)
(
)
1
(
)
1
(
2 2
k
k
k
k
ij
k
k
A
q
n
n
A
A
a
A
A
σ
σ
σ
σ
σ
=
−
−
≤
−
=
+
,

178
где
)
1
(
2 1
−
−
=
n
n
q
.
Очевидно, что
1 0
<
≤
q

при
2
≥
n
.
Поэтому

(
)
( )
( )
)
0
(
1
)
(
)
1
(
A
q
A
q
A
k
k
k
+
+
≤
≤
≤ σ
σ
.

Итак

( )
( )
)
0
(
)
(
A
q
A
k
k
σ
σ
≤
.

Так как
const
A
=
)
(
σ
,
то

( )
0
lim
)
(
→
∞
→
k
k
A
σ
,

т.е. процесс сходится.
В методе вращений (методе Якоби) решения полной проблемы
собственных значений симметричной матрицы A подобными
преобразованиями исходную матрицу A сводим к диагональной,
у которой на главной диагонали стоят собственные числа
матрицы A . На каждом шаге наших преобразований зануляем
максимальный по модулю недиагональный элемент матрицы
)
( k
A .
Число
недиагональных
элементов
конечно,
но
вычислительная погрешость заставляет рассматривать этот
метод как итерационный. Расчеты ведут до тех пор, пока станет
справедливым неравенство

179

( )
ε
σ
≤
)
( k
A
.
10.2. QR –
метод решения полной проблемы собственных значений

Наиболее универсальным (не требующим никакой предварительной информации относительно собственных чисел матрицы) в настоящее время является QR-метод.
QR-метод был предложен почти одновременно российским математиком В.Н. Кублановской (1961 г.) и системным программистом англичанином Д. Фрэнсисом (1962 г.). Суть метода заключается в следующем.
Пусть A – вещественная
n
n
×
матрица. При
,...
2
,
1
,
0
=
k
, начиная с
A
A
=
0
, строят последовательность матриц
)
(
k
A по формулам
k
k
k
R
Q
A
=
,
k
k
k
Q
R
A
=
+
1
,
(10.12)
первая из которых означает разложение матрицы
k
A в произведение матриц: ортогональной
k
Q и правой треугольной
k
R (такое разложение существует для любой квадратной матрицы), а вторая – перемножение полученных в результате факторизации матрицы
k
A матриц
k
Q и
k
R в обратном порядке.
На основе свойства ортогональных матриц
1
−
=
k
T
k
Q
Q
в соответствии с
(10.12) можно записать представление данной матрицы
A в виде
T
T
T
k
T
k
k
k
k
Q
Q
Q
Q
A
Q
Q
Q
Q
A
0 1
1 1
1 1
0
...
−
+
−
=
, или, иначе,
1 1
1 0
1 1
1 0
)
(
)
(
−
−
+
−
=
k
k
k
k
k
Q
Q
Q
Q
A
Q
Q
Q
Q
A

Значит, любая из этих матриц последовательности
)
(
k
A ортогонально подобна матрице A. Как известно, подобные матрицы имеют одинаковые собственные числа.
При определенных ограничениях, одним из которых выступает требование, чтобы матрица A не имела равных по модулю собственных

180
значений, генерируемая процессом (10.12) последовательность матриц
)
(
k
A сходится к правой треугольной матрице с диагональю из искомых
собственных чисел (это ограничение ввели, чтобы не усложнять изложение материала). Скорость обнуления поддиагональных частей матриц
k
A линейна и зависит от отношений
|
|
/
|
|
j
i
λ
λ
при
j
i
>
(считаем, что
|
|
|
|
|
|
2 1
n
λ
λ
λ
>
>
>
).
Наличие комплексно-сопряженных пар собственных чисел у данной вещественной матрицы A не является препятствием для применения QR- алгоритма, просто в этом случае предельной матрицей для последовательности
)
(
k
A будет матрица квазитреугольного (иначе, блочно – треугольного) вида. Каждой комплексной паре собственных чисел в такой матрице будет соответствовать
2 2
×
– блок, причем сходимость здесь наблюдается по форме матрицы, а не поэлементно (т.е. элементы внутри этих блоков могут изменяться без видимой зависимости от k при сохранении неизменными их собственных чисел).
Обычно QR-алгоритм применяют не к исходной матрице A, а к подобной ей правой почти треугольной матрице B , называемой также матрицей Хессенберга, вида




















=
−
−
−
−
−
−
−
nn
n
n
n
n
n
n
n
n
n
n
n
n
b
b
b
b
b
b
b
b
b
b
b
b
b
b
b
B
1
,
,
1 1
,
1 3
1
,
3 32 2
1
,
2 22 21 1
1
,
1 12 11 0
0 0
0
...
...
...
0
...

Предварительное приведение матрицы к верхней почти треугольной форме связано вот с чем. Разложение QR на каждом этапе в (10.12) требует
)
(
3
n
O
операций. Если же матрица является верхней почти треугольной, то
QR -разложение осуществляется лишь за
)
(
2
n
O
действий. Важно, что верхняя треугольная форма сохраняется при преобразованиях QR -алгоритма.
В основе преобразования матрицы A к матрице B (определяемому условием
0
=
ij
b
при
1
−
<
i
j
) лежит преобразование Хаусхолдера или, иначе, преобразование отражения
T
E
H
ω
ω
r r
2
−
=
, где
ω
r
– произвольный вектор- столбец, но такой, что его евклидова норма равна единице. Легко проверить, что
E
HH
T
=
, т.е. матрица отражения ортогональна, а значит, матрицы A и
B связаны соотношением
HAH
B
=
(
)
1
−
=
=
HAH
HAH
T
и являются подобными.
Теперь нужно распорядиться свободой задания элементов векторов
ω
r

181
при построении матриц отражения, чтобы за конечное число шагов преобразований Хаусхолдера произвольно заданную матрицу A привести к форме Хессенберга B .
А именно, можно показать, что начатым с
A
B
=
1
процессом
m
m
m
m
H
B
H
B
=
+
1
,
2
,...,
2
,
1
−
=
n
m
,
(10.13)
где
T
m
m
m
E
H
ω
ω
r r
2
−
=
, данная
n
n
×
матрица
A за
2
−
n
шага будет приведена к виду B , т.е. матрица
1
−
=
n
B
B
подобна
A , если задающие матрицы
Хаусхолдера
m
H векторы
m
ω
r по данной матрице
A строить следующим образом.
При
1
=
m
вектор
1
ω
r определяется равенством
),
;...;
;
;
0
(
1 31 1
21 1
1
n
T
a
a
s
a
−
= µ
ω
r
(10.14) где
∑
=
⋅
−
=
n
i
i
a
a
sign
s
2 2
1 21 1
)
(
,
)
(
2 1
21 1
1 1
a
s
s
−
=
µ
. Такое задание
1
ω
r обеспечивает ортогональность симметричной матрицы
T
E
H
1 1
1 2
ω
ω
r r
−
=
и одновременно получение с ее помощью нужных
2
−
n
нулей в первом столбце матрицы
1 1
1 2
H
B
H
B
=
(
)
1 1
AH
H
=
Вектор
2
ω
r по матрице
2
B строят совершенно аналогично, только фиксируют нулевыми не одну, а две первые его координаты, и определяющую роль играют теперь не первый, а второй столбец матрицы
2
B и его третий элемент. При этом у матрицы
2 2
2 3
H
B
H
B
=
окажется
3
−
n
нулевых элемента во втором столбце и сохранятся полученные на предыдущем шаге нули в первом столбце.
Этот процесс очевидным образом может быть продолжен до исчерпания и без особого труда описан общими формулами типа формул
(10.13) и (10.14), для чего нужно лишь ввести обозначения для элементов последовательности матриц
m
B (например, с помощью верхних индексов
m
).

Скачать 1.29 Mb.

Поделитесь с Вашими друзьями: