Wyznaczanie wartości i wektorów własnych macierzy

(problem własny)

Plan wykładu:

1. Pojęcia podstawowe, definicje

2. Metoda Kryłowa poszukiwania pierwiastków równania

charakterystycznego

3. Lokalizacja (szacowanie) wartości własnych – tw. Gerschgorina
4. Proste metody iteracyjne

a)

Metoda potęgowa dla wartości własnej o największym module

b)

przyśpieszanie zbieżności

c) redukcja macierzy metodą Hottelinga
d) redukcja macierzy metodą Wielandta

5. Sprowadzanie macierzy symetrycznych/hermitowskich do postaci

trójdiagonalnej

a) metoda Hauseholdera

b) Metoda Lanczosa

6. Sprowadzanie macierzy symetrycznych/hermitowskich do postaci

macierzy Hessenberga

7. Wyznaczanie wartości i wektorów własnych macierzy trójdiagonalnych

i macierzy Hessenberga

a) Metoda bisekcji
b) Rozkład LR i QR, rozkład QR metodą Hauseholdera

8.Uogólniony problem własny

2

Pojęcia podstawowe

Często przy tworzeniu modeli matematycznych
wykorzystywanych do symulacji zjawisk
fizycznych czy zachowania się układu, zachodzi
potrzeba rozwiązania tzw.

problemu własnego

(np. rów. Schrodingera):

- A jest macierzą kwadartową o nxn
- x

k

jest

wektorem własnym

macierzy

odpowiadającym

wartości własnej



k

Nie zawsze układ równań, którego chcemy
znaleźć rozwiązanie, przyjmuje tak prostą
postać. Nierzadko mamy do czynienia z tzw.

uogólnionym problemem własnym

:

Jeśli macierz B jest nieosobliwa to problem
uogólniony można przekształcić
do postaci:

Ax

k

= ¸

k

x

k

a

ij

; ¸

k

; x

(k)

m

2 C

A

¢ x = ¸B ¢ x

B

¡1

Ax = ¸x

Liczbę

 nazywamy wartością własną macierzy jeśli

istnieje taki niezerowy wektor x dla którego zachodzi:

Wektor x nazywamy (prawostronnym) wektorem
własnym przynależnym do wartości własnej

. Ciąg

wszystkich wartości własnych nazywamy widmem
macierzy A i oznaczamy: Sp(A).
Z powyższej definicji wynika:

Macierz (A-

I) jest osobliwa, więc:

Wyznacznik ten jest wielomianem stopnia n zmiennej
:

Każda wartość własna



k

jest pierwiastkiem

wielomianu charakterystycznego

macierzy A.

Ax = ¸x

(A

¡ ¸I)x = 0

det(A

¡ ¸I) = 0

w(¸)

=

det(A

¡ ¸I) =

=

(

¡1)

n

(¸

n

+ a

n

¡1

¸

n

¡1

+ : : : + a

0

)

A = [a

ii

]

3

Def. Wartości i wektory własne macierzy
transponowanej A

T

nazywamy

lewostronnymi

wartościami

i

lewostronnymi wektorami

własnymi

macierzy A.

Wyznacznik macierzy po jej transponowaniu nie
ulega zmianie. Dlatego widmo macierzy A jest
równe widmu lewostronnemu.

Tw. Jeżeli



p

jest wartością własną macierzy, a



l

jest jej lewostronną wartością własną
oraz gdy

Wówczas wektor własny x

l

jest ortogonalny do

lewostronnego wektora własnego x

p

.

Dla macierzy symetrycznej A=A

T

wektory

własne są zarazem wektorami

lewostronnymi.

Jeżeli więc wektory własne przynależą do różnych
wartości własnych to są do siebie

ortogonalne

.

¸

p

6= ¸

l

(¸

p

6= ¸

l

)

Def. Macierze A i B są podobne jeśli istnieje
nieosobliwa macierz podobieństwa P, że:

Tw. Jeżeli macierze A i B są podobne to mają
identyczne widmo wartości własnych.

Tw. Macierz Q

mxn

(m

≥ n) nazywamy ortogonalną

jeśli:

Tw. Jeżeli macierz Q

nxn

jest ortogonalna to:

Tw. Macierz symetryczna A jest ortogonalnie
podobna do macierzy diagonalnej D:

Tw. Wartości własne macierzy symetrycznej są
rzeczywiste.

Def. Macierz o elementach zespolonych i własności

nazywamy macierzą hermitowską.
Wartości własne macierzy hermitowskiej są
rzeczywiste a wektory własne są do siebie
ortogonalne.

P

¡1

AP = B

Q

T

Q = I

n

£n

QQ

T

= I

n

£n

Q

T

AQ = D

x

x

x

T

l

x

x

x

p

= 0

x

x

x

T

l

Ax

x

x

p

=

¡

A

T

x

x

x

l

¢

T

x

x

x

p

= ¸

l

x

x

x

T

l

x

x

x

p

(¸

l

¡ ¸

p

)x

x

x

l

x

x

x

p

= 0

x

x

x

T

l

Ax

x

x

p

= ¸

p

x

x

x

T

l

x

x

x

p

A = (A

T

)

¤

) a

ii

2 R

4

Tw. Dla dowolnej macierzy A istnieje macierz
nieosobliwa P, która może mieć elementy
zespolone i zachodzi pomiędzy nimi związek

Powyższa macierz definiuje

postać kanoniczną

Jordana

J

k

jest macierzą zdefiniowaną następująco



i

jest wartością własną macierzy A i może wystąpić

w wielu macierzach J

k

.

k = 1; 2; : : : ; K

· n

J

k

=

2
6

6

6

6

4

¸

i

1

0

: : :

0

0

0

¸

i

1

: : :

0

0

: : :

: : :

: : :

: : :

: : :

: : :

0

0

0

: : :

¸

i

1

0

0

0

: : :

0

¸

i

3
7

7

7

7

5

Wyznaczniki

są

dzielnikami elementarnymi

macierzy A.

Liczba m jest stopniem macierzy J

k

. Dla m=1

macierz J

k

stanowi

dzielnik liniowy

.

Jeśli wszystkie wartości własne macierzy A są różne
to wszystkie dzielniki elementarne są liniowe.

Iloczyn skalarny (wewnętrzny)

Iloczyn zewnętrzny

det(J

k

¡ ¸I) = (¸

i

¡ ¸)

m

P

¡1

AP =

2

6

6

4

J

1

0

: : :

0

0

J

2

: : :

0

: : :

: : :

: : :

: : :

0

0

: : :

J

k

3

7

7

5

haaa; bbbi = aaa

T

bbb = [®

1

; ®

2

; : : : ; ®

n

]

2
6

6

6

4

¯

1

¯

2

..

.

¯

n

3
7

7

7

5

a

a

a

ihbbb = aaabbb

T

=

2

6

6

6

4

®

1

®

2

..

.

®

n

3

7

7

7

5

[¯

1

; ¯

2

; : : : ; ¯

n

]

=

2
6

6

4

®

1

¯

1

®

1

¯

2

: : :

®

1

¯

n

®

2

¯

1

®

2

¯

2

: : :

®

2

¯

n

: : :

: : :

: : :

: : :

®

n

¯

1

®

n

¯

2

: : :

®

n

¯

n

3
7

7

5

5

Tw. (Schura) Suma kwadratów modułów wartości
własnych jest ograniczona od góry przez kwadrat
normy euklidesowej:

Tw. Widmo macierzy ulega przesunięciu po
dodaniu do niej macierzy jednostkowej
pomnożonej przez liczbę:

widmo:

zostaje zastąpione przez:

Tw. (Cayleya-Hamiltona)
Jeśli

jest równaniem charakterystycznym macierzy A
to

n

X

i=1

j¸j

2

· kAk

2

E

Sp(A + cI) = Sp(A) + c

Sp(A) =

f¸

1

; ¸

2

; : : :

g

Sp(A + cI) =

f¸

1

+ c; ¸

2

+ c; : : :

g

Metoda Kryłowa poszukiwania zer równania
charakterystycznego

Korzystając z tw. CH

Co dla dowolnego wektora y daje

układ n równań na n niewiadomych

Do jego utworzenia potrzeba jednak
n

3

obliczeń, oraz n

3

/3 aby go rozwiązać.

Uwaga:
Wyznaczenie zer wielomianu charakterystycznego
może być źle uwarunkowane.

w(¸) = det(A

¡ ¸I) = 0

w(A) = 0

w(¸) = ¸

n

+

n

¡1

X

i=0

b

i

¸

i

= 0

w(A) = A

n

+

n

¡1

X

i=0

b

i

A

i

= 0

b

0

; b

1

; b

2

; : : : ; b

n

¡1

A

n

yyy +

n

¡1

X

i=0

b

i

A

i

yyy = 0

6

Lokalizacja wartości własnych

Tw. (Gershgorina) Niech C

i

oznaczają koła

domknięte na płaszczyźnie zespolonej o
środkach w punktach a

ii

(elementy diagonalne

macierzy A) i promieniach równych sumie
modułów elementów z danego (i-tego) wiersza:

Wówczas:

Jeżeli k kół C

i

tworzy zbiór rozłączny z

pozostałymi kołami, to w zbiorze tym leży
dokładnie k wartości wlasnych macierzy A.

Wnioski:
1. Jeżeli macierz jest symetryczna i diagonalnie
dominująca o nieujemnych elementach na
diagonali, to jest nieujemnie określona, a jeśli
jest ona dodatkowo nieosobliwa to jest dodatnio
określona. Macierz symetryczna silnie
diagonalnie dominująca jest dodatnio określona
wtedy i tylko wtedy gdy elementy na diagonali
są dodatnie.
2. W każdym kole rozłącznym z pozostałymi
leży dokładnie jedna wartość własna.

C

i

=

fz : jz ¡ a

ii

j ·

n

X

j=1

j

6=i

ja

ij

jg

Sp(A)

½

n

[

i=1

C

i

Wartości własne macierzy A leżą na płaszczyźnie
zespolonej i zawarte są w kole o środku w 0 i promieniu
równym promieniowi spektralnemu tej macierzy.
Ponieważ:

więc można przyjać że:

Widma wartości własnych i lewostronnych wartości
własnych są identyczne.

Aby otrzymać lepsze

oszacowanie położenia wartości własnych można więc
zastosować twierdzenie Geshgorina dla A

T

. Koła

zawierające wartości własne mają środki w a

ii

i promienie

równe sumie modułów pozostałych elementów

w i-tych kulmnach.

Przykład. Podać lokalizację wartości własnych macierzy

½(A)

· kAk

p

;

p = 1; 2;

1; E

¸

i

2 fz : jzj · kAk

p

g

A =

2

4

¡2 ¡1 0

2

0

0

0

0

2

3

5

7

a) najgorsze oszacowanie – lokalizacja w kole o promieniu 4
b) tw. Gershgorina – lokalizacja w kole o promieniu 3
c) tw. Gershgorina dla A

T

– jedno z kół jest odseparowane (C'

3

) i zdegenerowane

znajduje się w nim dokładnie jedna wartość własna (

=2) - najlepsze oszacowanie

8

Metody wyznaczania wartości

i wektorów własnych

macierzy

Proste metody

iteracyjne

Metoda potęgowa

I jej modyfikacje

Macierze o elementach

I wartościach

własnych rzeczywistych

Macierze symetryczne

(w tym

hermitowskie

)

Macierze niesymetryczne

Redukcja macierzy

do postaci

trójdiagonalnej

Redukcja macierzy do

postaci Hessenberga

Metoda

Hauseholdera

Redukcja do postaci

diagonalnej

metodą Jacobiego

Metoda

Givensa

Rozkład QR

Metoda Lanczosa

(

macierze rzadkie

)

Wartości i wektory

własne

9

Metoda potęgowa wyznaczania
pojedynczych wartości własnych
i wektorów własnych.

Załóżmy że istnieje n liniowo niezależnych
wektorów własnych macierzy A, stanowią
bazę przestrzeni liniowej

Wówczas dla dowolnego wektora v

0

Jeśli



i

stanowią wartości własne macierzy

Zakładamy że wartości własne tworzą ciąg

Jeśli



1

jest dominującą wartością własną,

oraz wektor v

0

ma składową w kierunku x

1

to

wówczas zachodzi

vvv

0

=

n

X

i=1

a

i

x

x

x

i

x

x

x

1

; x

x

x

2

; x

x

x

3

; : : : ; x

x

x

n

Avvv

0

=

n

X

i=1

a

i

¸

i

x

x

x

i

vvv

m

= A

m

vvv

0

=

n

X

i=1

a

i

¸

m

i

x

x

x

i

j¸

1

j ¸ j¸

2

j ¸ j¸

3

j ¸ : : : ¸ j¸

n

j

Z czego można wysnuć wniosek że wartość własną
można obliczyć następująco

Dla dowolnego wektora y nieortogonalnego do x

1

.

Zazwyczaj y ma 1 na pozycji elementu o największym
module w v

m+1

a na pozostałych 0.

Jaka jest zbieżność metody?

Zależy od (



i

/



1

)

m

ale również od współczynników a

i

czyli od wyboru v

0

. Jeśli wartość własna o największym

module jest zespolona to ciąg nie jest zbieżny.

Jak wyznaczyć wektor własny x

1

?

Ponieważ

więc unormowany wektor własny będzie miał postać

lim

m

!1

A

m

vvv

0

¸

m

1

= a

1

x

x

x

1

¸

1

= lim

m

!1

yyy

T

vvv

m+1

yyy

T

vvv

m

x

x

x

1

=

vvv

m

jvvv

m

j

vvv

m

= ¸

m

1

"

a

1

x

x

x

1

+

n

X

i=2

µ

¸

i

¸

1

¶

m

a

i

x

x

x

i

#

vvv

m

¼ ¸

m

1

a

1

x

x

x

1

10

Jeśli wartość własna jest pierwiastkiem
wielokrotnym rówanania charakterystycznego to
metoda jest zbieżna bo składnik z



1

dominuje

Uwaga: problem pojawia się gdy



1

=-



j

tj.

identyczne moduły generują oscylacje (wtedy
wybieramy ciąg wektorów v

2k

)

Przyśpieszanie zbieżności

1. Proces



2

Po podzieleniu licznika i mianownika przez b

1



1

m

a

następnie rozwinięciu M w szereg potęgowy
dostajemy

Co dla dużych m można przybliżyć

vvv

m

= A

m

vvv

0

= ¸

m

1

k

X

i=1

a

i

x

x

x

i

+

n

X

i=k+1

¸

m

i

a

i

x

x

x

i

r =

¸

2

¸

1

¸

1

¼ R

m+2

¡

(¢R

m+1

)

2

¢

2

R

m

Wykorzystując różnice progresywne 1 i 2 rzędu można
znaleźć lepsze przybliżenie

2. Iloraz Rayleigha

Jeśli macierz jest symetryczna to jej wektory własne
są ortogonalne. Załóżmy że są również ortonormalne

Wtedy

Oraz

co daje lepszą zbieżność niż wariant podstawowy
metody

vvv

T

m

Avvv

m

= vvv

T

m

vvv

m+1

=

n

X

i=1

a

2

i

¸

2m+1

i

x

x

x

T

i

x

x

x

j

= ±

ij

vvv

T

m

vvv

m

=

n

X

i=1

a

2

i

¸

2m

i

vvv

T

m

vvv

m+1

vvv

T

m

vvv

m

= ¸

1

+ O

"µ

¸

i

¸

1

¶

2m

#

¸

1

¼ R

m

¡ ¯

2

r

m

eee

T

j

vvv

m

=

n

X

i=1

a

i

¸

m

i

eee

T

j

x

x

x

i

=

n

X

i=1

b

i

¸

m

i

eee

T

j

vvv

m+1

eee

T

j

vvv

m

=

b

1

¸

m+1

1

+

P

n
i=2

b

i

¸

m+1
i

b

1

¸

m

1

+

P

n
i=2

b

i

¸

m

i

eee

T

j

vvv

m+1

eee

T

j

vvv

m

= ¸

1

+ ¯

2

µ

¸

2

¸

1

¶

m

+ ¯

3

µ

¸

3

¸

1

¶

m

+ : : :

R

m

=

eee

T

j

vvv

m+1

eee

T

j

vvv

m

11

Wyznaczanie pozostałych wartości własnych

a) metoda redukcji wektora

Jeśli znamy wartość własną o największym
module to możemy wykorzystać ten fakt przy
wyznaczaniu kolejnej największej co do modułu
wartości własnej tj.



2

Ponieważ wektory v

m+1

oraz



1

v

m

są bliskie więc

metoda może być w pewnych przypadkach
nieużyteczna.

b) metoda zerowania składowej

Znając



1

możemy zdefiniować wektor

Czyli z v

0

usuwamy składową w kierunku x

1

.

Wektor w

0

nie ma składowej w kierunku x

1

więc

ciąg w

m

powinien być zbieżny do



2

oraz x

2

. Ze

względu na błędy zaokrągleń jednak usuwa się co
pewną ilość iteracji składową x

1

tj.

w

w

w

m

= vvv

m+1

¡ ¸

1

vvv

m

=

n

X

i=2

a

i

(¸

i

¡ ¸

1

)¸

m

i

x

x

x

i

w

w

w

0

= (A

¡ ¸

1

I)vvv

0

zzz

m

= (A

¡ ¸

1

I)w

w

w

m

Aby wyznaczyć kolejne wartości własne korzystamy z
wektora

Taki proces staje się mało wydajny dla kolejnych
wartości własnych.

Uwaga:

Metoda redukcji składowej jest odpowiednikiem

metody metody czasu urojonego z
ortogonalizacją Schmidta

stosowaną w mechanice

kwantowej do rozwiązania równ. Schrodingera.

Zaletami są: prostota oraz fakt że nie trzeba
pamietać postaci macierzy (elementy mogą być
wyznaczane na bieżąco).

w

w

w

0

= (A

¡ ¸

1

I)(A

¡ ¸

2

I)vvv

0

12

Redukcja macierzy

Tw. Jeśli



1

jest wartością własną macierzy A i x

1

odpowiadającym jej wektorem własnym oraz dla
dowolnego wektora v o własności

Macierz zredukowana

Ma te same wartości co macierz A oprócz



1

która

jest zerem.

Pełny dowód – A. Ralston „Wstęp do analizy...”

Zakładamy

Definiujemy wektor u

(e

i

- i -ty wersor w n-wymiarowej przestrzeni

kartezjańskiej, lub i-ta kolumna macierzy
jednostkowej)
oraz macierz T

W

1

= A

¡ ¸

1

x

x

x

1

vvv

T

x

(1)
1

= 1

u

u

u = x

x

x

1

¡ eee

1

T = I + u

u

ueee

T

1

T T

¡1

= I

) T

¡1

= I

¡ uuueee

T

1

Przy pomocy T konstruujemy macierz podobną do A

pierwsza kolumna C jest równa
(C

10

– pierwszy wiersz macierzy C)

Ponadto definiujemy macierz D

w której wiersze od 2 do n są równe 0, natomiast
wiersz pierwszy ma postać

Wobec powyższego można rozpatrywać macierz

która ma te same wartości własne co A za wyjątkiem


1

która jest równa 0.

C = T

¡1

AT

C

01

= ¸

1

eee

1

D = T

¡1

x

x

x

1

vvv

T

T

vvv

T

x

x

x

1

= 1

C

¡ ¸

1

D = T

¡1

(A

¡ ¸

1

x

x

x

1

vvv

T

)T = T

¡1

W

1

T

D

10

= vvv

T

+ (1

¡ v

(1)

)eee

T

1

13

Związek wektorów własnych macierzy W i A

Algorytm wykorzystujący redukcję macierzy

1. wyznaczamy metodą potęgową



1

oraz x

1

2. konstruujemy macierz W

1

i znajdujemy



2

oraz w

2

3. konstruujemy macierz

i szukamy



3

oraz w

3

Uwaga:
jeśli wartosć własna jest k-krotna to
zastosowanie powyższej procedury k-krotnie da
tę samą wartość własną oraz k różnych
wektorów.

x

x

x

i

= (¸

i

¡ ¸

1

)w

w

w

i

+ ¸

1

(vvv

T

w

w

w

i

)x

x

x

1

¸

i

6= ¸

1

x

x

x

i

= w

w

w

i

¸

i

= ¸

1

;

(A

¡ ¸

1

I)w

w

w

i

= 0

W

2

= W

1

¡ ¸

2

w

w

w

2

vvv

T

1

vvv

T

1

w

w

w

2

= 1

Metoda redukcji Hotellinga

Za wektor v przyjmujemy lewy wektor własny
przynależny do wartości własnej



1

. Ale na ogół nie

znamy lewych wektorów.

Metoda jest więc skutecza tylko w przypadku macierzy
symetrycznych, wtedy lewe wektory są identyczne z
prawymi

lub rekurencyjnie

vvv = x

x

x

1

W

1

= A

¡ ¸

1

x

x

x

1

x

x

x

T

1

W

0

= A

W

i

= W

i

¡1

¡ ¸

1

x

x

x

i

x

x

x

T

i

i = 1; 2; : : : ; n

¡ 1

14

Metoda redukcji Wielandta

Wektor v definiujemy następująco

Gdzie: x

1

(j)

jest j-tą składową wektora x

1

,

A

j0

jest j-tym wierszem macierzy A

Z powyższego wzoru wynika że j-ty wiersz macierzy
W

1

jest równy

vvv

T

=

A

j0

¸

1

x

(j)
1

vvv

T

x

x

x

1

=

A

j0

x

x

x

1

¸

1

x

(j)
1

=

¸

1

x

(j)
1

¸

1

x

(j)
1

= 1

W

1

= A

¡ ¸

1

x

x

x

1

vvv

T

= A

¡

x

x

x

1

A

j0

x

(j)
1

(W

1

)

j0

= A

j0

¡

x

(j)
1

A

j0

x

(j)
1

= 0

Wnioski:

a) wektory własne odpowiadające niezerowym

wartościom własnym mają zerową j-tą
składową

b) wartości własne W

1

można obliczać (

metodą

potęgową

) skreślając uprzednio j-ty wiersz i j-

tą kolumnę, powstała macierz jest stopnia n-1

c) wartość j dobieramy tak by odpowiadała

pozycji elementu o największym module w x

1

d) w każdym kolejnym kroku metody mamy do

czynienia z macierzami coraz niższego stopnia

15

Redukcja macierzy gęstej do prostszej
postaci - macierzy trójdiagonalnej lub
Hessenberga

Macierz pierwotną przekształcamy
iteracyjnie

tak aby końcowa macierz B była macierzą
podobną do A

0

a następnie w łatwy sposób wyznaczamy
wartości i wektory własne macierzy B

Uwagi:

a) nakład obliczeń potrzebnych do
wyznaczenia x

i

oraz



i

macierzy B powinien

być jak najmniejszy

A = A

0

! A

1

! A

2

! : : : ! A

m

b) algorytm numerycznego rozwiązania problemu
własnego macierzy B nie może być gorzej
uwarunkowany niż dla A

Dla

problem własny B będzie gorzej uwarunkowany niż
dla A. Ale ponieważ

więc

uwarunkowanie algorytmu zależy od postaci

P

i

A

i

= P

¡1

i

A

i

¡1

P

i

;

i = 1; 2; : : : ; m

B

=

A

m

= P

¡1

AP

P

=

P

1

P

2

: : : P

m

B = P

¡1

AP

¢A = P ¢BP

¡1

cond(P ) >> 1

cond(P )

=

cond(P

1

: : : P

m

)

· cond(P

1

) : : : cond(P

m

)

jBj · jP

¡1

j ¢ jAj ¢ jP j = cond(P )jAj

j¢Aj · cond(P )j¢Bj

j¢Aj

jAj

· (cond(P ))

2

j¢Bj

jBj

B + ¢B = P

¡1

(A + ¢A)P

Byyy

=

¸yyy

x

x

x

=

Pyyy = P

1

: : : P

m

yyy

Ax

x

x

=

¸x

x

x

16

Redukcja macierzy hermitowskiej do
postaci trójdiagonalnej metodą
Hauseholdera

Pierwotna macierz hermitowska

Transformujemy

Przy pomocy

macierzy Hauseholdera

Dokonujemy transformacji A

i-1

do A

i

. A

i-1

ma

postać

A

H

= A = A

0

A

H

= (A

T

)

¤

a

a

a

i

=

2
6

6

6

6

4

®

i+1;i

¢
¢
¢

®

ni

3
7

7

7

7

5

Zauważmy że w kolejnym kroku należy dokonać
transformacji
a) (n-i) elementowego wektora a

i

b) macierzy kwadratowej

stopnia (n-i)

~

A

i

¡1

2

4

J

i

¡1

ccc

ccc

H

±

i

3

5 =

2
6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

4

±

1

¹

°

2

¢

¢

¢

0

0

°

2

. .. ...

¢

¢

¢

. .. ... ...

¢

¢

¢

. .. ... ...

¢

¢

¢

. .. ... ¹°

i

¡1

0

0

¢

¢

¢

°

i

¡1

±

i

¡1

¹

°

i

0

¢

¢

¢

0

°

i

±

i

3
7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

5

A

i

= P

¡1

i

A

i

¡1

P

i

P

i

= P

H

i

= P

¡1

i

= I

¡ ¯

i

u

u

u

i

u

u

u

H

i

A

i

¡1

=

2
6

6

6

4

J

i

¡1

ccc

0

ccc

H

±

a

a

a

H

i

0

a

a

a

i

~

A

i

¡1

3
7

7

7

5

=

[®

jk

]

17

Trzeba utworzyć macierz Hauseholdera
stopnia (n-i)

taką aby transformowała wektor a

i

Macierz Hauseholdera konstruujemy
następująco

¯ =

8

<

:

1

¾(¾+

j®

i+1;i

j)

¾

6= 0

0

¾ = 0

¾ =

kaaa

i

k

2

=

v

u

u

t

n

X

j=i+1

j®

ji

j

2

k =

¡¾e

i'

() ®

i+1;i

= e

i'

j®

i+1;i

j

u

u

u =

2
6

6

6

6

6

6

4

e

i'

(¾ +

j®

i+1;i

j)

®

i+2;i

..

.

®

n;i

3
7

7

7

7

7

7

5

Dokonujemy transformacji

~

P

i

~

P

i

a

a

a

i

= k

¢ eee

1

~

P

i

= I

¡ ¯uu

uu

uu

H

P

¡1

i

A

i

¡1

P

i

= P

i

A

i

¡1

P

i

=

=

2
6

6

6

4

J

i

¡1

ccc

0

ccc

H

±

i

a

a

a

H

i

~

P

i

0

~

P

i

a

a

a

i

~

P

i

~

A

i

¡1

~

P

i

3
7

7

7

5

=

2
6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

4

±

1

¹

°

2

0

0

°

2

. .. ...

..

.

0

. .. ... ¹°

i

¡1

0

0

°

i

¡1

±

i

¡1

¹

°

i

0

: : :

0

°

i

±

i

~

°

i+1

0 : : : 0

°

i+1

0

0

..

.

~

P

i

~

A

i

¡1

~

P

i

0

3
7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

5

P

i

=

2
6

6

6

4

I

0

0

0

1

0

0

0

~

P

i

3
7

7

7

5

gi ¡ 1

18

Jak dokonać efektywnego mnożenia macierzy
(zwykłe jest nieekonomiczne)?

Definiujemy dwa wektory pomocnicze

ponieważ

oraz z faktu, że A

i-1

jest hermitowska wynika

°

i+1

= k

p

pp = ¯ ~

A

i

¡1

u

u

u

qqq = pp

p

¡

¯

2

(ppp

H

u

u

u)u

u

u

¯

¸ 0

p

pp

H

u

u

u = ¯u

u

u

H

~

A

i

¡1

u

u

u = (ppp

H

u

u

u)

H

Przeprowadzając proces przekształcania macierzy
do końca dostaniemy macierz B

Metoda Hauseholdera jest stabilna numerycznie.

B = A

n

¡2

=

2
6

6

6

6

4

±

1

~

°

2

¢ ¢ ¢

0

°

2

. .. ... ...

..

.

. .. ... ~°

n

0

¢ ¢ ¢

°

n

±

n

3
7

7

7

7

5

~

P

i

= I

¡ ¯uu

uu

uu

H

~

P

i

~

A

i

¡1

~

P

i

=

(I

¡ ¯uu

uu

uu

H

) ~

A

i

¡1

(I

¡ ¯uu

uu

uu

H

)

=

~

A

i

¡1

¡ ¯ ~

A

i

¡1

uu

uu

uu

H

¡ ¯uu

uu

uu

H

~

A

i

¡1

+

¯

2

uu

uu

uu

H

~

A

i

¡1

uu

uu

uu

H

~

P

i

~

A

i

¡1

~

P

i

= ~

A

i

¡1

¡ pu

pu

pu

H

¡ up

up

up

H

+ ¯up

up

up

H

uu

uu

uu

H

= ~

A

i

¡1

¡ uuu

µ

ppp

¡

¯

2

(ppp

H

u

u

u)u

u

u

¶

H

¡

µ

ppp

¡

¯

2

(ppp

H

u

u

u)u

u

u

¶

u

u

u

H

= ~

A

i

¡1

¡ uq

uq

uq

H

¡ qu

qu

qu

H

19

Redukcja macierzy (rzadkiej) hermitowskiej
do postaci trójdiagonalnej metodą Lanczosa

Naszym zadaniem jest znalezienie wartości i
wektorów własnych macierzy stopnia n

Jeśli jednak n jest bardzo duże (np.rzędu ~10

5

) a

nas interesuje jedynie mały wycinek widma
wartości (wektorów) własnych (np. m=50) to

wtedy należy zredukować macierz do postaci
trójdiagonalnej stopnia m

.

W metodzie Lanczosa wykorzystujemy do tego
celu

podprzestrzeń Kryłowa

Zakładamy że wektory rozpinające podprzestrzeń

są

ortonormalne

¸

1

; ¸

2

; : : : ; ¸

n

2 R

x

x

x

1

; x

x

x

2

; : : : ; x

x

x

n

2 C

n

qqq

2 C

n

K

m

(qqq; A)

=

span[qqq; Aqqq; : : : ; A

m

¡1

qqq]

m

¸ 1

K

0

(qqq; A)

=

f000g

dimK

m

=

m

K

m

= span[qqq

1

; qqq

2

; : : : ; qqq

m

]

qqq

H

i

qqq

j

= ±

ij

=

½

1

i = j

0

i

6= j

Metoda Lanczosa jest metodą iteracyjną – w każdej
iteracji poszukujemy nowego wektora bazowego
(wymiar podprzestrzeni zwiększa się o 1)

Sposób postępowania.
Wybieramy dowolny wektor startowy

zakładamy taże

a następnie wykonujemy rekurencyjnie ciąg obliczeń

Procedura zatrzyma się gdy

Wówczas wymiar wygenerowanej podprzestrzeni

qqq

1

2 C

n

;

qqq

1

6= 000

°

1

qqq

0

= 0

Aqqq

i

= °

i

qqq

i

¡1

+ ±

i

qqq

i

+ °

i+1

qqq

i+1

i

¸ 1

±

i

= qqq

H

i

Aqqq

i

°

i+1

=

krrr

i

k

rrr

i

= Aqqq

i

¡ ±

i

qqq

i

¡ °

i

qqq

i

¡1

qqq

i+1

=

rrr

i

°

i+1

, °

i+1

6= 0

±

i

; °

i

2 R

°

i+1

= 0

m = max

i

dim K

i

(qqq; A)

20

Schemat rekurencyjny można zapisać
macierzowo

Gdy warunek

jest spełniony wówczas zachodzi

Q

i

= [qqq

1

; : : : ; qqq

i

]

J

i

=

2
6

6

6

6

6

6

4

±

1

°

2

0

°

2

±

2

. ..

. .. ... °

i

0

°

i

±

i

3
7

7

7

7

7

7

5

AQ

i

=

Q

i

J

i

+ [0; : : : ; 0; °

i+1

q

i+1

]

=

Q

i

J

i

+ °

i+1

qqq

i+1

eee

T

i

i

=

1; 2; : : : ; m

eee

i

= [0; 0; : : : ; 1]

T

2 R

i

°

i+1

= 0

Q

H

i

Q

i

= I

AQ

m

= Q

m

J

m

Wartości własne J

m

są wartościami własnymi A

W skrajnym przypadku metoda nie zatrzyma się sama
i wówczas

Macierz Q

n

jest macierzą unitarną, a macierz J

n

jest

macierzą trójdiagonalną podobną do A.

Uwagi praktyczne:

1)można zredukować liczbę operacji wprowadzając
wektor pomocniczy

wtedy

ponieważ

J

m

zzz = ¸zzz;

zzz

6= 000

m = n

rrr

i

= u

u

u

i

¡ ±

i

qqq

i

±

i

= qqq

H

i

Aqqq

i

= qqq

H

i

u

u

u

i

qqq

H

i

qqq

i

¡1

= 0

u

u

u

i

= Aqqq

i

¡ °

i

qqq

i

¡1

x

x

x = Q

m

zzz

6= 0

Ax

x

x = AQ

m

zzz = Q

m

J

m

zzz = ¸Q

m

zzz = ¸x

x

x

21

2. Jeśli nie interesują nas wektory własne A to
można nie przechowywać wektorów q

i

W każdej iteracji (i->i+1) procedura wymaga obliczenia 5
iloczynów skalarnych oraz 1 przemnożenia wektora przez
macierz A.

3. Jeśli interesuje nas jedynie

Wartości (wektorów) własnych macierzy A to wówczas nie
czekamy aż



i+1

=0 ale przerywamy procedurę wcześniej.

Wykorzystujemy tu fakt, że zbieżność metody dla skrajnych
wartości własnych jest bardzo szybka.

4. Ze względu na błędy zaokrągleń wektory rozpinające
podprzestrzeń Kryłowa przestają być ortogonalne (gdy
rośnie wymiar podprzestrzeni). Konieczne staje się
przeprowadzenie tzw.

reortogonalizacji

nowego wektora

bazy

do już wyznaczonych

Przeprowadzenie reortogonalizacji jest kosztowne.

m << n

~

qqq

i+1

qqq

i+1

= ~

qqq

i+1

¡

i

X

j=1

(qqq

H

j

~

qqq

i+1

)qqq

j

22

Redukcja macierzy do postaci macierzy
Hessenberga (górnej) metodą
eliminacji Gaussa

Naszym zadaniem jest przekształcenie w
sposób rekurencyjny macierzy A

do postaci Hessenberga

T – macierz trójdiagonalna
U – macierz trójkątna górna

Uwagi:

1) każda macierz kwadratowa jest

ortogonalnie podobna do macierzy
Hessenberga

2) ponieważ macierz A jest

przekształcana przez podobieństwo
więc przekształcenie macierzy

hermitowskiej

do postaci

Hessenberga prowadzi do uzyskania
macierzy trójdiagonalnej ze względu
na symetrię A. Macierz trójdiagonalna
jest hermitowska.

A = A

0

! A

1

! : : : ! A

n

¡2

= B

H

A

i

= P

¡1

i

A

i

¡1

P

i

W każdej iteracji macierz przekształcenia P

i

konstruujemy z dwóch macierzy:

a) macierzy permutacji

macierz do niej odwrotna

¼

rs

=

2
6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

4

1

0

. ..

1

0

1

1

. ..

1

1

0

1

. ..

1

3
7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

5

à r

à s

¼

¡1

rs

= ¼

rs

B

H

=

2
6

6

6

6

4

¤ ¤ ¤ ¤ ¤
¤ ¤ ¤ ¤ ¤
0

¤ ¤ ¤ ¤

0

0

¤ ¤ ¤

0

0

0

¤ ¤

3
7

7

7

7

5

= T + U

23

b) macierzy eliminacji

z warunkiem

Macierz odwrotna G

j

-1

G

j

=

2
6

6

6

6

6

6

6

6

6

6

6

6

4

1

. ..

1

l

j+1;j

1

..

.

. ..

l

nj

1

3
7

7

7

7

7

7

7

7

7

7

7

7

5

l

ij

· 1

G

¡1

j

=

2
6

6

6

6

6

6

6

6

6

6

6

6

4

1

. ..

1

¡l

j+1;j

1

..

.

. ..

¡l

nj

1

3
7

7

7

7

7

7

7

7

7

7

7

7

5

Jakie są skutki mnożenia macierzy G

j

oraz



rs

z A?

1. - zamienia wiersze r i s w A

2. - zamienia kolumny r i s w A

2. - odejmuje wiersz j-ty
przemnożony przez l

rj

od

wiersza r w A

3. - przemnaża kolumnę r-tą przez
l

rj

a następnie dodaje do

kolumny j-tej w A

¼

¡1

rs

A

G

¡1

j

A

AG

j

A¼

rs

24

Załóżmy, że w macierzy A

i-1

(i-1) pierwszych

kolumn posiada postać Hessenberga

A

i

¡1

=

2
6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

6

4

¤ ¢ ¢ ¢ ¢ ¤ ¤ ¢ ¢ ¢ ¤

¤ ¢

¢

¢

¢

¢ ¢

¢

¢

¢

¢ ¢

¢

¢

¢

¢ ¢

¢

¢

¢

0

¤ ¤ ¤ ¢ ¢ ¢ ¤

0

¢ ¢ ¢ 0 ¤ ¤ ¢ ¢ ¢ ¤

¢

¢

¢

¢

¢

¢

0

¤ ¢ ¢ ¢ ¤

3
7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

7

5

=

2
6

6

6

4

B

i

¡1

d

d

d

b

A

i

¡1

ccc

±

i

bbb

0

a

a

a

~

A

i

¡1

3
7

7

7

5

a

a

a =

2
6

6

6

4

®

i+1;i

..

.

®

ni

3
7

7

7

5

25

Dokonujemy przekształcenia

macierz przekształcenia jest zdefiniowana w
poniższy sposób

Macierz przekształcenia P

i

nie zmienia postaci

Hessenberga w (i-1) pierwszych kolumnach.
Zmienia natomiast kolumnę i-tą tj. zeruje
elementy w tej kolumnie od i+2 do n.

A

i

= P

¡1

i

A

i

¡1

P

i

P

¡1

i

= G

¡1

i+1

¼

¡1

r;i+1

P

i

= ¼

r;i+1

G

i+1

r

¸ i + 1

¹

a

a

a =

2

6

6

6

4

¤
0

..

.

0

3

7

7

7

5

Aby określić które wiersze należy zamienić wierszami,
najpierw trzeba określić element o największym module
w a

Po wyznaczeniu wartości r zamieniamy miejscami
wiersze (r ↔ i+1) oraz kolumny (r ↔ i+1) w A

i-1

Dzięki temu element o największym module w wektorze
a znajduje się teraz na jego pierwszym miejscu.

Jak określić G

j+1

?

Warunek: trzeba wyzerować wszystkie elementy w a
poza jego pierwszym elementem

Przeprowadzając drugi etap przekształcenia

Dostajemy macierz Hessenberga w i kolumnach
macierzy A

i

.

j®

ri

j = max

i+1

·j·n

j®

ji

j;

r

¸ i + 1

P

¡1

i

¢

2

6

6

6

4

d

d

d

±

i

a

a

a

3

7

7

7

5

=

2

6

6

6

4

d

d

d

±

i

¹

a

a

a

3

7

7

7

5

A

0

= ¼

¡1

r;i+1

A

i

¡1

¼

r;i+1

= [®

0

jk

]

l

j;i+1

=

8

<

:

®

0
ji

®

0
i+1;i

®

i+1;i

6= 0

0

®

i+1;i

= 0

j

=

i + 2; i + 3; : : : ; n

A

i

= G

¡1

i+1

A

0

G

i+1

= P

¡1

i

A

¡1

P

i

26

Uwagi:

1) Po n-2 iteracjach macierz A zostaje

przekształcona do postaci Hessenberga
(górnej).

2) Metoda eliminacji Gaussa jest stabilna

numerycznie i wymaga wykonania

operacji mnożenia

3) Redukcję macierzy do postaci Hessenberga

można przeprowadzić także przy użyciu
metody Hauseholdera lub Givensa
W metodzie Givensa trzeba wykonać

a w metodzie Hausholdera

operacji mnożenia
Ze względu na błędy numeryczne, w obu
metodach macierz Hessenberga jest
macierzą podobną do A+E

M =

2
3

n

3

+ O(n

2

)

M =

10

3

n

3

+ O(n

2

)

M =

5
3

n

3

+ O(n

2

)

kE

Givens

k · k

1

"n

3=2

kAk;

k

1

» 1

kE

Hauseholder

k · k

2

"n

2

kAk;

k

2

» 1

27

Wyznaczanie wartości własnych
macierzy

Wyznaczanie wartości własnych
macierzy trójdiagonalnej hermitowskiej
metodą bisekcji

Po zredukowaniu macierzy A do postaci

chcemy znaleźć sposób na wyznaczenie
wartości własnych J.
Jeśli warunek

to macierz jest

nieredukowalna

.

W przeciwnym wypadku można ją zapisać w
postaci

J =

2
6

6

6

6

6

6

4

±

1

¹

°

2

0

°

2

. .. ...
. .. ... ¹°

n

0

°

n

±

n

3
7

7

7

7

7

7

5

°

i

6= 0; i = 2; : : : ; n

Macierze

są nieredukowalne. Ponadto widmo wartości własnych J
pokrywa się z widmem macierzy nieredukowalnych J

i

,

można więc badać je osobno.
Szukamy wielomianu charakterystycznego J

i

: W(

)

rozwijając wyznacznik względem kolejnych kolumn
macierzy

Macierz jest hermitowska więc

Do znalezienia wartości własnych można wykorzystać
metodę bisekcji.

J

i

; k = 1; : : : ; k

J =

2
6

6

6

6

6

6

4

J

1

0

J

2

. ..

0

J

k

3
7

7

7

7

7

7

5

±

i

;

j°

i

j

2

2 R

!

i

(¸) = det(J

i

¡ ¸I)

!

0

(¸)

=

1

!

1

(¸)

=

±

1

¡ ¸

: : :

: : :

: : :

!

i

(¸)

=

(±

i

¡ ¸)!

i

¡1

(¸)

¡ j°

2

i

j!

i

¡2

(¸)

i

=

2; 3; : : : ; n

W (¸)

=

!

n

(¸)

28

Sposób postępowania:

1) wybieramy dowolną liczbę

 i obliczamy

wartość wielomianu charakterystycznego
rekurencyjnie

2) następnie korzystamy z poniższych twierdzeń:

Tw. Jeżeli elementy





,



3

, ...,



n

(pozadiagonalne) są

niezerowe, to wartości własne macierzy J są
pojedyncze.

Tw. Jeżeli elementy



2

,



3

, ...,



n

(pozadiagonalne) są

niezerowe, to ciąg wartości









 

…

n



spełnia warunki:

a) jeżeli



i

 dla pewnego i<n, to



i-1



i+1



b) jeżeli



n

 jest różne od 0, to liczba zmian

znaków sąsiednich liczb



0



1



n



jest równa liczbie wartości własnych macierzy J
mniejszych od



c) Jeżeli



n

, to  jest wartością własną macierzy

J, a ponadto jest tyle wartości własnych
mniejszych niż

 ile nastąpiło zmian znaków w

ciągu



0



1



n-1



Metoda bisekcji jest bardzo dokładna. Wadą jest
uzyskiwanie dużych wartości ciągu:


0



1



n

 jeśli  znacznie różni się od

wartości własnych J. Zaletą natomiast możliwość
obliczenia wartości własnej o określonym indeksie
k. Liczba iteracji potrzebna do wyznaczenia



k

wynosi:









 

–

przedział poszukiwań wartości własnej



 

–

dokładność wyznaczenia wartości własnej

Wektory własne
Znając wartość własną macierzy J wektor własny
wyznaczamy według wzorów:

gdzie:



i

– elementy diagonalne J



i

– elementy pozadiagonalne J

IT = log

2

¯

0

¡ ®

0

½

x

1

= 1

x

2

=

¸

¡ ±

1

°

2

x

i+1

=

(¸

¡ ±

i

)x

i

¡ °

i

x

i

¡1

°

i+1

;

i = 2; 3; : : : ; n

¡ 1

29

Algorytm poszukiwania

k-tej

wartości własnych metodą

bisekcji

¯

0

=

kJk

®

0

=

¡kJk

30

Wyznaczanie wartości własnych
metodą LR

W metodzie tej iteracyjnie przekształcamy
macierz A uzyskując ciąg

W którym ostatni element stanowi macierz
trójkątna górna (Hessenberga). Elementy
diagonalne macierzy A

m

stanowią natomiast

ciąg wartości własnych macierzy A czyli

W każdej iteracji wyznaczamy rozkład A

i

na

iloczyn macierzy trójkątnej dolnej L z
jedynkami na diagonali oraz macierzy
trójkątnej górnej R:

Przekształcamy macierz w następujący sposób

Macierze A

i

oraz A

i+1

są podobne

A = A

1

! A

2

! A

3

! : : : ! A

m

A

i

= L

i

R

i

A

i+1

= L

i+1

R

i+1

= R

i

L

i

lim

i

!1

a

(i)
jj

= ¸

j

A

i

= L

i

R

i

;

L

¡1

i

A

i

= R

i

;

A

i

R

¡1

i

= L

i

A

i+1

= R

i

L

i

= L

¡1

i

AL

i

= R

i

AR

¡1

i

Wadą jest to że rozkład L

i

R

i

może nie istnieć lub może

istnieć ale jego znalezienie jest źle uwarunkowane. W
obu przypadkach algorytm przerywa działanie.

Wyznaczanie wartości własnych
metodą QR

Metoda wywodzi się z metody LR, przy czym macierz L
zastąpiono

macierzą ortogonalną

Q przez co metoda

jest stabilna numerycznie.

gdzie: R

i

jest macierzą trójkątna górną, a Q

i

jest

macierzą ortogonalną
W metodzie QR otrzymujemy ciąg macierzy

Macierze te są do siebie podobne więc mają te same
wartości własne. Jeśli m jest duże wówczas
spodziewamy się że na diagonali A

m

będą znajdować się

wartości własne A.

Tw. Jeżeli macierz A jest symetryczna i dodatnio
określona, to algorytm QR generuje ciąg macierzy A

(1

)

,A

(2)

,... zbieżny do macierzy diagonalnej.

Wadą metody QR jest wolna zbieżność dla macierzy
pełnych. Metoda jest szybkozbieżna dla
macierzy trójdiagonalnych i macierzy Hessenberga.

Q

H

Q = I

A

i

= Q

i

R

i

A

1

= A

A = A

1

! A

2

! A

3

! : : : ! A

m

A

i+1

= R

i

Q

i

31

Algorytm QR dla macierzy
Hessenberga.

Macierzą (górną) Hessenberga jest macierz,
którą można zapisać w postaci:

Każda macierz kwadratowa jest
ortogonalnie podobna do macierzy
Hessenberga, więc ma ona to samo widmo
wartości własnych co macierz H.

Wyznaczamy ciąg macierzy:

Wszystkie H

i

są macierzami Hessenberga

wobec czego wszystkie macierze H

i

,

i=1,2,3,.. są podobne. Elementy na
diagonali kolejnych H

i

dążą do wartości

własnych macierzy H

1

=H.

Jeżeli macierz H ma pojedyncze wartości
własne takie, że są pomiędzy nimi pary
sprzężone o identycznych modułach, to nie
wszystkie elementy poddiagonalne dążą
do 0. W granicy otrzymuje się macierz:

H = T + U

H

i

= Q

i

R

i

H

i+1

= R

i

Q

i

;

i = 1; 2; 3; : : :

H

i+1

= Q

T

i

H

i

Q

i

H

1

=

2
6

6

6

6

4

¤ ¤ ¢

¢

¢

¤ ¤ ¢

¢

¢

0

0

¤ ¢

¢

0

0

0

¤ ¤

0

0

0

¤ ¤

3
7

7

7

7

5

gwiazdki - elementy ustalone
kropki – elementy, których wartości mogą się
zmieniać.

Wartości własne leżą wtedy bezpośrednio na
diagonali lub są wartościami własnymi macierzy 2x2
leżących na diagonali.

Wadą metody może być powolna zbieżność.
Zwiększenie wydajności uzyskuje się dokonując
przesunięcia widma wartości własnych

Wartość k

i

wybiera się jako równe otrzymanym już

wartościom własnym

lub wartości k

i

oraz k

i+1

jako równe wartościom

własnym macierzy 2x2 znajdujących się w prawym
dolnym rogu macierzy H

i

.

H

i

¡ k

i

I = Q

i

R

i

H

i+1

= R

i

Q

i

+ k

i

I

k

i

= a

i

nn

32

Metoda Hauseholdera rozkładu QR

Definiujemy macierz Hauseholdera w
postaci

która ma własność

Macierz Hausholdera posłuży do znalezienia
rozkładu QR.
Określamy ciąg macierzy P

(1)

,P

(2)

,P

(3)

,...P

(n-1)

przy pomocy których definiujemy macierz R
(górną trójkatną):

¿ =

kzk

2

(

kzk

2

¡ ®)

P

(n

¡1)

P

(n

¡2)

: : : P

(1)

A = R

Zakładamy

Macierz H

(1)

jest macierzą Hauseholdera sprowadzającą

pierwsza kolumnę macierzy A (a

(1)

1

)do postaci:

Przez P

(2)

oznaczamy

Macierz H

(2)

sprowadza pierwszą kolumnę macierzy o

wymiarze (n-1)x(n-1) utworzonej z wierszy i kolumn o
numerach 2,3,4,...,n macierzy P

(1)

A (a

(2)

1

) do postaci

Postępujemy w ten sposób otrzymując kolejno P

(3)

,P

(4)

,

itd.

P

(1)

= H

(1)

®

(1)

ka

(1)
1

k

2

[1; 0; 0; : : : ; 0]

T

1

£n

P

(2)

=

2

6

6

4

1

..

.

0

: : :

: : :

: : :

0

..

.

H

(2)

3

7

7

5

®

2

ka

(2)
1

k

2

[1; 0; 0; : : : ; 0]

T

1

£(n¡1)

H = I

¡

1

¿

u

u

uu

u

u

H

u

u

u = zzz

¡ ®kzzzk

2

eee

1

zzz = [z

1

; z

2

; : : : ; z

n

]

T

eee

1

= [1; 0; : : : ; 0]

T

Hzzz = ®

kzzzk

2

[1; 0; 0; : : : ; 0]

T

® =

§1 = ¡sign(z

1

)

33

Macierz P

(n-1)

ma postać

A macierz H

(n-1)

ma wymiar 2x2.

Po wyznaczeniu wszystkich macierzy P

(i)

,

rozkład A=QR wyznaczamy według
wzorów

Liczba operacji potrzebnych do uzyskania
rozkładu Hauseholdera wynosi

a) mnożenie

b) dodawanie

c) pierwiastkowanie

P

(n

¡1)

=

2

6

6

4

1

1

1

H

(n

¡1)

3

7

7

5

Q = P

(1)

P

(2)

P

(3)

: : : P

(n

¡1)

R = Q

T

A

M =

2
3

n

3

+ 2n

2

¡

2
3

n

¡ 2

D =

2
3

n

3

+

10

3

n

¡ 4

p = n

¡ 1

Uwagi:

1) inne metody pozwalające uzyskać rozkład QR to

metoda Grama-Schmidta i Givensa

2) jeśli macierz A jest symetryczna to rozkład QR

zachowuje jej symetrię natomiast rozkład LR nie

3) jeśli macierz A jest symetryczna i dodatniookreślona

to można stosować rozkład LL

T

wówczas algorytm LR

zachowuje symetrię macierzy A

4) metoda szybka dla macierzy Hessenberga i

trójdiagonalnej

Wyznaczanie wektorów własnych dla rozkładu QR

Uwaga:

nie zakładamy postaci macierzy pierwotnej A

W metodzie QR dostajemy ciąg przekształceń

co można uogólnić

A

i

= Q

i

R

i

! Q

¡1

i

A

i

= R

i

A

i+1

= R

i

Q

i

= Q

¡1

i

AQ

i

= Q

i+1

R

i+1

A

i+2

= R

i+1

Q

i+1

= Q

¡1

i+1

Q

¡1

i

AQ

i

Q

i+1

A

k+1

= Q

¡1

k

Q

¡1

k

¡1

: : : Q

¡1

1

AQ

1

Q

2

: : : Q

k

Q

¡1

i+1

Q

¡1

i

AQ

i

= R

i+1

34

Oznaczmy

Wówczas można bezpośrednio pokazać że

Macierz A

k+1

jest macierzą trójkątną górną z

wartościami własnymi na diagonali.
Wektor własny x

(i)

tej macierzy (

H

)

przynależny do wartości własnej

wyznacza się według poniższych wzorów

P = P

k

= Q

1

Q

2

: : : Q

k

P

¡1

= P

¡1

k

= Q

¡1

k

Q

¡1

k

¡1

: : : Q

¡1

1

P

¡1

AP = A

k+1

= H

¸

i

= h

ii

x

(i)
j

=

0;

j = n; n

¡ 1; : : : ; i + 1

x

(i)
i

=

1

x

(i)
j

=

¡

P

i
k=j+1

h

jk

x

(i)
k

h

jj

¡ h

ii

;

j = i

¡ 1; i ¡ 2; : : : ; 1

Mając wyznaczone wektory własne macierzy H możemy
obliczyć wektory własne macierzy pierwotnej A

P

¡1

AP = H

Hx

x

x = ¸x

x

x

P

¡1

APx

x

x = Hx

x

x = ¸x

x

x

A(Px

x

x) = ¸Px

x

x

yyy = Px

x

x

Ayyy = ¸yyy

35

Uogólniony problem własny

Uogólniony problem własny definiujemy
następująco

gdzie: A i B są macierzami kwadratowymi.

Najprościej byłoby przekształcić powyższe
równanie tak aby przeprowadzić je do zwykłego
problemu własnego tj.

Problemem pozostaje jednak jak znaleźć B

-1

?

W przypadku, gdy B oraz A są macierzami
symetrycznymi możemy posłużyć się rozkładem
LL

T

(w ogólnym przypadku można skorzystać z

rozkładu LU)

wówczas wykorzystując rozkład LL

T

można

znaleźć macierz podobną do B

-1

A

B = LL

T

BB

¡1

= I = LL

T

(L

T

)

¡1

L

¡1

B

¡1

= (L

T

)

¡1

L

¡1

Dzięki temu przekształceniu, macierz G jest
symetryczna jak A i posiada identyczne widmo
wartości własnych (ale inne wektory własne).

Jak znaleźć G?
Najpierw należy znaleźć macierz F

rozwiązując układ równań

a następnie wyznaczamy G

rozwiązując układ równań

Rozkład LL

T

wymaga wykonania n

3

/6 mnożeń a

wyznaczenie macierzy G (2/3)n

3

. Macierz G jest

symetryczna więc w celu wyznaczenia jej wartości i
wektorów własnych korzystamy z metod
przeznaczonych dla tej klasy macierzy.

Wektory własne macierzy A wyznaczamy
przekształcając wektory macierzy G lub
rozwiązując układ

G = L

¡1

F

LG = F

Ax

x

x = ¸Bx

x

x

B

¡1

Ax

x

x = Cx

x

x = ¸x

x

x

L

T

(B

¡1

A)(L

T

)

¡1

=

L

T

(L

T

)

¡1

L

¡1

A(L

T

)

¡1

=

L

¡1

A(L

T

)

¡1

=

G

F = A(L

T

)

¡1

F L

T

= A =

) LF

T

= A

T

= A

Gyyy = ¸yyy

L

T

x

x

x = yyy