Een kwart eeuw iteratieve methoden

1

116

NAW 5/10 nr. 2 juni 2009

Een kwart eeuw iteratieve methoden

Kees Vuik

Kees Vuik Technische Universiteit Delft Faculteit Elektrotechniek, Wiskunde en Informatica Delft Institute of Applied Mathematics Mekelweg 4, 2628 CD Delft [email protected]

Onderzoek

Een kwart eeuw iteratieve methoden Lineaire algebra vormt het wiskundige fundament waarop bijvoorbeeld in computer-animaties de meest fantastische voorstellingen worden gebouwd. Voor de weergave van deze beelden in films of games is een ontzaglijke partij rekenwerk nodig, en er is dan ook doorlopend behoefte aan snellere hardware en slimmere rekenmethoden. Kees Vuik geeft een overzicht van de ontwikkelingen in de afgelopen vijfentwintig jaar in de iteratieve methoden, waarbij Nederlandse onderzoekers een vooraanstaande plaats innemen. In dit artikel worden iteratieve methoden beschreven, die de afgelopen kwart eeuw ontwikkeld zijn voor het oplossen van lineaire stelsels van de vorm Ax = b,

R n×n

waarbij A ∈ een niet-singuliere matrix is en de oplossing x en de rechterlidvector b beide element zijn van R n . In dit overzicht veronderstellen we dat de gebruikte matrices veel nullen bevatten: A is een ijle matrix. Veelal zijn de matrices afkomstig van gediscretiseerde partiële differentiaalvergelijkingen. De matrices zijn groot tot zeer groot: stelsels met miljarden vergelijkingen en miljarden onbekenden zijn geen uitzondering. Dit betekent vooral dat methoden die optimaal gebruik maken van het computergeheugen en efficiënt zijn in rekentijd, erg in de aandacht staan. Als een methode gebruikt wordt door niet-experts dan is het van groot belang dat de methode robuust is, dat wil zeggen dat de methode altijd met een betrouwbaar antwoord komt, weinig of geen parameters bevat en tenminste sneller is dan de methoden die tot nu toe gebruikt worden. Een aardige uitzondering op het gebruik van grote stelsels is de toepassing van het oplossen van stelsels binnen een computer-

spel. Hierbij komen stelsels van 80 vergelijkingen en 80 onbekenden voor. Echter omdat de oplossing binnen een milliseconde gevonden moet worden, maken de spelletjes toch gebruik van geavanceerde parallelle iteratieve methoden [30]. We starten met een kort overzicht van de ‘state of the art’ methoden in 1980. Daarna worden de nieuwe ontwikkelingen besproken. Eerst worden Krylov-deelruimtemethoden besproken. Deze methoden kunnen alleen zinvol gebruikt worden als ze gecombineerd worden met geschikte preconditioneringen. Een andere belangrijke ontwikkeling bij het oplossen van partiële differentiaalvergelijkingen zijn domein decompositie-methoden. Deze methoden leiden automatisch tot iteratieve methoden. Multigridmethoden kunnen zowel als oplosmethode en als preconditionering gebruikt worden. Als laatste maken we een uitstapje naar speciale computers. In de conclusies geven we een korte samenvatting en een vooruitblik van de toekomst van iteratieve methoden. Als illustratie beschouwen we het 2 × 2 stelsel: 2x + y = 50, x + 2y = 55.

De exacte oplossing wordt gegeven door x = 15 en y = 20. Als we dit stelsel oplossen met de Gauss Jacobi iteratieve methode dan gaat dit als volgt. Kies een startbenadering, bijvoorbeeld x0 = 0 en y0 = 0. Om de volgende benadering (x1 , y1 ) te vinden, lossen we x1 en y1 op uit: 2x1 + y0 = 50, x0 + 2y1 = 55.

De oplossing hiervan wordt gegeven door x1 = 25 en y1 = 27.5. Zo doorgaande vinden we: x2 = 11.25 en y2 = 15, x3 = 17.5 en y3 = 21.875, x4 = 14.0625 en y4 = 18.75, . . . en we zien dat er snelle convergentie optreedt naar de exacte oplossing x = 15 en y = 20. Nederlandse onderzoekers nemen een vooraanstaande plaats in bij de ontwikkeling van snelle en robuuste iteratieve oplosmethoden. Omdat de hoeveelheid literatuur over iteratieve oplosmethoden erg groot is, beperken we ons tot de belangrijkste verwijzingen. In elke paragraaf zal verwezen worden naar boeken [2, 25, 39, 63] of overzichtsartikelen [43]. waarin verdere verwijzingen te vinden zijn. Iteratieve methoden tot aan 1980 Boeken waarin een overzicht gegeven wordt van de literatuur uit deze periode zijn [67, 72, 75]. Opvallend hierbij is dat er in de jaren 1970-1980 weinig of geen boeken verschenen zijn over iteratieve methoden. Een mogelijke oorzaak zou kunnen zijn dat de ontwikkeling van standaard iteratieve methoden min

1

2


NAW 5/10 nr. 2 juni 2009

117

Illustratie: Ryu Tajiri

Kees Vuik

2

3

118

NAW 5/10 nr. 2 juni 2009

Krylov-deelruimtemethoden De meest eenvoudige methode om het stelsel Ax = b op te lossen is de Richardson methode. Kies een startbenadering x0 . Om de nieuwe iterand te bepalen wordt de volgende iteratie gebruikt: xj+1 = xj + (b − Axj ) = xj + rj .

Het is eenvoudig om in te zien dat xj een lineaire combinatie is van x0 , r0 , Ar0 , ...Aj−1 r0 .

Dit blijkt voor veel iteratieve methoden het geval te zijn. De ruimte opgespannen door r0 , Ar0 , ...Aj−1 r0 wordt de Krylovdeelruimte genoemd. Een Krylov-deelruimtemethode is een methode waarbij de iterand xj element is van de Krylovdeelruimte en waarbij of een minimalisatie eigenschap geldt (CG, GMRES), of het residu (bi)-orthogonaal staat op een andere deelruimte met dimensie j .

of meer voltooid was, terwijl de stormachtige ontwikkeling van Krylovmethoden net startte in deze periode. Onder standaard iteratieve methoden verstaan we methoden die gebaseerd zijn op een splitsing van de matrix, zoals Gauss-Jacobi, Gauss-Seidel en Successieve Overrelaxatie (SOR). Veel verschillende varianten (blok versies, Symmetrische SOR, SOR met lokale overrelaxatie enz.) zijn bekend aan het einde van de jaren zeventig. Verder is een goede analyse van de convergentie gemaakt en is er theorie ontwikkeld om schattingen te maken van de optimale (lokale) overrelaxatie parameter [7]. In 1968 is de SIP (Strongly Implicit Procedure) beschreven door Stone [52]. De methode lijkt op een incomplete Choleski ontbinding, met dit verschil: de ongewenste fill wordt niet weggelaten maar wordt zodanig opgenomen in de andere niet-nulelementen van de Choleski factor, dat de fout in de incomplete ontbinding dezelfde orde van grootte heeft als de afbreekfout. Lange tijd is dit één van de snelste oplosmethoden. We zien bijvoorbeeld dat de ICCG methode pas doorbreekt als blijkt dat deze methode sneller is dan de SIP methode. Als laatste noemen we de TDMA (tridiagonal matrix algorithm) methode. Deze is door Patankar [36] gebruikt om gediscretiseerde Navier-Stokes vergelijkingen op te lossen. De methode is nog steeds heel populair bij deze toepassingen. De Chebychev semi-iteratieve methode [23–24] is eerst voor symmetrisch positief


definiete (SPD) matrices ontwikkeld en geanalyseerd en later pas voor algemene matrices. Een nadeel van deze methode is dat er informatie over de eigenwaarden bekend moet zijn. Een voordeel van deze methode is dat er geen inproducten bepaald behoeven te worden, wat de methode heel geschikt maakt voor massaal parallelle computers. Verder wordt deze methode gebruikt om de convergentie(snelheid) van Krylovdeelruimtemethoden theoretisch te schatten. In 1952 is voor het eerst de Conjugate Gradient methode beschreven voor SPD matrices [29]. De methode had een slechte start als directe oplosmethode. Pas toen in 1977 de methode gecombineerd werd door Meijerink en Van der Vorst met een preconditionering gebaseerd op een incomplete Choleski ontbinding ontstond de zeer succesvolle en robuuste ICCG methode [32]. Uit de 872 verwijzingen binnen Google Scholar blijkt dat deze in Nederland ontwikkelde methode zijn weg gevonden heeft bij allerlei toepassingen. Het generaliseren naar algemene matrices is in deze periode nog niet erg succesvol. Een andere belangrijke oplosmethode, die ontwikkeld is vanaf 1960 is de geometrische multigrid methode. Deze is in eerste instantie toepasbaar voor het snel oplossen van elliptische partiële differentiaalvergelijkingen op gestructureerde roosters. De methode wordt in Nederland voor het eerst toegepast op de gediscretiseerde Navier Stokes vergelijkingen. Theoretische analyse laat zien, dat de methode een optimale complexiteit heeft. Bij een vergelijking van de ICCG en de geometrische multigrid methode blijkt dat de laatste voor Poisson-achtige problemen altijd sneller is, als het aantal onbekenden heel groot is. Het toepassingsgebied van de geometrische multigrid methode is in deze periode nogal beperkt. De ICCG methode heeft een veel hoger ‘black box’ gehalte. Domein decompositie methoden zijn al vanaf 1870 bekend. Onder de naam Schwarz wordt een hele grote verscheidenheid van methoden geplaatst. Deze methoden kunnen gebruikt worden om met behulp van analytische technieken partiële differentiaalvergelijkingen op te lossen in gecompliceerde gebieden. Bij een andere toepassing worden er op verschillende gebieden verschillende modellen opgelost. In zekere zin kunnen de blok standaard iteratieve methoden ook gezien worden als domein decompositie methoden. De domein decompositie methoden hebben altijd een iteratieve component in zich. Soms worden ze gebruikt als preconditionering voor een Krylovmethode. Na de opkomst van de

Kees Vuik

Figuur 1 A.N. Krylov (1863–1945)

parallelle computer na 1980 staan deze methoden erg in de aandacht. We komen hier later op terug. In de periode van 1940 tot 1980 zien we dat de digitale computer een enorme ontwikkeling heeft doorgemaakt. Zijn in de beginjaren stelsels met n = 80 al zeer groot, in de jaren tachtig is het al mogelijk om vectoren en (ijle) matrices op te slaan voor n = 106 . Het blijkt dat er een interactie is tussen de groei van de computersnelheid en de ontwikkeling van de iteratieve methode. Een voorbeeld hiervan is dat iteratieve oplosmethoden sneller worden dan directe oplosmethoden als n groter genomen kan worden. Alhoewel er omstreeks 1980 speciale computerarchitecturen op de markt komen (vector, transputer, RISC, etc.) vinden de grote ontwikkelingen op dit gebied pas na 1980 plaats. Krylov-deelruimtemethoden Zoals uit het voorgaande blijkt waren de Krylovmethoden voor symmetrische positief definiete matrices al voor 1980 uitgekristalliseerd. De situatie is geheel anders voor nietsymmetrische matrices. Rondom 1980 wordt zelfs bewezen [18, 69] dat het onmogelijk is om een Krylovmethode te vinden, die alle prettige eigenschappen van de Conjugate Gradient-methode bezit en die gebruikt kan worden voor algemene matrices. Eén van de weinige ontwikkelingen die wel plaatsvindt in de tachtiger jaren is een goed begrip van de superlineaire convergentie van de Conjugate Gradient methode [59–60]. Voor de Conjugate Gradient methode gelden de volgende eigenschappen: • de benadering xk is een element van de Krylov-deelruimte gebaseerd op A en het beginresidu r0 = b − Ax0 ,

3

4

Kees Vuik


• bij de Conjugate Gradient methode wordt

de fout geminimaliseerd in de A-norm, • de Conjugate Gradient methode maakt gebruik van een korte recurrenties (2 of 3 terms recurrentie, afhankelijk van de implementatie). In [18, 69] wordt aangegeven voor welke klasse van matrices Krylovmethoden ontworpen kunnen worden met bovenstaande eigenschappen. Voor algemene matrices zal er dus één van de eigenschappen niet gelden. Dit leidt tot de volgende drie klassen van methoden: • De Conjugate Gradient methode toegepast op de normaalvergelijkingen. Hierbij is niet voldaan aan de eerste eigenschap. • BiCG methoden. Deze voldoen niet aan een optimalisatie eigenschap. • GMRES methoden. Bij deze methoden worden lange recurrenties gebruikt. Conjugate Gradient toegepast op de normaalvergelijkingen Dit is een voor de hand liggend idee: pas Conjugate Gradient toe op AT Ax = AT b.

Als de matrix A niet singulier is, dan is de matrix AT A SPD. De beste methode binnen deze klasse is de LSQR methode [35]. Een sterk punt bij deze methode is dat er goede foutschattingen gemaakt worden; deze maken de methode robuust en zeer betrouwbaar. Een nadeel is dat de convergentie bepaald wordt door κ(AT A) = κ(A2 ). In het algemeen leidt dit tot trage convergentie. Ook is het vinden van een geschikte preconditionering geen eenvoudige zaak. BiCG methoden Bij de BiCG methode, die beschreven is in [19], wordt de CG methode toegepast op ˜ ˜ = b. Ax = b en AT x ˜ − AT x ˜k het schaduwHierbij wordt r˜k = b residu genoemd. Het blijkt dat in deze methode de zoekrichtingen en de residuen biorthogonaal zijn, dat wil zeggen: riT r˜j = 0 als i 6= j.

De methode convergeert snel als A een kleine afwijking heeft van een SPD matrix. Nadelen zijn dat er per iteratie twee matrix-vector producten nodig zijn, verder is AT nodig, die niet altijd beschikbaar is, er kan ‘serious break down’ optreden en er is geen convergentie analyse mogelijk. Het aandeel nadelen wordt minder bij de CGS (CG Squared) methode [50]. Bewerkingen met AT zijn niet langer nodig. Er

zijn nog wel twee matrix-vector producten nodig per iteratie, echter de convergentie gaat ook twee keer zo snel als bij de BiCG methode. De methode vulde een gat in de markt en wordt heel vaak gebruikt en geciteerd (651 verwijzingen binnen Google Scholar). Een nadeel van CGS is dat lokaal het residu enorm kan toenemen. Een oplossing hiervan is de Bi-CGSTAB methode ontwikkeld door Henk van der Vorst [62]. Deze methode is nog steeds één van de favoriete Krylov oplosmethoden. Naast deze methode zijn nog veel andere varianten ontwikkeld. De meeste daarvan worden niet zo veel gebruikt. Een uitzondering hierop is de QMR (Quasi Minimal Residual) methode [21]. Deze methode lijkt op de Minimal Residual methode, echter, het residu wordt niet geminimaliseerd in een echte norm. Voor een bepaalde klasse van matrices is dit een hele geschikte methode. Het idee van Sonneveld om van BiCG over te stappen op CGS is later breed toegepast. Voor bijna elke methode werden er na verloop van tijd Transpose Free varianten ontwikkeld zoals TFQMR en Transpose Free Lanczos [10, 20]. De ontwikkelingen komen vanaf 1995 bijna tot stilstand bij het ontwerpen van nieuwe Krylovmethoden. Onlangs is een nieuwe ontwikkeling gestart binnen de wereld van de Krylovmethoden. Sonneveld en Van Gijzen hebben de 30 jaar oude IDR (een voorloper van CGS) nieuw leven ingeblazen. Een nieuwe variant IDR(s ) met meerdere schaduw residuen blijkt verrassend goed te convergeren [51]. Zij hebben bewezen, dat IDR(1) gelijk is aan Bi-CGSTAB zonder rekening te houden met het effect van afrondfouten. Voor waarden van s groter dan 1 ontstaan er methoden (IDR(4) is een goede keuze), die korte recurrenties gebruiken, terwijl de convergentie veel dichter bij full GMRES ligt dan voor de Bi-CGSTAB methode. De IDR(s ) methode is ook vergeleken met BICGSTAB2 [28], Bi-CGSTAB(l) [47] en ML(k)BiCGSTAB [74]. Voor de details verwijzen we naar [51]. GMRES methoden Als laatste bespreken we de klasse van GMRES methoden. Dit zijn methoden, die een benadering vinden binnen de Krylovdeelruimte en waarvan het residu geminimaliseerd wordt. Als gevolg van de stelling kan dit alleen maar door lange recurrenties te gebruiken. De GMRES en FOM methoden zijn ontwikkeld door Saad en Schultz [40] in 1986. Als het aantal iteraties klein is dan is GMRES de beste Krylovmethode. Als er veel iteraties nodig zijn, dan moet GMRES herstart worden.

NAW 5/10 nr. 2 juni 2009

119

Figuur 2 De niet-nul structuur van een ijle matrix

De methoden GCR [16] en Orthomin [68] leveren dan betere prestaties, omdat deze ook getrunkeerd kunnen worden en zelfs gebruik kunnen maken van variabele preconditioneringen. Door de optimalisatie eigenschap kunnen voor deze methoden ook convergentie resultaten bewezen worden. Eén daarvan is eenvoudig te gebruiken als het spectrum omsloten wordt door een cirkel in het complexe vlak, die de oorsprong niet bevat [40]. In deze stelling wordt gebruik gemaakt van het conditiegetal van de matrix van eigenvectoren. Bovendien is er voor deze methoden bewezen, dat er superlineaire convergentie optreedt [64]. Omdat de klasse van niet-symmetrische matrices ook niet normale matrices bevat kan het conditiegetal van de matrix van eigenvectoren heel groot (of ∞) zijn. Een gevolg hiervan is dat dan het spectrum geen bruikbare informatie oplevert voor de convergentie [26]. Een wat ongebruikelijke methode is de ENmethode [15]. De methode is niet erg bekend maar berust wel op een origineel idee, namelijk het oplossen van een lineair stelsel met een Newton Raphson achtige methode. Deze methode is de motivatie geweest voor de ontwikkeling van de GMRESR methode [65]. Afsluitende opmerkingen We zien dat elke klasse van methoden bepaalde voor- en nadelen heeft. Dit heeft ertoe geleid dat er ‘hybride’ methoden ontwikkeld zijn om de voordelen te behouden en de nadelen te verminderen. Voorbeelden zijn: GMRESR [65], een combinatie van een GCR buitenloop en GMRES binnenloop, Bi-CGSTAB(l) [47], waarbij het BiCG polynoom gecombineerd wordt met een GMRES(l) polynoom en Krylovmethoden gecombineerd met een polynoom preconditionering. Na verloop van tijd bleek het aantrekke-

4

5

120

NAW 5/10 nr. 2 juni 2009


Kees Vuik

lijk te zijn om Krylovmethoden toe te passen, die gebruik kunnen maken van variabele preconditioneringen en/of benaderende matrix vector vermenigvuldigingen. Dit heeft geleid tot flexibele Krylov varianten. GCR is standaard geschikt voor een variabele preconditionering [65]. In 1992 is de F(lexible)GMRES methode [38] ontwikkeld. Hierbij is een kleine aanpassing van de GMRES methode nodig. Onlangs is de F(lexible)CG beschreven in [33]. Voor recente resultaten over benaderende matrix vector vermenigvuldigingen verwijzen we naar [45–46, 58]. Het kiezen van een geschikte Krylovmethode is niet altijd eenvoudig. Belangrijk is het aantal benodigde/verwachte iteraties mg en de verhouding f tussen de rekentijd van een matrix vector product en een vector update. Als mg klein is en f groot dan is GMRES de beste methode. Als mg groot is en f klein dan is Bi-CGSTAB de beste methode terwijl voor tussenliggende waarden van mg en f GMRESR de beste methode is. Figuur 3 Het convergentiegedrag van IDR(s) vergeleken met andere Krylovmethoden

Preconditioneringen Voor veel stelsels is de convergentie van Krylovmethodes traag. Dit kan verbeterd worden door een geschikte preconditionering toe te passen [3]. Hierbij wordt het linker- en rechterlid van het stelsel vermenigvuldigd met een matrix zodat de oplossing niet verandert maar de eigenschappen van de gepreconditioneerde matrix veel gunstiger zijn voor de convergentie van een Krylovmethode. We hebben in het overzicht tot 1980 gezien dat voor SPD matrices een incomplete Choleski methode tot een goede verbetering leidt. Deze methode is gegeneraliseerd voor algemene matrices en staat bekend als een ILU preconditionering. Bij deze preconditioneringen wordt een benaderende ontbinding gemaakt van A. Het product van een onderdriehoeksmatrix L en ˆ een bovendriehoeksmatrix U is gelijk aan A ˆ en de waarde van kA − Ak is klein. Door de ‘fill in’ te verwaarlozen kan de geheugenruimte beperkt blijven. In veel gevallen zijn deze preconditioneringen snel en robuust. We geven nu een aantal observaties over deze klasse van preconditioneringen. • Het blijkt dat de ordening van de onbekenden van groot belang is voor de kwaliteit van de ILU preconditioneringen. Methoden om de bandbreedte te beperken (zoals het reversed Cuthill-McKee algoritme (RCM) [11]) of het profiel te minimaliseren (zoals het Sloan algoritme [48] en het Approximate Minimum Degree algoritme (AMD) [1]) blijken ook goed te werken voor ILU preconditioneringen.

• Voor Poisson problemen blijkt dat het aan-

tal iteraties toeneemt als de grootte van het rekenrooster toeneemt. Dit is een nadeel van ILU preconditioneringen. • Voor sommige herordeningen (MRILU [8] en ARMS [41]) blijft het aantal iteraties min of meer constant ook als de grootte van het rekenrooster toeneemt. • Een nadeel van ILU preconditioneringen is dat ze in het algemeen lastig te gebruiken zijn op parallelle computers. Voor sommige speciale matrices zijn parallelle varianten beschikbaar [13–14, 39, 71]. • Het toelaten van extra ‘fill in’ leidt vaak tot minder iteraties en soms tot minder rekentijd. De laatste opmerking zullen we wat verder uitwerken. Er zijn twee belangrijke methoden om ‘fill in’ toe te laten: ‘level fill in’ ILU(k) [32] gebaseerd op de niet-nul structuur van de matrix en ‘threshold fill in’ waarbij elementen onder een bepaalde drempel weggegooid worden zoals ILUT(). Meestal blijkt de laatste methode erg gevoelig te zijn voor de keuze van . Als te klein is, dan zijn de matrices L en U bijna vol en als te groot is, dan treedt er trage convergentie op. Op dit moment is ILUpack [6] één van de beste preconditionerings pakketten gebaseerd op ‘threshold fill in’. Een manier om ILU te generaliseren is door gebruik te maken van blok varianten. Hiervoor zijn verschillende toepassingen mogelijk. De eerste is dat er per roosterpunt verschillende onbekenden zijn zoals bij stromingsproble-

men: drie snelheden, de druk en de temperatuur. Het kan aantrekkelijk zijn om de matrix die behoort bij deze onbekenden als één 5×5 blok te zien. Voordeel is dat de convergentie vaak sterk verbeterd. Een nadeel is dat zo’n ontbinding meer geheugenruimte kost [2, 31, 66]. Een andere blokversie is een opdeling van de matrix in deelmatrices. De buitendiagonaalblokken worden weggelaten bij het bepalen van de blok ILU ontbinding. Een voordeel is dat deze BILU preconditionering goed parallelliseerbaar is. Een nadeel is dat de convergentie slechter wordt als het aantal blokken toeneemt [70]. Het is ook mogelijk om als preconditionering gebruik te maken van één of meerdere iteraties van een standaard iteratieve methode. Kandidaten hiervoor zijn: SOR, SSOR en multigrid methoden. De laatste preconditionering is heel geschikt voor Poisson-achtige problemen. Het blijkt dat het aantal iteraties voor zulke problemen onafhankelijk is van de grootte van het rekenrooster. Verder is deze combinatie erg robuust. Veranderingen in het rekenrooster, rechterlid, coëfficiënten en of randvoorwaarden beïnvloeden de convergentie nauwelijks [56]. Een andere klasse van preconditioneringen is gebaseerd op ijle benaderingen van de inverse, de zogenaamde SPAI (Sparse Approximate Inverse) methoden. Hierbij wordt een matrix B geconstrueerd zodanig dat B een gekozen niet-nul patroon heeft. De resterende

5

6

Kees Vuik


NAW 5/10 nr. 2 juni 2009

121

de convergentie niet meer. De combinatie van traditionele preconditioneringen met een ‘second level’ preconditionering leidt vaak tot robuuste en snelle oplosmethoden [9, 42]. Een recent overzicht is te vinden in [54]. Andere ontwikkelingen We sluiten af met recente ontwikkelingen op het gebied van domein decompositie, multigrid methoden en speciale computers.

Figuur 4 Roosters gebruikt bij een V -cycle van multigrid

elementen van B worden dan bepaald door BA − I

in één of andere norm te minimaliseren. Een groot voordeel is dat een SPAI preconditionering eenvoudig te parallelliseren is. Dit verklaart ook de grote populariteit van deze methode. Een groot nadeel is dat de kwaliteit van de preconditionering erg afhangt van de gekozen norm en het gekozen niet-nul patroon [27]. Momenteel is het efficiënt oplossen van zadelpunt problemen een belangrijk onderwerp. Dit type problemen treedt op bij het oplossen van optimaliseringsproblemen en het oplossen van de gediscretiseerde incompressibele Navier Stokes vergelijkingen. Hiervoor zijn veel verschillende preconditioneringen voorgesteld. Het blijkt dat de kwaliteit van deze preconditioneringen erg afhankelijk is van het op te lossen probleem. Op dit moment is er geen preconditionering, die het goed doet voor alle zadelpunt problemen. Verder is er nog geen preconditionering bekend van de gediscretiseerde incompressibele Navier Stokes vergelijkingen, waarbij de gebruikte discretisatie cellen een hoge aspect ratio hebben [4, 12, 57]. Als laatste voorbeeld van een preconditionering noemen we ‘operator based’ preconditioneringen. Het idee is als volgt: het op te lossen stelsel is afkomstig van een gediscretiseerde operator en is moeilijk oplosbaar met een Krylovmethode. We kunnen nu een operator bepalen die veel lijkt op de oorspronkelijke operator, maar na discretisatie is het resulterende stelsel wel efficiënt op te lossen. De-

ze gediscretiseerde operator kan dan gebruikt worden als preconditionering. Als voorbeeld beschouwen we de Helmholtz vergelijking −∆p − k2 p = 0,

met geschikte randvoorwaarden. Na discretisatie levert dit een stelsel Ax = b op. Dit stelsel is moeilijk oplosbaar voor grote waarden van k omdat dan de eigenwaarden van A positieve en negatieve waarden hebben. Een goede preconditionering is de discretisatie van −∆p + k2 p = 0.

De resulterende matrix geven we aan met B . We kunnen dan een Krylovmethode toepassen op B −1 Ax = B −1 b.

De vermenigvuldiging B −1 r = v wordt uitgevoerd door het stelsel Bv = r benaderend op te lossen. Andere mogelijkheden zijn: de preconditionering baseren op ‘eenvoudige randvoorwaarden’, op constante coëfficiënten, op een lagere orde discretisatie van de operator of op het (anti-)symmetrische deel van de operator. [17]. Voor moeilijke problemen blijkt dat een preconditionering wel helpt, maar dat er vaak één of meerdere eigenwaarden de convergentie enorm vertragen. Om dit te verbeteren zijn de zogenaamde ‘second level’ preconditioneringen ontwikkeld. Nadat er een schatting gemaakt is van de ‘slechte’ eigenvectoren, wordt er een projectie operator gemaakt die er voor zorgt dat deze eigenvectorcomponenten niet meer voorkomen in het residu. Als gevolg hiervan beïnvloeden deze eigenwaarden

Domein decompositie De domein decompositie methoden zijn enorm ontwikkeld in de afgelopen 25 jaar. Het toepassingsgebied is uitgebreid van eenvoudige Poisson-achtige problemen naar problemen op het gebied van stromingsleer, elektromagnetische problemen en mechanica. Generalisering naar niet-symmetrische matrices heeft plaatsgevonden. Verder is de combinatie van een domein decompositie methode en een parallelle computer heel goed. Net als bij blok preconditioneringen zien we vaak dat de convergentie slechter wordt als het aantal domeinen toeneemt en dat er problemen zijn bij discontinue coëfficiënten. Beide problemen kunnen verholpen worden door een coarse grid acceleratie toe te voegen. Hierbij is een grote overeenkomst met de ‘second level’ preconditioneringen beschreven in de vorige paragraaf. Er is veel onderzoek gedaan naar optimale koppelingsvoorwaarden [22, 53]. Domein decompositie methoden worden vaak gecombineerd met Krylovmethoden. Als de subdomein problemen benaderend opgelost worden en een Krylovmethode wordt gebruikt, dan blijkt vaak dat eenvoudige koppelingsvoorwaarden voldoende zijn voor snelle convergentie. Voor verdere details verwijzen we naar [49, 55]. Multigrid methoden Ook bij de multigrid methoden zien we een ontwikkeling naar algemene toepassingsgebieden. Multigrid methoden worden breed toegepast bij problemen in stromingsleer die gediscretiseerd zijn op gestructureerde rekenroosters. De belangrijkste ontwikkeling is het verschijnen van Algebraïsche Multigrid Methoden [37]. Het voordeel hiervan is dat deze methoden een hoog black box gehalte hebben. Een nadeel is dat er weinig theorie beschikbaar is. Bovendien zijn voor vergelijkbare problemen Algebraïsche Multigrid Methoden een orde duurder in rekentijd dan de Geometrische Multigrid Methoden. Multigrid methoden sluiten goed aan bij de toenemende interesse voor het oplossen van ‘multiscale’ problemen. Voor verdere details verwijzen we

6

7

122

NAW 5/10 nr. 2 juni 2009

naar [56, 73]. Speciale computers Ontwikkelingen binnen de numerieke lineaire algebra gaan vaak hand in hand met ontwikkelingen op het gebied van de beschikbare computers. Een eerste gevolg hiervan is dat de problemen tegenwoordig veel groter kunnen zijn dan de problemen 25 jaar geleden. Dit heeft belangrijke gevolgen voor de methoden die het best gebruikt kunnen worden. In de jaren 80 zien we de opkomst van vectorcomputers, parallelle computers en risc processoren. Al deze computers hebben vaak speciale eigenschappen die gebruikt moeten worden om ze efficiënt in te zetten bij numerieke lineaire algebra methoden. Eén van de belangrijkste punten is dat deze computers prima werken op regelmatige vectoren en heel traag werken op vectoren met indirecte adressering. We hebben al gezien dat veel preconditioneringen lastig parallel te maken zijn [14, 61]. Meestal zien we dat er eerst een stap terug gedaan wordt naar een eenvoudige preconditionering en pas later worden de complexe preconditioneringen ook geschikt gemaakt voor de nieuwe typen computers.


In de jaren 90 komen de Beowulf clusters op. Dit zijn clusters van PC’s gekoppeld via ethernet die gebruikt worden als een poor man’s parallelle computer. Bij de massief parallelle computers is ethernet vervangen door zeer snelle verbindingen. De opkomst van PVM, MPI, BSP [5] en OpenMP heeft het gebruik van deze computers erg vereenvoudigd. In de 21-eeuw zien we de opkomst van PC’s met twee tot zestien centrale processor eenheden binnen één computer. Dit geeft opnieuw mogelijkheden om de methoden te parallelliseren. Op dit moment worden er methoden ontwikkeld om te gebruiken op de zeer snelle videokaarten. Deze ontwikkeling staat nog in de kinderschoenen, maar heeft de potentie om uit te groeien tot een nieuwe standaard voor high performance computing. Voor verdere details verwijzen we naar [13, 34, 44]. Conclusies We zien dat er een grote ontwikkeling geweest is op het gebied van iteratieve methoden. De toepassingsgebieden zijn enorm toegenomen. Er blijven echter nog steeds gebieden waar geschikte iteratieve methoden ontbreken. De ontwikkeling van algemene Krylov-

Kees Vuik

methoden lijkt tot stilstand gekomen te zijn, uitgezonderd de recente ontwikkeling van de IDR(s) methode. De verwachting is dat er voor ingewikkelde problemen steeds meer gebruik gemaakt gaat worden van combinaties van verschillende numerieke lineaire algebra methoden. Nog steeds staat het ontwikkelen van geschikte preconditioneringen enorm in de aandacht. Ook het beschikbaar komen van nieuwe computer-architecturen zal leiden tot nieuwe methoden. De grenzen van het rekenen met ‘double precision’-getallen komen in zicht. De reden hiervoor is dat bij de discretisatie de gebruikte stapgrootte steeds kleiner gekozen wordt. Dit betekent dat de problemen steeds groter worden. Het resultaat hiervan is dat nκ(A)µ (n is de dimensie van de oplosvector, κ(A) is het conditiegetal en µ de machineprecisie, die gelijk is aan 10−16 voor ‘double precision’-getallen) in de buurt van 1 komt te liggen. Een gevolg hiervan is dat de nauwkeurigheid van de berekende oplossing niet meer gegarandeerd kan worden. De verwachting is dat er overgestapt moet worden naar machines met een relatieve machine precisie van 10−32 . k

Referenties 1

Patrick R. Amestoy, Timothy A. Davis, and Iain S. Duff. An approximate minimum degree ordering algorithm. SIAM J. Matrix Anal. Appl., 17(4):886–905, 1996.

2

O. Axelsson. Iterative solution methods. Cambridge University Press, Cambridge, 1994.

3

M. Benzi. Preconditioning techniques for large linear systems: A survey. Journal of Computational Physics, 182:418–477, 2002.

4

Michele Benzi, Gene H. Golub, and Jörg Liesen. Numerical solution of saddle point problems. Acta Numer., 14:1–137, 2005.

5

R. Bisseling. Parallel Scientific Computation: A Structured Approach Using BSP and MPI. Oxford University Press, Oxford, 2004.

6

M. Bollhöfer and Y. Saad. Multilevel preconditioners constructed from inverse-based ILUs. SIAM J. Sci. Comput., 27:1627–1650, 2006.

7

E.F.F. Botta and A.E.P. Veldman. On localrelaxation methods and their application to

convection diffusion equations. J. Comp. Phys., 48:127–149, 1982. 8

E.F.F. Botta and F.W. Wubs. Matrix renumbering ILU: An effective algebraic multilevel ILU preconditioner for sparse matrices. SIAM. J. Matrix Anal. and Appl., 20:1007–1026, 1999.

9

B. Carpentieri, L. Giraud, and S. Gratton. Additive and multiplicative spectral two-level preconditioners for general linear systems. SIAM Journal on Scientific Computing, 29:1593–1612, 2007.

10 T.F. Chan, L. de Pillis, and H.A. van der Vorst. Transpose-free formulations of Lanczos-type methods for nonsymmetric linear systems. Numerical Algorithms, 17:51–66, 1998. 11

E. Cuthill and J. McKee. Reducing the bandwidth of sparse symmetric matrices. In Proceedings of the 1969 24th national conference, pages 157– 172, New York, NY, USA, 1969. ACM.

12 H. Sue Dollar, Nicholas I. M. Gould, Wil H. A. Schilders, and Andrew J. Wathen. Implicit-factorization preconditioning and iterative solvers for regularized saddle-point systems. SIAM J. Matrix Anal. Appl., 28(1):170–189, 2006. 13

J.J. Dongarra, L.S. Duff, D.C. Sorensen, and H.A. van der Vorst. Solving Linear Systems on Vector and Shared Memory Computers. SIAM, Philadelphia, 1991.

14 I.S. Duff and H.A. van der Vorst. Developments and trends in the parallel solution of linear systems. Parallel Computing, 25:1931–1970, 1999. 15

T. Eirola and O. Nevanlinna. Accelerating with rank-one updates. L.A.A., 121:511–520, 1989.

16 S.C. Eisenstat, H.C. Elman, and M.H. Schultz. Variational iterative methods for nonsymmetric systems of linear equations. SIAM J. Num. Anal., 20:345–357, 1983.

7

8

Kees Vuik

17

Y.A. Erlangga, C.W. Oosterlee, and C. Vuik. A novel multigrid based preconditioner for heterogeneous Helmholtz problems. SIAM J. Sci. Comput., 27:1471–1492, 2006.

18 V. Faber and T. Manteuffel. Necessary and sufficient conditions for the existence of a conjugate gradient method. SIAM J. Numer. Anal., 21:352– 362, 1984. 19 R. Fletcher. Conjugate gradient methods for indefinite systems. Lecture notes in Mathematics, 506:73–89, 1976. Springer-Verlag, Berlin, Heidelberg. 20 R.W. Freund. A transpose-free quasi-minimal residual algorithm for non-Hermitian linear systems. SIAM Journal on Scientific Computing, 14:470–482, 1993. 21 R.W. Freund and N.M. Nachtigal. QMR: a quasiminimal residual method for non-Hermitian linear systems. Numerische Mathematik, 60:315– 339, 1991. 22 M. Genseberger. Domain decomposition in the Jacobi-Davidson method for eigenproblems. PhD thesis, University Utrecht, Utrecht, 2001. 23 G.H. Golub and R.S. Varga. Chebyshev semiiterative methods, successive over-relaxation iterative methods, and second order Richardson iterative methods. I. Numer. Math., 3:147–156, 1961. 24 G.H. Golub and R.S. Varga. Chebyshev semiiterative methods, successive over-relaxation iterative methods, and second order Richardson iterative methods. II. Numer. Math., 3:157–168, 1961. 25 A. Greenbaum. Iterative methods for solving linear systems. Frontiers in applied mathmatics 17. SIAM, Philadelphia, 1997. 26 A. Greenbaum, V. Ptak, and Z. Strakos. Any nonincreasing convergence curve is possible for GMRES. SIAM J. Matrix Anal. Appl., 17:465–469, 1996. 27 M.J. Grote and T. Huckle. Parallel preconditioning with sparse approximate inverses. SIAM J. Sci. Comp., 18:838–853, 1997. 28 Martin H. Gutknecht. Variants of BICGSTAB for matrices with complex spectrum. SIAM Journal on Scientific Computing, 14:1020–1033, 1993. 29 M.R. Hestenes and E. Stiefel. Methods of Conjugate Gradients for Solving Linear Systems. J. Res. Nat. Bur. Stand., 49:409–436, 1952. 30 C. Lacoursière. A parallel block iterative method for interactive contacting rigid multibody simulations on multicore pcs. In F. G. Gustavson, B. Kagstrom, J.Dongarra, F. Wolf, A.C. Elster, D. Kressner, M. Sala, X. Cai, A. Laaksonen, and M. Gerndt, editors, Applied Parallel Computing. State of the Art in Scientific Computing, 8th International Workshop, PARA 2006, Umea, Sweden, June 18-21, pages 956–965, Berlin, 2007. Springer. 31

J.A. Meijerink. Iterative methods for the solution of linear equations based on incomplete block factorization of the matrix. Proc. 1983 Reservoir Simulation Symposium, San Francisco, CA, 1983 (paper SPE 12262):297–303, 1983.

32 J.A. Meijerink and H.A. Van der Vorst. An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. Math. Comp., 31:148–162, 1977. 33 Y. Notay. Flexible conjugate gradients. SIAM Journal on Scientific Computing, 22:1444–1460, 2000. 34 J.M. Ortega. Introduction to parallel and vector solution of linear systems. Plenum Press, New York, 1988. 35 C.C. Paige and M.A. Saunders. LSQR : An algorithm for sparse linear equations and sparse


NAW 5/10 nr. 2 juni 2009

123

least squares. A.C.M. Trans. Math. Softw., 8:43– 71, 1982.

56 U. Trottenberg, C.W. Oosterlee, and A. Schüller. Multigrid. Academic Press, London, 2000.

36 S.V. Patankar. Numerical heat transfer and fluid flow. McGraw-Hill, New York, 1980.

57 M. ur Rehman, C. Vuik, and G. Segal. A comparison of preconditioners for incompressible Navier-Stokes solvers. Int. J. Num. Meth. in Fluids, 57:1731–1751, 2008.

37 J.W. Ruge and K. Stüben. Algebraic multigrid. In S.F. McCormick, editor, Multigrid Methods, chapter 4, pages 73–130. SIAM, Philadelphia, Pennsylvania, 1987. 38 Y. Saad. A flexible inner-outer preconditioned GMRES algorithm. SIAM J. Sci. Comput., 14:461– 469, 1993. 39 Y. Saad. Iterative methods for sparse linear systems, Second Edition. SIAM, Philadelphia, 2003. 40 Y. Saad and M.H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput., 7:856–869, 1986. 41 Y. Saad and B. Suchomel. ARMS: an algebraic recursive multilevel solver for general sparse linear systems. Num. Lin. Alg. Appl., 9:359– 378, 2002. 42 Y. Saad, M. Yeung, J. Ehrel, and F. Guyomarc’h. A deflated version of the conjugate gradient algorithm. SIAM Journal on Scientific Computing, 21:1909–1926, 2000. 43 Yousef Saad and Henk A. Van Der Vorst. Iterative solution of linear systems in the 20th century. Journal of Computational and Applied Mathematics, 123:1–33, 2000. 44 R.W. Shonkwiler and L. Lefton. An introduction to parallel and vector scientific computation. Cambridge University Press, Cambridge, 2006. 45 V. Simoncini and D.B. Szyld. Theory of inexact Krylov subspace methods and applications to scientific computing. SIAM Journal on Scientific Computing, 25:454–477, 2003. 46 V. Simoncini and D.B. Szyld. On the occurrence of superlinear convergence of exact and inexact Krylov subspace methods. SIAM Review, 47:247–272, 2005. 47 G.L.G. Sleijpen and D.R. Fokkema. BiCGstab(l) for linear equations involving unsymmetric matrices with complex spectrum. ETNA, 1:11–32, 1993. 48 S.W. Sloan. An algorithm for profile and wavefront reduction of sparse matrices. International Journal for Numerical Methods in Engineering, 23:239–251, 1986. 49 B. Smith, P. Bjorstad, and W. Gropp. Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations. Cambridge University Press, Cambridge, 1996. 50 P. Sonneveld. CGS: a fast Lanczos-type solver for nonsymmetric linear systems. SIAM J. Sci. Stat. Comp., 10:36–52, 1989. 51

Peter Sonneveld and Martin B. van Gijzen. IDR(s ): A family of simple and fast algorithms for solving large nonsymmetric systems of linear equations. SIAM Journal on Scientific Computing, 31:1035–1062, 2008.

52 H.L. Stone. Iterative solution of implicit approximations of multidimensional partial differential equations. SIAM J. Numer. Anal., 5:530–558, 1968. 53 K.H. Tan. Local Coupling in Domain Decomposition. PhD thesis, University Utrecht, Utrecht, 1995.

58 J. van den Eshof and G.L.G. Sleijpen. Inexact Krylov subspace methods for linear systems. SIAM Journal on Matrix Analysis and Applications, 26:125–153, 2004. 59 A. van der Sluis and H.A. van der Vorst. The rate of convergence of conjugate gradients. Numer. Math., 48:543–560, 1986. 60 A. van der Sluis and H.A. van der Vorst. The convergence behaviour of Ritz values in the presence of close eigenvalues. Lin. Alg. Appl., 88/89:651–694, 1987. 61 H.A. van der Vorst. High performance preconditioning. SIAM J. Sci. Stat. Comput., 10:1174– 1185, 1989. 62 H.A. van der Vorst. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for solution of non-symmetric linear systems. SIAM J. Sci. Stat. Comp., 13:631–644, 1992. 63 H.A. van der Vorst. Iterative Krylov Methods for Large Linear Systems. Cambridge Monographs on Applied and Computational Mathematics, 13. Cambridge University Press, Cambridge, 2003. 64 H.A. van der Vorst and C. Vuik. The superlinear convergence behaviour of GMRES. J. Comp. Appl. Math., 48:327–341, 1993. 65 H.A. van der Vorst and C. Vuik. GMRESR: a family of nested GMRES methods. Num. Lin. Alg. Appl., 1:369–386, 1994. 66 S. van Veldhuizen, C. Vuik, and C.R. Kleijn. A class of projected Newton methods to solve laminar reacting flow problems. Report 08-03, Delft University of Technology, Delft Institute of Applied Mathematics, Delft, 2008. 67 R.S. Varga. Matrix iterative analysis. PrenticeHall, London, 1962. 68 P.K.W. Vinsome. ORTHOMIN: an iterative method solving sparse sets of simultaneous linear equations, Proceedings of the Fourth Symposium on Reservoir Simulation. Society of Petroleum Engineers of AIME, Richardson, 1976. pp. 149–159. 69 V. V. Voevodin. The problem of nonselfadjoint extension of the conjugate gradient method is closed. U.S.S.R. Comput. Math. and Math. Phys., 23:143–144, 1983. 70 C. Vuik and J. Frank. Coarse grid acceleration of a parallel block preconditioner. Future Generation Computer Systems, 17:933–940, 2001. 71

C. Vuik, R.R.P. van Nooyen, and P. Wesseling. Parallelism in ILU-preconditioned GMRES. Paral. Comp., 24:1927–1946, 1998.

72 E.L. Wachspress. Iterative solution of elliptic systems. Prentice-Hall, Englewood Cliffs, 1966. 73 P. Wesseling. An introduction to multigrid methods. Wiley, Chichester, 1992. 74 Man-Chung Yeung and Tony F. Chan. ML(k)BiCGSTAB: A BiCGSTAB variant based on multiple Lanczos starting vectors. SIAM Journal on Scientific Computing, 21:1263–1290, 1999. 75 D.M Young. Iterative solution of large linear systems. Academic Press, New York, 1971.

54 J.M. Tang. Two-level preconditioned conjugate gradient methods with applications to bubbly flow problems. PhD thesis, Delft University of Technology, Delft, 2008. 55 A. Toselli and W. Widlund. Domain Decomposition Methods. Computational Mathematics, Vol. 34. Springer, Berlin, 2005.

8

Een kwart eeuw iteratieve methoden

Recommend Documents