2. Domácí úkol: Mocninná metoda
Úvod
Tento úkol se zabývá vlastními čísly a vektory matic. Nebudeme se pro zadanou matici snažit počítat všechna vlastní čísla, ale podíváme se na algoritmus umožňující nalézt tzv. dominantní vlastní číslo (největší v absolutní hodnotě) a pomocí jednoduché modifikace i nejmenší v absolutní hodnotě, resp. nejbližší k zadané komplexní hodnotě.
Konkrétně si zde ozkoušíme tzv. mocninnou metodu (nebo mocninnou iteraci; power method resp. power iteration). Důvod k tomuto pojmenování bude patrný z popisu níže. V aplikacích na tuto metodu můžete narazit například v Google PageRank algoritmu, při hledání základních stavů kvantových systémů, nebo při vylepšování aproximace ke konkrétnímu hledanému vlastnímu číslu matice.
Teoretický background
Hledání dominantního vlastního čísla
Mějme matici , mající dominantní vlastní číslo jednoduché násobnosti. Tj. pro libovolné další vlastní číslo matice platí .
Pro jednoduchost v tento okamžik dále předpokládejme, že matice je diagonalizovatelná s vlastními čísly , pro která platí nerovnosti , a vlastními vektory tvořícími ortonormální bázi prostoru .
Zvolme nyní vektor , typicky náhodně (nebo pomocí nějaké heuristiky plynoucí z daného problému; to v obecné metodě ale neřešíme). Pak ho lze vyjádřit ve tvaru
kde jsou jisté komplexní konstanty a typicky (toto při náhodné volbě vysoce pravděpodobně bude splněno, případně efektivně následně vynuceno numerickými chybami ve výpočtech).
Nyní pojďme zkoumat vlastnosti posloupnosti zadané rekurentně vztahem
kde označuje standardní Euklidovskou normu vektoru v , tj.
Platí tedy vztah (odtud je vidět důvod pro označení "mocninná metoda")
kde působení mocnin matice lze vyjádřit explicitně s využitím linearity a předpokladů výše takto:
Po dosazení a vytknutí dominantní vlastní hodnoty v čitateli i jmenovateli dostaneme
Všimněme si, že členy této posloupnosti se limitně blíží násobkům vlastního vektoru příslušnému vlastnímu číslu , tedy výrazu tvaru . Odtud můžeme získat aproximaci odpovídajícího vlastního čísla pomocí tzv. Rayleighova podílu
kde označuje standardní skalární součin v prostoru .
Pro zastavení výpočtu použijeme dvě podmínky (v konjunkci):
- abslutní hodnota rozdílu dvou posledních aproximací vlastních čísel je menší než zadaná tolerance ,
- norma rezidua, tj. , kde je poslední aproximace vlastního čísla s přibližným vlastním vektorem , je menší než zadaná tolerance .
Nejmenší a nejbližší vlastní číslo
Nyní využijme dvě následující jednoduchá pozorování plynoucí přímo z definice:
- je vlastním číslem matice , právě když je vlastním číslem matice .
- Máme-li regulární matici , pak je vlastním číslem (nutně nenulovým) matice , právě když je vlastním číslem matice .
V obou případech jsou vždy vlastní vektory shodné. Z toho plyne, že aplikujeme-li mocninnou metodu popsanou výše na matici , pro zadané , pak najdeme vlastní číslo matice , které je nejblíže . Ovšem za předpokladu, že takovéto vlastní číslo existuje právě jedno.
Implementační pokyny
V souboru main.jl
doplňte implementaci funkcí maxeig
a mineig
.
První metoda se snaží nalézt dominantní vlastní číslo.
Druhá metoda pak hledá vlastní číslo nejblíže k zadanému číslu.
Níže obě metody podrobněji rozebereme.
Obě metody v případě úspěchu vrací dvojici mající v první složce aproximaci vlastního čísla a v druhé aproximaci příslušného vlastního vektoru.
maxeig
Funkce Na vstupu máme jeden poziční argument
- čtvercovou matici
A
s prvky typuComplex{T}
, kdeT
je podtyp typuAbstractFloat
(pokud není čtvercová, pak vyvoláme výjimkuArgumentError
)
a dále tři parametry
epsilon
: absolutní tolerace pro porovnání aproximace vlastních hodnot (výchozí hodnotasqrt(eps(T))
),norm_epsilon
: absolutní tolerace pro porovnání nulovosti normy rezidua (výchozí hodnotaepsilon * length(A)
),max_iter
: maximální povolený počet iterací (výchozí hodnota je1_000
).
Všechny tyto parametry musí být nezáporné, jinak vyvoláme výjimku ArgumentError
.
Algoritmus nejprve vygeneruje náhodný vektor a postupně napočítává normalizované mocniny A
na tomto vektoru.
Iterace je zastavena, pokud je absolutní hodnota rozdílu dvou posledních aproximací vlastních čísel menší než epsilon
a současně norma rezidua (tj. vektoru , kde je poslední aproximace vlastního čísla a vlastního vektoru) je menší než norm_epsilon
.
V tomto případě metoda vrací dvojici (tuple) složenou z aproximace vlastního čísla a vlastního vektoru.
Pokud k situaci popsané v předchozím odstavci nedojde do max_iter
iterací, pak vyvoláme výjimku ErrorException
.
Pro pohodlí uživatele nechť je dále k dispozici metoda maxeig
akceptující na vstupu matici typu Matrix{T}
, kde T
je podtyp typu Real
.
V tomto případě matici přetypujeme na typ Complex{Float64}
a dále vše běží jako výše.
mineig
Funkce Tato funkce se od maxeig
liší pouze druhým nepovinným pozičním argumentem mu
(výchozí hodnota 0
), který představuje hodnotu, kolem které hledáme nejbližší vlastní číslo matice A
.
Funkce tedy implementuje mocninnou metodu pro matici (A - mu * I)^(-1)
a vrací (aproximaci) vlastního čísla matice A
, které je nejblíže číslu mu
, spolu s odpovídajícím přibližným vlastním vektorem (nebo opět vyvolá výjimku, pokud narazíme na singulární matici či nedosáhneme konvergence).
Kritérium pro zastavení je stejné jako v případě základní mocninné metody, reziduum počítáme vzhledem k původní matici A
a jejím vlastním číslům/vektorů.
Testy
Vaši implementaci prověříte spuštěním testů (automaticky se spouštějí i na Gitlabu).
$ julia test/runtests.jl
┌ Info: Running tests in test_maxeig.jl...
└ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Test Summary: | Pass Total Time
Exceptions. | 5 5 0.9s
Test Summary: | Pass Total Time
Matrices with complex (floating point) coefficients. | 3 3 3.1s
Test Summary: | Pass Total Time
Matrices with real (floating point) coefficients. | 2 2 0.2s
Test Summary: | Pass Total Time
Matrices with rational coefficients. | 2 2 0.3s
Test Summary: | Pass Total Time
Matrices with integer coefficients. | 2 2 0.1s
Test Summary: | Pass Total Time
Matrices with big integer coefficients. | 2 2 0.9s
Test Summary: | Pass Total Time
Matrices with big float coefficients. | 2 2 0.2s
Test Summary: | Pass Total Time
Matrices with non-unique dominant eigenvalue. | 2 2 0.0s
Test Summary: | Pass Total Time
Matrices with degenerate dominant eigenvalue. | 2 2 0.4s
Test Summary: | Pass Total Time
An example of a larger matrix. | 1 1 0.1s
┌ Info: Running tests in test_mineig.jl...
└ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Test Summary: | Pass Total Time
Exceptions. | 7 7 1.5s
Test Summary: | Pass Total Time
Matrices with complex (floating point) coefficients. | 4 4 1.7s
Test Summary: | Pass Total Time
Matrices with real (floating point) coefficients. | 2 2 0.0s
Test Summary: | Pass Total Time
Matrices with rational coefficients. | 1 1 0.2s
Test Summary: | Pass Total Time
Matrices with integer coefficients. | 2 2 0.1s
Test Summary: | Pass Total Time
Matrices with big integer coefficients. | 2 2 0.0s
Test Summary: | Pass Total Time
Matrices with big float coefficients. | 2 2 0.1s
Test Summary: | Pass Total Time
Matrices with non-unique 'smallest dominant' eigenvalue. | 2 2 0.0s
Test Summary: | Pass Total Time
An example of a larger matrix. | 3 3 0.1s
Řešení a odevzdání
Opět vytvořte větev odvozenou z větve assignment/02-power
a nezvěte ji solution/02-power
.
Do solution/02-power
vložte své řešení editací souboru main.jl
případně přidáním testů do složky test
.
Až budete se svým řešením spokojeni, vytvořte MR (to můžete i dříve, aspoň uvidíte výsledek testů, pokud je nespouštíte lokálně) a přiřaďte mě k němu jako assignee
.
Tímto aktem úkol odevzdáte.