Matematické Fórum

MariPecha · 07. 03. 2014 23:19

Dobrý den,

obracím se na Vás s prosbou. Už delší dobou mi vrtá v hlavě otázka, jak může matlab spočítat dostatečně rychle a dostatečné přesně nejmenší vlastní vektory příslušné k nejmenším vlastním číslům obecné matice.

Díval jsem se do kódu matlabovské funkce eigs a zjistil jsem, že se počítají pomocí Shift-Invert Implicity Restarted Arnoldi Method (výpočet matlab provádí pomocí knihovny ARPACK ). Narazil jsem však na problém, kterému nerozumím.

1. Shift-invert mód řešíče arpack je mód, při kterém se posouvá spektru matice:

( A - \lambda*I )^(-1) , kde \lambda je číslo, ke kterému mají být požadovaná vlastní čísla blízké.

se stávají největšímu, čili díky velkému eigen-gapu je zaručeno, že problém bude rychleji konvergovat. Ale teď jsem si položil otázku. Co když chci spočítat k-nejmenších vlastních čísel Laplaceovy matice L, která je čtvercová symetrická positivně semi definitní a SINGULÁRNÍ, nejmenší vlastní číslo je 0.

Takže sigma nastavím na 0 a dostanu:

(L)^(-1)

A teď přichází část které nerozumím.

Ale vím, že matice je singulární a k takové neexistuje inverze. Díval jsem se, že tento problém matlab řeší pomocí LU faktoru matice. Ale stejný problém LU rozklad je definován jenom pro regulární matice a L je singulární. Na internetu jsem četl něco o pseudo-inverzích, které jsou definovány pro obecné matice (regulární singulární, nesymetrické, nečtvercové apod. ). Ale pořád nějak nemohu pochopit, jak to funguje.

Děkuji za odpovědi,
Mari

pospik · 08. 03. 2014 00:30

Zdravím a pozdravuji,

Matlabovským suprřěšičům moc nerozumím, ale o pseudoinverzích něco málo vím.

Představte si, že úkolem je řešit soustavu
Ax = b, kde A je symetrická positivně semidefinitní matice (nějaké to nulové vlastní číslo, čili netriviální jádro) a vektor b je z oboru hodnot matice A.

Soustava má očividně řešení (protože vektor b je z oboru hodnot matice A, tedy existuje takové x, pro které Ax = b - viz definice oboru hodnot). Těchto řešení je nekonečně mnoho a dají se vyjádřit ve tvaru
x = A^{+} b + d
kde d je libovolný prvek jádra matice A (čili Ad = 0) a A^{+} je pseudoinverze matice A, což je matice taková, že
A A^{+} A = A
A^{+} A A^{+} = A^{+}
(doporučuji googlit pod heslem "Moore-Penrose pseudoinverse")
Skutečně dosazením Ax = A (A^{+} b + d) = A A^{+} b + A d = A A^{+} b + 0 = A A^{+} b = .. (předpokládám, že b je z oboru hodnot, čili existuje nějaké y takové, že Ay = b) ... = A A^{+} A y = A y = b

Toliko k motivaci. Doporučuji porovnat s klasickou inverzí.

Symetrické matice mají jednu úžasnou vlastnost - maticové operace lze převést na spektrum, tj. pokud A = U^T D U (spektrální rozklad), pak pro libovolnou maticovou funkci se symetrickou maticí
f(A) = f(U^T D U) = U^T f(D) U

Pokud máte regulární matici, pak jsou vlastní čísla nenulová, tedy inverzní matice bude mít ve spektru převrácené hodnoty vlastních čísel. Se singulárními maticemi je problém v tom, že nulou se nedělí.

Nicméně, u pseudoinverzí se to dá udělat obdobně
A^{+} = U^T D^{+} U
kde aplikuji "pseudoinverzi" na diagonální prvky předpisem
[D^{+}]_ii = ...
... 1/D_ii pokud D_ii > 0
... 0 pokud D_ii = 0 (tím se docílí nedělení nulou)

A ještě jedno pozorování. Pokud hledáte vlastní čísla příslušné nulovým vlastním číslům, tak de-facto řešíte rovnici
A v = lambda v
A v = 0 v
A v = 0
Vlastní vektory jsou lineárně nezávislé.... nehledáte náhodou bázi jádra matice A?

MariPecha · 08. 03. 2014 13:35

↑ pospik:

Ahoj, Lukáši,

čekal jsem, že ty budeš mezi prvními (nakonec jsi jediný), kteří budou reagovat. To co jsi mi napsal je zajímavé, až bude chvilku čas, tak si to budu muset promyslet, ale v podstatě to neřeší moji otázku, resp. aspoň nevidím souvislosti s tím, co jsem včera našel.

Než se pustím do popisování, toho, co jsem našel, nejdříve budu reagovat na Tvou poslední otázku.

A ještě jedno pozorování. Pokud hledáte vlastní čísla příslušné nulovým vlastním číslům, tak de-facto řešíte rovnici
A v = lambda v
A v = 0 v
A v = 0

Ano přesně tohle řeším. Souhlasím s tím, že jsou lineárně nezávislé, čili jak správně píšeš hledám bázi jádra matice. A pokud řeším vlastní problém Laplaceovy matice, tak dokonce vím, že pokud je algebraická násobnost nulového vlastního čísla 1, tak vlastní vektor je konstantní, pokud > 1, tak jsou vlastní vektory po částech konstantní, ale tohle není předmětem mojeho dotazu.

Teď dál.

>> Pokud máte regulární matici, pak jsou vlastní čísla nenulová, tedy inverzní matice bude mít ve spektru převrácené hodnoty vlastních čísel. Se singulárními maticemi je problém v tom, že nulou se nedělí.

V podstatě s Tebou souhlasím, ale v úvaze je jedna malá chybička. Tím, že se otočí spektrum tak mi solver bude vracet největší čísla a u těch přece můžu udělat převrácenou hodnotu. Teoreticky pokud je největší vlastní číslo inf, tak nejmenší je identická 0. :)

Problém je tedy jinde a to je jak získat inverzní matici k singulární. Jelikož jsem fanda Laplaceových matic, zůstaňme u ní. A problém trochu zkontretizuji. Jak získat inverzi k singulární čtvercové. A teď vím pro obecnou čtvecovou matici platí:

P A = L U , P je permutační matice
pro takové řešení je vhodné např. http://www.gnu.org/software/gsl/manual/ … ition.html

A jelikož už mám LU faktor, tak můžu v pohodě spočítat inverzi A. Takže tím je vyřešený první problém. Ještě to úplně nevidím, ale předpokládám, že to co jsem popsal souvisí s tím, co jsi mi popsal ty.

Matlab uvažuje stabilnější rozklad:

P A Q = L U , k tomu používá knihovnu UMFPACK. :)

viz. http://www.mathworks.com/help/matlab/ref/lu.html.

Takže tím mám vyřešeno, jak může být matlabovská funkce eigs tak rychlá.

V nápovědě k lu je ještě jeden argument, kterému nerozumím, co je "row-scaling matrix"?

Matematické Fórum

#1 07. 03. 2014 23:19

Inverze obecné matice pomocí LU faktoru při řešení vlastního problému

#2 08. 03. 2014 00:30

Re: Inverze obecné matice pomocí LU faktoru při řešení vlastního problému

#3 08. 03. 2014 13:35

Re: Inverze obecné matice pomocí LU faktoru při řešení vlastního problému

#4 15. 06. 2018 22:08 Příspěvek uživatele stucja01 byl skryt uživatelem stucja01. Důvod: chyně zadáno

Zápatí