Zastosowania warto ´sci i wektorów własnych macierzy

Wydział Matematyki i Nauk Informacyjnych Politechnika Warszawska

Metody numeryczne algebry - wykłady

Alicja Smoktunowicz

Metody numeryczne algebry

Zagadnienie własne (wprowadzenie)

Wydział Matematyki i Nauk Informacyjnych

Politechnika Warszawska

Alicja Smoktunowicz

Zastosowania warto ´sci i wektorów własnych macierzy

W wyszukiwarkach internetowych (np. w Google wyznaczenie PageRank– ustalenie wa˙zno´sci strony internetowej).

Przy badaniu bezpiecze ´nstwa samolotów, okr ˛etów, pociagów, mostów etc.

W statystyce matematycznej.

W teorii sterowania i równaniach ró˙zniczkowych przy badaniu stabilno´sci układów.

Przy badaniu zbie˙zno´sci metod iteracyjnych rozwi ˛azywania wielkich układów równa ´n liniowych.

W genetyce.

Ciekawe zastosowania praktyczne podaj ˛a:

[1] B. N. Datta, Numerical Linear Algebra and Applications, SIAM, 2010,

[2] P. Krzy˙zanowski, Obliczenia in˙zynierskie i naukowe. Szybkie, skuteczne, efektowne, PWN 2011.

Macierze symetryczne

Warto´sci własne λ₁, λ₂, . . . , λ_n macierzy symetrycznej A ∈ R^n×n s ˛a rzeczywiste.

Ponadto istnieje baza ortonormalna {q₁, q₂, . . . , q_n} przestrzeni Rⁿ zło˙zona z wektorów własnych A.

Wtedy Aq_i = λ_iq_i dla i = 1, . . . , n, q^T_i q_i = 1 oraz q^T_i q_j = 0 dla i 6= j.

Rozkład spektralny macierzy symetrycznej A:

A = QDQ^T, D = diag(λ₁, . . . , λ_n), gdzie Q = (q₁, . . . , q_n) jest macierz ˛a ortogonaln ˛a.

W pakiecie MATLAB rozkład spektralny macierzy A mo˙zna

Macierz stowarzyszona i pierwiastki wielomianu

Dany jest wielomian w(t) = p₁tⁿ + p₂tⁿ⁻¹ + . . . + p_nt + p_n+1, gdzie p₁ 6= 0.

Ciekawy fakt: pierwiastki t₁, t₂, . . . , t_n wielomianu w(t) s ˛a warto´sciami własnymi macierzy stowarzyszonej C postaci

C =





















−^p_p²

1 −^p_p³

1 −^p_p⁴

1 . . . −^pn+1_p

1

1 0 0 . . . 0

0 1 0 . . . 0

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 0 0 . . . 1





















W pakiecie MATLAB warto´sci własne t₁, t₂, . . . , t_n macierzy C mog ˛a by´c obliczane poleceniem t=eig(C), a pierwiastki

t₁, t₂, . . . , t_n wielomianu w(t) poleceniem t=roots(p).

Jak nie powinno si ˛e oblicza ´c warto ´sci własnych macierzy?

Dla danej macierzy A ∈ C^n×n szukamy liczby λ ∈ C oraz niezerowego wektora x ∈ Cⁿ, dla których zachodzi rowno´s´c Ax = λx. Liczb ˛e λ nazywamy warto´sci ˛a własn ˛a, a x– wektorem własnym macierzy A.

Liczba λ jest warto´sci ˛a własn ˛a macierzy A wtedy i tylko wtedy, gdy λ jest pierwiastkiem wielomianu charakterystycznego

macierzy A, tj.

w(λ) = det(A − λI) = 0.

Istniej ˛a zatem liczby p_i ∈ C takie, ˙ze

w(t) = det(A − tI) = p₁tⁿ + p₂tⁿ⁻¹ + . . . + p_nt + p_n+1. Czy ta posta´c wielomianu w(t) zawsze nadaje si ˛e do

numerycznego obliczania warto´sci własnych macierzy?

Ciekawy przykład (1)

Chcemy obliczy´c warto´sci własne λ_i macierzy diagonalnej A = diag(1, 2, . . . , n).

Jest to bardzo proste zadanie, gdy˙z λ_i = i dla i = 1, . . . , n.

Zobaczmy teraz, jakie wyniki numeryczne otrzymamy, je´sli

najpierw wyznaczymy współczynniki p_i postaci jawnej wielomianu charakterystycznego

w(t) = det(A − tI) = p₁tⁿ + p₂tⁿ⁻¹ + . . . + p_nt + p_n+1, a nast ˛epnie obliczymy numerycznie pierwiastki wielomianu w(t).

Zauwa˙zmy, ˙ze współczynniki p_i s ˛a liczbami całkowitymi i mog ˛a by´c bardzo du˙ze. Na przykład, |p_n+1| = n!.

W literaturze w_n(t) = (t − 1)(t − 2) . . . (t − n) = (−1)ⁿw(t) nazywa si ˛e wielomianem Wilkinsona (lub wielomianem

perfidnym). Czy zasługuje na t ˛e nazw ˛e? Spójrzmy na wyniki testów numerycznych!

Testy dla n=23

Aby obliczy´c warto´sci własne macierzy A = diag(1, 2, . . . , 23), skorzystamy z nast ˛epuj ˛acego kodu w pakiecie MATLAB:

n=23;

z=1:n; %z=[1,2,...,n]

p=poly(z); %p=[p(1),p(2),...,p(n+1)]

r=roots(p);

Precyzja oblicze ´n: ε_M ≈ 2.2e−16.

Oto obliczone przybli˙zenia kilku warto´sci własnych macierzy A:

22.9599 22.2648

20.9035 + 0.9420i 20.9035 - 0.9420i 18.8006 + 1.6438i 18.8006 - 1.6438i

Ciekawy przykład (2)

Chcemy obliczy´c warto´sci własne macierzy jednostkowej I(10 × 10).

Macierz I ma jedn ˛a warto´s´c własn ˛a λ = 1 o krotno´sci 10.

Korzystamy z nast ˛epuj ˛acego kodu w pakiecie MATLAB:

n=10;

z=ones(1,n); %z=[1,1,...,1]

p=poly(z); %p=[p(1),p(2),...,p(n+1)]

r=roots(p);

Precyzja oblicze ´n: ε_M ≈ 2.2e−16.

Oto obliczone 4 przybli˙zenia warto´sci własnej λ = 1:

1.0562 + 0.0187i 1.0562 - 0.0187i 1.0336 + 0.0482i 1.0336 - 0.0482i

Metody numeryczne algebry

Plan wykładów

Wydział Matematyki i Nauk Informacyjnych Politechnika

Warszawska

Alicja Smoktunowicz

Plan

Tematy wykładów.

Literatura w j ˛ezyku polskim.

Literatura w j ˛ezyku angielskim.

Regulamin zaliczenie przedmiotu.

Tematy wykładów

1 Zagadnienie własne (wprowadzenie).

2 Uogólnione zagadnienie własne.

3 Lokalizacja warto´sci własnych macierzy.

4 Metody wyznacznikowe.

5 Metoda pot ˛egowa i jej warianty.

6 Klasyczna metoda Grama–Schmidta.

7 Zmodyfikowana metoda Grama–Schmidta.

8 Kwadratury Gaussa, wielomiany ortogonalne i warto´sci własne macierzy.

9 Uogólnione odwrotno´sci macierzy.

10 Tematy wykładów.

Literatura w j ˛ezyku polskim

1. G. Dahlquist and Å. Björck, Metody numeryczne, PWN, Warszawa, 1987 (wyd. 2).

2. J. i M. Jankowscy (M.Dryja), Przegl ˛ad metod i algorytmów numerycznych, cz. 1 i 2, WNT, Warszawa, 1988 (wyd. 2).

3. A. Ma´ckiewicz, Algorytmy algebry liniowej. Metody bezpo´srednie, Wyd. Politechniki Pozna ´nskiej, Pozna ´n 2002.

4. P. Krzy˙zanowski, Obliczenia in˙zynierskie i naukowe. Szybkie, skuteczne, efektowne, PWN 2011.

5. Praca zbiorowa pod red. J.W ˛asowskiego: ´Cwiczenia

laboratoryjne z metod numerycznych, OWPW, Warszawa 2002.

6. D. Kincaid, W. Cheney: Analiza numeryczna, WNT, Warszawa 2005.

Literatura w j ˛ezyku angielskim

1. B. N. Datta, Numerical Linear Algebra and Applications, PA, SIAM, 2010.

2. Å. Björck, Numerical Methods for Least Squares Problems, SIAM, PA, USA, 1996.

3. J. W. Demmel, Applied Numerical Linear Algebra, SIAM, PA, USA, 1997.

4. N. J. Higham, Accuracy and stability of numerical algorithms, SIAM, PA, USA, 1996.

5. Ilse C. F. Ipsen, Numerical Matrix Analysis: Linear Systems and Least Squares, SIAM, PA, USA, 2009.

6. G. H. Golub and Ch. F. Van Loan, Matrix Computations, Third Edition, The Johns Hopkins University Press, Baltimore and London, 1996.

7. R. S. Varga, Geršgorin and his circles, Springer-Verlag, Berlin,

Regulamin zaliczenia przedmiotu

Warunkiem zaliczenia przedmiotu jest uzyskanie co najmniej 51

punktów na 100 mo˙zliwych.

Z ´cwicze ´n mo˙zna otrzyma´c 50 punktów (2 sprawdziany oceniane po 25 punktów).

Z zaj ˛e´c laboratoryjnych mo˙zna zdoby´c tak˙ze 50 punktów

(aktywno´s´c na zaj ˛eciach laboratoryjnych - 2 pkt., sprawozdanie z projektów w LaTeX-u- po 2 pkt. z ka˙zdego projektu, czyli w sumie 4 pkt. oraz 2 projekty po 22 punkty).

Ostateczna ocena z przedmiotu wynika z sumy punktów uzyskanych z ´cwicze ´n i zaj ˛e´c laboratoryjnych:

51-60pkt. =⇒ dostateczny,

61-70pkt. =⇒ trzy i pół,

71-80pkt. =⇒ dobry,

81-90pkt. =⇒ cztery i pół,

od 91pkt. =⇒ bardzo dobry.

Metody numeryczne algebry

Uogólnione zagadnienie własne

Wydział Matematyki i Nauk Informacyjnych

Politechnika Warszawska

Alicja Smoktunowicz

Podstawy teoretyczne

Podamy podstawowe informacje o uogólnionym zagadnieniu własnym.

Dla danych macierzy A, B ∈ C^n×n szukamy liczby λ ∈ C i niezerowego wektora x ∈ Cⁿ, dla których zachodzi rowno´s´c

Ax = λBx (λ jest uogólnion ˛a warto´sci ˛a własn ˛a, a x– uogólnionym wektorem własnym). taki, ˙ze Ax = λBx.

Takie zagadnienia pojawiaj ˛a si ˛e w wielu in˙zynierskich i

naukowych problemach. Na ogół warto´sci własne λ reprezentuj ˛a pewien rodzaj drga ´n.

W praktyce obliczeniowej (np. w analizie stanów nieustalonych w polu elektromagnetycznym, w projektowaniu rezonatorów

kwarcowych) bardzo istotn ˛a cech ˛a jest konstrukcja algorytmów o dobrych własno´sciach numerycznych.

Je´sli B jest macierz ˛a jednostkow ˛a (B = I), to mamy do czynienia z klasycznym zagadnieniem własnym: Ax = λx.

Uogólnione warto ´sci własne

Zbiór macierzy postaci A − tB, przy dowolnych warto´sciach t, gdzie t ∈ C, nosi nazw ˛e wi ˛azki macierzy przyporz ˛adkowanej parze (A, B) (ang. pencil). Dla uproszczenia równie˙z par ˛e uporz ˛adkowan ˛a (A, B) b ˛edziemy nazywa´c wi ˛azk ˛a macierzy.

Niech spect(A, B) oznacza zbiór wszystkich warto´sci własnych wi ˛azki (A, B), tzn.

spect(A, B) = {λ ∈ C : det(A − λB) = 0}.

Je´sli detB 6= 0, to zadanie uogólnione Ax = λBx jest równowa˙zne klasycznemu zagadnieniu własnemu B⁻¹Ax = λx. W praktyce obliczeniowej nie nale˙zy jednak stosowa´c tego sposobu, je´sli macierz B jest ´zle uwarunkowana.

Je´sli detB = 0, to spect(A, B) mo˙ze by´c zbiorem pustym, sko ´nczonym lub zawiera´c niesko ´nczenie wiele elementów.

Przykład 1

Niech

A = 2 1 0 1

!

, B = 1 0 0 0

! .

Wtedy

w(λ) = det(A − λB) = 2 − λ, zatem

spect(A, B) = {2}.

Wyniki obliczone w pakiecie MATLAB poleceniem e=eig(A,B): e(1) = 2, e(2) = Inf.

Przykład 2

Niech

A = 2 1 0 0

!

, B = 1 0 0 0

! .

Wtedy

w(λ) = det(A − λB) = 0, zatem

spect(A, B) = C.

Wyniki obliczone w pakiecie MATLAB poleceniem e=eig(A,B): e(1) = 2, e(2) = NaN.

Przykład 3

Niech

A = 1 1 0 1

!

, B = 0 1 0 0

! .

Wtedy

w(λ) = det(A − λB) = 1, zatem

spect(A, B) = Φ.

Wyniki obliczone w pakiecie MATLAB poleceniem e=eig(A,B): e(1) = Inf, e(2) = Inf.

Algorytm QZ

Uogólnione warto´sci własne mo˙zna obliczy´c, stosuj ˛ac algorytm QZ.

Algorytm QZ stanowi uogólnienie algorytmu qr dla klasycznego zagadnienia własnego Ax = λx.

Podstaw ˛a tego algorytmu jest nast ˛epuj ˛ace twierdzenie, stanowi ˛ace uogólnienie twierdzenia Schura.

Twierdzenie 1. Dla dowolnych macierzy A, B ∈ C^{n×n istniej ˛a macierze}

unitarne Q i Z takie, ˙ze

Q^HAZ = S oraz Q^HBZ = T,

gdzie S i T s ˛a macierzami trójk ˛atnymi górnymi.

Wówczas

spect(A, B) = {s /t : t 6= 0}.

Symetria nie wystarcza!

Okazuje si ˛e, ˙ze je´sli macierze rzeczywiste A i B s ˛a symetryczne, to uogólnione warto´sci własne dla zadania Ax = λBx nie musz ˛a by´c rzeczywiste!

Przykład 4.

Niech

A = 1 0

0 −1

!

, B = 0 1 1 0

! . Wtedy

w(λ) = det(A − λB) = −1 − λ², zatem

spect(A, B) = {−i, i}.

Wi ˛ecej o symetrii...

Okazuje si ˛e, ˙ze je´sli zało˙zymy dodatkowo, ˙ze macierz B jest dodatnio okre´slona, to uogólnione warto´sci własne dla zadania Ax = λBx b ˛ed ˛a rzeczywiste.

(Rozkład Cholesky’ego–Banachiewicza).

Je´sli B ∈ R^n×n jest symetryczna i dodatnio okre´slon ˛a, to istnieje macierz trójk ˛atna dolna L ∈ R^n×n z dodatnimi elementami na przek ˛atnej, ˙ze

B = LL^T.

Wtedy zagadnienie własne Ax = λBx jest równowa˙zne

klasycznemu zadaniu własnemu Cy = λy, gdzie y = L^Tx, a C jest macierz ˛a symetryczn ˛a postaci

C = L⁻¹A(L⁻¹)^T.

Twierdzenie o symetrii

Twierdzenie 2. Je´sli

A, B ∈ R^n×n s ˛a symetryczne oraz

B jest macierz ˛a dodatnio okre´slon ˛a, to

uogólnione warto´sci własne dla wi ˛azki (A, B) s ˛a rzeczywiste oraz istnieje n liniowo niezale˙znych wektorów własnych {u₁, u₂, . . . , u_n}

tworz ˛acych układ B–ortonormalny w Rⁿ_{, tzn.}

Au_i = λ_iBu_i,

u^T_i Bu_j = δ_i,j.

Metody numeryczne algebry

Lokalizacja warto ´sci własnych macierzy

Wydział Matematyki i Nauk Informacyjnych

Politechnika Warszawska

Alicja Smoktunowicz

O twierdzeniu Gerszgorina

S. A. Gerszgorin (1901–1933) w roku 1931 opublikował wa˙zne twierdzenie o lokalizacji warto´sci własnych macierzy.

Twierdzenie Gerszgorina ma znaczenie teoretyczne (por.

Przykład 4) oraz praktyczne.

Olga Taussky wykorzystywała je w roku 1947 przy badaniu bezpiecze ´nstwa lotów (przy falowaniu skrzydeł samolotów)!

Rys. 1. S.A. Gerszgorin.

Twierdzenie Gerszgorina (1931)

Twierdzenie 1. Je´sli λ∈ C jest warto´sci ˛a własn ˛a macierzy A ∈ C^n×n, to istnieje

i ∈ {1, . . . , n} takie, ˙ze

λ ∈ G_i = {z ∈ C : |z − a_ii| ≤

Xn j=1,j6=i

|a_ij|},

gdzie G_i jest i-tym kołem Gerschgorina. Mamy wi ˛ec spect(A) ⊂

n

[

i=1

G_i = G(A). (1)

Ponadto, je´sli s kół tworzy zbiór rozł ˛aczny z innymi kołami, to w zbiorze tym znajduje si ˛e dokładnie s warto´sci własnych macierzy A.

Przykład 1

Pewne koła Gerszgorina mog ˛a nie zawiera´c ˙zadnych warto´sci własnych danej macierzy.

Oto ciekawy przyklad takiej macierzy (por. R. S. Varga, Geršgorin and his circles, Springer-Verlag, Berlin, Germany, 2004).

Warto´sci własne macierzy s ˛a zaznaczone symbolem "x".

A₁ =

























0 4 0 0 0 0 0

1 2 0 0 0 0 0

0 1 −2 0 0 0 0

0 0 1/8 −i 1/8 0 0

0 0 0 1/4 −2 ∗ i 1/4 0

0 0 0 0 0 9/2 1/2

0 0 0 0 0 1/2 −9/2

























Przykład 1

−5 −4 −3 −2 −1 0 1 2 3 4 5

−5

−4

−3

−2

−1 0 1 2 3 4 5

Rys. 2. Koła Gerszgorina dla macierzy A₁.

Przykład 2 (macierz magiczna)

A₂ = magic(4) =











16 2 3 13

5 11 10 8

9 7 6 12

4 14 15 1











−30 −20 −10 0 10 20 30 40

−30

−20

−10 0 10 20 30

Rys. 3. Koła Gerszgorina macierzy A₂.

Przykład 3

A₃ =







7 0 1

i/2 5 −i/2

0 1 1







0 1 2 3 4 5 6 7 8

−4

−3

−2

−1 0 1 2 3 4

Rys. 3. Koła Gerszgorina macierzy A₃.

Wska´znik uwarunkowania macierzy

Twierdzenie Gerszgorina mo˙ze by´c u˙zyteczne przy szacowaniu wska´znika uwarunkowania macierzy A ∈ R^n×n.

Jak wiadomo, wska´znik uwarunkowania macierzy (cond(A)) stanowi miar ˛e wra˙zliwo´sci rozwi ˛azania układu równa ´n liniowych Ax = b na małe wzgl ˛edne zaburzenia macierzy A.

Je´sli A jest symetryczna i dodatnio okre´slona, to

cond(A) = kAk₂ kA⁻¹k₂ = λ_max/λ_min, gdzie λ_max jest najwi ˛eksz ˛a warto´sci ˛a własn ˛a A, a λ_min- najmniejsz ˛a warto´sci ˛a własn ˛a.

Podamy przykład takiej macierzy A(n × n).

Przykład 4

Niech

A =











n 1 . . . 1 1 n 1 . . . 1 . . . . . . . . . . . .

1 1 . . . n









 .

Z twierdzenie Gerszgorina wynika, ˙ze wszystkie koła Gerszgorina s ˛a takie same, zatem je´sli λ jest warto´sci ˛a własn ˛a A, to

|λ − n| ≤ n − 1.

Macierz A jest symetryczna, wi ˛ec jej warto´sci własne s ˛a rzeczywiste, zatem 1 ≤ λ ≤ 2n − 1, a wi ˛ec A jest dodatnio okre´slona. Mamy wi ˛ec

cond(A) = λ_max/λ_min ≤ 2n − 1.

Kryterium residualne

Twierdzenie 2. Je´sli A ∈ R^n×n jest macierz ˛a symetryczn ˛a o warto´sciach własnych

λ₁, . . . ,λ_n, wtedy dla dowolnej liczby α ∈ R i dowolnego wektora x ∈ Rⁿ takiego,

˙ze kxk₂ = 1, mamy

i=1,...,nmin |λ_i − α| ≤ kAx − αxk₂.

Przykład 4. Je´sli

A =







4 1 0 1 4 1 0 1 4





 , α = 5, x =

√3

3 (1, 1, 1)^T,

to kAx − αxk₂ =

√ 3

3 ≈ 0.5774. Przybli˙zone warto´sci własne A, obliczone w pakiecie MATLAB, wynosz ˛a 2.5858, 4.0000, 5.4142, zatem

min_i=1,...,n |λ_i − α| ≈ 0.4142.

Metody numeryczne algebry

Metody wyznacznikowe

Wydział Matematyki i Nauk Informacyjnych

Politechnika Warszawska

Alicja Smoktunowicz

Idea metod wyznacznikowych

Nie nale˙zy (na ogół) oblicza´c numerycznie warto´sci własnych macierzy A(n × n) z postaci naturalnej wielomianu

charakterystycznego

w(t) = det(A − tI) = p₁tⁿ + p₂tⁿ⁻¹ + . . . + p_nt + p_n+1. Idea metod wyznacznikowych polega na zastosowaniu metod iteracyjnych wyznaczania zer funkcji w(t) = det(A − tI) (np.

metody siecznych, Newtona, Halleya). Konstruuje si ˛e ci ˛ag {t_k}, który przy pewnych dodatkowych zało˙zeniach jest zbie˙zny do pewnej warto´sci własnej macierzy A. Przy obliczaniu warto´sci w(t_k) w kolejnych krokach iteracyjnych stosuje si ˛e metody na obliczenie wyznaczników macierzy A − t_kI.

Ze wzgl ˛edu na koszt oblicze ´n metody wyznacznikowe stosuje si ˛e głównie dla macierzy strukturalnych, np. dla macierzy

trójdiagonalnych i prawie trójk ˛atnych (Hessenberga).

Metoda siecznych dla w(t) = det(A − tI)

Dane s ˛ a dwa przybli˙zenia t

ⁱ t

^.

Zakładamy, ˙ze w(t

) 6= w(t

) . Dla k = 1, 2, . . . oblicz

t

= t

− w(t

)(t

− t

)

w(t

) − w(t

) .

Metoda stycznych (Newtona) dla

Dane jest przybli˙zenie t

^.

Zakładamy, ˙ze w

(t

) 6= 0 ^.

Dla k = 1, 2, . . . oblicz

t

= t

− w(t

)

w

(t

) .

Algorytm dla macierzy trójdiagonalnej

Jak obliczy ´c warto ´sci w(t) = det(A − tI) i w^′(t) dla danego t? Zakładamy, ˙ze A(n × n) jest macierz ˛a trójdiagonaln ˛a postaci

A =















a₁ b₁ 0 . . . 0 c₁ a₂ b₂ . . . 0 . . . . . . . . . . . . 0

0 0 c_n−2 a_n−1 b_n−1

0 0 0 c_n−1 a_n













 .

Niech w_i = w_i(t) b ˛edzie i. minorem głównym macierzy A − tI. Algorytm: w₀ = 1, w₁ = a₁ − t,

w_i = (a_i − t)w_i−1 − b_i−1c_i−1w_i−1, i = 2, 3, . . . , n,

Algorytm dla macierzy trójdiagonalnej

Jak obliczy ´c warto ´s ´c pochodnej w^′(t)?

Niech w_i = w_i(t) b ˛edzie i. minorem głównym macierzy A − tI. Niech p_i = p_i(t) b ˛edzie pochodn ˛a wielomianu w_i(t) (p_i = w_i^′).

w₀ = 1, w₁ = a₁ − t, p₀ = 0, p₁ = −1, Dla i = 2, 3, . . . , n oblicz

w_i = (a_i − t)w_i−1 − b_i−1c_i−1w_i−1,

p_i = (a_i − t)p_i−1 − w_i−1 − b_i−1c_i−1p_i−1. w(t) = det(A − tI) = w_n, w^′(t) = p_n.

Koszt: O(n) operacji arytmetycznych.

Algorytm Hymana dla macierzy Hessenberga

Omówimy teraz algorytm Hymana obliczenia warto´sci wielomianu charakterystycznego w(t) = det(A − tI) i jego pochodnej w^′(t).

Zakładamy, ˙ze A(n × n) jest macierz ˛a Hessenberga postaci

A =





















a₁₁ a₁₂ . . . a_1,n−1 a_1n b₂ a₂₂ . . . a_2,n−1 a_2n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 0 b_n−1 a_n−1,n−1 a_n−1,n

0 0 0 b_n a_nn



















 .

Koszt algorytmu Hymana: O(n) operacji arytmetycznych.

Dekompozycja zadania

Zakładamy, ˙ze b_i 6= 0 dla i = 2, . . . , n.

Je´sli ten warunek nie jest spełniony, to

A − tI = A₁ − tI A₂ 0 A₃ − tI

! , gdzie A₁ i A₃ s ˛a macierzami Hessenberga.

Wtedy

spect(A) = spect(A₁) ∪ spect(A₃), gdy˙z

det(A − tI) = det(A₁ − tI)det(A₃ − tI).

Idea algorytmu Hymana dla n = 3

Wyznacznik macierzy A − tI nie ulegnie zmianie, gdy do ostatniej kolumny A − tI dodamy kombinacj ˛e liniow ˛a pozostałych kolumn.

Po tym przekształceniu mamy w(t) = det(A − tI) = det(A − tI),^˜ gdzie

A − tI =˜ x₁







a₁₁ − t b₂

0





 + x₂







a₁₂ a₂₂ − t

b₃





 +







a₁₃ a₂₃ a₃₃ − t





.

Chcemy tak dobra´c współczynniki x₁ i x₂, aby ostatnia kolumna macierzy A − tI^˜ miała posta´c





 K

0





.

Idea algorytmu Hymana dla n = 3

Wtedy

A − tI =˜







a₁₁ − t a₁₂ K b₂ a₂₂ − t 0

0 b₃ 0





,

zatem

w(t) = det(A − tI) = (−1)^{˜ 3+1}K(b₂b₃).

Idea algorytmu Hymana dla n = 3

St ˛ad dostajemy zale˙zno´sci

a₂₃ + b₂x₁ + x₂(a₂₂ − t) = 0, (a₃₃ − t) + x₂b₃ = 0.

Poniewa˙z b₂ 6= 0 6= b₃, zatem

x₂ = −(a₃₃ − t) b₃ ,

x₁ = −(a₂₃ + x₂(a₂₂ − t))

b₂ .

Widzimy, ˙ze

K = a₁₃ + x₁(a₁₁ − t) + x₂a₁₂,

Metoda Hymana dla dowolnego n

W tym przypadku macierz A − tI jest macierz ˛a Hessenberga postaci

A − tI =





















a₁₁ − t a₁₂ . . . a_1,n−1 a_1n b₂ a₂₂ − t . . . a_2,n−1 a_2n . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

0 0 b_n−1 a_n−1,n−1 − t a_n−1,n

0 0 0 b_n a_nn − t



















 .

Zakładamy, ˙ze b_i 6= 0 dla i = 2, . . . , n.

Metoda Hymana obliczenia

Mamy

x_n−1 = −(a_nn − t) b_n

Dla i = n − 1, n − 2, . . . , 2 wyznaczamy

x_i−1 = −(a_in + x_i(a_ii − t) + P_n

j=i+1 a_ijx_j) b_n

Obliczamy

K = a_1n + x₁(a₁₁ − t) + x₂a₁₂ + . . . + x_n−1a_1,n−1. St ˛ad

w(t) = (−1)ⁿ⁺¹K(b b . . .b ).

Metoda Hymana obliczenia

Chcemy obliczy´c warto´s´c w^′(t), gdzie w(t) = det(A − tI).

Niech y_i = x_i^′ dla i = 1, . . . , n − 1 (pochodne wzgl ˛edem t).

y_n−1 = 1 b_n

Dla i = n − 1, n − 2, . . . , 2 wyznaczamy

y_i−1 = −(y_i(a_ii − t) − x_i + P_n

j=i+1 a_ijy_j) b_n

Obliczamy

K^′ = y₁(a₁₁ − t) − x₁ + y₂a₁₂ + . . . + y_n−1a_1,n−1, a nast ˛epnie

w^′(t) = (−1)ⁿ⁺¹K^′(b₂b₃ . . . b_n).

Metody numeryczne algebry

Metoda pot ˛egowa i jej warianty

Wydział Matematyki i Nauk Informacyjnych

Politechnika Warszawska

Alicja Smoktunowicz

Idea metody pot ˛egowej

Za twórc ˛e metody pot ˛egowej (Power Iteration) jest uwa˙zany Richard Edler Von Mises (1883-1953).

Chcemy wyznaczy´c dominuj ˛ac ˛a warto´s´c własn ˛a λ₁ macierzy A ∈ C^n×n i odpowiadaj ˛acy jej wektor własny p₁.

Niech A ∈ C^n×n ma warto´sci własne λ₁, λ₂, . . . , λ_n.

Mówimy, ˙ze λ₁ jest dominuj ˛ac ˛a warto´sci ˛a własn ˛a A, je´sli

|λ₁| > |λ₂| ≥ . . . ≥ |λ_n|. (1)

Zakładamy, ˙ze A jest macierz ˛a diagonalizowaln ˛a (podobn ˛a do macierzy diagonalnej), tj. istnieje macierz nieosobliwa P taka, ˙ze

A = PDP⁻¹, D = diag(λ₁, λ₂, . . . , λ_n). (2)

Kolumny p₁, p₂, . . . , p_n macierzy P s ˛a zatem wektorami własnymi macierzy A i s ˛a liniowo niezale˙zne.

Idea metody pot ˛egowej

Dany jest wektor pocz ˛atkowy z₀ ∈ Cⁿ. Tworzymy ci ˛ag wektorów {z_k}:

z_k = A^kz₀, k = 1, 2, . . . . Co mo˙zna powiedzie´c o zbie˙zno´sci ci ˛agu {z_k}?

Przedstawmy z₀ w bazie wektorów własnych {p₁, p₂, . . . , p_n}:

z₀ = c₁p₁ + c₂p₂ + . . . + c_np_n. Wtedy dla k = 1, 2, . . . mamy

z_k = c₁A^kp₁ + c₂A^kp₂ + . . . + c_nA^kp_n.

Zakładamy dodatkowo, ˙ze wektor z₀ nie posiada niezerowej

Idea metody pot ˛egowej

Poniewa˙z Ap_i = λ_ip_i, wi ˛ec A^kp_i = λ^k_i p_i, a zatem

z_k = c₁λ₁^kp₁ + c₂λ₂^kp₂ + . . . + c_nλ_n^kp_n, z_k

λ₁^k = c₁p₁ + c₂ λ₂ λ₁

^k

p₂ + . . . + c_n λ_n λ₁

^k p_n.

Poniewa˙z |_λ^λi

1| < 1 dla i = 2, . . . , n, wi ˛ec

k→ ∞lim

 λ_i λ₁

^k

= 0, i = 2, . . . , n.

St ˛ad dostajemy

k→ ∞lim

z_k

λ₁^k = c₁p₁.

Iloraz Rayleigha

W praktyce obliczeniowej normuje si ˛e wektory z_k, aby unikn ˛a´c nadmiaru.

Jak obliczy´c przybli˙zenie warto´sci własnej λ₁?

Zdefiniujmy iloraz Rayleigha dla dowolnego wektora z ∈ Cⁿ takiego, ˙ze kzk₂ = 1:

R(z) = z^∗Az.

Zauwa˙zmy, ˙ze je´sli p₁ jest wektorem własnym A dla λ₁, tj.

Ap₁ = λ₁p₁, kp₁k₂ = 1, to

Praktyczny algorytm metody pot ˛egowej

Dany jest wektor q₀ ∈ Cⁿ taki, ˙ze kq₀k₂ = 1.

Dla k = 1, 2, . . .

oblicz z_k = Aq_k−1, oblicz q_k = z_k/kz_kk₂,

oblicz α_k = R(q_k) = q^∗_kAq_k.

WARUNKI STOPU:

|α_k+1 − α_k| ≤ tol,

k ≤ N- ograniczenie na liczb ˛e iteracji.

Zbie˙zno ´s ´c metody pot ˛egowej

Twierdzenie.

Zakładamy, ˙ze A ∈ C^n×n jest diagonalizowalna (A = PDP⁻¹, D = diag(λ_i)) oraz posiada dominuj ˛ac ˛a warto´s´c własn ˛a λ₁:

|λ₁| > |λ₂| ≥ . . . ≥ |λ_n|. (3)

Zakładamy, ˙ze wektor pocz ˛atkowy q₀ posiada niezerow ˛a

współrz˛edn ˛a c₁ w kierunku wektora własnego p₁ (Ap₁ = λ₁p₁).

Wówczas metoda pot ˛egowa jest zbie˙zna, tj.

k→ ∞lim α_k = λ₁ (4)

oraz

dist{span{q }, span{p }} = O |λ₂|k

.