Ominaisarvoista ja -vektoreista

[viesti Survo-keskustelupalstalla (2001-2013)]

Olen viime viikkoina käyttänyt verrattain paljon aikaa selvittääkseni,
mikä on nopea tapa laskea isoista symmetrisistä  matriiseista
suurimmat ominaisarvot ja niitä vastaavat ominaisvektorit.

Tämä tehtävä tulee vastaan mm. pääkomponenttianalyysissa (lähtökohtana
kovarianssi- tai korrelaatiomatriisi) ja klassisessa moniulotteisessa
skaalauksessa (lähtökohtana etäisyysmatriisin keskistetty muoto).

Tilastollisissa sovelluksissa riittää yleensä saada lasketuksi
kohtalaisella tarkkuudella po. matriisin muutamat suurimmat ominaisarvot
ja niitä vastaavat ominaisvektorit. Esim. moniulotteisessa skaalauksessa
tyydytään usein kahteen ensimmäiseen. Pääkomponenttianalyysissa ja
pääakselifaktoroinnissa kiinnostus rajoittuu tyypillisesti alle
kymmeneen.

Viime vuosikymmeninä numeerisessa matematiikassa vallalla ovat olleet
menetelmät, jotka tuottavat samalla kertaa kaikki ominaisarvot ja
-vektorit.
Survossa olen 1970-luvulla ottanut käyttöön algoritmin, joka perustuu
symmetrisen matriisin muuntamiseen tridiagonaaliseksi (Householderin
kierroin) ja ominaisarvojen (sekä -vektorien) laskentaan tästä
supistetusta muodosta QL-menetelmällä. Tietojeni lähteenä on ollut
alan klassikko, Wilkinsonin ja Reinschin "Handbook for Automatic
Computation", Volume 2 vuodelta 1971.
Ymmärtääkseni mitään ratkaisevasti parempaa ei ole tullut esiin
viime vuosikymmeninä. Rinnakkaislaskennassa (mihin PC:t eivät pysty)
tosin esim. vanhasta Jacobin menetelmästä on tullut uudelleen
varteenotettava vaihtoehto (kts. esim. Golub, van Loan: Matrix
Computations, 2nd edition, 1989, ss. 444-).

Tilastollisten monimuttujamenetelmien edellyttämissä ominaisarvo-
tehtävissä, kuten alussa totesin, riittää laskea vain muutamat
suurimmat ominaisarvot ja vastaavat ominaisvektorit.
Tällöin em. menetelmät eivät suurilla dimensioilla voi olla nopeimpia,
koska ne tavallaan laskevat turhan paljon.

Yksinkertaisin keino tähän rajoitettuun tehtävään on vanha potenssi-
menetelmä, mutta se on melko hidas ja epätarkka varsinkin silloin, kun
suurimmat ominaisarvot ovat lähellä toisiaan. Paljon tehokkaampi ja
luotettavampi on Lanczosin menetelmä.

Kornel Lanczos (1893-1974) oli unkarilainen fyysikko ja matemaatikko,
joka nuorena tutki suhteellisuusteoriaa ja toimi yhteistyössä Albert
Einsteinin kanssa. Myöhemmin hän kunnostautui myös numeerisessa
matematiikassa ja keksi mm. nopean Fourier-muunnoksen (25 vuotta ennen
John Tukeyta, jonka nimiin tämä saavutus nykyisin luetaan).

Hänen ominaisarvotehtävää koskeva tärkeä artikkelinsa "An Iteration
Method for the Solution of the Eigenvalue Problem of Linear Differential
and Integral Operators" ilmestyi vuonna 1950.
Lanczosin menetelmän keskeiset ominaisuudet on kuvattu em. Golubin
ja van Loanin kirjassa (luku 9: Lanczos Methods).

Lanczosin idea on luoda alkuperäisestä symmetrisestä n*n-matriisista A
ja jostain yksikkövektorin mittaisesta lähtövektorista v
olennaisesti muotoa A*v olevalla k peräkkäisellä kertolaskulla
k*k-tridiagonaalinen matriisi T, jolla on likimain samat
ominaisarvot kuin mitä ovat A:n suurimmat ominaisarvot.
On hämmästyttävää, että, vaikka n olisi tuhansia, tyypillisesti jo
arvolla k=10 saadaan hyviä likiarvoja suurimmille ominaisarvoille
laskemalla ne em. Wilkinsonin ja Reinschin menetelmällä alkuperäiseen
matriisiin verrattuna erittäin pienestä k*k-matriisista T.

Vastaavat ominaisvektorit saadaan em. prosessin sivutuotteena, mutta
valitettavasti ne menettävät käytännössä ortogonaalisuuttaan, mikä
ei johdu ainoastaan pyöristysvirheistä. Lanczosin algoritmi antaa
myös helposti uudelleen samoja ominaisarvoja ja jopa perättömiä arvoja.
Näiden erottaminen "oikeista" ja monikertaisten ominaisarvojen
tunnistaminen ovat ongelmia, joihin on esitetty erilaisia ratkaisuja.

Koska en ole löytänyt mitään ehdottoman varmaa reseptiä po. ongelmien
selvittämiseen, olen päätynyt menettelyyn, joka ei ole ehkä nopeuden
suhteen aivan paras mahdollinen, mutta luotettava ja riittävän tarkka.

Kuvaan menettelyni pääpiirteet:

Tehtävänä on siis ratkaista symmetrisen n*n-matriisin A
ominaisarvotehtävä suurimpien ominaisarvojen l_1>=l_2>=,...,>=l_m
ja ominaisvektorien v_1,v_2,...,v_m osalta, jotka toteuttavat
yhtälöt
              A*v_i = l_i*v_i, i=1,2,...,m.

Sovelletaan Lanczosin iteraatiota matriisiin A aluksi 10 kertaa ja
lasketaan saadusta tridiagonaalisesta 10*10-matriisista T ominaisarvot
ja -vektorit sekä näistä A:n suurimman ominaisarvon l_1 ja ominais-
vektorin v_1 likiarvot. Tutkitaan ratkaisun tarkkuutta laskemalla
erotusmatriisi
              A*v_1-l_1*v_1

ja tämän matriisin alkioitten neliösumma s2. Jos sqrt(s2/n) alittaa
tarkkuusvakion eps (oletusarvoksi olen valinnut 1e-13), ratkaisu
L_1,v_1 kelpuutetaan. Jos ei, jatketaan Lanzos-iteraatiota 20
askeleella ja tehdään uusi tarkistus. Tätä prosessia jatketaan edelleen
lisäämällä iteraatioiden lukumäärää askelmäärin 40,60,80,...
kunnes haluttu tarkkuus on saavutettu.
Seuraava pari l_2,v_2 määrätään tämän jälkeen redusoimalla (samaan
tapaan kuin potenssimenetelmässä) matriisista ensimmäisen parin l_1,v_1
osuus laskemalla matriisi
              A2 = A - l_1*v_1*v_1'
ja toistetaan edellä esitetty menettely matriisin A2 suhteen.
Vastaavalla tavalla "kuoritaan" sitten A2:sta parin l_2,v_2 osuus
laskemalla
              A3 = A2 - l_2*v_2*v_2'
ja tätä toistetaan, kunnes on saatu lasketuksi l_m,V_m.

Tässä menettelyssä siis kussakin vaiheessa tarkkaillaan vain yhtä
ominaisarvo-vektori-paria ja jätetään käyttämättä kaikki muu informaatio
em. ortogonaalisuus- ja valeominaisarvo-ongelmien vuoksi.

Matriisikomento
MAT #EIGLAN(A,m,S,L)
laskee symmetrisen matriisin A m suurinta ominaisarvoa (L) ja niitä
vastaavat ominaisvektorit (S). Se on yleensä nopeampi kuin
MAT SPECTRAL DECOMPOSITION OF A TO S,L
(joka laskee kaikki ominaisarvot ja -vektorit) vasta kun matriisin
A dimensio n on 300 tai suurempi. Huomattakoon, että nykyisillä
koneilla dimensiolla n=300 po. laskenta sujuu alle puolessa sekunnissa
eli menetelmällä ei ole paljon merkitystä.
Kun n on tuhansia, mikä on täysin mahdollinen tilanne moniulotteisessa
skaalauksessa, MAT #EIGLAN on huomattavasti nopeampi. Tätä kuvaa
seuraava esimerkki.

Tehdään symmetrinen 3000*3000-matriisi A, jonka suurimmat ominaisarvot
ovat 10, 10, 9, 8.999, 8.998, 7, 6 ja loput ovat "satunnaisia arvoja"
välillä (0,1). Tämä on hiukan hankala tilanne, sillä kaksi suurinta
on samoja ja seuraavat kolme ovat varsin lähellä toisiaan.

Matriisi A luodaan muodossa Q*D*Q', missä Q on "satunnainen" ortogonaa-
linen matriisi ja D ominaisarvojen muodostama lävistäjämatriisi.
A syntyy seuraavalla laskentakaaviolla:

MATRIX K ///
10 10 9 8.999 8.998 7 6

n=3000
MAT A=ZER(n,n)
MAT #TRANSFORM A BY rand(2007)  / satunnainen neliömatriisi A
MAT A=MTM(A)                    / A'*A on symmetrinen
MAT QR OF A TO Q,R       / QR-hajotelmasta saadaan ortogonaalinen Q.
MAT D=ZER(n,1)           / Ominaisarvojen lävistäjämatriisin laskenta
MAT #TRANSFORM D BY rand(2005)
MAT SAVE K
MAT D(1,1)=K'
MAT D=DV(D)
MAT A=Q*D*Q'             / Lopullinen A-matriisi

MAT #EIGLAN laskee 2 ensimmäistä ominaisarvoa ja vektoria 5 sekunnissa:

TIME COUNT START
MAT #EIGLAN(A,2,S,L)
TIME COUNT END   4.987

ja 7 ensimmäistä noin 20 sekunnissa

TIME COUNT START
MAT #EIGLAN(A,7,S,L)
TIME COUNT END   20.249

Lasketut ominaisarvot:
MAT LOAD L,1234.123456789012345,CUR+1
MATRIX L
Eigenvalues
///                  eigenval
eigen1     10.000000000000050
eigen2      9.999999999999988
eigen3      8.999999999999986
eigen4      8.998999999999898
eigen5      8.997999999999966
eigen6      7.000000000000017
eigen7      6.000000000000077

Vastaava laskenta täyden spektraalihajotelman kautta ottaa lähes 800
sekuntia:
TIME COUNT START
MAT SPECT DECOMP of A TO S2,L2
TIME COUNT END   797.206

ja alkupään ominaisarvot ovat:
MAT L2=L2(1:7)
MAT LOAD L2,1234.123456789012345,CUR+1
MATRIX L2
L2(1:7,1)
///                  eigenval
ev1        10.000000000000057
ev2         9.999999999999925
ev3         8.999999999999995
ev4         8.998999999999871
ev5         8.997999999999944
ev6         7.000000000000027
ev7         6.000000000000060

Olen (osittain opetuksellisista syistä) laatinut myös potenssimenetelmää
käyttävän MAT #EIGFEW -komennon, joka selvästi häviää sekä ajoajan että
tulosten tarkkuuden osalta MAT #EIGLANille:

TIME COUNT START
MAT #EIGFEW(A,7,S0,L0,1e-16,300)
TIME COUNT END   69.900

MAT LOAD L0,1234.123456789012345,CUR+1
MATRIX L0
Eigenvalues
///                  eigenval
eigen1      9.999999999999943
eigen2     10.000000000000032
eigen3      8.999816976439480
eigen4      8.998966658162855
eigen5      8.998216345468908
eigen6      7.000000000000023
eigen7      6.000000000000064

Tulosten tarkkuutta kuvaavat seuraavat laskelmat.
Lasketaan kunkin menetelmän osalta, millä tarkkuudella ominais-
arvoyhtälöt toteutuvat erotusmatriisien alkioiden neliösummina
jaettuna alkioiden lukumäärällä:

MAT #EIGLAN:
MAT E=A*S-S*DV(L)
MAT F=SUM(SUM(E,2)')
MAT_F(1,1)/7/2000=5.1581745805117e-027

Spektraalihajotelma antaa tarkemman tuloksen:
MAT E=A*S2(*,1:7)-S2(*,1:7)*DV(L2(1:7))
MAT F=SUM(SUM(E,2)')
MAT_F(1,1)/7/2000=4.7009490616147e-031

Potenssimenetelmän tulos on selvästi edellisiä heikompi:
MAT E=A*S0(*,1:7)-S0(*,1:7)*DV(L0(1:7))
MAT F=SUM(SUM(E,2)')
MAT_F(1,1)/7/2000=0.00000000005124

Ominaisvektorien ortogonaalisuutta voi tarkastella samalla tavalla:

MAT #EIGLAN:
MAT E=MTM(S)-IDN(7,7)  / Verrataan matriisia S'*S yksikkömatriisiin.
MAT F=SUM(SUM(E,2)')
MAT_F(1,1)/7/7=4.0144282174305e-028

Spektraalihajotelma:
MAT E=MTM(S2(*,1:7))-IDN(7,7)
MAT F=SUM(SUM(E,2)')
MAT_F(1,1)/7/7=1.1562424801766e-030

Potenssimenetelmä:
MAT E=MTM(S0)-IDN(7,7)
MAT F=SUM(SUM(E,2)')
MAT_F(1,1)/7/7=0.00000000018077

MAT #EIGLAN antaa tilastollisia sovelluksia ajatellen kohtuullisen
tarkat tulokset vaikka häviääkin tässä tapauksessa spektraali-
hajotelmalle.
Nopeudeltaan se on kuitenkin tässä 40-kertainen tai jopa 160-kertainen
jos tyydytään vain kahteen ensimmäiseen ominaisarvoon.

MAT #EIGLAN tulee mukaan SURVO MM:n versiossa 2.32.

- Seppo Mustonen

Luettelo keskustelunaiheista

Survo-keskustelupalstan (2001-2013) viestit arkistoitiin aika ajoin sukrolla, joka automaattisesti rakensi viesteistä (yli 1600 kpl) HTML-muotoisen sivukokonaisuuden. Vuoden 2013 alusta Survo-keskustelua on jatkettu entistäkin aktiivisemmin osoitteessa forum.survo.fi. Tervetuloa mukaan!