Wyznaczanie wartości i wektorów własnych macierzy (problem własny)

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, przesuwanie widma 9. Rozkład SVD

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

= ¸

x

a

; ¸

; x

2 C

A ¢ x = ¸B ¢ x

B

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)

(¸

+ a

¸

+ : : : + a

)

A = [a

]

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 prawostronną 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.

¸

6= ¸

(¸

6= ¸

)

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

AP = B

Q

Q = I

QQ

= I

Q

AQ = D x x x

x x x

= 0

x

x x

Ax x x

= ¡

A

x x x

¢

x x x

= ¸

x x x

x x x

(¸

¡ ¸

)x x x

x x x

= 0 x

x x

Ax x x

= ¸

x x x

x x x

A = (A

)

) a

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

= 2 6 6 6 6 4

¸

1 0 : : : 0 0 0 ¸

1 : : : 0 0 : : : : : : : : : : : : : : : : : :

0 0 0 : : : ¸

1 0 0 0 : : : 0 ¸

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

¡ ¸I) = (¸

¡ ¸)

P

AP = 2 6 6 4

J

0 : : : 0 0 J

: : : 0 : : : : : : : : : : : :

0 0 : : : J

3 7 7 5

haaa; bbbi = aaa

bbb = [®

; ®

; : : : ; ®

] 2 6 6 6 4

¯

¯

.. .

¯

3 7 7 7 5

a

a a ihbbb = aaabbb

= 2 6 6 6 4

®

®

.. .

®

3 7 7 7 5

[¯

; ¯

; : : : ; ¯

]

= 2 6 6 4

®

¯

®

¯

: : : ®

¯

®

¯

®

¯

: : : ®

¯

: : : : : : : : : : : :

®

¯

®

¯

: : : ®

¯

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

X

j¸j

· kAk

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

; ¸

; : : : g

Sp(A + cI) = f¸

+ c; ¸

+ 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³ obliczeń, oraz n³/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

X

b

¸

= 0

w(A) = A

+

n

X

b

A

= 0

b

; b

; b

; : : : ; b

A

yyy +

n

X

b

A

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

= fz : jz ¡ a

j · X

j=1 j6=i

ja

jg

Sp(A) ½ [

C

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 = 1; 2; 1; E

¸

2 fz : jzj · kAk

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'₃) 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₀

Jeśli_i stanowią wartości własne macierzy

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

Jeśli ₁ jest dominującą wartością własną, oraz wektor v₀ ma składową w kierunku x₁ to wówczas zachodzi

vvv

= X

i=1

a

x x x

x x x

; x x x

; x x x

; : : : ; x x x

Avvv

= X

i=1

a

¸

x x x

vvv

= A

vvv

= X

i=1

a

¸

x x x

j¸

j ¸ j¸

j ¸ j¸

j ¸ : : : ¸ j¸

j

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

Dla dowolnego wektora y nieortogonalnego do x₁. 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/₁)^m ale również od współczynników a_i czyli od wyboru v₀. Jeśli wartość własna o największym module jest zespolona to ciąg nie jest zbieżny.

Jak wyznaczyć wektor własny x₁? Ponieważ

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

m

lim

A

vvv

¸

= a

x x x

¸

= lim

m!1

yyy

vvv

yyy

vvv

x x x

= vvv

jvvv

j vvv

= ¸

"

a

x x x

+ X

i=2

µ ¸

¸

¶

a

x x x

#

vvv

¼ ¸

a

x x x

10 Jeśli wartość własna jest pierwiastkiem

wielokrotnym rówanania charakterystycznego to metoda jest zbieżna bo składnik z ₁ dominuje

Uwaga: problem pojawia się gdy ₁=-_j tj.

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

Przyśpieszanie zbieżności - 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

= A

vvv

= ¸

X

i=1

a

x x x

+

X

¸

a

x x x

vvv

Avvv

= vvv

vvv

= X

i=1

a

¸

x

x x

x x x

= ±

vvv

vvv

= X

i=1

a

¸

vvv

vvv

vvv

vvv

= ¸

+ O

"µ

¸

¸

¶

#

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. ₂

Ponieważ wektory v_m+1 oraz ₁v_m są bliskie więc metoda może być w pewnych przypadkach nieużyteczna.

b) metoda zerowania składowej Znając ₁ możemy zdefiniować wektor

Czyli z v₀ usuwamy składową w kierunku x₁. Wektor w₀ nie ma składowej w kierunku x₁ więc ciąg w_m powinien być zbieżny do ₂ oraz x₂. Ze względu na błędy zaokrągleń jednak usuwa się co pewną ilość iteracji składową x₁ tj.

w w w

= vvv

¡ ¸

vvv

=

X

a

(¸

¡ ¸

)¸

x x x

w w w

= (A ¡ ¸

I)vvv

zzz

= (A ¡ ¸

I)w w w

11 Aby wyznaczyć kolejne wartości własne

korzystamy z wektora

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

c) redukcja macierzy

Tw. Jeśli ₁ jest wartością własną macierzy A i x₁ 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 ₁ która jest zerem.

w

w w

= (A ¡ ¸

I)(A ¡ ¸

I)vvv

W

= A ¡ ¸

x x x

vvv

vvv

x x x

= 1

przykład: redukcja Hotellinga

Za wektor v przyjmujemy lewy wektor własny

przynależny do wartości własnej ₁. 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

W

= A ¡ ¸

x x x

x x x

W

= A

i = 1; 2; : : : ; n ¡ 1

W

= W

¡ ¸

x x x

x x x

12 Redukcja macierzy gęstej do prostszej

postaci - macierzy trójdiagonalnej lub macierzy Hessenberga

Macierz pierwotną przekształcamy iteracyjnie

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

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

! A

! A

! : : : ! A

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

= P

A

P

; i = 1; 2; : : : ; m

B = A

= P

AP P = P

P

: : : P

B = P

AP

¢A = P ¢BP

cond(P ) >> 1

cond(P ) = cond(P

: : : P

)

· cond(P

) : : : cond(P

) jBj · jP

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

j¢Aj · cond(P )j¢Bj j¢Aj

jAj · (cond(P ))

j¢Bj jBj

B + ¢B = P

(A + ¢A)P

Byyy = ¸yyy

x x x = Pyyy = P

: : : P

yyy

Ax x x = ¸x x x

13 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

= A = A

A

= (A

)

a a a

=

2 6 6 6 6 4

®

¢

¢

¢

®

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

rzędu (n-i)

A ~

2

4 J

ccc ccc

±

3 5 =

2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4

±

° ¹

¢ ¢ ¢ 0 0

°

. .. ... ¢ ¢

¢ . .. ... ... ¢ ¢

¢ . .. ... ... ¢ ¢

¢ . .. ... ¹°

0 0 ¢ ¢ ¢ °

±

° ¹

0 ¢ ¢ ¢ 0 °

±

3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 A

= P

A

P

P

= P

= P

= I ¡ ¯

u u u

u u u

A

= 2 6 6 6 4

J

ccc 0 ccc

± a a a

0 a a a

A ~

3 7 7 7 5

= [®

]

14 Trzeba utworzyć macierz Hauseholdera

rzędu (n-i)

taką aby transformowała wektor a_i

Macierz Hauseholdera konstruujemy następująco

¯ = 8 <

:

1

¾(¾+j®i+1;ij)

¾ 6= 0

0 ¾ = 0

¾ = kaaa

k

= v u u t

X

j®

j

k = ¡¾e

() ®

= e

j®

j

u u u =

2 6 6 6 6 6 6 4

e

(¾ + j®

j)

®

.. .

®

3 7 7 7 7 7 7 5

Dokonujemy transformacji

P ~

P ~

a a a

= k ¢ eee

P ~

= I ¡ ¯uu uu uu

P

A

P

= P

A

P

=

= 2 6 6 6 4

J

ccc 0 ccc

±

a a a

P ~

0 P ~

a a a

P ~

A ~

P ~

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

±

° ¹

0 0

°

. .. ... .. . 0

. .. ... ¹°

0 0 °

±

° ¹

0 : : : 0 °

±

° ~

0 : : : 0

°

0

0 .. . P ~

A ~

P ~

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

=

2 6 6 6 4

I 0 0 0 1 0 0 0 P ~

3 7 7 7 5

gi ¡ 1

15 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

°

= k

p pp = ¯ ~ A

u u u qqq = pp p ¡ ¯

2 (ppp

u u u)u u u

¯ ¸ 0

p pp

u u u = ¯u u u

A ~

u u u = (ppp

u u u)

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

Metoda Hauseholdera jest stabilna numerycznie.

B = A

= 2 6 6 6 6 4

±

° ~

¢ ¢ ¢ 0

°

. .. ... ...

.. . . .. ... ~°

0 ¢ ¢ ¢ °

±

3 7 7 7 7 5 P ~

= I ¡ ¯uu uu uu

P ~

A ~

P ~

= (I ¡ ¯uu uu uu

) ~ A

(I ¡ ¯uu uu uu

)

= A ~

¡ ¯ ~ A

uu uu uu

¡ ¯uu uu uu

A ~

+ ¯

uu uu uu

A ~

uu uu uu

P ~

A ~

P ~

= ~ A

¡ pu pu pu

¡ up up up

+ ¯up up up

uu uu uu

= ~ A

¡ uuu µ

ppp ¡ ¯

2 (ppp

u u u)u u u

¶

¡ µ

ppp ¡ ¯

2 (ppp

u u u)u u u

¶ u u u

= ~ A

¡ uq uq uq

¡ qu qu qu

16 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⁵ ) 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

¸

; ¸

; : : : ; ¸

2 R x x x

; x x x

; : : : ; x x x

2 C

qqq 2 C

K

(qqq; A) = span[qqq; Aqqq; : : : ; A

qqq]

m ¸ 1

K

(qqq; A) = f000g dimK

= m

K

= span[qqq

; qqq

; : : : ; qqq

] qqq

qqq

= ±

=

½ 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ń korzystając z relacji:

Procedura zatrzyma się gdy

Wówczas wymiar wygenerowanej podprzestrzeni

qqq

2 C

; qqq

6= 000

°

qqq

= 0

Aqqq

= °

qqq

+ ±

qqq

+ °

qqq

i ¸ 1

±

= qqq

Aqqq

°

= krrr

k rrr

= Aqqq

¡ ±

qqq

¡ °

qqq

qqq

= rrr

°

, °

6= 0

±

; °

2 R

°

= 0

m = max

i

dim K

(qqq; A)

17 Schemat rekurencyjny można zapisać

macierzowo

Gdy warunek

jest spełniony wówczas zachodzi

Q

= [qqq

; : : : ; qqq

]

J

= 2 6 6 6 6 6 6 4

±

°

0

°

±

. ..

. .. ... °

0 °

±

3 7 7 7 7 7 7 5

AQ

= Q

J

+ [0; : : : ; 0; °

q

]

= Q

J

+ °

qqq

eee

i = 1; 2; : : : ; m

eee

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

2 R

°

= 0

Q

Q

= I

AQ

= Q

J

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

zzz = ¸zzz; zzz 6= 000

m = n

rrr

= u u u

¡ ±

qqq

±

= qqq

Aqqq

= qqq

u u u

qqq

qqq

= 0 u

u u

= Aqqq

¡ °

qqq

x x x = Q

zzz 6= 0

Ax x x = AQ

zzz = Q

J

zzz = ¸Q

zzz = ¸x x x

18 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. ortogonalizacji nowego wektora bazy

do już wyznaczonych

Przeprowadzenie reortogonalizacji jest kosztowne.

m << n

qqq ~

qqq

= ~ qqq

¡ X

(qqq

qqq ~

)qqq

19 Przesuwanie widma (shift-invert)

Ax = ¸Bx

(A ¡ ¾B)x = (¸ ¡ ¾)Bx

y = Bx

(A ¡ ¾B)

¢ = (A ¡ ¾B)B

y = (¸ ¡ ¾)y

B ¢ = 1

¸ ¡ ¾ y = B(A ¡ ¾B)

y

Cy = ~ ¸y

C = B(A ¡ ¾B)

¸ = ~ 1

¸ ¡ ¾

np. zastosowanie w metodzie Lanczosa Iteracyjnie wyznaczamy kolejne wektory q⁽ⁱ⁾ czyli:

Najpierw szukamy wektora w⁽ⁱ⁾

W tym celu należy rozwiązać układ równań

np. przy użyciu CG lub ILU, wtedy

Rozwiązanie układu w każdej iteracji jest tanie dzieki czemu zastosowanie modyfikacji

„shift-invert” pozwala zmniejszyć ilość iteracji około 10-krotnie.

Cq

= q

B(A ¡ ¾B)

q

= q

(A ¡ ¾B)

q

= w

q

= (A ¡ ¾B)w

Bw

= q

wprowadzamy oznaczenie

i do rozwiązania pozostaje nam standardowy problem

¾ 6= ¸ ¡ shift/przesuniecie

= : (¸ ¡ ¾)

20 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

! A

! : : : ! A

= B

A

= P

A

P

W każdej iteracji macierz przekształcenia P_i konstruujemy z dwóch macierzy:

a) macierzy permutacji

macierz do niej odwrotna

¼

= 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

¼

= ¼

B

=

2 6 6 6 6 4

¤ ¤ ¤ ¤ ¤

¤ ¤ ¤ ¤ ¤ 0 ¤ ¤ ¤ ¤ 0 0 ¤ ¤ ¤ 0 0 0 ¤ ¤

3 7 7 7 7

5 = T + U

21 b) macierzy eliminacji

z warunkiem

Macierz odwrotna G_j^-1

G

= 2 6 6 6 6 6 6 6 6 6 6 6 6 4

1

. ..

1

l

1

.. . . ..

l

1

3 7 7 7 7 7 7 7 7 7 7 7 7 5

l

· 1

G

= 2 6 6 6 6 6 6 6 6 6 6 6 6 4

1

. ..

1

¡l

1

.. . . ..

¡l

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

¼

A

G

A

AG

A¼

22 Załóżmy, że w macierzy A_i-1 (i-1) pierwszych

kolumn posiada postać Hessenberga

A

= 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

d d d A b

ccc ±

bbb 0 a a a A ~

3 7 7

7 5 a a a = 2 6 6 6 4

®

.. .

®

3

7 7

7 5

23 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

= P

A

P

P

= G

¼

P

= ¼

G

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®

j = max

i+1·j·n

j®

j; r ¸ i + 1

P

¢ 2 6 6 6 4

d d d

±

a a a

3 7 7 7 5 =

2 6 6 6 4

d d d

±

a ¹ a a

3 7 7 7 5

A

= ¼

A

¼

= [®

]

l

=

8 <

:

®⁰_ji

®⁰_i+1;i

®

6= 0 0 ®

= 0 j = i + 2; i + 3; : : : ; n

A

= G

A

G

= P

A

P

24 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

+ O(n

)

M = 10

3 n

+ O(n

)

M = 5

3 n

+ O(n

)

kE

k · k

"n

kAk; k

» 1

kE

k · k

"n

kAk; k

» 1

25 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

±

° ¹

0

°

. .. ...

. .. ... ¹°

0 °

±

3 7 7 7 7 7 7 5

°

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

; k = 1; : : : ; k

J = 2 6 6 6 6 6 6 4

J

0

J

. ..

0 J

3 7 7 7 7 7 7 5

±

; j°

j

2 R

!

(¸) = det(J

¡ ¸I)

!

(¸) = 1

!

(¸) = ±

¡ ¸ : : : : : : : : :

!

(¸) = (±

¡ ¸)!

(¸) ¡ j°

j!

(¸) i = 2; 3; : : : ; n

W (¸) = !

(¸)

26 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 _, ₃, ..., _n (pozadiagonalne) są niezerowe, to wartości własne macierzy J są

pojedyncze.

Tw. Jeżeli elementy ₂, ₃, ..., _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

₀₁_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 ₀₁_n-1

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

₀₁_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

¯

¡ ®

½

x

= 1

x

= ¸ ¡ ±

°

x

= (¸ ¡ ±

)x

¡ °

x

°

; i = 2; 3; : : : ; n ¡ 1

27 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

! A

! A

! : : : ! A

A

= L

R

A

= L

R

= R

L

lim

a

= ¸

A

= L

R

; L

A

= R

; A

R

= L

A

= R

L

= L

AL

= R

AR

Wadą jest to że rozkład L_iR_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⁽¹⁾,A⁽²⁾,... 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

Q = I A

= Q

R

A

= A

A = A

! A

! A

! : : : ! A

A

= R

Q

28 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₁=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

= Q

R

H

= R

Q

; i = 1; 2; 3; : : :

H

= Q

H

Q

H

= 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

¡ k

I = Q

R

H

= R

Q

+ k

I

k

= a

29 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⁽¹⁾,P⁽²⁾,P⁽³⁾,...P^(n-1) przy pomocy których definiujemy macierz R (górną trójkatną):

¿ = kzk

( kzk

¡ ®)

P

P

: : : P

A = R

Zakładamy

Macierz H⁽¹⁾ jest macierzą Hauseholdera sprowadzającą pierwsza kolumnę macierzy A (a⁽¹⁾₁)do postaci:

Przez P⁽²⁾ oznaczamy

Macierz H⁽²⁾ sprowadza pierwszą kolumnę macierzy o wymiarze (n-1)x(n-1) utworzonej z wierszy i kolumn o numerach 2,3,4,...,n macierzy P⁽¹⁾A (a⁽²⁾₁) do postaci

Postępujemy w ten sposób otrzymując kolejno P⁽³⁾,P⁽⁴⁾, itd.

P

= H

®

ka

k

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

P

= 2 6 6 4

1 .. . 0 : : : : : : : : :

0 .. . H

3 7 7 5

®

ka

k

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

H = I ¡ 1

¿ u u uu u u

u u u = zzz ¡ ®kzzzk

eee

zzz = [z

; z

; : : : ; z

]

eee

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

Hzzz = ® kzzzk

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

® = §1 = ¡sign(z

)

30 Macierz P^(n-1) ma postać

A macierz H^(n-1) ma wymiar 2x2.

Po wyznaczeniu wszystkich macierzy P⁽ⁱ⁾, 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

= 2 6 6 4

1 1

1

H

3 7 7 5

Q = P

P

P

: : : P

R = Q

A

M = 2

3 n

+ 2n

¡ 2

3 n ¡ 2 D = 2

3 n

+ 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

= Q

R

! Q

A

= R

A

= R

Q

= Q

AQ

= Q

R

A

= R

Q

= Q

Q

AQ

Q

A

= Q

Q

: : : Q

AQ

Q

: : : Q

Q

Q

AQ

= R