Eigenwaarde-algoritmes voor Hamiltoniaanse
Matrices
Harald Monsuur
13 juli 2018
Bachelorscriptie Wiskunde Begeleiding: dr. Jan Brandts
Korteweg-de Vries Instituut voor Wiskunde
Samenvatting
In deze tekst zullen we de stabiliteit van verschillende algoritmes van de SR-decompositie onderzoeken evenals de condities van de bijbehorende factoren. Vervolgens passen we deze algoritmes toe om eigenwaardes van Hamiltoniaanse matrices te berekenen. Dit laatste al-goritme maakt gebruik van algemene stelling over het zogeheten GR-alal-goritme. Voordat we dit alles doen geven we eerst een introductie over bilineaire vormen. We defini¨eren hier ook de Lie-algebra van een bilineaire vorm en de automorfismegroep. Vervolgens zoomen we in op het symplectische product. We proberen een beter begrip te krijgen van de symplectische matrices die afkomstig zijn van de automorfismegroep horend bij het symplectische product. Deze voorkennis zal aanleiding geven tot de te onderzoeken SR-decompositie. Omdat deze decompositie niet uniek is zullen we kijken naar de mogelijk te behalen numerieke winst.
Titel: Eigenwaarde-algoritmes voor Hamiltoniaanse Matrices Auteur: Harald Monsuur,
Begeleiding: dr. Jan Brandts,
Tweede beoordelaar: dr. Chris Stolk, Einddatum: 13 juli 2018
Korteweg-de Vries Instituut voor Wiskunde Universiteit van Amsterdam
Science Park 904, 1098 XH Amsterdam http://www.kdvi.uva.nl
Inhoudsopgave
1 Introductie 4
2 Bilineaire vormen 6
2.1 Algemene introductie bilineaire vormen . . . 6
2.2 Automorfismegroep . . . 8
2.3 Lie-algebra en Jordan-algebra . . . 9
2.4 Equivalentie van bilineaire vormen . . . 11
3 Symplectische matrices 14 3.1 G-reflectoren . . . 14
3.1.1 The mapping problem . . . 15
3.1.2 Opspansel van G-reflectoren . . . 16
3.2 Ortho-symplectische matrices . . . 18 3.2.1 Givens-rotaties . . . 19 3.2.2 Symplectische Householder-transformaties . . . 19 3.2.3 Gauss-transformaties . . . 20 4 De SR-decompositie 22 4.1 Symplectische Gram-Schmidt . . . 22
4.2 Onderzoek naar existentie en uniciteit . . . 24
4.3 Symplectische triangularisatie met G-reflectoren . . . 28
4.3.1 Aanpak . . . 28
4.3.2 Conditiegetallen G-reflectoren . . . 30
4.4 SR-decompositie met ortho-symplectische matrices . . . 31
4.5 Optimale conditie S en R . . . 33
5 Eigenwaardenproblemen 36 5.1 Het SR-algoritme . . . 36
5.1.1 Schets bewijs . . . 37
5.2 Keuze van polynoom en deflatie . . . 38
5.3 Hessenberg-matrices . . . 40 5.3.1 Impliciete SR-algoritme . . . 41 6 Testresultaten 43 7 Conclusie 47 8 Appendix 48 9 Populaire samenvatting 49 10 Bibliografie 51 11 Matlabcodes 52
1 Introductie
Dit is een scriptie die gaat over numerieke wiskunde. Op het eerste gezicht lijkt numerieke wiskunde een beetje saai. De wiskunde heeft in dit geval wel een concreet doel. Namelijk het oplossingen aanbieden van problemen die mensen in de praktijk ook echt willen oplossen. Maar daar lijkt de pret wel op te houden. Er is een probleem waarvan je de uitkomst niet goed meer met de hand kan uitrekenen. Dus gebruik je de computer er maar voor. Het zal dan nog even een kwestie van goed programmeren zijn, maar dan ben je er ook wel. Zou je denken. De ervaring leert dat er nog best wat wiskunde naar voren komt.
Dat er interessante wiskunde naar voren komt ligt aan het volgende. Een computer kan maar met eindige nauwkeurigheid rekenen. Dit heeft een aantal vervelende gevolgen. Er zijn namelijk problemen die zeer slecht reageren op rekenfouten. Dit heten slecht-geconditioneerde problemen. Voor een numeriek wiskundige is dit slecht nieuws, er kan namelijk niets aan verandert worden.
Er zijn echter ook algoritmes die niet goed reageren op rekenfouten. In dit geval ligt de oorzaak niet in het wiskundig probleem, maar aan het algoritme zelf. Een klassiek voorbeeld is het Gram-Schmidt algoritme. Dit algoritme transformeert een normale basis naar een orthonormale basis. Wiskundig klopt het algoritme, het algoritme gaat er alleen vanuit dat de reeds geproduceerde kolommen orthonormaal zijn. Door de rekenfouten zal dit in de praktijk niet zo zijn. Als je daar geen rekening mee houdt worden deze rekenfouten steeds verder uitvergroot. Dit is in dit geval niet te wijten aan de eindige nauwkeurigheid van de computer zelf, maar eerder aan de manier waarop het algoritme reageert op deze rekenfouten. In dit geval kunnen we het algoritme gelukkig aanpassen om deze uitvergroting van fouten te verminderen. U merkt al dat er in de numerieke wiskunde veel aandacht geschonken moet worden aan wat het algoritme doet. Niet alleen op zuiver wiskundig gebied, maar ook in de minder sprookjesachtige wereld vol met rekenfouten moet het algoritme acceptabel zijn.
In deze scriptie komt er ook een Gram-Schmidt proces voor in een soortgelijke context. Ook in dit geval is het algoritme erg na¨ıef, wat we ook terugzien in de resultaten. Zie hiervoor het laatste hoofdstuk. Dit geeft aanleiding tot interessante wiskunde die het probleem op een andere manier oplost.
Eigenlijk is deze introductie nog maar beperkt, er is enorm veel te beleven als numerieke wiskundige. Bekijk de volgende situatie. Stel dat er vrije keuze in een stel parameters is, hoe kies je die dan? Een normale wiskundige zal hier niet heel erg moeilijk over doen en gewoon maar wat kiezen. Een numerieke wiskundige ziet hier potentie. De vrije parameters kiest hij om de algoritmes stabieler te laten werken. Het doel zal dan ook zijn om een methode van kiezen te zoeken en te bewijzen dat die keuze ook echt goed is om zekerheid te verschaffen.
Deze scriptie heeft als belangrijkste doel om de SR-decompositie stabiel uit te voeren en de conditionering van de factoren acceptabel te krijgen. Als toepassing hiervan zullen we ook een eigenwaarde algoritme bekijken.
Om de context van de SR-decompositie te begrijpen geven we eerst een introductie over bilinieaire vormen. Dit zal vrij algebra¨ısch van aard zijn. In dat hoofdstuk zullen we zien dat bilineaire vormen aanleiding geven tot definities van bepaalde klassen van matrices. Over deze
matrices kunnen we met behulp van de bilineaire vorm een aantal bekende resultaten bewijzen. We zullen in dat hoofdstuk ook het symplectische product defini¨eren wat een bilineaire vorm is. Met dit product zullen we de sciptie voortzetten. Dit product induceert zoals gezegd een aantal klassen van matrices, waaronder de groep van symplectische matrices en de Hamiltoniaanse matrices. In de hoofdstuk 3 bekijken we symplectische matrices uit [5] en [2] die we vervolgens in hoofdstuk 4 gebruiken voor algoritmes voor het berekenen van de SR-decompositie van een matrix. We gebruiken hier algoritmes uit [1] en [7]. We geven een alternatief bewijs dat de algoritmen niet vastlopen. We onderzoeken in hoofdstuk 4 ook het bestaan en de uniciteit van zo’n decompositie van een matrix. We bewijzen het bekende resultaat dat de decompositie niet uniek is. Dit levert, zoals daarnet al benoemd, een extra vraag op over stabiliteit van de algoritmes. Op deze vraag gaan we verder in. In het laatste hoofdstuk geven we een toepassing van de SR-decompositie, namelijk een eigenwaarde-algoritme uit [1] Dankzij de structuurbehoudende werking van die decompositie op Hamiltoniaanse matrices, zal dit een relatief snel algoritme zijn. Tevens is besteden we kort aandacht waarom dit eigenwaarde-algoritme werkt. Hiervoor gebruiken we het bewijs uit [8]. Als laatste geven we een aantal testresultaten van de verschillende algoritmen.
Dankwoord
Mijn dank gaat uit naar mijn begeleider dr. Jan Brandts. Hij heeft mij tijdens zijn hoorcolleges weten te enthousiasmeren een onderwerp te kiezen in de numerieke wiskunde. Zijn begeleiding was ontspannen en leerzaam.
2 Bilineaire vormen
In dit hoofdstuk zullen we allereerst de definitie geven van een bilineaire vorm. Uit de lineaire algebra kennen we al een bilineaire vorm, namelijk het inproduct. Met behulp van dit inpro-duct konden we een aantal andere zaken defini¨eren. We kennen de orthogonale matrices en de symmetrische matrices die symmetrisch blijven onder conjugatie met orthogonale matrices. Ook konden we een vectorruimte laten corresponderen met zijn duale via het inproduct. Een bilineaire vorm geeft ook aanleiding tot dit soort matrices en symmetrie¨en. We zullen dit in hoofdstuk ook speciale aandacht geven aan de zogeheten symplectische vorm met zijn struc-turen. Dit hoofdstuk is gebaseerd op [4]. We hebben enkele bewijzen herschreven zodat er geen lineaire afbeeldingen naar de duale worden gebruikt. Daarnaast geven we definities van de verschillende verzamelingen lineaire afbeeldingen in de paragrafen 2.2 en 2.3, deze komen uit [2]. We geven hierover een aantal bewijzen van algemeen bekende resultaten.
2.1 Algemene introductie bilineaire vormen
Definitie 2.1.1. Laat V een eindig-dimensionale re¨ele of complexe vectorruimte zijn. Een bilineaire vorm is een afbeelding B : V × V → F, die lineair is in beide argumenten:
B(av + w, q) = aB(v, q) + B(w, q) en B(q, av + w) = aB(q, v) + B(q, w) ∀v, w, q ∈ V, a ∈ F. In deze definitie staat F voor de re¨ele getallen of de complexe getallen. Dit zullen we voortaan zo blijven doen.
Voortaan zullen we een bilineaire vorm vaak gewoon een vorm noemen. Er komen namelijk alleen bilineaire vormen voor in dit artikel. Eerst een aantal voorbeelden van bilineaire vormen. Voorbeeld 2.1.2. Wanneer V = Rn en M ∈ Rn×n dan is de afbeelding
B : V × V → R (v, w) 7→ vTM w een bilineaire vorm.
Voorbeeld 2.1.3. Een ander voorbeeld is de determinant met V = R2: B : R2× R2→ R; (v, w) 7→ det(v|w).
Deze vorm is duidelijk niet symmetrisch. Hij is zelfs anti-symmetrisch zoals in de komende definitie. Merk meteen op dat we ook diezelfde B kunnen defini¨eren door
B(v, w) = vT 0 1 −1 0
Voorbeeld 2.1.4. Een voorbeeld uit [4], waarbij V = C[0, 1] en dus oneindigdimensionaal is. Kies hiervoor een continue (of Lebesgue integreerbare) k : [0, 1]2→ R en definieer:
Bk(f, g) =
Z
[0,1]2
f (x)g(y)k(x, y)dxdy. Dat B bilineair is volgt uit de lineariteit van de integraal.
Definitie 2.1.5. Een bilineaire vorm noemen we
• Symmetrisch wanneer B(v, w) = B(w, v) ∀v, w ∈ V ;
• Scheef-symmetrisch (ook wel anti-symmetrisch) als B(v, w) = −B(w, v) ∀v, w ∈ V ; Voor een bilineaire vorm gaan we nu een notie van orthogonaliteit defini¨eren, namelijk v ⊥B w als B(v, w) = 0. Dit zal echter in het algemeen niet een symmetrische relatie zijn,
het kan zo zijn dat B(w, v) 6= 0 terwijl B(v, w) = 0. Het is zelfs zo dat orthogonaliteit alleen maar een symmetrische relatie is als B symmetrisch is of alternerend. Het bewijs voor het volgende lemma kunt u vinden in [4].
Lemma 2.1.6. Zij V een vectorruimte met B een bilineaire vorm. Dan geldt dat orthogona-liteit een equivalentierelatie is dan en slechts dan als B symmetrisch of alternerend is.
Voor een deelruimte van een vectorruimte met vorm B kunnen we het orthogonale comple-ment defini¨eren. Het orthogonale complement van U ⊂ V is
U⊥B = {v ∈ V : B(v, u) = 0 ∀u ∈ U }.
Merk hierbij wel op dat U ∩ U⊥B niet alleen maar nul hoeft te bevatten. Voor een
anti-symmetrische vorm S geldt bijvoorbeeld dat v ∈ span{v}⊥S.
Een andere eigenschap die een bilineaire vorm kan hebben is niet-gedegenereerdheid. Definitie 2.1.7. Een bilineaire vorm heet gedegenereerd wanneer er voor elke niet-triviale vector een links-en rechtsproduct ongelijk aan nul is:
∀v ∈ V − {0} ∃w1, w2 zodat B(v, w1) 6= 0 6= B(w2, v).
Zoals we zagen in Voorbeeld 2.1.2 en Voorbeeld 2.1.3, kunnen we veel vormen beschrijven met behulp van een matrix. Dit blijkt vaker waar te zijn. Voor eindigdimensionale ruimten kunnen we elke bilineaire vorm defini¨eren met behulp van een matrix op eenduidige wijze. Dit lemma komt uit [4].
Lemma 2.1.8. Zij B een bilineaire vorm op een eindigdimensionale re¨ele of complexe vec-torruimte V zijn met basis β = {v1, . . . vn}. Noem coβ de afbeelding zijn die een vector naar
zijn co¨ordinaten stuurt. Dan bestaat er een unieke matrix M ∈ Rn×n zodat voor alle v, w ∈ V geldt dat
B(v, w) = coβ(v)TM coβ(w).
Deze matrix M noemen we de matrix horend bij B ten opzichte van de basis β. Bewijs. Kies M = [B(vi, vj)]ni,j=1.
Opmerking 2.1.9. Wanneer B een vorm is op V = Rn dan zullen we de matrix M , zodat voor alle v, w ∈ V geldt dat
B(v, w) = vTM w,
de matrix horend bij B noemen als er geen basis wordt genoemd. De matrix M is in dit geval gewoon de matrix horend bij B ten opzichte van de standaardbasis.
Het niet-gedegenereerd zijn van een bilineaire vorm kunnen we ook beschrijven met behulp van de bijbehorende matrix t.o.v. een basis.
Lemma 2.1.10. Zij B een bilineaire vorm op een eindigdimensionale re¨ele of complexe vec-torruimte V zijn met een basis β = {v1, . . . vn}. Laat Mβ de bijbehordende matrix zijn t.o.v.
de basis β. Dan is B niet gedegenereerd dan en slecht dan als Mβ inverteerbaar is
Bewijs. Alternatief. We bewijzen de stelling alleen maar voor ´e´en van de links-en rechtspro-ducten. Het bewijs kan gemakkelijk afgemaakt worden.
=⇒ Stel dat B niet gedegenereerd is, maar Mβ is inverteerbaar. Dan is er een v ∈ Rn met
M v = 0, dus dan geldt B(co−1β (w), co−1β (v)) = wTM v = 0 voor elke w ∈ V . Dit is een tegenspraak.
⇐= Stel dat M inverteerbaar is. Zij v ∈ V dan geldt dat (coβ(v)TMT)M coβ(v) 6= 0, want
M coβ(v) 6= 0. Dus B(v, co−1β (M coβ(v))) 6= 0.
Een andere eigenschap van bilineaire vormen is dat ze elementen uit de duale vectorruimte V∗ defini¨eren. De functie
B(·, w) : V → F; v 7→ B(v, w) is lineair en dus een elemenent van V∗.
Wanneer B niet gedegenereerd is hebben we het volgende lemma uit [4].
Lemma 2.1.11. Zij B een niet-gedegenereerde bilineaire vorm op een eindigdimensionale re¨ele of complexe vectorruimte V . Dan geldt dat de lineaire functionalen B(·, w) heel V∗ omvatten:
{B(·, v) : V → R; w 7→ B(w, v) met v ∈ V } = V∗.
Bewijs. Alternatief. Het is duidelijk dat de functie B(·, v) met v ∈ V in V∗ zit. Nu het omgekeerde. Laat β een basis zijn van V en M de matrix horend bij B zijn t.o.v. β. Definieer ook een inproduct op V door middel van hv, wi = coβ(v)Tcoβ(w). Laat f ∈ V∗. We weten
dan dat f (·) gelijk is aan coβ(·)Tcoβ(v) voor een zekere v ∈ V door middel van een stelling
van Riesz. Omdat B niet gedegenereerd is volgt dat M inverteerbaar is. Kies dus nu w ∈ Rn
zodat M w = coβ(v) en definieer u = co−1β (w). Dan vinden we dat f (·) = B(·, u).
2.2 Automorfismegroep
Bij een inproductruimte waren er matrices/lineaire afbeeldingen die het inproduct intact hiel-den: de orthogonale matrices. Voor bilineare vormen kunnen we ook zo’n groep van lineaire operatoren defini¨eren. In deze paragraaf geven we een aantal bekende resultaten.
Definitie 2.2.1. Laat V een re¨ele of complexe vectorruimte zijn met B een niet-gedegenereerde bilineaire vorm op V . Dan is
een verzameling van automorfismen horend bij B. Wanneer V eindig-dimensionaal is dan heet G de automorfismegroep horend bij B.
Wanneer V eindig-dimensionaal is en B niet gedegenereerd is dan is G ook echt een groep door de volgende observaties. Laat A, U ∈ G:
B(Av, Aw) = B(v, w) = B(U v, U w) ∀v, w ∈ V
Merk op, samen met de niet-gedegenereerdheid van B, dat A en U alleen nul in de kern hebben. Ze zijn dus inverteerbare lineaire operatoren. Nu hebben we
B(A−1U v, A−1U w) = B(AA−1U v, AA−1U w) = B(U v, U w) = B(v, w) ∀v, w ∈ V Dus A−1U is een element van G. Dus we concluderen dat G een groep is. Dit bewijs kunnen we niet overhevelen naar de oneindige dimensie. Voor een vectorruimte met oneindige dimensie is het niet altijd waar dat een element in G inverteerbaar is, wanneer zijn kern alleen uit nul bestaat.
Opmerking 2.2.2. Wanneer V = Fn kunnen we de automorfismegroep ook anders aandui-den. In dit geval geldt dat er een M is zodat voor alle v, w ∈ Fn geldt dat B(v, w) = vTM w. Dus B(Gv, Gw) = vTGTM Gw. Als G ∈ G, dan moet dus gelden dat GTM G = G. Dat laatste geeft ons een equivalente definitie voor elementen in de automorfismegroep als V = Fn.
Het volgende lemma gaat over de eigenwaarden van elementen uit de automorfisme groep. We bewijzen dat het spectrum van elementen in de autormorfismegroep een bepaalde structuur hebben.
Lemma 2.2.3. Laat V met dim(V ) < ∞ en B een niet-gedegenereerde vorm. Dan geldt voor elke A ∈ G, dat als λ ∈ σ(A), dan ook 1/λ ∈ σ(A).
Bewijs. Alternatief. We bewijzen het voor V = Fn. Het bewijs voor ene willekeurige V met
dim(V ) < ∞ gaat dat via de matrix horend bij de vorm.
Stel dat S ∈ G en Sv = λv, S is in dit geval een matrix. Dan geldt ook STM λv = STM Sv = M v.
Dus M v 6= 0 is een eigenvector van ST met eigenwaarde 1λ. Dus S heeft ook die eigenwaarde, want S en ST hebben hetzelfde spectrum.
Opmerking 2.2.4. Neem weer V = Rn. Omdat voor G ∈ G geldt dat det(M ) = det(GTM G) = det(GT) det(M ) det(G) zien we dat er moet gelden dat det(G) = ±1. Omdat det(G) = Q
λ∈σ(G)λ kunnen we het voorgaande lemma sterker maken door te zeggen dat er evenveel
eigenwaarden gelijk aan λ zijn als aan 1/λ.
2.3 Lie-algebra en Jordan-algebra
Met bilineaire vormen kunnen we ook andere klasses van matrices defini¨eren. We geven hier enkele bekende resultaten.
Definitie 2.3.1 (Geadjungeerde). Zij V een vectorruimte met B een bilineaire vorm daarop. Laat T een begrensde lineaire operator (een continue lineaire afbeelding) zijn, dan noemen we de lineaire operator T∗: V → V die voldoet aan
B(T v, w) = B(v, T∗w) ∀v, w ∈ V de geadjungeerde van T .
Opmerking 2.3.2. Merk op dat wanneer V eindigdimensionaal is en B niet gedegenereerd dat dan de geadjungeeerde van een G ∈ G gelijk is aan zijn inverse. Dit volgt door:
B(v, w) = B(Gv, Gw) = B(v, G∗Gw) ∀v, w ∈ V
en dus B(v, G∗Gw − w) = 0 ∀v ∈ V waaruit volgt dat G∗Gw = w voor elke w wanneer B niet-gedegenereerd is.
We defini¨eren nu de twee klassen lineaire operatoren wiens structuur bepaald worden door middel van de bijbehorende bilineaire vorm.
Definitie 2.3.3. Zij V een vectorruimte met B een bilineaire vorm op V . De Jordan-algebra is een verzameling lineaire operatoren:
J = {S : S een lineaire operator en S∗ = S}
Definitie 2.3.4. Zij V een vectorruimte met B een bilineaire vorm op V . De Lie-algebra is een verzameling lineaire operatoren:
L = {S : S een lineaire operator en S∗ = −S}
Opmerking 2.3.5. Wanneer V = Fn kunnen we de Lie-en Jordan-algebra ook anders de-fini¨eren. Laat hierbij M de matrix zijn zodat B(v, w) = vTM W voor alle v, w ∈ V (zie Opmerking 2.1.9).
J = {S ∈ Fn×n| STM = M S} en
L = {S ∈ Fn×n| STM = −M S}.
Een handige eigenschap is dat de Jordan en Lie-algebra gesloten zijn onder conjugatie met elementen van G:
Lemma 2.3.6. Zij V een vectorruimte met een niet-gedegenereerde vorm B. Dan geldt voor elke G ∈ G, J ∈ J en L ∈ L het volgende:
G−1J G ∈ J en G−1LG ∈ L.
Bewijs. Alternatief. We gaan de definitie van een element uit J na voor G−1J G, het bewijs voor de Lie-algebra gaat hetzelfde. We hebben dat
B(G−1J Gv, w) = B(J Gv, Gw) = B(Gv, J Gw) = B(v, G−1J Gw),
waarbij we twee keer Opmerking 2.3.2 gebruik voor G en G−1. Dus (G−1J G)∗ = G−1J G ∈ J.
Door dit lemma is het mogelijk om in algoritmes de structuur van de matrices te behouden. We kunnen elementen uit de Lie-algebra namelijk conjugeren met automorfismen van B, wat soms ook echt nuttig is. In de eigenwaardes van lineaire afbeeldingen uit de Lie-algebra kunnen we bijvoorbeeld een bepaalde structuur vinden. Dit staat samengevat in het volgende lemma: Lemma 2.3.7. Zij V een vectorruimte met een niet-gedegenereerde vorm B. Laat L ∈ L. Dan λ ∈ σ(L) ⇐⇒ −λ ∈ σ(L).
Bewijs. Alternatief. We kiezen weer eerst V = Fn. Laat M de matrix horend bij B zijn t.o.v.
de standaarbasis. Laat L ∈ L met Lv = λv. Dan hebben we door B(w, Lv) = B(−Lw, v) het volgende:
M λv = M Lv = −LTM v.
Dus we zien dat M v 6= 0 een eigenvector is voor L met eigenwaarde −λ. Doordat de spectra van M en MT gelijk zijn hebben we −λ ∈ σ(M ).
Merk op dat we dezelfde bewijsstrategie kunnen proberen op een element J uit de Jordan-algebra, dat levert ons op dat λ ∈ σ(J ). Dat is natuurlijk geen nieuwe informatie.
2.4 Equivalentie van bilineaire vormen
Twee bilineaire vormen B1 en B2 op eindig-dimensionale vectorruimten kan je als equivalent
beschouwen wanneer ze hetzelfde zijn, alleen dan ten opzichte van verschillende bases. Dit betekent concreet dat er twee bases β1 en β2 zijn voor V zodat de matrixrepresentaties van
de bilineaire vormen ten opzichte van de twee bases overeenkomen. Dus Mβ1,B1 = Mβ2,B2.
Het volgende lemma uit [4] geeft in matrixtaal aan wanneer twee vormen equivalent zijn. Het bewijs laten we achterwege.
Lemma 2.4.1. Laat V eindigdimensionaal en re¨eel of complex zijn en B1 en B2 twee
vor-men op V . Dan zijn de twee vorvor-men equivalent wanneer voor een basis β geldt dat er een inverteerbare matrix C is zodat
CTMβ,B1C = Mβ,B2.
We kunnen bewijzen dat elke niet-gedegenereerde ani-symmetrische vorm equivalent is. Tevens kan zo’n vorm alleen voorkomen in een even dimensionale ruimte. Om dit te bewijzen hebben we eerst het volgende gereedschap uit [4] nodig:
Lemma 2.4.2. Zij V een eindig-dimensionale vectorruimte met een anti-symmetrische of symmetrische niet-gedegenereerde vorm B daarop. Als U een lineaire deelruimte is van V zodat BU: U × U → R een niet-gedegenereerde vorm is, dan geldt
V = U ⊕ U⊥B.
Bewijs. Alternatief. Omdat BU een niet-gedegenereerde vorm is op U , geldt dat U∗ =
{B(u, ·) : U → F : u ∈ U}. Neem een v ∈ V en stel dat B(v, q) 6= 0 voor een zekere q ∈ U . Dan weten we dat er een u0 ∈ U is zodat B(v, u) = B(u0, u) voor alle u ∈ U . Maar
dan geldt dus dat v − u0 ∈ U⊥B. Dus v = u0+ (v − u0) en we concluderen V = U⊥B + U .
Ook geldt dat U ∩ U⊥B = {0}, want als u in de doorsnede zit dan is B
U gedegenereerd. Dus
de som is ook uniek, en dus hebben we een directe som. Dat BU⊥B niet gedegenereerd is
volgt doordat U ⊂ (U⊥B)⊥B en doordat B niet gedegenereerd is op V . Stel namelijk dat voor
q ∈ U⊥B geldt dat B(q, w) = 0 ∀w ∈ U⊥B, dan moet gelden dat B(q, u) 6= 0 voor een u ∈ U .
Maar dan q 6∈ U⊥B, wat natuurlijk een tegenspraak is.
Dit lemma komt uit [4].
Lemma 2.4.3. Zij V een eindig-dimensionale vectorruimte met een niet-gedegenereerde anti-symmetrische bilineaire vorm B. Dan is de dimensie van V deelbaar door twee. Tevens zijn alle inverteerbare anti-symmetrische vormen horend bij de vectorruimte V equivalent.
Bewijs. We zullen gebruikmaken van Lemma 2.1.11. In dit bewijs gaan we een basis {e1, . . . , en, f1. . . , fn} van V construeren zodat de bijbehorende matrix M gelijk is aan:
J = 0 In −In 0 .
Oftewel B(ei, fj) = δij en B(ei, ej) = B(fi, fj) = 0 voor elke i en j. Omdat het equivalent
zijn van vormen een equivalentierelatie is, zal het bewijs dan rond zijn. We gaan inductie doen. Stel dus dat voor elke dimensie lager dan n er een basis is zodat de matrix horend bij een anti-symmetrische vorm gegeven wordt door J . Nu kies e1 willekeurig en kies f1 zodat
B(e1, f1) = 1, dit kan want B is niet-gedegenereerd. Nu is B een niet gedegenereerde vorm op
U1 = Span{e1, f1}, dus V = U1L U1⊥B, waarbij BU⊥B
1
weer niet gedegenereerd is. Op U⊥B
1
is B zelfs anti-symmetrisch, dus we kunnen een basis {e2, . . . , en, f2. . . , fn} vinden zodat de
matrix horend bij BU⊥B 1
van het type J is. De matrix horend bij {e1, . . . , en, f1. . . , fn} is
vervolgens ook van het type J . Dus het bewijs is rond.
Nu alle niet-gedegenereerde anti-symmetrische vormen equivalent zijn hoeven we eigenlijk alleen maar ´een van die vormen te bestuderen. Toch zullen we twee van deze vormen bestu-deren. De Lie-algebra van de eerste vorm blijkt makkelijker te bestuderen, maar bepaalde matrix-decomposities zijn weer overzichtelijker in de tweede vorm. Beide vormen defini¨eren we hieronder:
Definitie 2.4.4. Zij V = Rn dan heet de vorm
B(v, w) = vTJ w, waarbij J = 0 In −In 0
Definitie 2.4.5. Zij V = Rn dan heet de vorm B(v, w) = vTJ w,ˆ waar ˆ J = 0 1 0 · · · 0 −1 0 0 · · · 0 0 0 . .. .. . ... 0 1 0 0 −1 0 het quasi-symplectische product.
Zoals net gezegd, de Lie-algebra van het symplectische product is ’goed’ te bestuderen. Het volgende lemma komt uit [5]
Lemma 2.4.6. Wanneer B het symplectische product is op R2n×2n dan geldt voor A ∈ R2n het volgende. A = A11 A12 A21 A22
∈ L, met Aij ∈ Rn×n ⇐⇒ A12 en A21 symmetrisch en A22= −AT11.
Bewijs. Alternatief. We moeten hebben dat ATJ = −J A. Dit uitwerken geeft precies het resultaat wat we willen.
Elementen uit de automorfismegroep horend bij het (quasi-)symplectische product noemen we voortaan (quasi-)symplectische matrices. Evenzo, elementen van de Lie-algebra horend bij het (quasi-)symplectische product noemen we voortaan (quasi-)Hamiltoniaanse matrices. In hoofstuk 5 zullen we eigenwaarde-algoritmes bestuderen voor deze Hamiltoniaanse matrices. Om resultaten over te hevelen van het quasi-symplectische product naar het symplectische product maken we gebruik van de volgende permutatiematrix
P = [e1|e3|e5| . . . |e2n−1|e2|e4| . . . |e2n]. Voor n = 3 hebben we dus
P = 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 1 .
In het volgende lemma staat precies hoe de twee vormen in verband staan met elkaar. Het bewijs laten we achterwege
Lemma 2.4.7. Zij V = R2n en laat J en ˆJ zijn zoals in de definities hierboven. Er gelden
de volgende uitspraken: • PTJ P = J .ˆ
• Als S symplectisch is dan is P SPT quasi-symplectisch. Omgekeerd, als S
quasi-symplectischis dan is PTSP symplectisch.
• Als H Hamiltoniaans is dan is P HPT Hamiltoniaans. Omgekeerd, als H
3 Symplectische matrices
In dit hoofdstuk gaan we de automorfismegroep van het symplectische product onderzoeken. Het doel van dit hoofdstuk is om verschillende soorten symplectische matrices te vinden en hun eigenschappen. In dit hoofdstuk vinden we twee verschillende verzamelingen die de symplectische groep voortbrengen, beide met hun eigen voordelen. We zullen in het begin rang 1 verstoorde identiteitmatrices onderzoeken, dit zijn matrices waar we ook makkelijk mee kunnen werken. Daarna zullen we symplectische matrices zoeken die ook orthogonaal zijn. Orthogonaliteit van matrices is een belangrijke eigenschap, want die matrices hebben een conditiegetal t.o.v. de 2-norm gelijk aan 1. Na dit hoofdstuk zullen we genoeg gereedschap hebben om de zogeheten SR-decompositie in hoofdstuk 4 stabiel uit te voeren. Met deze algoritme zullen we een eigenwaarde-algoritme maken die we in hoofstuk 5 gaan bekijken.
3.1 G-reflectoren
In deze paragraaf zullen we eerst rang 1 verstoorde identiteitmatrices zoeken die in de auto-morfisme groep G horend bij een bilineaire vorm op de Fn zitten. Dit zijn dus matrices van
de vorm
I + β · uvT
die in G zitten voor het bijbehorende product B. Dit soort matrices noemen we G-reflectoren. Dit soort matrices zijn handig. De inverse is makkelijk op te schrijven en dus bruikbaar in algoritmes. Dan is het al helemaal handig als deze matrices de hele automorfismegroep voortbrengen. Dit is bijna het geval bij orthogonale matrices. Deze worden daarom gebruikt in een algoritme voor de QR-decompositie. Dit is een algoritme die met behulp van verme-nigvuldigingen met orthogonale matrices het Gram-Schmidt proces stabieler uitvoert. Later zullen we G-reflectoren horend bij het symplectische product gebruiken voor een soortgelijk algoritme met symplectische matrices.
In [2] is een criterium ontwikkeld voor wanneer matrices een G-reflector zijn horend bij een vorm B.
Stelling 3.1.1. Laat Fn uitgerust zijn met een niet-gedegeneerde bilineaire vorm B en laat G zijn automorfisme groep zijn.
• Als G een G-reflector is dan kan G geschreven worden als G = I + βuuTM waarbij M de matrix is horend bij de bilineaire vorm B.
• Een matrix G van deze vorm zit in G en is dus een G-reflector wanneer aan de volgende voorwaarde is voldaan:
Merk hier gelijk op dat als B anti-symmetrisch is dat de voorwaarde in het tweede gedeeldte van de stelling triviaal waar is. Dus alle matrices van de vorm G = I+βuuTJ zijn G-reflectoren horend bij het symplectische product. Merk op dat we dit ook direct kunnen controleren.
3.1.1 The mapping problem
We weten nu wanneer een rang 1 verstoorde matrix een G-reflector. We willen nu nog on-derzoeken wat we precies met deze matrices kunnen doen. Om dit soort matrices te kunnen gebruiken in algoritmen is het handig om te weten wanneer het mogelijk is om een G-reflector te vinden die een vector x naar een andere vector y stuurt. Dit wordt in het Engels ’the mapping problem’ genoemd. In [2] wordt hiervoor een oplossing gegeven in het geval we te maken hebben met een symmetrische of een anti-symmetrische vorm.
Hier volgen we grotendeels het bewijs uit [2]. Laat B een niet-gedegenereerde vorm op Fn zijn met M de bijbehorende matrix. Merk hiervoor op dat als Gx = (I + βuuTM )x = y we het volgende kunnen zeggen:
Gx = (I + βuuTM )x = x + βuB(u, x) = y. Dus we vinden dat
y − x = βB(u, x)u. Deze vergelijking heeft als oplossing
u = α(y − x) en
β = 1
αB(α(y − x), x).
Dus we zien hier al dat B(y − x, x) mag niet nul zijn. Alles bij elkaar zien we dat G = 1
B(y − x, x)(y − x)(y − x)
TM.
We weten dus nu dat een G-reflector die x op y afbeeldt van deze vorm moet zijn. De vraag is alleen of deze G ook daadwerkelijk een G-reflector is. In het anti-symmetrische geval was dit altijd waar. Voor het symmetrische geval moeten we de voorwaarde in Stelling 3.1 nagaan:
(2 + 1
B(y − x, x)B(y − x, y − x))M (y − x) = 0. Omdat y − x 6= 0 en M inverteerbaar is, moeten we nagaan of
2 + 1
B(y − x, x)B(y − x, y − x) = 0. Merk hiervoor eerst op dat moet gelden dat
B(y, y) = B(Gx, Gx) = B(x, x), want G ∈ G. Nu kunnen we de gelijkheid nagaan:
2 + 1
B(y − x, x)B(y − x, y − x) = 2 +
B(y − x, y) − B(y − x, x) B(y − x, x) =
1 + B(y − x, y) B(y − x, x) = 1 −
B(x, y) − B(y, y) B(x, y) − B(x, x) = 0 Het voorgaande kunnen we samenvatten in de volgende stelling.
Stelling 3.1.2 (Uit [2]). Laat Fn uitgerust zijn met een (anti)symetrische bilineaire vorm.
Dan bestaat er voor elke x en y, met B(y − x, x) 6= 0 en B(y, y) = B(x, x) een G-reflector G zodat Gx=y. Deze G is uniek en kan geschreven worden als
G = I + 1
B(y − x, x)(y − x)(y − x)
TM.
In het geval B anti-symmetrisch is, hebben we dat B(y − x, x) = B(y, x).
3.1.2 Opspansel van G-reflectoren
We zullen nu bewijzen dat de elke (quasi-)symplectische matrix te schrijven is als product van G-reflectoren. Stellingen over (quasi-)symplectische matrices kunnen dan namelijk bewezen worden met behulp van deze ’makkelijk’ te begrijpen matrices. Ook hebben we dan zekerheid dat we in algoritmes die gebruikmaken van (quasi-)symplectische matrices alleen maar de G-reflectoren hoeven te gebruiken. Hier volgen we het bewijs uit [3], maar dan voor het in dit geval makkelijker te begrijpen quasi-symplectische product. Tevens veralgemeniseren we het resultaat naar alle niet-gedegenereerde anti-symmetrische vormen.
We hebben eerst een aantal lemma’s nodig. De rest van deze paragraaf zullen we werken met het quasi-symplectische product, zie Definitie 2.4.5. De G-reflectoren horen dus ook bij het quasi-symplectische product. De volgende twee lemma’s komen uit [3].
Lemma 3.1.3. Voor elke x ∈ Fn− {0} bestaan er twee (of ´een) G-reflectoren zodat G2G1x =
e1.
Bewijs. Als B(e1− x, x) = B(e1, x) 6= 0 dan weten we uit de voorgaande stelling dat er zo’n
G-reflector bestaat. Dus we moeten ervan uitgaan dat B(e1, x) = 0, dus x2 = 0. We weten
ook dat x 6= 0, dus er is een kleinste entry ongelijk aan nul, zeg xi 6= 0. Dus x heeft een
quasi-symplectisch product ongelijk aan nul met ei+1 of ei−1 (dit ligt eraan of i deelbaar is
door 2 of niet). We kunnen x dus afbeelden op y = e2+ ei+1of y = e2+ ei−1, want B(x, y) is
dan ongelijk aan nul. Vervolgens geldt dat B(y, e1) 6= 0, dus vinden we een weer met behulp
van Stelling 3.1.2 een G-reflector die y op e1 afbeeldt.
Lemma 3.1.4. Stel G is een G-reflector zodat Gx = y. Dan geldt dat Ge1 = e1 dan en slechts
dan als x2 = y2.
Bewijs. =⇒ Stel Ge1= Ge1, dan x2 = B(e1, x) = B(Ge1, Gx) = B(e1, y) = y2.
⇐= We weten dat de G-reflector die x op y afbeeldt gegeven wordt door G = I + 1
B(y, x)(y − x)(y − x)
TJ .ˆ
Dus als x2 = y2 dan vinden we dat Ge1 = e1+ B(y,x)(y−x)B(y − x, e1) = e1, want B(y − x, e1) =
y2− x2.
Opmerking 3.1.5. Merk op dat we het voorgaande lemma ook mogen veralgemeniseren. Er geldt dat Ge2i−1= e2i−1 dan en slechts dan als x2i= y2i, voor een G-reflector zodat Gx = y.
Stelling 3.1.6 (Uit [3]). Elke symplectische matrix is te schrijven als het product van eindig veel (maximaal 4n) G-reflectoren.
Bewijs. Dit bewijs zal constructief zijn. Laat S een symplectische matrix. We zullen S links vermenigvuldigen totdat we de identiteit hebben. In Lemma 3.1.3 zagen we al dat we met twee G-reflectoren de eerste kolom kunnen afbeelden op e1. Noem deze nieuwe matrix S2
Omdat het resultaat nog steeds symplectisch is weten we dat (S2)2,2 = B((S2)1, (S2)2) = 1.
Nu zijn er twee gevallen
• (S2)1,2 6= 0. In dit geval kunnen we gelijk een G-reflector vinden die de tweede kolom
op e2 afbeeldt. Merk op dat de eerste kolom dan met Lemma 3.1.4 behouden blijft.
• (S2)1,2 = 0. In dit geval kunnen we eerst de tweede kolom op e1+ e2 afbeelden om in
het eerste geval te belanden, want B((S2)2, e1+ e2) = −1.
Noem deze nieuwe matrix nu S3. We hebben nu dit:
S3= 1 0 ~v1 0 1 ~v2 ~0 ~0 S˜
Merk op dat ~v1 en ~v2 allebei nul moeten zijn, want het quasi-symplectische product van die
kolommen moet nul zijn met de eerste twee kolommen. Dus S3= 1 0 ~0 0 1 ~0 ~0 ~0 ˜S
We kunnen nu met inductie verder gaan met ˜S totdat onze matrix gelijk is aan de identiteit. Hiervoor hebben we wel nodig dat een matrix van de vorm
I2 0
0 G
met G een G-reflector ook weer een G-reflector is. Dit is het geval. Stel namelijk dat G = I + βuuTJ . Dan is de bovenstaande matrix gelijk aan I + β ˆuˆuTJ , met ˆu = (0, 0, uT)T. Nu
hebben we dus dat Qm
i=1GiS = I. De inverse van een G-reflector is ook een G-reflector. Dus
S =Q1
i=mG −1
i laat zien dat G opgespannen wordt door G-reflectoren.
Gevolg 3.1.7. Elke automorfismegroep horend bij een anti-symmetrische bilineaire vorm op R2n×2n wordt opgespannen door zijn G-reflectoren.
Bewijs. Dit volgt uit het feit dat elke anti-symmetrische vorm B equivalent is met de quasi-symplectische vorm. Stel dus dat M de matrixrepresentatie is voor B. Dan weten we dat er een matrix C is zodat CTM C = ˆJ . Zo zien we dus als G een G-reflector is voor het quasi-symplectische product, dat dan CGC−1 een G-reflector is voor B, want CGC−1 is een rang 1 verstoorde identiteitmatrix. Een element Q uit de automorfismegroep kunnen we ook transformeren naar S := C−1QC, een quasi-symplectische matrix. Dus S = G1· · · Gn, met
Gi een G-reflector voor elke i. Dus we hebben
Q = CG1· · · GnC−1 = (CG1C−1)(CG2C−1) · · · (CGnC−1),
Hier een veralgemenisering van een stelling uit [3].
Gevolg 3.1.8. Elke matrix uit een antisymmetrische automorfismegroep heeft determinant 1.
Bewijs. Alternatief. Om dit te bewijzen hoeven we alleen maar aan te tonen dat een quasi-sympectische G-reflector determinant 1 heeft. Voor de andere vormen hebben we dan det(C−1GC) = 1, wat het bewijs oplevert.
We kijken naar de eigenwaardes van een G-reflector G: G = I + βuuTJ .ˆ
Allereerst geldt dat er een lineaire deelruimte van dimensie 2n − 1 is loodrecht op de vector u. Voor een v hierin geldt dus dat
Gv = v + βu(uTJ v) = v + 0 = v.ˆ
Deze vectoren zijn daarom allemaal eigenvectoren van G met eigenwaarde 1. Daarnaast geldt dat als λ een eigenwaarde is van G dat λ1 ook een eigenwaarde is voor G. Dus de laatste eigenwaarde kan alleen maar 1 zijn. De determinant is het product van alle eigenwaarden en is dus gelijk aan 1.
3.2 Ortho-symplectische matrices
In deze paragraaf gaan we eerst zoeken naar symplectische matrices die ook orthogonaal zijn, oftewel ortho-symplectische matrices. Omdat niet elke symplectische matrix orthogonaal is zullen we ook zoeken naar de overige symplectische matrices die niet orthogonaal zijn. We zullen in dit hoofdstuk ook aandacht schenken aan wat we met deze matrices kunnen doen. Dit blijkt allemaal genoeg te zijn om de algoritmes in de volgende hoofdstukken te kunnen uit-voeren. Deze paragraaf is gebaseerd op een gedeelte uit [5]. Daar worden de lemma’s van deze paragraaf gepresenteerd als een algoritme zonder bewijs. Wij presenteren hier wel de bewijzen. In deze paragraaf zullen we werken met het gewone symplectische product. Allereerst het volgende lemma over de vorm van een ortho-symplectische matrix. Dit lemma uit [5] is ook meteen de reden dat we werken met symplectische matrices in deze paragraaf, de vorm van ortho-symplectische matrices is namelijk overzichtelijk.
Lemma 3.2.1. Voor een matrix Q =
Q11 Q12
Q21 Q22
met Qij ∈ Rn×n geldt dat
Q is ortho-symplectisch ⇐⇒ Q11= Q22, Q12= −Q21 en QTQ = I
Bewijs. Alternatief. =⇒ We weten dat (−J QTJ )Q = I, dus als Q orthogonaal is dan moet gelden dat QT = −J QTJ . De rechterkant kunnen we uitschrijven
−J QTJ = −J−Q T 21 QT11 −QT 22 QT12 = Q T 22 −QT12 −QT 21 QT11 . Dus we zien precies onze voorwaarden terugkomen.
We zullen nu de verschillende soorten ortho-symplectische matrices defini¨eren. Met dit lemma kunnen we telkens bewijzen dat de matrices ook daadwerkelijk ortho-symplectisch zijn.
3.2.1 Givens-rotaties
Een symplectische Givensrotatie is van de vorm G(k, θ) = C S
−S C
,
waarbij C = diag(1, . . . , 1, cos θ, 1, . . . , 1) en S = diag(0, . . . , 0, sin θ, 0, . . . , 0). De cos θ en sin θ zitten beide op de (k, k)-de entry van C resp. S. Duidelijk is uit Lemma 3.2.1 dat G(k, θ) ook echt ortho-symplectisch is.
De Givensrotaties kunnen we net als de G-reflectoren gebruiken om een vector op een andere vector af te beelden. Dit zal alleen in minder gevallen kunnen. Een Givensrotatie bevindt zich namelijk in een twee-dimensionaal vlak. Het volgende lemma beschrijft hoe we deze matrices later zullen gebruiken.
Lemma 3.2.2. Voor x ∈ R2n geldt dat er een symplectische Givensrotatie is die het n + k-de entry van x op nul afbeeldt, met k een geheel getal tussen 1 en n. Namelijk G(k, θ), met cos θ = xk
||(xk,xn+k)||.
Bewijs. Merk op dat geldt 1 ||(xk, xn+k)|| xk xn+k −xn+k xk xk xn+k =||(xk, xn+k)|| 0 .
De entries van deze Givensmatrix moeten we nu op de juiste plek zetten om een symplectische Givensrotatie te krijgen we vinden precies G(k, θ), met cos θ = xk
||(xk,xn+k)||.
3.2.2 Symplectische Householder-transformaties
Een symplectische Householder-transformatie is afgeleid van de ’normale’ Householder-reflectie. Een Householder-reflectie is een orthogonale matrix van de vorm P = I − 2wwT/wTw. Een Householder-reflectie is een matrix die vectoren in de Rn spiegelt op een n − 1-dimensionaal vlak. Stel namelijk dat v ⊥ w, dan hebben we dat P v = v − (2/wTw)wwTv = v, dus de n − 1-dimensionale deelruimte hwi⊥ wordt op zichzelf afgebeeld. De vector w wordt op door P op −w afgebeeld.
Een symplectische Householder-transformatie is gegeven door
H(k, w) = Ik−1 0 0 P 0 0 Ik−1 0 0 P ,
waarbij P = I − 2wwT/wTw een Householder-reflectie is. Weer is het duidelijk dat deze matrix aan de voorwaarden van Lemma 3.2.1 voldoet.
Lemma 3.2.3. Voor x ∈ R2n is er een Householder-reflectie die de k + 1 tot en met de n-de entries op nul afbeeldt, met k een geheel getal tussen 1 en n. Namelijk H(k, w) waarbij w1= xk+ sign(xk)α, α = ||(xk, . . . , xn||, en wi= xk+i−1 voor i = 2 . . . n − k. In de expressie
voor w1 kiezen we voor sign(xk) voor een betere stabiliteit in algoritmen.
Bewijs. We willen het product H(k, w)x uitrekenen. We zijn alleen ge¨ınteresseerd in de entries k t/m n. Dus we hoeven alleen P x uit te rekenenen. Definieer y = (xk, . . . xn) Merk op dat
x = 1/2xk+ 1/2||y|| 1/2x2 .. . 1/2xn + 1/2xk− 1/2||y|| 1/2x2 .. . 1/2xn
waarbij de eerste vector in een veelvoud is van w en de tweede vector loodrecht staat op w. Met de eerdere opmerkingen over de Householder-refecties zien we nu in dat de stelling waar is.
Meer ortho-symplectische matrices kunnen we niet vinden De symplectische Householder-reflecties en de symplectische givens-rotaties brengen alle ortho-symplectische matrices voor. Zie hiervoor [5]. In dat bewijs wordt gebruikgemaakt van een symplectische singuliere waarden decompositie.
3.2.3 Gauss-transformaties
Niet alle symplectische matrices zijn ook orthogonaal. Daarom is het van belang om een derde verzameling aan voortbrengers voor de symplectische groep te vinden. Deze matrices heten Gauss-transformaties.
Een Gauss-transformatie is gegeven door
G(k, v) =D V 0 D−1
,
waar V alleen ongelijk is aan nul in de entries (k, k − 1) en (k − 1, k), die entries zijn gelijk aan
v (1 + v2)1/4,
en D is een diagonaalmatrix gegeven met diagonaalelementen gelijk aan 1 behalve in de entries (k − 1, k − 1) en (k, k), die entries zijn gelijk aan
1 (1 + v2)1/4.
De Gauss-reflecties gaan we op de manier als in het volgende lemma gebruiken.
Lemma 3.2.4. Zij x ∈ R2n met xn+k−1 6= 0 dan is een Gauss-transformatie die het k-de
Bewijs. We rekenen G(k, v)x uit. Omdat G(k, v) grotendeels de identiteit is zullen alleen de entries k − 1, k, n + k − 1 en n + k veranderen. We zijn ge¨ınteresseerd in het k-de entrie. Er geldt dat (G(k, v)x)k= (1+vxk2)1/4 +
vxn+k−1
(1+v2)1/4 =
xk−xk
(1+v2)1/4 = 0.
Een terechte opmerking is dat er ook andere matrices zijn die hetzelfde kunnen bewerk-stelligen als de Gauss-transformatie. Bekijk hiervoor wel het algoritme in hoofdstuk 4, de Gauss-transformatie tast daar de reeds gecre¨erde nullen niet aan. Een voordeel van de Gauss-transformatie is dat de inverse gemakkelijk op te schrijven is. We hadden mede door tijdgebrek niet in andere artikelen gevonden of dit de meest stabiele matrix is die dit werk op deze wijze verricht.
4 De SR-decompositie
Zoals we al zagen zijn er voor de R2n zogeheten symplectische bases {e1, . . . , en, f1, . . . , fn}.
Het blijkt zo te zijn dat we meestal, met behulp van een proces dat lijkt op het Gram-Schmidt proces, een symplectische basis kunnen maken vanuit een ’gewone’ basis. In dit proces ontstaat er een symplectische matrix (de basis) S en als bijproduct een bovendriehoeksmatrix R. In dit hoofdstuk zullen we dit eerst wat preciezer maken. Daarnaast zullen we ook onderzoeken hoe dit proces stabiel stabiel gemaakt kan worden voor computeralgoritmes. Dit heeft onder andere te maken met de vrije keuze van een aantal parameters.
4.1 Symplectische Gram-Schmidt
In dit paragraaf zal B de quasi-symplectische vorm zijn. Met ˆJ bedoelen we de matrix
ˆ J = 0 1 0 · · · 0 −1 0 0 · · · 0 0 0 . .. .. . ... 0 1 0 0 −1 0 ,
en we zullen werken in de vectorruimte R2n. Het algoritme wat we zullen onderzoeken is het quasi-symplectische Gram-Schmidt algoritme. Hier zullen we van de basis/matrix A (zie dit als de basis α = {A1, . . . , A2n}) een quasi-symplectische basis/matrix S ∈ G maken. We zullen
nu in elke stap ervan uitgaan dat het algoritme niet vastloopt. Later zullen we terugkomen op de voorwaarden voor het niet vastlopen van het algoritme. Allereerst de definitie van de SR-decompositie, het gewenste resultaat van ons algoritme.
Definitie 4.1.1 (SR-decompositie; quasi-symplectisch). Zij A ∈ R2n een matrix. Dan heet A = SR, met S quasi-symplectisch en R ∈ R2n bovendriehoeks een (niet unieke) SR-decompositie van A.
Opmerking 4.1.2. Voor elke i ≤ n geldt dat span{A1, . . . , Ai} = span{S1, . . . , Si}, waarbij
Aien Side i-de kolom is van A resp. S. Dus we kunnen concluderen dat ImA(span{e1, . . . , ei})
= ImS(span{e1, . . . , ei}) voor elke i ≤ n. Het is ook duidelijk dat we de kolommen van S
maken met behulp van eerdere kolommen van A.
Opmerking 4.1.3. De SR decompositie in het symplectische geval is afgeleid van het quasi-sympletische geval. Stel dat A = SR met S quasi-symplectisch en R bovendriehoeks. Dan geldt ook dat A = (PTSP )(PTRP ) een decompositie van A met PTSP symplectisch en waarbij PTRP eruit ziet als.
PTRP = ∗ ∗ 0 ∗ ∗ ∗ 0 ∗ 0 ∗ 0 0 ∗ ∗ 0 ∗ .
De matrix PTRP wordt ook wel J -bovendriehoeks genoemd. We hebben dus de volgende definitie.
Definitie 4.1.4 (SR-decompositie; symplectisch). Zij A ∈ R2neen matrix. Dan heet A = SR, met S symplectisch en R ∈ R2n J -bovendriehoeks een (niet unieke) SR-decompositie van A.
We gaan nu kijken hoe we in het quasi-symplectische geval een SR-decompositie kunnen vinden van een 4 × 4 matrix. Het algoritme dat we hier zo zien is het quasi-symplectische Gram-Schmidt algoritme. In het 4 × 4 zullen we eigenlijk alles al zien wat van belang is voor dit algoritme. Kies dus een A ∈ R4×4. We gaan nu een quasi-symplectische matrix S construeren die dus voldoet aan de definitie. Tegelijk construeren we de bijbehorende R.
Eerste vector S1. Hier kiezen de eerste vector van A, namelijk A1. We mogen deze vector
vermenigvuldigen met een willekeurige constante ongelijk aan nul, er liggen namelijk geen op de lengte van deze vector. Dus S1 = A1/R11.
Tweede vector S2. Deze vector moet bestaan uit de eerste twee vectoren van A en er
moet gelden dat B(S1, S2) = 1. Dus we kunnen S2 = −R12RS221+A2 kiezen, waarbij R22 =
B(S1, −R12S1+ A2) = B(S1, A2). Merk hierbij op dat we R12 kunnen kiezen zoals we willen.
Het algoritme loopt hier vast als B(S1, A2) = 0.
Derde vector S3. Deze vector moet orthogonaal zijn met de eerste twee kolommen van S.
We zullen dus beginnen door een ˜S3te defini¨eren die orthogonaal is op S1 en S2. Definieer dus
S3 = 1/R33(A3− B(S1, A3)S2+ B(S2, A3)S1), de lezer kan controleren dat deze S3 volstaat.
De entries van R zijn dan R13 = −B(S2, A3), R23 = B(S1, A3) en R33 willekeurig ongelijk
aan nul. Hier kan het algoritme duidelijk niet vastlopen als A inverteerbaar is.
Vierde vector S4. Ook deze vector moet orthogonaal zijn met de eerste twee vectoren van A.
Ook moet gelden dat B(S3, S4) = 1. Kies daarom eerst ˜S4 = A4− B(S1, A4)S2+ B(S2, A4)S1.
Vervolgens zorgen we voor B(S3, S4) = 1. Dit doen we door S4 = −CSR344+ ˜S4 te kiezen. De
entries van R zijn dan R14 = −B(S2, A4), R24 = B(S1, A4) en R34 = C met C
willekeu-rig en R44 = B(S3, ˜S4) 6= 0. Hier merken we weer op dat het algoritme vastloopt wanneer
B(S3, ˜S4) = 0.
Hieronder nog een keer een beschrijving van het algoritme zonder de vorming van R voor algemene 2n × 2n matrices.
Algoritme 4.1.5. Input: matrix A ∈ R2n, output:S quasi-symplectisch en een matrix R ∈
R2n bovendriehoeks.
Iteratie 1: Laat S1 = C1A1 zijn. Dan S2= B(SA2+C1,A2S21).
Iteratie 2: Laat S3 = C3(A3− B(S1, A3)S2+ B(S2, A3)S1) en ˜S4= A4+ C4S3− B(S1, A4+
C4S3)S2+ B(S2, A4+ C4S3)S1. Dan S4 = ˜ S4
B(S3, ˜S4.
Iteratie i: Laat S2i−1= C2i−1(A2i−1−Pi−1j=1B(S2j−1, A2i−1)S2j+Pi−1j=1B(S2j, A2i−1)S2j−1)
en ˜S2i= A2i+C2iS2i−1+Pi−1j=1B(S2j−1, A2i +C2iS2i−1)S2j+Pi−1j=1B(S2j, A2i +C2iS2i−1)S2j−1. Vervolgens S2i= ˜ S2i B(S2i−1, ˜S2i).
4.2 Onderzoek naar existentie en uniciteit
In deze paragraaf zullen we ingaan op de vraag wanneer er een SR-decompositie bestaat en of die dan uniek is. We zullen zien dat er geen uniciteit is van de decompositie. In het laatste paragraaf van dit hoofdstuk zullen we ingaan op de gevolgen hiervan. Allereerst een opmerking over het Gramm-Schmidt algoritme van de vorige paragraaf.
Opmerking 4.2.1. Merk op dat het algoritme vastloopt dan en slechts dan als geldt dat B(S2i−1, ˜S2i) = 0. We moeten dan namelijk door nul delen. Deze voorwaarde is onafhankelijk
van de gekozen parameters, door het komende lemma. De voorwaarde is op het eerste gezicht moeilijk vooraf te controleren. De stelling na het lemma biedt ons uitkomst door te laten dat een voorwaarde op de determinant van [ATJ A]ˆ 2ii,j=1 equivalent is met wat we hier zien. Deze stelling zullen we later nog vaker gaan gebruiken.
Lemma 4.2.2. In het quasi-symplectische Gram-Schmidt algoritme geldt dat B(S2i−1, ˜S2i) =
0 voor de parameters C2i= 0 en C2i−1= 1 met i = 1 . . . n dan en slechts dan B(S2i−1, ˜S2i) = 0
voor alle mogelijke keuzes van de parameters Ci i = 1 . . . 2n.
Bewijs. Allereerst bewijzen we dat A2i+Pi−1j=1B(S2j−1, A2i)S2j+Pi−1j=1B(S2j, A2i)S2j−1
al-tijd gelijk is onafhankelijk van de gekozen parameters. Dit volgt uit Lemma 2.4.2, want op span{S1. . . S2(i−1)} = span{A1. . . A2(i−1)} is de quasi-symplectische vorm niet-gedegenereerd,
dus is er een unieke directe som van A2i. Nu kunnen we het lemma bewijzen. Dus neem een
stel parameters en bereken B(S2i−1, ˜S2i) en laat ˆS2i−1 = 1/CS2i−1 horen bij de eerste
para-meterkeuze in ons lemma.
B(S2i−1, ˜S2i) = B(C ˆS2i−1, A2i+ C2iC ˆS2i−1+ i−1 X j=1 B(S2j−1, A2i)S2j + i−1 X j=1 B(S2j, A2i)S2j−1) = CB( ˆS2i−1, A2i+ i−1 X j=1 B(S2j−1, A2i)S2j+ i−1 X j=1 B(S2j, A2i)S2j−1).
Dus B(S2i−1, ˜S2i) is op een constante ongelijk aan nul na gelijk voor elke keuze van parameters
Ci i = 1 . . . 2n.
Het eerste gedeelte van de volgende stelling komen onder andere voor in [1] en [10]. Wij geven hier ook het bewijs.
Stelling 4.2.3. Voor een inverteerbare matrix A ∈ R2n bestaat er een SR-decompositie dan en slechts dan als [ATJ A]ˆ 2mi,j=1inverteerbaar is voor elke 1 ≤ m ≤ n. De keuze van parameters is niet van belang in algoritmes.
Bewijs. Alternatief. We zagen in het vorige lemma al dat de gekozen parameters in het algoritme niet uitmaken. =⇒ Stel dat A = SR met S ∈ G en R bovendriehoeks. Dan geldt dat
[ATJ A]ˆ 2mi,j=1= [RTSTJ SR]ˆ 2mi,j=1= [RTJ R]ˆ 2mi,j=1= [RT]2mi,j=1[ ˆJ ]2mi,j=1[R]2mi,j=1,
want STJ S = ˆˆ J en voor de laatste stap gebruiken we dat R bovendriehoeks is. Nu zijn [R]2mi,j=1, [RT]2mi,j=1en [ ˆJ ]2mi,j=1inverteerbaar dus hun product ook. Dus [ATJ A]ˆ 2mi,j=1is inverteerbaar voor elke 1 ≤ m ≤ n.
⇐= Deze kant is moeilijker en langer. We gaan bewijzen dat er een bovendriehoeksmatrix R is zodat ATJ A = Rˆ TJ R. Dan weten namelijk dat R inverteerbaar is (Rˆ TJ R is inverteerbaar)ˆ en AR−1∈ G, wat komt door
(AR−1)TJ ARˆ −1= R−TATJ ARˆ −1 = R−TRTJ RRˆ −1= ˆJ . Dus A = (AR−1)R is dan een SR-decompositie van A.
We gaan nu dus voor een A zijn een bovendriehoekse R construeren die voldoet aan ATJ A =ˆ RTJ R. We gaan R stapsgewijs ’bouwen’. We zullen de kolommen van R van links naar rechtsˆ afwerken. Elke kolom zelf zullen we van boven naar beneden construeren. Dit zal goed gaan, door in te zien dat we hebben dat [ATJ A]ˆ ij = RTi J Rˆ j = B(Ri, Rj) voor de gewenste R. We
maken voor de derde en de daaropvolgende kolommen telkens gebruik van het feit dat R geen diagonaalentry’s ongelijk aan nul mag hebben omdat hij inverteerbaar is. De lezer kan nagaan dat dit allemaal goed gaat.
Opmerking 4.2.4. Ter verduidelijking. De matrix ATJ A kunnen we schrijven alsˆ B(A1, A1) . . . B(A1, A2n) .. . . .. ... B(A2n, A1) . . . B(A2n, A2n) .
We zien hier ook, zoals we al wisten zagen in de discussie van het algoritme, dat een voorwaarde op [ATJ A]ˆ 2mi,j=1 zegt dus alleen iets over de eerste 2m kolommen. Ook zijn de diagonaalentries van deze matrix gelijk aan nul.
Voorbeeld 4.2.5. Eerst voorbeeld wat we later zullen gebruiken van een matrix waarvoor de SR-decompositie niet bestaat:
A = a11 a12 a13 a14 a15 a16 0 a22 a23 a24 a25 a26 0 0 a33 a34 a35 a36 0 0 0 0 a45 a46 0 0 0 a54 a55 a56 0 0 0 a64 a65 a66 .
Om in te zien dat voor deze matrix er geen SR-decompositie is kunnen we naar de determi-nanten van [ATJ A]ˆ 2mi,j=1 kijken. In dit geval geldt dat [ATJ A]ˆ 4i,j=1 niet inverteerbaar is. Dit volgt doordat [ATJ A]ˆ 4i,j=1= BTJ B, metˆ
B = a11 a12 a13 a14 0 a22 a23 a24 0 0 a33 a34 0 0 0 0 .
Omdat B niet inverteerbaar is kan BTJ B ook niet inverteerbaar zijn (alternatief: merk op datˆ de matrix bestaande uit de eerste vier kolommen van A getransponeerd een niet-nul vector in de kern heeft).
In [1] en [10] wordt niet ingegaan op de situatie wanneer A singulier is. We kunnen hier wel iets zeggen over zeggen over het bestaan van een SR-decompositie voor zo’n A. Dit is echter ingewikkelder doordat voor deze matrix ATJ A altijd singulier is en dus is het moeilijker teˆ controleren of we in ons algoritme door nul moeten delen. Het volgende analogon met Stelling 4.2.3 kunnen we opstellen, maar dan zonder een bi-implicatie.
Lemma 4.2.6. Zij A ∈ R2n×2n singulier. Als A een quasi-symplectische SR-decompositie
heeft dan geldt dat ATJ T = Rˆ TJ R voor een bovendriehoekse matrix R. Voor deze R geldt datˆ de rang van de eerste i kolomen van R en A gelijk zijn voor elke i.
Met dit lemma kunnen we voor de volgende matrix bewijzen dat die geen SR-decompositie heeft.
Voorbeeld 4.2.7. Zij A ∈ R4×4 singulier
A = 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 2 0 0 0 0 0 0 1 .
We kunnen nu de eerste vier kolommen van ATJ A berekenenˆ
[ATJ A]ˆ 4i,j=1= 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 .
We kunnen nu proberen een R te construeren die voldoet aan de eisen van het lemma hier-boven. We komen na wat stappen uit op
R = C1 C1 ∗ 0 0 0 0 0 0 C2 0 0 0 .
De constantes zijn ongelijk aan nul en de sterretjes mogen willekeurig zijn. Het lukt nu echter niet om een vierde kolom toe te voegen zodat R voldoet aan de eisen van het lemma, namelijk dat de rang 3 is en dat RTJ R = Aˆ TJ A. Dus A heeft geen SR-decompositie.ˆ
Een omkering van Lemma 4.2.6 is moeilijker te onderzoeken. Het volgende voorbeeld illu-streert wat we wel kunnen zeggen
Voorbeeld 4.2.8. Zij A ∈ R6×6 singulier en laat R ∈ R6×6 zodat RTJ R = Aˆ TJ A. Stel datˆ
R = 1 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 2 0 0 0 0 0 0 1 .
We kunnen dus de lineair afhankelijke kolom R2 vervangen en daarmee R singulier maken.
Dit doen we, definieer dus
ˆ R = 1 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 2 0 0 0 0 0 0 1 .
Dan kunnen we bewijzen dat we de tweede kolom van A kunnen vervangen zodat ˆATJ ˆˆA = ˆ RTJ ˆˆR. ˆ RTJ ˆˆR = 0 1 0 1 0 0 −1 0 0 −1 0 0 0 0 0 1 1 0 −1 1 −1 0 0 0 0 0 −1 0 0 2 0 0 0 0 −2 0 .
We willen namelijk een kolom ˆA2 vinden zodat AT2J A = −1ˆ ∗ 0 −1 0 0. Deze kolom
bestaat, dit kan je inzien door weer de tweede kolom van A te vervangen zodat de matrix inverteerbaar is. Dan houden we een gewone matrix-vector vergelijking over. Voor ˆA bestaat er dan vervolgens een SR-decompositie. Vanuit deze SR-decompositie kunnen we vervolgens door ˆR met R te vervangen een SR-decompositie voor A extraheren.
Als R niet van deze vorm zou zijn, kunnen we nog iets proberen. We kunnen een quasi-symplectische matrix D zoeken zodat DR wel van de gewenste vorm is. Merk op dat dan RTDTJ DR = Rˆ TJ R = Aˆ TJ A. Maar dan moet die wel bestaan. Als deze altijd bestaatˆ kunnen we de omkering van Lemma 4.2.6 opstellen.
Opmerking 4.2.9. Merk op dat als we de lineair afhankelijke kolom van A kunnen vervangen zodat ˜A een SR-decompositie heeft dat A dan ook een SR-decompositie heeft. Als we de bi-implicatie van daarnet kunnen bewijzen, zal dit ook een bi-bi-implicatie zijn.
Opmerking 4.2.10. Het is terecht om op te merken dat het controleren van de voorwaarden voor het bestaan van de SR-decompositie erg duur is voor een hogerdimensionale A. In de praktijk mag je er echter er wel vanuit gaan dat de voorwaarden voldaan zijn. Er is namelijk een kans van nul dat een determinant gelijk aan nul is. Ook in het singuliere geval is de kans erg klein dat er geen SR-decompositie bestaat. Wel is het natuurlijk verstandig om te kijken naar de context van de matrix, het kan zo zijn dat de logische uitkomst van een probleem een matrix is waarvoor geen SR-decompositie bestaat. We zullen ons hier niet verder mee bezig houden.
In ons algoritme zagen we ook dat er veel vrije parameters waren. Een andere keuze voor een parameter zal logischerwijs een andere SR-decompositie opleveren. De volgende stelling met bewijs uit [1] geeft precies aan hoe de verschillende SR-decomposities aan elkaar gerelateerd zijn.
Stelling 4.2.11. Laat A een matrix zijn zodat A = SR een quasi-symplectische SR-decompositie is. Dan geldt dat S2R2een quasi-symplectische SR-decompositie is dan en slechts
Bewijs. =⇒ . Stel dat A = S2R2 een SR-decompositie is. We hebben dan dat SR = S2R2,
wat we kunnen omschrijven tot S = S2R2R−1 en R = S−1S2R2. Het eerste impliceert dat
D een bovendriehoeksmatrix is en het tweede dat hij ook quasi-symplectisch is. Dat D ook bidiagonaal is volgt uit het feit dat een bovendriehoekse matrix die quasi-symplectisch is, noodzakelijk bidiagonaal is. De laatste kolom moet bijvoorbeel orthogonaal zijn met de eerste kolom, dit kan alleen als de tweede entry van de laatste kolom nul is. Met inductie kunnen de nullen boven de tweede diagonaal aangetoond worden.
⇐= . Als L bidiagonaal en quasi-symplectisch is dan is SL quasi-symplectisch en L−1R bovendriehoeks, dus (SL)(L−1R) is ook een SR-decompositie.
Opmerking 4.2.12. Een bidiagonale quasi-symplectische matrix is altijd van de vorm λ1 b1 . . . 0 0 0 1/λ1 0 . . . 0 0 λ2 b2 ... ... 0 1/λ2 .. . ... . .. 0 0 . . . λn bn 0 0 . . . 0 1/λn .
De vrijheid in de keuze van de λi’s komt hier overeen met de schaling van de eerste vector
in elke paar. De keuze van de bi’s komt overeen met de keuze van het optellen van de eerste
vector in een paar bij de tweede.
4.3 Symplectische triangularisatie met G-reflectoren
In het bewijs dat de determinant van een quasi-symplectische matrix een is maakten we ge-bruik van G-reflectoren om een quasi-symplectische matrix tot de identiteit te transformeren. Net zo kunnen we een ’gewone’ matrix met behulp van G-reflectoren transformeren naar een bovendriehoeksmatrix (mits er voor die matrix een SR-decompositie bestaat). In dit para-graaf zullen we bekijken hoe dat precies werkt, met aandacht voor het niet vastlopen van de algoritmes. We zullen ook een poging doen om dit SR-algoritme redelijk stabiel te laten werken. In [7] wordt ook dit algoritme gepresenteerd, maar dan voor het symplectische geval. Wij bewijzen hier dat het algoritme niet vastloopt wanneer er voor de betreffende matrix een SR-decompositie bestaat.
4.3.1 Aanpak
Neem aan dat onze matrix A = A(1) voldoet aan de voorwaarde voor een SR-decompositie, namelijk dat [ATJ A]ˆ 2mi,j=1 inverteerbaar is voor elke 1 ≤ m ≤ n.
We willen allereerst proberen de eerste kolom op een veelvoud van e1 te krijgen. Door
Lemma 3.1.3 weten we dat er G-reflectoren zijn die dit werk voor ons verrichten, noem hun product S en noem A(2) = SA. Als we dit gedaan hebben kunnen we het volgende bewijzen voor de tweede kolom van A(2). Namelijk dat (A(2))22 6= 0. Dit komt doordat
[ATSTJ SA]ˆ 2i,j=1 = [ATJ A]ˆ 2i,j=1 inverteerbaar is, dus A(2) = SA heeft een SR-decompositie. Met behulp van Voorbeeld 4.2.5 kunnen we nu concluderen dat (A(2))22 6= 0. Met deze
informatie gaan we nu verder. We hebben de volgende situatie: A(2) = A(2)11 A(2)12 0 A(2)22 C ~0 v B
Om nu de vector v op nul te krijgen kunnen we een G-reflector vinden die de tweede kolom op C A(2)22 0
afbeeldt, met C een willekeurige constante ongelijk aan nul, die aan de voorwaarde in Stelling 3.1.2 is voldoet (alleen voor C = A(2)12 gaat het mis). Merk op dat de eerste kolom behouden blijft door Lemma 3.1.4.
We hebben nu de eerste twee kolommen bovendriehoeks gekregen. We kunnen nu inductief verder gaan met de rest van de kolommen. Om dit te verantwoorden hebben we het volgende lemma. Dit bewijst dat het algoritme uit [7] niet vastloopt.
Lemma 4.3.1. Voor inverteerbare matrices van de vorm A = A11 A12 0 A22 C 0 0 B geldt het volgende:
A heeft een SR-decompositie ⇐⇒ B een SR-decompositie heeft. Wat gedetailleerder: er geldt dat als B = SR dan is
A = I2 0 0 S A11 A12 0 A22 C 0 0 R
een SR-decompositie voor A. Alle SR-decomposities van A zijn van deze vorm.
Bewijs. We gaan de laaste zin van de stelling bewijzen. Stel dus dat A = SR en deel de matrices op in blokmatrices: A11 A12 0 A22 C 0 0 B = S11 S12 S21 S22 R1 M 0 R2 .
Dan vinden we dus dat
B = S21M + S22R2.
Ook zien we dat
Omdat R1 inverteerbaar is volgt dan dat S21= 0. Omdat S quasi-symplectisch is en S21= 0
moet gelden dat S12= 0. Hieruit volgt dan dat S22quasi-symplectisch is. Dus we zien
B = S21M + S22R2= S22R2,
wat een SR-decompositie is voor B.
We zien dus dat we in onze huidige situatie net zo goed verder kunnen werken met de matrix B. Ook hiervan kunnen we de eerste twee kolommen bovendriehoeks krijgen. Zo zien we met behulp van inductie dus dat we met G-reflectoren een SR-decompositie vinden. Over het algemeen zal deze methode beter werken dan het quasi-symplectische Gram-Schmidt algoritme van eerder in dit hoofdstuk. Zie hiervoor de testresultaten.
4.3.2 Conditiegetallen G-reflectoren
In het algoritme hierboven maken we gebruik van vermenigvuldigingen met G-reflectoren. Om dit zo stabiel mogelijk te doen moeten we meer te weten komen over de conditiegetallen van (quasi-)symplectische matrices. Dit gedeelte is gebaseerd op een gedeelte uit [7].
Lemma 4.3.2. Voor een (quasi-)symplectische matrix S ∈ R2n×2n geldt dat κ2(S) = ||S||22.
Bewijs. Alternatief. We bewijzen dit alleen voor symplectische matrices. Allereerst bewijzen we dat ST symplectisch is. We weten dat S−1 symplectisch is, dus
S−TJ S−1= J. Als we nu de inverse nemen aan beide kanten krijgen we
−SJ ST = −J,
dus ST is ook symplectisch.
Om nu het conditiegetal van S te berekenen gaan we de grootste en kleinste singuliere waarde van S berekenen. Dit doen we door de eigenwaarden van STS te bepalen. De grootste eigenwaarde van STS is de grootste singuliere waarde van S in het kwadraat, evenzo is de kleinste eigenwaarde de kleinste singuliere waarde in het kwadraat. Omdat STS symplectisch is weten we dat de eigenwaardes in paren van {λ, 1/λ} komen. Dus we kunnen concluderen dat σ2n = 1/σ1= 1/||S||2, dus κ2(S) = ||S||22.
Met behulp van dit lemma kunnen we bepalen welke G-reflectoren het kleinste conditiegetal hebben in elke stap van ons algoritme. Dit zal hopelijk een redelijk optimaal algoritme opleveren. Omdat ons algoritme in essentie maar twee stappen heeft, zals het niet al te veel werk zijn. Allereerst een lemma uit [7] die een bruikbare conditie geeft wanneer een G-reflector een minimale conditiegetal heeft.
Lemma 4.3.3. Laat T = I + βuuTJ een G-reflector zijn voor de (quasi-)symplectische groep.
Dan is het conditiegetal ten opzichte van de 2-norm gegeven door κ2(T ) =
2 + β2||u||4
2+pβ4||u||82+ 4β2||u||42
2 .
Bewijs. Idee. We zagen al dat κ(T ) = ||T ||2. Dus we berekenen de eigenwaarden van TTT . Vervolgens kunnen we een slimme orthogonale gelijkvormigheidstranformatie toepassen om de eigenwaarden van TTT bloot te leggen.
We kunnen nu constantes vinden zodat de G-reflectoren in elke stap minimale conditiege-tallen hebben. We zullen in het eerste lemma er vanuit gaan dat we de beide stappen in het algoritme in ´e´en keer kunnen doen. Dit betekent dat we maar ´e´en G-reflector uit Lemma 3.1.3 gebruiken. De volgende twee lemma’s komen uit [7].
Lemma 4.3.4. De G-reflector T1 voor de eerste stap heeft minimale conditietgetal t.o.v. de
2-norm wanneer hij de eerste kolom van A afbeeldt op sign(a11)||A1||2e1. De transformatie
wordt dan gegeven door
T1 = I +
1 B(ρe1, A1)
(ρe1− A1)(ρe1− A1)TJ ,ˆ
waar ρ = sign(a11)||A1||2.
Lemma 4.3.5. De G-reflector T1 voor de tweede stap heeft minimale conditietgetal t.o.v. de
2-norm wanneer hij de tweede kolom van A afbeeldt op v = (A12± η)e1+ A22e2, met µ de
2-norm van de vector (a32, a42, . . . , a2n,2). De transformatie wordt dan gegeven door
T1 = I +
1
±µa22(v − A2)(v − A2)
TJ .ˆ
Beide lemma’s worden bewezen door een minimum te vinden voor β2||u||4
2 door af te leiden.
Vervolgens werk je de bijbehorende G-reflector uit. Het volledige bewijs kunt u vinden in [7]. We hebben nog niet de uitzondering behandeld van de eerste stap. Het kan namelijk zo zijn dat B(A1, e1) = 0, dan moesten we eerst een andere G-refector toepassen die A1 op een
vector afbeeldde die wel een symplectisch product had met e1. Deze G-refector is het echter
veel moeilijker optimaal te kiezen. Dit komt doordat de keuzevrijheid voor deze G-reflector ongeveer 2n-dimensionaal is. We hoeven namelijk alleen maar A1naar een vector af te beelden
die een quasi-symplectisch product met e1 heeft. Hier is erg veel keuzevrijheid in. Wel zouden
we nadat we voldoende extra eisen hebben gesteld aan deze vector gemakkelijk een optimale G-reflector kunnen vinden die dit werk doet. Het is overigens natuurlijk ook van belang om die eisen zo te kiezen zodat de volgende G-reflector ook stabiel gekozen kan worden.
4.4 SR-decompositie met ortho-symplectische matrices
We kunnen ook met behulp van de ortho-symplectische matrices een algoritme maken voor de SR-decompositie. Dit algoritme werkt dus hetzelfde als het algoritme van hierboven, maar nu maken we zoveel als mogelijk gebruik van de ortho-symplectische matrices. Deze matrices hebben de eigenschap dat ze gemaakte fouten niet opblazen, maar juist gelijk houden. Wel moeten we ook een aantal Gauss-transformaties gebruiken die niet deze eigenschap hebben. Allereerst een beschrijving van het algoritme met behulp van een 6 × 6 matrix. Dit algoritme wordt ook op deze manier gepresenteerd in [1]. We geven hier ook weer de redenen waarom dit algoritme niet vastloopt.
Merk eerst op dat een de factor R in een symplectische SR-decompositie van de vorm R11 R12
R21 R22
is waarbij Rij ∈ Rn × n bovendriehoeks zijn met R21 zelfs strikt bovendriehoeks.
Laat nu A ∈ R6×6 een matrix zijn waarvoor een symplectische SR-decompositie bestaat
A = x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x
We kunnen met Givens-rotaties achtereenvolgens de entries (6, 1), (5, 1) en (4, 1) gelijk aan nul maken. Merk daarbij op dat bijvoorbeeld G(5, θ) (deze beeldt (5, 1) af op nul) alleen de tweede en vijfde kolommen en rijen aanpast. Dus de nullen die eerder gecree¨erd waren blijven behouden. Dus nu A = x x x x x x x x x x x x x x x x x x 0 x x x x x 0 x x x x x 0 x x x x x .
Om nu verder te gaan gebruiken we een symplectische Householder-transformatie om de entries (2, 1) en (3, 1) op nul af te beelden, zie Lemma 3.2.3. Weer blijven de nullen gecree¨erd met de Givens-rotaties behouden door de vorm van de Householder-transformatie. Nu hebben we dus A = x x x x x x 0 x x x x x 0 x x x x x 0 x x x x x 0 x x x x x 0 x x x x x .
Vervolgens kunnen we met Givens-rotaties nullen cre¨eren op de entries (6, 4) en (5, 4), het entrie (5, 3) kunnen we niet op nul krijgen zonder de nul in (6, 1) te verliezen. Daarna kunnen we weer een Householder-reflectie toepassen om entrie (3, 4) op nul te krijgen (weer kunnen we niet meer doen dan dit). We hebben dan
A = x x x x x x 0 x x x x x 0 x x 0 x x 0 x x x x x 0 x x 0 x x 0 x x 0 x x .
Als laatste passen we een Gauss-transformatie toe om de entrie (2, 4) nul te krijgen, Lemma 3.2.4. Hiervoor moeten we wel hebben dat a446= 0. Maar als a246= 0 en a44= 0, dan kunnen
laaste stap mogen we toepassen, we krijgen A = x x x x x x 0 x x 0 x x 0 x x 0 x x 0 x x x x x 0 x x 0 x x 0 x x 0 x x .
We kunnen nu ons probleem reduceren door verder te werken met de matrix ˆA die de eerste en de vierde kolom mist en tevens ook de eerste en vierde rij. We hebben namelijk te maken met de symplectische situatie van Stelling 4.3.1.
Dit proces kunnen we generaliseren naar matrices met de afmetingen 2n × 2n voor n ∈ N. De code hiervoor staat in de bijlage.
4.5 Optimale conditie S en R
In dit gedeelte zullen we ingaan op een vraag die ontstaat door Stelling 4.2.11. Doordat de SR-decompositie niet uniek is kunnen we ons afvragen of we op een gemakkelijke wijze ons algoritme zo kunnen inrichten dat ´e´en of beide factoren S en R goed geconditioneerd zijn. Deze vraag is bijvoorbeeld van belang in het SR-algoritme van het volgende hoofdstuk. Hier conjugeren we een matrix met symplectische matrices, waarbij het behouden van het spectrum van essentieel belang is. Hier is het dus fijn als de symplectische matrices goed geconditio-neerd zijn.
Voor een soortgelijk probleem is al een respectabele oplossing gevonden. In dat probleem mogen we een matrix A vermenigvuldigen met een diagonaalmatrix D (in ons probleem is D bidiagonaal en symplectisch). De vraag is dan hoe we het conditiegetal van A klein kunnen krijgen met vermenigvuldiging van zo’n diagonaalmatrix. Een goede manier is om de kolommen van A op gelijke grootte te schalen. De matrix D wordt dus in Matlabtaal D = diag(1/||A1||2, 1/||A2||2, . . . , 1/||An||2). Er wordt in [6] onder andere bewezen dat voor
deze D geldt dat
κ2(AD) ≤
√ n min
D κ2(AD),
We hebben hier dus een relatief makkelijk te construeren diagonaalmatrix, waarbij we ook niet al te ver van het minimum van het conditiegetal zitten.
Voor ons probleem kunnen we soortgelijke bewijsmethodes gebruiken als in [6]. Het idee van deze oplossing is dat als we het maximum van de normen van de kolommen zo klein mogelijk krijgen, dat dan de norm van de hele matrix ook niet al te groot meer is. Het conditiegetal van een symplectische matrix zal dan ook niet al te groot zijn.
Bekijk nu eerst het volgende lemma uit [6].
Lemma 4.5.1. Voor elke matrix A ∈ Rn×n gelden de volgende ongelijkheden. max
i ||Ai||2≤ ||A||2 ≤
√ n max
i ||Ai||2.
Bewijs. Alternatief. De eerste ongelijkheid. Stel de kolom Ai is maximaal, dan hebben we
