Algoritmus vlastní hodnoty Jacobiho - Jacobi eigenvalue algorithm

v numerická lineární algebra, Algoritmus vlastní hodnoty Jacobiho je iterační metoda pro výpočet vlastní čísla a vlastní vektory a nemovitý symetrická matice (proces známý jako diagonalizace ). Je pojmenován po Carl Gustav Jacob Jacobi, který tuto metodu poprvé navrhl v roce 1846,^[1] ale široce se začal používat až v padesátých letech s příchodem počítačů.^[2]

Popis

Nechat ${displaystyle S}$ být symetrická matice a ${displaystyle G = G (i, j, heta)}$ být Dává rotační matici. Pak:

{displaystyle S '= GSG ^ {op},}

je symetrický a podobný na ${displaystyle S}$ .

Dále ${displaystyle S ^ {prime}}$ má záznamy:

{displaystyle {egin {aligned} S '_ {ii} & = c ^ {2}, S_ {ii} -2, sc, S_ {ij} + s ^ {2}, S_ {jj} S' _ { jj} & = s ^ {2}, S_ {ii} + 2sc, S_ {ij} + c ^ {2}, S_ {jj} S '_ {ij} & = S' _ {ji} = (c ^ {2} -s ^ {2}), S_ {ij} + sc, (S_ {ii} -S_ {jj}) S '_ {ik} & = S' _ {ki} = c, S_ { ik} -s, S_ {jk} & keq i, j S '_ {jk} & = S' _ {kj} = s, S_ {ik} + c, S_ {jk} & keq i, j S'_ {kl} & = S_ {kl} & k, leq i, jend {zarovnáno}}}

kde ${displaystyle s = sin (heta)}$ a ${displaystyle c = cos (heta)}$ .

Od té doby ${displaystyle G}$ je kolmý, ${displaystyle S}$ a ${displaystyle S ^ {prime}}$ mít stejné Frobeniova norma ${displaystyle || cdot || _ {F}}$ (druhá odmocnina součet druhých mocnin všech složek), ale můžeme si vybrat ${displaystyle heta}$ takhle ${displaystyle S_ {ij} ^ {prime} = 0}$ , v jakém případě ${displaystyle S ^ {prime}}$ má větší součet čtverců na úhlopříčce:

{displaystyle S '_ {ij} = cos (2 heta) S_ {ij} + {frac {1} {2}} sin (2 heta) (S_ {ii} -S_ {jj})}

Nastavte tuto hodnotu na 0 a přeskupte:

{displaystyle an (2 heta) = {frac {2S_ {ij}} {S_ {jj} -S_ {ii}}}}

-li ${displaystyle S_ {jj} = S_ {ii}}$

{displaystyle heta = {frac {pi} {4}}}

Za účelem optimalizace tohoto efektu S_ij by měl být mimo diagonální prvek s největší absolutní hodnotou, která se nazývá pivot.

Metoda vlastní hodnoty Jacobi opakovaně provádí rotace dokud se matice nestane téměř diagonální. Pak jsou prvky v úhlopříčce aproximacemi (skutečných) vlastních čísel S.

Konvergence

Li ${displaystyle p = S_ {kl}}$ je otočným prvkem, tedy podle definice ${displaystyle | S_ {ij} | leq | p |}$ pro ${displaystyle 1leq i, jleq n, ieq j}$ . Nechat ${displaystyle Gamma (S) ^ {2}}$ označte součet čtverců všech off-diagonálních záznamů z ${displaystyle S}$ . Od té doby ${displaystyle S}$ má přesně ${displaystyle 2N: = n (n-1)}$ mimo diagonální prvky, máme ${displaystyle p ^ {2} leq gama (S) ^ {2} leq 2Np ^ {2}}$ nebo ${displaystyle 2p ^ {2} geq Gamma (S) ^ {2} / N}$ . Nyní ${displaystyle Gamma (S ^ {J}) ^ {2} = Gamma (S) ^ {2} -2p ^ {2}}$ . Z toho vyplývá ${displaystyle Gamma (S ^ {J}) ^ {2} leq (1-1 / N) Gamma (S) ^ {2}}$ nebo ${displaystyle Gamma (S ^ {J}) leq (1-1 / N) ^ {1/2} Gamma (S)}$ ,tj. sled Jacobiho rotací konverguje alespoň lineárně o faktor ${displaystyle (1-1 / N) ^ {1/2}}$ na diagonální matici.

Počet ${displaystyle N}$ Jacobiho rotace se nazývá zametání; nechat ${displaystyle S ^ {sigma}}$ označit výsledek. Předchozí odhad přináší výnosy

{displaystyle Gamma (S ^ {sigma}) leq left (1- {frac {1} {N}} ight) ^ {N / 2} Gamma (S)}

,

tj. posloupnost rozmítání konverguje alespoň lineárně s faktorem ≈ ${displaystyle e ^ {1/2}}$ .

Následující výsledek však Schönhage^[3] poskytuje lokálně kvadratickou konvergenci. Za tímto účelem S mít m odlišné vlastní hodnoty ${displaystyle lambda _ {1}, ..., lambda _ {m}}$ s multiplicitami ${displaystyle u _ {1}, ..., u _ {m}}$ a nechte d > 0 je nejmenší vzdálenost dvou různých vlastních čísel. Zavoláme na číslo

{displaystyle N_ {S}: = {frac {n (n-1)} {2}} - součet _ {mu = 1} ^ {m} {frac {1} {2}} u _ {mu} (u _ {mu} -1) leq N}

Jacobi otáčí Schönhage. Li ${displaystyle S ^ {s}}$ pak označuje výsledek

{displaystyle Gamma (S ^ {s}) leq {sqrt {{frac {n} {2}} - 1}} left ({frac {gamma ^ {2}} {d-2gamma}} ight), quad gamma: = Gamma (S)}

.

Konvergence se tak stane kvadratickou, jakmile ${displaystyle Gamma (S) <{frac {d} {2+ {sqrt {{frac {n} {2}} - 1}}}}}$

Náklady

Každou Jacobiho rotaci lze provést v O (n) kroky, když je otočný prvek str je známo. Nicméně hledání str vyžaduje kontrolu všech N ≈ ½ n² mimo diagonální prvky. Můžeme to snížit na O (n) složitost, pokud zavedeme další pole indexů ${displaystyle m_ {1} ,, tečky ,,, m_ {n-1}}$ s majetkem, který ${displaystyle m_ {i}}$ je index největšího prvku v řádku i, (i = 1, …, n - 1) proudu S. Pak indexy pivot (k, l) musí být jedním z dvojic ${displaystyle (i, m_ {i})}$ . Aktualizaci indexového pole lze provést také v O (n) průměrná složitost: Nejprve maximální položka v aktualizovaných řádcích k a l lze najít v O (n) kroky. V ostatních řádcích i, pouze položky ve sloupcích k a l změna. Smyčka přes tyto řádky, pokud ${displaystyle m_ {i}}$ není ani jedno k ani l, stačí porovnat staré maximum v ${displaystyle m_ {i}}$ k novým položkám a aktualizaci ${displaystyle m_ {i}}$ Pokud je potřeba. Li ${displaystyle m_ {i}}$ by se měl rovnat k nebo l a odpovídající položka se během aktualizace snížila, maximum přes řádek i musí být nalezen od nuly v O (n) složitost. K tomu však dojde v průměru pouze jednou za rotaci. Každá rotace má tedy O (n) a jedno zametání O (n³) průměrná složitost, která odpovídá jednomu násobení matic. Navíc ${displaystyle m_ {i}}$ musí být inicializován před zahájením procesu, což lze provést v n² kroky.

Jacobiho metoda obvykle konverguje v rámci numerické přesnosti po malém počtu tažení. Všimněte si, že několik vlastních čísel od té doby snižuje počet iterací ${displaystyle N_ {S}$ .

Algoritmus

Následující algoritmus je popisem Jacobiho metody v matematické notaci. Vypočítá vektor E který obsahuje vlastní čísla a matici E který obsahuje odpovídající vlastní vektory, tj. ${displaystyle e_ {i}}$ je vlastní číslo a sloupec ${displaystyle E_ {i}}$ orthonormální vlastní vektor pro ${displaystyle e_ {i}}$ , i = 1, …, n.

postup jacobi (S ∈ R^n×n; ven E ∈ Rⁿ; ven E ∈ R^n×n)  var    i, k, l, m, Stát ∈ N    s, C, t, str, y, d, r ∈ R    ind ∈ Nⁿ    změněno ∈ Lⁿ  funkce maxind (k ∈ N) ∈ N ! index největšího mimo diagonálního prvku v řádku k    m := k+1    pro i := k+2 na n dělat      -li │S_ki│ > │S_km│ pak m := i endif    konec    vrátit se m  endfunc  postup Aktualizace(k ∈ N; t ∈ R) ! aktualizace e_k a jeho stav    y := E_k; E_k := y+t    -li změněno_k a (y=E_k) pak změněno_k : = false; Stát := Stát−1    elsif (ne změněno_k) a (y≠E_k) pak změněno_k : = true; Stát := Stát+1    endif  endproc  postup točit se(k,l,i,j ∈ N) ! provést rotaci S._ij, S._kl    ┌ ┐    ┌     ┐┌ ┐    │S_kl│    │C  −s││S_kl│    │ │ := │     ││ │    │S_ij│    │s   C││S_ij│    └ ┘    └     ┘└ ┘  endproc  ! init e, E a pole ind se změnily  E := Já; Stát := n  pro k := 1 na n dělat ind_k : = maxind (k); E_k := S_kk; změněno_k : = pravda konec  zatímco Stát≠0 dělat ! další rotace    m := 1 ! najít index (k, l) kontingenčního p    pro k := 2 na n−1 dělat      -li │S_k ind_{_k}│ > │S_m ind_{_m}│ pak m := k endif    konec    k := m; l := ind_m; str := S_kl    ! vypočítat c = cos φ, s = sin φ    y := (E_l−E_k)/2; d := │y│+√(str²+y²)    r := √(str²+d²); C := d/r; s := str/r; t := str²/d    -li y<0 pak s := −s; t := −t endif    S_kl : = 0,0; Aktualizace(k,−t); Aktualizace(l,t)    ! otáčet řádky a sloupce k a l    pro i := 1 na k−1 dělat točit se(i,k,i,l) konec    pro i := k+1 na l−1 dělat točit se(k,i,i,l) konec    pro i := l+1 na n dělat točit se(k,i,l,i) konec    ! otáčet vlastní vektory    pro i := 1 na n dělat      ┌ ┐    ┌     ┐┌ ┐      │E_ik│    │C  −s││E_ik│      │ │ := │     ││ │      │E_il│    │s   C││E_il│      └ ┘    └     ┘└ ┘    konec    ! řádky k, jsem změnil, aktualizovat řádky ind_k, ind_l    ind_k : = maxind (k); ind_l : = maxind (l)  smyčkaendproc

Poznámky

1. Logické pole změněno udržuje stav každého vlastního čísla. Pokud je číselná hodnota ${displaystyle e_ {k}}$ nebo ${displaystyle e_ {l}}$ změny během iterace, odpovídající složka změněno je nastaven na skutečný, jinak Nepravdivé. Celé číslo Stát počítá počet komponent změněno které mají hodnotu skutečný. Iterace se zastaví, jakmile Stát = 0. To znamená, že žádná z aproximací ${displaystyle e_ {1} ,, ... ,, e_ {n}}$ nedávno změnila svou hodnotu, a proto není příliš pravděpodobné, že k tomu dojde, pokud bude iterace pokračovat. Zde se předpokládá, že operace s plovoucí desetinnou čárkou jsou optimálně zaokrouhleny na nejbližší číslo s plovoucí desetinnou čárkou.

2. Horní trojúhelník matice S je zničen, zatímco spodní trojúhelník a úhlopříčka se nezmění. Je tedy možné obnovit S v případě potřeby podle

pro k := 1 na n−1 dělat ! obnovit matici S    pro l := k+1 na n dělat        S_kl := S_lk    koneckonec

3. Vlastní čísla nejsou nutně v sestupném pořadí. Toho lze dosáhnout jednoduchým třídicím algoritmem.

pro k := 1 na n−1 dělat    m := k    pro l := k+1 na n dělat        -li E_l > E_m pak            m := l        endif    konec    -li k ≠ m pak        vyměnit E_m,E_k        vyměnit E_m,E_k    endifkonec

4. Algoritmus se zapisuje pomocí maticového zápisu (1 pole namísto 0).

5. Při implementaci algoritmu musí být část zadaná pomocí maticového zápisu provedena současně.

6. Tato implementace správně nezohledňuje případ, kdy je jedna dimenze nezávislým podprostorem. Například pokud je dána úhlopříčná matice, výše uvedená implementace se nikdy nezastaví, protože se žádný z vlastních čísel nezmění. Proto v reálných implementacích musí být pro tento případ přidána další logika.

Příklad

Nechat ${displaystyle S = {egin {pmatrix} 4 & -30 & 60 & -35 -30 & 300 & -675 & 420 60 & -675 & 1620 & -1050 -35 & 420 & -1050 & 700end {pmatrix}}}$

Pak jacobi po 3 taženích (19 iterací) vytvoří následující vlastní čísla a vlastní vektory:

${displaystyle e_ {1} = 2585.25381092892231}$

${displaystyle E_ {1} = {egin {pmatrix} 0,0291933231647860588 -0,328712055763188997 0,791411145833126331 -0,514552749997152907end {pmatrix}}}$

${displaystyle e_ {2} = 37.1014913651276582}$

${displaystyle E_ {2} = {egin {pmatrix} -0,179186290535454826 0,741917790628453435 -0.100228136947192199 -0.638282528193614892end {pmatrix}}}$

${displaystyle e_ {3} = 1,4780548447781369}$

${displaystyle E_ {3} = {egin {pmatrix} -0,582075699497237650 0,370502185067093058 0,509578634501799626 0,514048272222164294end {pmatrix}}}$

${displaystyle e_ {4} = 0,1666428611718905}$

${displaystyle E_ {4} = {egin {pmatrix} 0,792608291163763585 0,451923120901599794 0,322416398581824992 0,252161169688241933end {pmatrix}}}$

Aplikace pro skutečné symetrické matice

Když jsou známa vlastní čísla (a vlastní vektory) symetrické matice, lze snadno vypočítat následující hodnoty.

Singulární hodnoty: Singulární hodnoty (čtvercové) matice A jsou druhé odmocniny (nezáporných) vlastních čísel z ${displaystyle A ^ {T} A}$ . V případě symetrické matice S máme ${displaystyle S ^ {T} S = S ^ {2}}$ , proto singulární hodnoty S jsou absolutní hodnoty vlastních čísel S

2-norma a spektrální poloměr: 2-norma matice A je norma založená na euklidovském vektorovém normálu, tj. největší hodnota ${displaystyle | Axe | _ {2}}$ když x prochází všemi vektory s ${displaystyle | x | _ {2} = 1}$ . Je to největší singulární hodnota A. V případě symetrické matice se jedná o největší absolutní hodnotu jejích vlastních vektorů a tedy rovnou jejímu spektrálnímu poloměru.

Číslo podmínky: Číslo podmínky nesingulární matice A je definován jako ${displaystyle {mbox {cond}} (A) = | A | _ {2} | A ^ {- 1} | _ {2}}$ . V případě symetrické matice se jedná o absolutní hodnotu kvocientu největšího a nejmenšího vlastního čísla. Matice s velkými čísly podmínek mohou způsobit numericky nestabilní výsledky: malá porucha může mít za následek velké chyby. Hilbertovy matice jsou nejznámější špatně podmíněné matice. Například Hilbertova matice čtvrtého řádu má podmínku 15514, zatímco pro řád 8 je 2,7 × 10⁸.

Hodnost: Matice A má hodnost r pokud ano r sloupce, které jsou lineárně nezávislé, zatímco zbývající sloupce jsou na nich lineárně závislé. Ekvivalentně r je rozměr rozsahuA. Dále je to počet nenulových singulárních hodnot.; V případě symetrické matice r je počet nenulových vlastních čísel. Bohužel kvůli chybám zaokrouhlování nemusí být numerická aproximace nulových vlastních čísel nulová (může se také stát, že numerická aproximace je nula, zatímco skutečná hodnota není). Lze tedy pouze vypočítat numerické pořadí rozhodnutím, které z vlastních čísel jsou dostatečně blízko nule.

Pseudo-inverzní: Pseudo inverze matice A je jedinečná matice ${displaystyle X = A ^ {+}}$ pro který SEKERA a XA jsou symetrické a pro které AXA = A, XAX = X drží. Li A je nesmyslné, pak ' ${displaystyle A ^ {+} = A ^ {- 1}}$ .; Když se volá procedura jacobi (S, e, E), pak relace ${displaystyle S = E ^ {T} {mbox {Diag}} (e) E}$ platí kde Diag (E) označuje diagonální matici s vektorem E na úhlopříčce. Nechat ${displaystyle e ^ {+}}$ označit vektor kde ${displaystyle e_ {i}}$ je nahrazen ${displaystyle 1 / e_ {i}}$ -li ${displaystyle e_ {i} leq 0}$ a o 0, pokud ${displaystyle e_ {i}}$ je (číselně blízký) nule. Protože matice E je ortogonální, z toho vyplývá, že pseudo-inverze S je dána vztahem ${displaystyle S ^ {+} = E ^ {T} {mbox {Diag}} (e ^ {+}) E}$ .

Řešení nejmenších čtverců: Pokud matice A nemá úplnou hodnost, nemusí existovat řešení lineárního systému Sekera = b. Lze však hledat vektor x, pro který ${displaystyle | Ax-b | _ {2}}$ je minimální. Řešení je ${displaystyle x = A ^ {+} b}$ . V případě symetrické matice S jako předtím, jeden má ${displaystyle x = S ^ {+} b = E ^ {T} {mbox {Diag}} (e ^ {+}) Eb}$ .

Exponenciální matice: Z ${displaystyle S = E ^ {T} {mbox {Diag}} (e) E}$ jeden najde ${displaystyle exp S = E ^ {T} {mbox {Diag}} (exp e) E}$ kde expE je vektor kde ${displaystyle e_ {i}}$ je nahrazen ${displaystyle exp e_ {i}}$ . Stejně, F(S) lze vypočítat zjevným způsobem pro jakoukoli (analytickou) funkci F.

Lineární diferenciální rovnice: Diferenciální rovnice X' = Sekera, X(0) = A má řešení X(t) = exp (t A) A. Pro symetrickou matici S, z toho vyplývá, že ${displaystyle x (t) = E ^ {T} {mbox {Diag}} (exp te) Ea}$ . Li ${displaystyle a = součet _ {i = 1} ^ {n} a_ {i} E_ {i}}$ je expanze A vlastními vektory S, pak ${displaystyle x (t) = součet _ {i = 1} ^ {n} a_ {i} exp (te_ {i}) E_ {i}}$ .; Nechat ${displaystyle W ^ {s}}$ být vektorový prostor překlenut vlastními vektory z S které odpovídají zápornému vlastnímu číslu a ${displaystyle W ^ {u}}$ analogicky pro kladná vlastní čísla. Li ${displaystyle ain W ^ {s}}$ pak ${displaystyle {mbox {lim}} _ {t infty} x (t) = 0}$ tj. rovnovážný bod 0 je atraktivní pro X(t). Li ${displaystyle ain W ^ {u}}$ pak ${displaystyle {mbox {lim}} _ {t infty} x (t) = infty}$ , tj. 0 je odpudivé pro X(t). ${displaystyle W ^ {s}}$ a ${displaystyle W ^ {u}}$ jsou nazývány stabilní a nestabilní rozdělovače pro S. Li A má komponenty v obou rozdělovačích, pak je přitahována jedna komponenta a jedna je odpuzována. Proto X(t) přístupy ${displaystyle W ^ {u}}$ tak jako ${displaystyle t o infty}$ .

Zobecnění

Jacobiho metoda byla zobecněna na složité hermitovské matice, obecné nesymetrické reálné a komplexní matice i blokové matice.

Vzhledem k tomu, že singulární hodnoty skutečné matice jsou druhou odmocninou vlastních čísel symetrické matice ${displaystyle S = A ^ {T} A}$ lze jej také použít pro výpočet těchto hodnot. V tomto případě je metoda upravena takovým způsobem, že S nesmí být výslovně vypočítáno, což snižuje nebezpečí chyby zaokrouhlování. Všimněte si, že ${displaystyle JSJ ^ {T} = JA ^ {T} AJ ^ {T} = JA ^ {T} J ^ {T} JAJ ^ {T} = B ^ {T} B}$ s ${displaystyle B,: = JAJ ^ {T}}$ .

Jacobiho metoda je také velmi vhodná pro paralelismus.

Reference

^ Jacobi, C.G.J. (1846). „Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen“. Crelle's Journal (v němčině). 1846 (30): 51–94. doi:10,1515 / crll.1846.30.51.
^ Golub, G.H .; van der Vorst, H.A. (2000). „Výpočet vlastních čísel ve 20. století“. Journal of Computational and Applied Mathematics. 123 (1–2): 35–65. doi:10.1016 / S0377-0427 (00) 00413-1.
^ Schönhage, A. (1964). „Zur quadratischen Konvergenz des Jacobi-Verfahrens“. Numerische Mathematik (v němčině). 6 (1): 410–412. doi:10.1007 / BF01386091. PAN 0174171.

Další čtení

Stiskněte, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007), „Oddíl 11.1. Jacobiho transformace symetrické matice“, Numerické recepty: Umění vědecké práce na počítači (3. vyd.), New York: Cambridge University Press, ISBN 978-0-521-88068-8
Rutishauser, H. (1966). „Příručka Lineární algebra: Jacobiho metoda pro skutečné symetrické matice“. Numerische Mathematik. 9 (1): 1–10. doi:10.1007 / BF02165223. PAN 1553948.
Sameh, A.H. (1971). „Na Jacobiho a Jacobiho podobných algoritmech pro paralelní počítač“. Matematika výpočtu. 25 (115): 579–590. doi:10.1090 / s0025-5718-1971-0297131-6. JSTOR 2005221. PAN 0297131.
Shroff, Gautam M. (1991). "Paralelní algoritmus pro vlastní hodnoty a vlastní vektory obecné komplexní matice". Numerische Mathematik. 58 (1): 779–805. CiteSeerX 10.1.1.134.3566. doi:10.1007 / BF01385654. PAN 1098865.
Veselić, K. (1979). „Na třídu Jacobiho postupů pro diagonalizaci libovolných reálných matic“. Numerische Mathematik. 33 (2): 157–172. doi:10.1007 / BF01399551. PAN 0549446.
Veselić, K .; Wenzel, H. J. (1979). "Kvadraticky konvergentní Jacobiho metoda pro skutečné matice se složitými vlastními hodnotami". Numerische Mathematik. 33 (4): 425–435. doi:10.1007 / BF01399324. PAN 0553351.

externí odkazy

[1] Jacobi, C.G.J. (1846). „Über ein leichtes Verfahren, die in der Theorie der Säkularstörungen vorkommenden Gleichungen numerisch aufzulösen“. Crelle's Journal (v němčině). 1846 (30): 51–94. doi:10,1515 / crll.1846.30.51.

[2] Golub, G.H .; van der Vorst, H.A. (2000). „Výpočet vlastních čísel ve 20. století“. Journal of Computational and Applied Mathematics. 123 (1–2): 35–65. doi:10.1016 / S0377-0427 (00) 00413-1.

[3] Schönhage, A. (1964). „Zur quadratischen Konvergenz des Jacobi-Verfahrens“. Numerische Mathematik (v němčině). 6 (1): 410–412. doi:10.1007 / BF01386091. PAN 0174171.

[1]

[2]

[3]

Numerická lineární algebra
Klíčové koncepty	Plovoucí bod Numerická stabilita
Problémy	Systém lineárních rovnic Maticové rozklady Násobení matic (algoritmy ) Rozdělení matice Řídké problémy
Hardware	Mezipaměť CPU TLB Algoritmus zapomínající na mezipaměť SIMD Multiprocesing
Software	MATLAB Základní podprogramy lineární algebry (BLAS) LAPACK Specializované knihovny Software pro všeobecné účely