Lanczos algoritm

För noll space-finding-algoritmen, se block Lanczos-algoritmen .

För interpolationsmetoden, se Lanczos resampling .

För approximationen av gammafunktionen, se Lanczos approximation .

Lanczos -algoritmen är en iterativ metod utarbetad av Cornelius Lanczos som är en anpassning av potensmetoder för att hitta de $m$ "mest användbara" (som tenderar mot extremt högsta/lägsta) egenvärden och egenvektorer för en $n\ gånger n$ Hermitisk matris , där $m$ ofta men inte nödvändigtvis är mycket mindre än $n$ . Även om den är beräkningseffektiv i princip, var metoden som den ursprungligen formulerades inte användbar på grund av dess numeriska instabilitet .

1970 visade Ojalvo och Newman hur man gör metoden numeriskt stabil och tillämpade den på lösningen av mycket stora tekniska strukturer utsatta för dynamisk belastning. Detta uppnåddes genom att använda en metod för att rena Lanczos-vektorerna (dvs. genom att upprepade gånger ortogonalisera varje nygenererad vektor med alla tidigare genererade) till vilken grad av noggrannhet som helst, som när den inte utfördes producerade en serie vektorer som var starkt kontaminerade av de som var associerade med de lägsta naturliga frekvenserna.

I sitt ursprungliga arbete föreslog dessa författare också hur man väljer en startvektor (dvs. använd en slumptalsgenerator för att välja varje element i startvektorn) och föreslog en empiriskt bestämd metod för att bestämma m {\displaystyle m} , $reducerade$ talet av vektorer (dvs. den bör väljas till att vara ungefär 1,5 gånger antalet önskade exakta egenvärden). Strax därefter följdes deras arbete av Paige, som också gav en felanalys. 1988 producerade Ojalvo en mer detaljerad historia av denna algoritm och ett effektivt egenvärdesfeltest.

Algoritmen

Mata in en hermitisk matris $A$ av storleken $n\times n$ , och valfritt ett antal iterationer $m$ (som standard, låt $m=n$ ).

• Strängt taget behöver algoritmen inte tillgång till den explicita matrisen, utan endast en funktion $v\mapsto Av$ som beräknar produkten av matrisen med en godtycklig vektor. Denna funktion kallas högst $m$ gånger.

Mata ut en $n\times m$ matris $V$ med ortonormala kolumner och en tridiagonal reell symmetrisk matris $T=V^{*}AV$ av storleken $m\times m$ . Om $m=n$ , då $V$ är enhetlig , och $A=VTV^{*}$ .

Varning Lanczos-iterationen är benägen till numerisk instabilitet. När de körs i icke-exakt aritmetik bör ytterligare åtgärder (som beskrivs i senare avsnitt) vidtas för att säkerställa resultatens giltighet.

1. Låt $v_{1}\in \mathbb {C} ^{n}$ vara en godtycklig vektor med euklidisk norm $1$ .
2. Förkortat inledande iterationssteg:
   1. Låt $w_{1}'=Av_{1}$ .
   2. Låt $\alpha _{1}=w_{1}'^{*}v_{1}$ .
   3. Låt $w_{1}=w_{1}'-\alpha _{1}v_{1}$ .
3. För $j=2,\dots ,m$ gör:
   1. Låt $\beta _{j}=\|w_{j-1}\|$ (även euklidisk norm ).
   2. Om $\beta _{j}\neq 0$ , låt $v_{j}=w_{j-1}/\beta _{j }$ ,

      annars välj som $v_{j}$ en godtycklig vektor med euklidisk norm $1$ som är ortogonal mot alla $v_{1},\dots ,v_{j-1}$ .

   3. Låt $w_{j}'=Av_{j}$ .
   4. Låt $\alpha _{j}=w_{j}'^{*}v_{j}$ .
   5. Låt $w_{j}=w_{j}'-\alpha _{j}v_{j}-\beta _{j} v_{j-1}$ .
4. Låt $V$ vara matrisen med kolumner $v_{1},\dots ,v_{m}$ . Låt $T={\begin{pmatrix}\alpha _{1}&\beta _{2}&&&&0\\\beta _{2}&\alpha _{2}&\beta _{ 3}&&&\\&\beta _{3}&\alpha _{3}&\ddots &&\\&&\ddots &\ddots &\beta _{m-1}&\\&&&&\beta _{m- 1}&\alpha _{m-1}&\beta _{m}\\0&&&&\beta _{m}&\alpha _{m}\\\end{pmatrix}}$ .

Notera ${\displaystyle Av_{j}=w_{j}'=\beta _{j+1 }v_{j+1}+\alpha _{j}v_{j}+\beta _{j}v_{j-1}} för$ 1 $\displaystyle 1<j<m}$ .

Det finns i princip fyra sätt att skriva iterationsproceduren. Paige och andra verk visar att ovanstående ordningsföljd är den mest numeriskt stabila. I praktiken kan initialvektorn $v_{1}$ tas som ett annat argument för proceduren, med $\beta _{j}=0$ och indikatorer på numerisk unrecision inkluderas som ytterligare loopavslutningsvillkor.

Utan att räkna matris-vektormultiplikationen, gör varje iteration $O(n)$ aritmetiska operationer. Matris-vektormultiplikationen kan göras i $O(dn)$ aritmetiska operationer där $d$ är det genomsnittliga antalet element som inte är noll i en rad. Den totala komplexiteten är alltså $O(dmn)$ , eller $O(dn^{2})$ om $m=n$ ; Lanczos-algoritmen kan vara mycket snabb för glesa matriser. Schema för att förbättra numerisk stabilitet bedöms vanligtvis mot denna höga prestanda.

Vektorerna $v_{j}$ kallas Lanczos-vektorer . Vektorn $w_{j}'$ används inte efter att $w_{j}$ har beräknats, och vektorn $w_{j}$ används inte efter $v_{j+1}$ beräknas. Därför kan man använda samma lagring för alla tre. På samma sätt, om bara den tridiagonala matrisen $T$ söks, då behöver den råa iterationen inte $v_{j-1}$ efter att ha beräknat $w_{j}$ , även om vissa scheman för att förbättra det numeriska stabilitet skulle behöva det senare. Ibland beräknas de efterföljande Lanczos-vektorerna om från $v_{1}$ vid behov.

Tillämpning på egenproblemet

Lanczos-algoritmen tas oftast upp i samband med att hitta egenvärdena och egenvektorerna för en matris, men medan en vanlig diagonalisering av en matris skulle göra egenvektorer och egenvärden uppenbara från inspektion, är detsamma inte sant för tridiagonaliseringen som utförs av Lanczos. algoritm; icke-triviala ytterligare steg behövs för att beräkna även ett enda egenvärde eller egenvektor. Icke desto mindre är tillämpningen av Lanczos-algoritmen ofta ett betydande steg framåt för att beräkna egennedbrytningen. Om $\lambda$ är ett egenvärde av $A$ , och om $Tx=\lambda x$ ( $x$ är en egenvektor till $T$ ) så är $y =Vx$ är motsvarande egenvektor för $A$ (eftersom $Ay=AVx=VTV^{*}Vx=VTIx=VTx=V(\lambda x)=\lambda Vx= \lambda y$ ). Sålunda transformerar Lanczos-algoritmen egennedbrytningsproblemet för $A$ till egennedbrytningsproblemet för $T$ .

1. För tridiagonala matriser finns det ett antal specialiserade algoritmer, ofta med bättre beräkningskomplexitet än generella algoritmer. Till exempel, om $T$ är en $m\times m$ tridiagonal symmetrisk matris då:
   • Den kontinuerliga rekursionen gör det möjligt att beräkna det karakteristiska polynomet i ${\displaystyle O(m^{2})}-$ operationer och utvärdera det vid en punkt i ${\displaystyle O(m)}-$ operationer.
   • Dela -och-erövra-egenvärdesalgoritmen kan användas för att beräkna hela egenuppdelningen av $T$ i ${\displaystyle O(m^{2})}-$ operationer.
   • Den snabba multipolmetoden kan beräkna alla egenvärden i bara $O(m\log m)$ operationer.
2. Vissa allmänna egennedbrytningsalgoritmer, särskilt QR-algoritmen , är kända för att konvergera snabbare för tridiagonala matriser än för allmänna matriser. Asymptotisk komplexitet för tridiagonal QR är $O(m^{2})$ precis som för dela-och-härska-algoritmen (även om konstantfaktorn kan vara annorlunda); eftersom egenvektorerna tillsammans har $m^{2}$ element är detta asymptotiskt optimalt.
3. Även algoritmer vars konvergenshastigheter är opåverkade av enhetliga transformationer, såsom effektmetoden och invers iteration , kan ha prestandafördelar på låg nivå genom att appliceras på den tridiagonala matrisen $T$ snarare än den ursprungliga matrisen $A$ . Eftersom $T$ är mycket gles med alla element som inte är noll i mycket förutsägbara positioner, tillåter den kompakt lagring med utmärkt prestanda gentemot caching . Likaså $T$ är en reell matris med alla egenvektorer och egenvärden reella, medan $A$ i allmänhet kan ha komplexa element och egenvektorer, så reell aritmetik är tillräcklig för att hitta egenvektorerna och egenvärdena för $T$ .
4. Om $n$ är mycket stor, då reducering av $m$ så att $T$ har en hanterbar storlek kommer fortfarande att göra det möjligt att hitta de mer extrema egenvärdena och egenvektorerna för $A$ ; i $m\ll n$ -regionen kan Lanczos-algoritmen ses som ett förlustkompressionsschema för hermitiska matriser, som betonar att bevara de extrema egenvärdena.

Kombinationen av bra prestanda för glesa matriser och förmågan att beräkna flera (utan att beräkna alla) egenvärden är de främsta anledningarna till att man väljer att använda Lanczos-algoritmen.

Applikation för tridiagonalisering

Även om egenproblemet ofta är motivet för att tillämpa Lanczos-algoritmen, är operationen som algoritmen primärt utför tridiagonalisering av en matris, för vilken numeriskt stabila hushållartransformationer har gynnats sedan 1950-talet. Under 1960-talet ignorerades Lanczos-algoritmen. Intresset för det föryngrades av Kaniel-Paiges konvergensteorin och utvecklingen av metoder för att förhindra numerisk instabilitet, men Lanczos-algoritmen förblir den alternativa algoritmen som man bara försöker om Householder inte är tillfredsställande.

Aspekter där de två algoritmerna skiljer sig är:

• Lanczos drar fördel av att $A$ är en gles matris, medan Householder inte gör det, och kommer att generera fill-in .
• Lanczos fungerar genomgående med den ursprungliga matrisen $A$ (och har inga problem med att den bara är känd implicit), medan raw Householder vill modifiera matrisen under beräkningen (även om det kan undvikas).
• Varje iteration av Lanczos-algoritmen producerar ytterligare en kolumn i den slutliga transformationsmatrisen $V$ , medan en iteration av Householder producerar en annan faktor i en enhetlig faktorisering $Q_{1}Q_{2 }\dots Q_{n}$ av $V$ . Varje faktor bestäms dock av en enda vektor, så lagringskraven är desamma för båda algoritmerna, och $V=Q_{1}Q_{2}\dots Q_{n}$ kan beräknas i $O(n^{3})$ tid.
• Hushållaren är numeriskt stabil, medan rå Lanczos inte är det.
• , med endast $O(n)$ synkroniseringspunkter (beräkningarna av $\alpha _{j}$ och $\beta _{j}$ ). Hushållaren är mindre parallell, med en sekvens av $O(n^{2})$ beräknade skalära kvantiteter som var och en beror på den föregående kvantiteten i sekvensen.

Härledning av algoritmen

Det finns flera resonemang som leder till Lanczos-algoritmen.

En mer försynt kraftmetod

Huvudartikel: Power iteration

Potensmetoden för att hitta egenvärdet av största storleken och en motsvarande egenvektor för en matris $A$ är ungefär

1. Välj en slumpmässig vektor $u_{1}\neq 0$ .
2. För $j\geqslant 1$ (tills riktningen för $u_{j}$ har konvergerat) gör:
   1. Låt $u_{j+1}'=Au_{j}.$
   2. Låt $u_{j+1}=u_{j+1}'/\|u_{j+1}'\|.$

• I den stora $j$ -gränsen närmar sig $u_{j}$ den normerade egenvektorn som motsvarar det största egenvärdet.

En kritik som kan riktas mot denna metod är att den är slösaktig: den lägger ner mycket arbete (matris-vektor-produkterna i steg 2.1) på att extrahera information från matrisen A {\displaystyle A} , men uppmärksammar $det$ mycket sista resultat; implementeringar använder vanligtvis samma variabel för alla vektorer ${\displaystyle u_{j}} ,$ där varje ny iteration skriver över resultaten från den föregående. Tänk om vi istället behöll alla mellanresultat och organiserade deras data?

En bit information som trivialt är tillgänglig från vektorerna $u_{j}$ är en kedja av Krylov-underrymden . Ett sätt att påstå det utan att införa set i algoritmen är att hävda att den beräknar

en delmängd $\{v_{j}\}_{j=1}^{m}$ av en bas av $\mathbb {C} ^{n}$ så att $Ax\in \operatorname {span} (v_{1},\dotsc ,v_{j+1})$ för varje $x\in \operatörsnamn {span} (v_{1},\dotsc ,v_{j})$ och alla $1\leqslant j<m;$

detta är trivialt tillfredsställt av ${ \displaystyle u_{1},\dotsc ,u_{j-1}}$$\displaystyle v_{j}=u_{j}}$ så länge som $u_{j}$ är linjärt oberoende av (och om det finns ett sådant beroende kan man fortsätta sekvensen genom att välja som $v_{j}$ en godtycklig vektor linjärt oberoende av ${\displaystyle u_{1},\dotsc ,u_{j-1}} )$ . En bas som innehåller ${\displaystyle u_{j}}-$ vektorerna är dock sannolikt numeriskt dåligt konditionerad , eftersom denna sekvens av vektorer är avsedd att konvergera till en egenvektor av $A$ . För att undvika det kan man kombinera kraftiterationen med en Gram–Schmidt-process , för att istället producera en ortonormal bas av dessa Krylov-delrum.

1. Välj en slumpmässig vektor $u_{1}$ av euklidisk norm $1$ . Låt $v_{1}=u_{1}$ .
2. För $j=1,\dotsc ,m-1$ gör:
   1. Låt $u_{j+1}'=Au_{j}$ .
   2. För alla $k=1,\dotsc ,j$ låt $g_{k,j}=v_{k}^ {*}u_{j+1}'$ . (Dessa är koordinaterna för $Au_{j}=u_{j+1}'$ med avseende på basvektorerna $v_{1},\dotsc ,v_{j}$ .)
   3. Låt $w_{j+1}=u_{j+1}'-\summa _{k=1}^ {j}g_{k,j}v_{k}$ . (Avbryt komponenten av $u_{j+1}'$ som är i ${\displaystyle \operatorname {span} (v_{1},\dotsc ,v_{j})} .$ )
   4. Om $w_{j+1}\neq 0$ så låt $u_{j+1}=u_ {j+1}'/\|u_{j+1}'\|$ och ${\displaystyle v_{j+1}=w_{j+1}/\|w_{j+1}\|} , annars välj som u j + 1 = v$ j

      + $} =v_{j+1}}$ en godtycklig vektor av euklidisk norm $1$ som är ortogonal mot alla $v_{1},\dotsc ,v_{j}$ .

Relationen mellan power iterationsvektorerna $u_{j}$ och de ortogonala vektorerna $v_{j}$ är att

$Au_{j}=\|u_{j+1}'\|u_{j+1}=u_{j+1}'=w_{j+1}+\ summa _{k=1}^{j}g_{k,j}v_{k}=\|w_{j+1}\|v_{j+1}+\summa _{k=1}^{j }g_{k,j}v_{k}$ .

Här kan det observeras att vi faktiskt inte behöver $u_{j}$ vektorerna för att beräkna dessa $v_{j}$ , eftersom ${\displaystyle u_{j}-v_{j}\in \operatorname {span} (v_{1},\dotsc ,v_{j-1})} och därför skillnaden$ mellan $u_{j+1}'=Au_{j}$ och $w_{j+1}'=Av_{j}$ är i ${\displaystyle \operatorname {span} (v_{1},\dotsc ,v_{j})} ,$ som avbryts av ortogonaliseringsprocessen. Således beräknas samma grund för kedjan av Krylov-delrum av

1. Välj en slumpmässig vektor $v_{1}$ av euklidisk norm $1$ .
2. För $j=1,\dotsc ,m-1$ gör:
   1. Låt $w_{j+1}'=Av_{j}$ .
   2. För alla $k=1,\dotsc ,j$ låt $h_{k,j}=v_{k}^ {*}w_{j+1}'$ .
   3. Låt $w_{j+1}=w_{j+1}'-\summa _{k=1}^ {j}h_{k,j}v_{k}$ .
   4. Låt $h_{j+1,j}=\|w_{j+1}\|$ .
   5. Om $h_{j+1,j}\neq 0$ så låt $v_{j+1} =w_{j+1}/h_{j+1,j}$ ,

      annars välj som $v_{j+1}$ en godtycklig vektor av euklidisk norm $1$ som är ortogonal till alla $v_{1},\dotsc ,v_{j}$ .

uppfyller koefficienterna $h_{k,j}$

$Av_{j}=\sum _{k=1}^{j+1}h_{k,j}v_{k}$ för alla $j<m$ ;

definitionen $h_{j+1,j}=\|w_{j+1}\|$ kan verka lite udda, men passar det allmänna mönstret $h_{k,j}=v_{k}^{*}w_{j+1}'$ sedan

$v_{j+1}^{*}w_{j+1}'=v_{j+1}^{*}w_{j+1}=\|w_{j+1}\|v_{ j+1}^{*}v_{j+1}=\|w_{j+1}\|.$

Eftersom power iterationsvektorerna $u_{j}$ som eliminerades från denna rekursion uppfyller $u_{j}\in \operatorname {span } (v_{1},\ldots ,v_{j}),$ vektorerna $\{v_{j}\}_{j=1}^{m}$ och koefficienter $h_{k,j}$ innehåller tillräckligt med information från $A$ för att alla $u_{1},\ldots ,u_{m}$ kan beräknas , så ingenting gick förlorat genom att byta vektorer. (Det visar sig faktiskt att data som samlas in här ger betydligt bättre approximationer av det största egenvärdet än vad man får från lika många iterationer i effektmetoden, även om det inte nödvändigtvis är uppenbart vid denna tidpunkt.)

Denna sista procedur är Arnoldi-iterationen . Lanczos-algoritmen uppstår sedan som förenklingen man får genom att eliminera beräkningssteg som visar sig vara triviala när $A$ är hermitisk — särskilt de flesta $h_{k,j}$ koefficienterna visa sig vara noll.

Elementärt, om $A$ är hermitisk då

$h_{k,j}=v_{k}^{*}w_{j+1}'=v_{k}^{*}Av_{j}=v_{k}^{*}A^{ *}v_{j}=(Av_{k})^{*}v_{j}.$

För $k<j-1$ vet vi att $Av_{k}\in \operatorname {span} ( v_{1},\ldots ,v_{j-1})$ och eftersom $v_{j}$ genom att konstruktionen är ortogonal mot detta delrum måste denna inre produkt vara noll. (Detta är i grunden också anledningen till att sekvenser av ortogonala polynom alltid kan ges en tre-term recurrencerelation .) För $k=j-1$ får man

$h_{j-1 ,j}=(Av_{j-1})^{*}v_{j}={\overline {v_{j}^{*}Av_{j-1}}}={\overline {h_{j, j-1}}}=h_{j,j-1}$

eftersom den senare är verklig på grund av att den är normen för en vektor. För $k=j$ får man

$h_{j,j}=(Av_{j})^{*}v_{ j}={\overline {v_{j}^{*}Av_{j}}}={\overline {h_{j,j}}},$

vilket betyder att detta också är sant.

Mer abstrakt, om $V$ är matrisen med kolumner $v_{1},\ldots ,v_{m}$ då talen $h_{ k,j}$ kan identifieras som element i matrisen $H=V^{*}AV$ , och $h_{k,j}=0$ för $k>j+1;$ matrisen $H$ är övre Hessenberg . Eftersom

$}=\left(V^{*}AV\right)^{*}=V ^{*}A^{*}V=V^{*}AV=H}$

matrisen $H$ är hermitisk. Detta innebär att $H$ också är lägre Hessenberg, så det måste i själva verket vara tridiagional. Eftersom den är hermitisk är dess huvuddiagonal verklig, och eftersom dess första subdiagonal är verklig genom konstruktion, gäller samma sak för dess första superdiagonal. Därför $H$ en verklig, symmetrisk matris - matrisen $T$ i Lanczos algoritmspecifikation.

Samtidig approximation av extrema egenvärden

Ett sätt att karakterisera egenvektorerna för en hermitisk matris $A$ är som stationära punkter i Rayleigh-kvotienten

$r(x)={\frac {x^{*}Ax}{x^{*}x}},\qquad x\in \mathbb {C} ^{n}.$

I synnerhet är det största egenvärdet $\lambda _{\max }$ det globala maximum av $r$ och det minsta egenvärdet $\lambda _{\min }$ är det globala minimum av $r$ .

Inom ett lågdimensionellt delrum ${\mathcal {L}}$ av $\mathbb {C} ^{n}$ kan det vara möjligt att lokalisera maximalt $x$ och minimum $y$ av $r$ . Att upprepa det för en ökande kedja ${\mathcal {L}}_{1}\subset {\mathcal {L}}_{2}\subset \cdots$ producerar två sekvenser av vektorer : $x_{1},x_{2},\ldots$ och $y_{1},y_{2},\dotsc$ så att $x_{j},y_{j}\in {\mathcal {L}}_{j}$ och

${\begin{aligned}r(x_{1})&\ leqslant r(x_{2})\leqslant \cdots \leqslant \lambda _{\max }\\r(y_{1})&\geqslant r(y_{2})\geqslant \cdots \geqslant \lambda _{ \min }\end{aligned}}$

Frågan uppstår då hur man väljer delrymden så att dessa sekvenser konvergerar med optimal hastighet.

Från ${\displaystyle x_{j}} är$ den optimala riktningen för att söka större värden på $r$ den för gradienten $\nabla r(x_{j} )$ , och likaså från ${\displaystyle y_{j}} är$ den optimala riktningen för att söka mindre värden på $r$ den negativa gradienten $-\nabla r(y_{j})$ . I allmänhet

$\nabla r(x)={\frac {2}{x^{*}x}}(Ax-r (x)x),$

så intresseriktningarna är lätta nog att beräkna i matrisaritmetik, men om man vill förbättra både $x_{j}$ och $y_{j}$ så finns det två nya riktningar för att ta hänsyn till: $Ax_{j}$ och $Ay_{j};$ eftersom $x_{j}$ och $y_{j}$ kan vara linjärt oberoende vektorer (är faktiskt nära ortogonala), man kan i allmänhet inte förvänta sig att $Ax_{j}$ och $Ay_{j}$ är parallella. Är det därför nödvändigt att öka dimensionen på ${\mathcal {L}}_{j}$ med $2$ för varje steg? Inte om $\{{\mathcal {L}}_{j}\}_{j=1}^{m}$ tas för att vara Krylov-underrymden, eftersom $Az\in {\mathcal {L}}_{j+1}$ för alla $z\in { \mathcal {L}}_{j},$ alltså särskilt för både $z=x_{j}$ och $z=y_{j}$ .

Med andra ord kan vi börja med någon godtycklig initialvektor $x_{1}=y_{1},$ konstruera vektorrymden

${\mathcal {L}}_{j}=\operatörsnamn {span} (x_{1},Ax_{ 1},\ldots ,A^{j-1}x_{1})$

och sök sedan $x_{j},y_{j}\i {\mathcal {L}}_{j}$ så att

$r(x_{j})=\max _{z\in {\mathcal {L}}_{j}}r(z)\qquad {\text{and}}\qquad r(y_{j })=\min _{z\in {\mathcal {L}}_{j}}r(z).$

Eftersom $j$ :te potensmetoden iterate $u_{j}$ tillhör ${\mathcal {L}}_{j},$ följer det att en iteration för att producera x $x_{j}$ och ${\displaystyle y_{j}} kan inte konvergera långsammare än det$ potensmetoden, och kommer att uppnå mer genom att approximera båda egenvärdesextrema. För underproblemet att optimera $r$ på vissa ${\mathcal {L}}_{j}$ , det är bekvämt att ha en ortonormal bas $\{v_{1},\ldots ,v_{ j}\}$ för detta vektorutrymme. Sålunda leds vi återigen till problemet med att iterativt beräkna en sådan grund för sekvensen av Krylov-underrum.

Konvergens och annan dynamik

När man analyserar dynamiken i algoritmen är det bekvämt att ta egenvärdena och egenvektorerna för $A$ som givna, även om de inte är explicit kända för användaren. För att fixa notation, låt $\lambda _{1}\geqslant \lambda _{2}\geqslant \dotsb \geqslant \lambda _{n}$ vara egenvärdena (dessa är kända för alla att vara verkliga och därmed möjliga att beställa) och låt $z_{1},\dotsc ,z_{n}$ vara en ortonormal uppsättning egenvektorer så att $Az_{k}=\lambda _{k}z_{k}$ för alla $k=1,\dotsc ,n$ .

Det är också bekvämt att fixa en notation för koefficienterna för den initiala Lanczos-vektorn $v_{1}$ med avseende på denna egenbas; låt $d_{k}=z_{k}^{*}v_{1}$ för alla $k=1,\dotsc ,n$ , så att $\textstyle v_{1}=\sum _{k=1}^{n}d_{k}z_{k}$ . En startvektor $v_{1}$ utarmad på någon egenkomponent kommer att fördröja konvergensen till motsvarande egenvärde, och även om detta bara kommer ut som en konstant faktor i felgränserna förblir utarmningen oönskad. En vanlig teknik för att undvika att konsekvent drabbas av det är att välja $v_{1}$ genom att först rita elementen slumpmässigt enligt samma normalfördelning med medelvärde $0$ och skala sedan om vektorn till norm $1$ . Före omskalningen gör detta att koefficienterna $d_{k}$ också är oberoende normalfördelade stokastiska variabler från samma normalfördelning (eftersom ändringen av koordinater är enhetlig), och efter omskalning av vektorn $(d_{1},\dotsc ,d_{n})$ kommer att ha en enhetlig fördelning på enhetssfären i $\mathbb {C} ^{n}$ . Detta gör det möjligt att binda sannolikheten för att till exempel $|d_{1}|<\varepsilon$ .

Det faktum att Lanczos-algoritmen är koordinatagnostisk – operationer tittar bara på inre produkter av vektorer, aldrig på enskilda element av vektorer – gör det enkelt att konstruera exempel med känd egenstruktur att köra algoritmen på: gör A {\displaystyle $}$ en diagonal matris med de önskade egenvärdena på diagonalen; så länge som startvektorn $v_{1}$ har tillräckligt många element som inte är noll, kommer algoritmen att mata ut en allmän tridiagonal symmetrisk matris som $T$ .

Kaniel–Paiges konvergensteori

Efter $m$ iterationssteg av Lanczos-algoritmen är $T$ en $m\times m$ reell symmetrisk matris, som på samma sätt som ovan har $m$ egenvärden $\theta _{1}\geqslant \theta _{2}\geqslant \dots \geqslant \theta _{m}.$ Med konvergens förstås i första hand konvergensen av $\theta _{1}$ till $\lambda _{1}$ (och den symmetriska konvergensen av $\theta _{m}$ till $\ lambda _{n}$ ) när $m$ växer, och sekundärt konvergensen av något område $\theta _{1},\ldots ,\theta _{k}$ av egenvärden för $T$ till sina motsvarigheter $\lambda _{1},\ldots ,\lambda _{k}$ av $A$ . Konvergensen för Lanczos-algoritmen är ofta storleksordningar snabbare än den för kraftiterationsalgoritmen.

Gränserna för $\theta _{1}$ kommer från ovanstående tolkning av egenvärden som extremvärden för Rayleigh-kvoten $r(x)$ . Eftersom $\lambda _{1}$ är a priori maximum av $r$ på hela $\mathbb {C} ^{n},$ medan ${ \displaystyle \theta _{1}}$ är bara maxvärdet på en $m$ -dimensionellt Krylov-underrum får vi trivialt $\lambda _{1}\geqslant \theta _{1}$ . Omvänt ger vilken punkt ${\displaystyle x} som helst$ i det Krylov-underrummet en nedre gräns $r(x)$ för $\theta _{1}$ , så om en punkt kan vara visas för vilken $\lambda _{1}-r(x)$ är liten så ger detta en snäv gräns på $\theta _{1}$ .

Dimensionen $m$ Krylovs delrum är

$\operatorname {span} \left\{v_{1},Av_{1},A^{ 2}v_{1},\ldots ,A^{m-1}v_{1}\right\},$

så vilket element som helst av det kan uttryckas som $p(A)v_{1}$ för något polynom $p$ av högst grad $m-1$ ; koefficienterna för det polynomet är helt enkelt koefficienterna i den linjära kombinationen $v_{1},Av_{1},A^{2}v_{1},\ldots ,A^{m-1}v_{1}$ . Polynomet vi vill ha kommer att visa sig ha reella koefficienter, men för tillfället bör vi även tillåta komplexa koefficienter, och vi kommer att skriva $p^{*}$ för polynomet som erhålls genom att komplext konjugera alla koefficienter för $p$ . I denna parametrisering av Krylov-underrummet har vi

$r(p(A)v_{1})={\frac {(p(A) )v_{1})^{*}Ap(A)v_{1}}{(p(A)v_{1})^{*}p(A)v_{1}}}={\frac {v_ {1}^{*}p(A)^{*}Ap(A)v_{1}}{v_{1}^{*}p(A)^{*}p(A)v_{1}} }={\frac {v_{1}^{*}p^{*}(A^{*})Ap(A)v_{1}}{v_{1}^{*}p^{*}( A^{*})p(A)v_{1}}}={\frac {v_{1}^{*}p^{*}(A)Ap(A)v_{1}}{v_{1 }^{*}p^{*}(A)p(A)v_{1}}}$

Genom att nu använda uttrycket för $v_{1}$ som en linjär kombination av egenvektorer får vi

$Av_{1}=A\summa _{k=1}^{n}d_{k }z_{k}=\summa _{k=1}^{n}d_{k}\lambda _{k}z_{k}$

och mer generellt

$q(A)v_{1}=\summa _{k=1}^{n}d_{k}q (\lambda _{k})z_{k}$

för vilket polynom som helst $q$ .

Således

$\lambda _{1}-r(p(A)v_{1})=\lambda _{1}-{\frac {v_{1}^{*}\sum _{k=1}^ {n}d_{k}p^{*}(\lambda _{k})\lambda _{k}p(\lambda _{k})z_{k}}{v_{1}^{*}\ summa _{k=1}^{n}d_{k}p^{*}(\lambda _{k})p(\lambda _{k})z_{k}}}=\lambda _{1} -{\frac {\sum _{k=1}^{n}|d_{k}|^{2}\lambda _{k}p(\lambda _{k})^{*}p(\lambda _{k})}{\sum _{k=1}^{n}|d_{k}|^{2}p(\lambda _{k})^{*}p(\lambda _{k} )}}={\frac {\summa _{k=1}^{n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})\left|p( \lambda _{k})\right|^{2}}{\sum _{k=1}^{n}|d_{k}|^{2}\left|p(\lambda _{k}) \right|^{2}}}.$

En viktig skillnad mellan täljare och nämnare här är att $k=1$ termen försvinner i täljaren, men inte i nämnaren. Så om man kan välja $p$ är stor vid $\lambda _{1}$ men liten vid alla andra egenvärden, får man en snäv gräns för felet ${\ displaystyle \lambda _{1}-\theta _{1}}$ .

Eftersom $A$ har många fler egenvärden än $p$ har koefficienter, kan detta tyckas vara en hög ordning, men ett sätt att möta det är att använda Chebyshev-polynom . Skriva $c_{k}$ för graden $k$ Chebyshev-polynom av det första slaget (det som uppfyller $c_{k}(\cos x)=\cos(kx)$ för alla $x$ ), har vi ett polynom som stannar i intervallet $[-1 ,1]$ på det kända intervallet $[-1,1]$ men växer snabbt utanför det. Med viss skalning av argumentet kan vi låta det mappa alla egenvärden utom $\lambda _{1}$ till $[-1,1]$ . Låta

$p(x)=c_{m-1}\left({\frac {2x-\lambda _ {2}-\lambda _{n}}{\lambda _{2}-\lambda _{n}}}\right)$

(om $\lambda _{2}=\lambda _{1}$ , använd istället det största egenvärdet strikt mindre än $\lambda _{1}$ ), sedan maximalt värde på $|p(\lambda _{k})|^{2}$ för $k\geqslant 2$ är $1$ och det minimala värdet är $0$ , alltså

$\lambda _{1}-\theta _{1}\leqslant \lambda _{1}-r(p(A)v_{1})={\frac {\sum _{k=2}^ {n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})|p(\lambda _{k})|^{2}}{\summa _{k =1}^{n}|d_{k}|^{2}|p(\lambda _{k})|^{2}}}\leqslant {\frac {\sum _{k=2}^{ n}|d_{k}|^{2}(\lambda _{1}-\lambda _{k})}{|d_{1}|^{2}|p(\lambda _{1})| ^{2}}}\leqslant {\frac {(\lambda _{1}-\lambda _{n})\summa _{k=2}^{n}|d_{k}|^{2}} {|p(\lambda _{1})|^{2}|d_{1}|^{2}}}.$

dessutom

$p(\lambda _{1})=c_{m-1}\left({\frac {2\lambda _{1}-\lambda _{2}-\lambda _{n}}{\ lambda _{2}-\lambda _{n}}}\right)=c_{m-1}\left(2{\frac {\lambda _{1}-\lambda _{2}}{\lambda _ {2}-\lambda _{n}}}+1\höger);$

kvantiteten

$\rho ={\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}-\lambda _{n} }}$

(dvs förhållandet mellan det första egengapet och diametern för resten av spektrumet ) är således av nyckelvikt för konvergenshastigheten här. Skriver också

$R=e^{\operatörsnamn {arcosh} (1+2\rho )}=1+2\rho +2{\sqrt {\rho ^{2}+\rho }},$

vi kan dra slutsatsen att

${\begin{aligned}\lambda _{1}-\theta _{1}&\leqslant {\frac {(\lambda _{1}-\lambda _{n})\left(1-| d_{1}|^{2}\höger)}{c_{m-1}(2\rho +1)^{2}|d_{1}|^{2}}}\\[6pt]&= {\frac {1-|d_{1}|^{2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {1 }{\cosh ^{2}((m-1)\operatörsnamn {arcosh} (1+2\rho ))}}\\[6pt]&={\frac {1-|d_{1}|^{ 2}}{|d_{1}|^{2}}}(\lambda _{1}-\lambda _{n}){\frac {4}{\left(R^{m-1}+R ^{-(m-1)}\höger)^{2}}}\\[6pt]&\leqslant 4{\frac {1-|d_{1}|^{2}}{|d_{1} |^{2}}}(\lambda _{1}-\lambda _{n})R^{-2(m-1)}\end{aligned}}$

Konvergenshastigheten styrs alltså huvudsakligen av $R$ , eftersom denna gräns krymper med en faktor $R^{-2}$ för varje extra iteration.

Som jämförelse kan man överväga hur konvergenshastigheten för potensmetoden beror på ${\displaystyle \rho } ,$ men eftersom potensmetoden i första hand är känslig för kvoten mellan absoluta värden på egenvärdena behöver vi $|\lambda _{n}|\leqslant |\lambda _{2}|$ för egengapet mellan $\lambda _{1}$ och $\lambda _{2}$ att vara den dominerande. Under den begränsningen är det fall som gynnar potensmetoden mest att ${\displaystyle \lambda _{n}=-\lambda _{2}} ,$ så tänk på det. Senare i effektmetoden, iterationsvektorn:

$u=(1-t^{2})^{1/2}z_{1}+ tz_{2}\approx z_{1}+tz_{2},$

där varje ny iteration effektivt multiplicerar $z_{2}$ -amplituden $t$ med

${\frac {\lambda _{2}}{\lambda _{1}}}={\frac {\lambda _{2}}{\lambda _{2}+(\lambda _{1} -\lambda _{2})}}={\frac {1}{1+{\frac {\lambda _{1}-\lambda _{2}}{\lambda _{2}}}}}= {\frac {1}{1+2\rho }}.$

Uppskattningen av det största egenvärdet är då

$u^{*}Au=(1-t^{2})\lambda _{1}+t^{2} \lambda _{2},$

så ovanstående gräns för Lanczos-algoritmens konvergenshastighet bör jämföras med

$\lambda _{1}-u^{*}Au=(\lambda _{1}-\lambda _{2}) t^{2},$

som krymper med en faktor på $(1+2\rho )^{-2}$ för varje iteration. Skillnaden kokar alltså ner till den mellan $1+2\rho$ och $R=1+2\rho +2{\sqrt { \rho ^{2}+\rho }}$ . I $\rho \gg 1$ region, den senare är mer som ${\displaystyle 1+4\rho } ,$ och presterar som effektmetoden med ett egengap dubbelt så stort; en märkbar förbättring. Det mer utmanande fallet är dock det för $\rho \ll 1,$ där $R\approx 1+2{\sqrt {\rho }}$ är en ännu större förbättring av egengapet; ρ $displaystyle \rho \gg 1}$ regionen är där Lanczos-algoritmen konvergensmässigt gör den minsta förbättringen av effektmetoden.

Numerisk stabilitet

Stabilitet betyder hur mycket algoritmen kommer att påverkas (dvs. kommer den att ge det ungefärliga resultatet nära det ursprungliga) om det introduceras och ackumuleras små numeriska fel. Numerisk stabilitet är det centrala kriteriet för att bedöma användbarheten av att implementera en algoritm på en dator med avrundning.

För Lanczos-algoritmen kan det bevisas att med exakt aritmetik , uppsättningen av vektorer $v_{1},v_{2},\cdots ,v_{m+ 1}$ konstruerar en ortonormal bas, och de lösta egenvärdena/vektorerna är goda approximationer till den ursprungliga matrisen. Men i praktiken (eftersom beräkningarna utförs i flyttalsaritmetik där inexakthet är oundviklig), går ortogonaliteten snabbt förlorad och i vissa fall kan den nya vektorn till och med vara linjärt beroende av den mängd som redan är konstruerad. Som ett resultat kan vissa av egenvärdena för den resulterande tridiagonala matrisen inte vara approximationer till den ursprungliga matrisen. Därför är Lanczos-algoritmen inte särskilt stabil.

Användare av denna algoritm måste kunna hitta och ta bort dessa "falska" egenvärden. Praktiska implementeringar av Lanczos-algoritmen går i tre riktningar för att bekämpa denna stabilitetsfråga:

1. Förhindra förlust av ortogonalitet,
2. Återställ ortogonaliteten efter att basen har genererats.
3. Efter att alla goda och "falska" egenvärden har identifierats, ta bort de falska.

Variationer

Variationer på Lanczos-algoritmen finns där de involverade vektorerna är höga, smala matriser istället för vektorer och de normaliserande konstanterna är små kvadratiska matriser. Dessa kallas "block" Lanczos-algoritmer och kan vara mycket snabbare på datorer med ett stort antal register och långa minneshämtningstider.

Många implementeringar av Lanczos-algoritmen startar om efter ett visst antal iterationer. En av de mest inflytelserika omstartade varianterna är den implicit omstartade Lanczos-metoden, som är implementerad i ARPACK . Detta har lett till ett antal andra omstartade varianter som omstartad Lanczos-bidiagonalisering. En annan framgångsrik omstartad variant är Thick-Restart Lanczos-metoden, som har implementerats i ett mjukvarupaket som heter TRLan.

Nullutrymme över ett ändligt fält

Huvudartikel: Block Lanczos algoritm

1995 publicerade Peter Montgomery en algoritm, baserad på Lanczos-algoritmen, för att hitta delar av nollutrymmet för en stor gles matris över GF(2) ; eftersom mängden personer som är intresserade av stora glesa matriser över ändliga fält och mängden personer som är intresserade av stora egenvärdeproblem knappast överlappar varandra, kallas detta ofta även för block Lanczos-algoritmen utan att orsaka orimlig förvirring.

Ansökningar

Lanczos algoritmer är mycket attraktiva eftersom multiplikationen med $A\,$ är den enda linjära operationen i stor skala. Eftersom vägda texthämtningsmotorer implementerar just denna operation, kan Lanczos-algoritmen tillämpas effektivt på textdokument (se latent semantisk indexering ) . Egenvektorer är också viktiga för storskaliga rankningsmetoder som HITS-algoritmen utvecklad av Jon Kleinberg , eller PageRank -algoritmen som används av Google.

Lanczos algoritmer används också inom den kondenserade materiens fysik som en metod för att lösa Hamiltonians av starkt korrelerade elektronsystem, samt i skalmodellkoder inom kärnfysik .

Genomföranden

NAG -biblioteket innehåller flera rutiner för lösning av storskaliga linjära system och egenproblem som använder Lanczos-algoritmen.

MATLAB och GNU Octave kommer med ARPACK inbyggt. Både lagrade och implicita matriser kan analyseras genom eigs() -funktionen ( Matlab / Octave ).

På liknande sätt, i Python , har SciPy - paketet scipy.sparse.linalg.eigsh som också är ett omslag för funktionerna SSEUPD och DSEUPD från ARPACK som använder den implicit omstartade Lanczos-metoden.

En Matlab-implementering av Lanczos-algoritmen (observera precisionsproblem) är tillgänglig som en del av Gaussian Belief Propagation Matlab-paketet . GraphLab - samarbetsfiltreringsbiblioteket innehåller en storskalig parallell implementering av Lanczos-algoritmen (i C++) för flerkärniga.

PRIMME - biblioteket implementerar också en Lanczos-liknande algoritm.

Vidare läsning

•   Golub, Gene H .; Van Loan, Charles F. (1996). "Lanczos metoder" . Matrisberäkningar . Baltimore: Johns Hopkins University Press. s. 470–507. ISBN 0-8018-5414-8 .
• Ng, Andrew Y .; Zheng, Alice X.; Jordan, Michael I. (2001). "Länkanalys, egenvektorer och stabilitet" (PDF) . IJCAI'01 Proceedings of the 17th International Joint Conference on Artificiell Intelligens . 2 : 903-910.