Iterative methods and preconditionersVictor [email protected], University of TennesseeCS 594, March 2003OverviewTheory of iterative methodsComputational kernelsPreconditionersPro/con iterative or direct methodsAgainst direct Variants of Gaussian elimination: Time and spacedemand can be huge.No use of analytical properties.Pro iterative Use analytical properties of systems.Less time and space. Possibly.Pro direct Predictable in advance. Solution up to round-off.Against iterative Unpredictable. Solution not guaranteed.Iterative methodsStationary methodsIterative refinement:• To solve: Ax = b.• Suppose M ≈ A and solve Mx1= b• Error: e1= x1− x; residual: r1= Ax1− b• Since x = A−1b = x1− A−1r1,– define x2= x1− M−1r1,– and x3= x2− M−1r2,– et ceteraAnalysis of stationary methodsBasic step:r2= Ax2− b = A(x1− M−1r1) − b = (I − AM−1)r1Inductively:rn= (I − AM−1)n−1r1so convergence if all|λ(I − AM−1)| < 1.Examples of stationary iterative methodsJacobi method:M = DA= diag(A)Gauss-Seidel method:M = DA+ LASOR method:M = DA+ ωLAThese methods converge for M matrices:• A positive definite• aij≤ 0 for i 6= j.Steepest descentIterative process for Ax = b:xi+1= xi+ M−1riSlight generalization:xi+1= xi+ αiM−1riline search;with optimal αi: steepest descentDerivation of Krylov space methodsStationary methods:rn= (I − AM−1)rn−1, rn= (I − AM−1)n−1r1Try more general (with M = I):rn= Pn−1(A)r1or: define the Krylov sequence:k1= r1, kn+1= Akn,and set:rn∈ span(k1, . . . , kn) = span(r1, Ar1, A2r1. . . , Anr1)or equivalentlyrn∈ span(r1, . . . , rn−1, Arn).AccelerationIn the definitionrn= Pn−1(A)r1choose the polynomial coefficients.Two choices for polynomial coefficients:• Get Pnas close to zero as possible.– Pn(0) = 1 consistency condition– knowledge of spectrum required– ⇒ Chebyshev polynomials• Make residuals mutually orthogonal under some inner product.OrthogonalityMaking error orthogonal to generated subspace⇔ projecting error on subspace⇔ minimising error on subspacePrerequisite: inner product⇒ matrix needs to be definite⇒ the further from elliptic you get, the worse an iterative method willperform(or change your inner product)Properties of orthogonalising methods• Theoretically: convergence in N steps (N matrix dimension).• In practice: error decreases, though not monotonically.• Theoretically: iterands are optimal under energy norm.• In practice: not optimal under L2norm.• Monotonic bound: #it ∼pκ(A).The matter of SymmetryQuestion: doesrn∈ span(r1, . . . , rn−1).imply that all rnhave to be retained?Answer:• If A is symmetric: no. (CG; 3 vectors)• If A is not symmetric and you want orthogonality: yes. (OrthoDir,OrthoRes)Practical solutions:– Restart– Truncate• If A is not symmetric and you can relax orthogonality: no. (BiCG;5 vectors)OptimalityIdea: take span(r1, . . . , rn−1) and construct a minimal-length affinecombination.• Theoretical optimum obtained exactly. (MinRes from CG, GMRESfrom OrthoRes)• Optimum obtained up to small factor. (QMR from BiCG)• Long sequences iff the original method has them.Speed of convergenceNumber of iterations bounded bypκ(A) orpκ(M−1A).Actual number is complicated function of spectrumPreconditioner: reduce condition number by constant or orderUse of transposeSome methods for non-symmetric matrices need the transposematrix-vector multiplication:• OrthoRes / OrthoDir and GMRES: no.• BiCG and QMR: yes.Problem for implicitly defined matrix-vector product.Sketch of CGFOR j = 1, 2, . . . , mzj= M−1rjρj= ztjrjpj= rj−1−ρjρj−1pj−1πj= ptjApjrj+1= rj+ρjπjApjxj+1= xj+ρjπjpjObservations about CG• preconditioner / matrix not applied as combo• inner products introduce synchronization points• limited storage needs• no coupling of xjand rjsequencesSketch of GMRESFOR j = 1, 2, . . . , mw = M−1AvjFOR i = 1, . . . , jhij= (w, vi)w = w − hijvihj+1,j= kwk2; vj+1= w/kwk2Solve Hmym= kr0ke1x = x0+ [v1, . . . , vm] ymif necessary repeatComputational aspects of GMRES• Matrix-vector product• Preconditioner solve• Inner products: linear in iteration #• Vector updatesComputational kernelsSparse matrix-vector productdo i=1,ns = 0do j=ptr(i),ptr(i+1)-1s = s + a(j) * x(idx(j))end dox(i) = send do3 mem refs / 2 ops ⇒ low efficiency, little cache reuseSparse matrix-vector product, distributedVector elements referenced in stencil fashion:border variableowned variableedge variable, ownedSparse matvec, distributed, cont’dProcessor subdomain references ‘ghost’ border region:subdomainghost regionSparse matvec, distributed, ccont’dMatrix-vector product has local and remote references:ABP0P1Sparse matrix-vector product, algorithmOverlapping communication/computation:Initiate send/recv of boundary dataPerform local matvecComplete recv of boundary dataPerform matvec with boundary data,and updateInner productsInner products are almost unavoidable:Chebyshev method needs spectrum informationInner products are bad news in parallel: synchronisation points ⇒ try tobundle them.Combine inner products:• Three-term methods: combine two inner products per iteration(Chronopoulos/Gear, Saad/Meurant, D’Azevedo/Eijkhout/Romine)⇒ can be stable• Arnoldi methods: use non-modified Gramm-Schmid, GS twice,block orthogonalisation on partial Krylov series.⇒ stability a problem.or hide them (vd Vorst)Advanced iterative methodsReduction GMRES-iterationsrestrict number of GMRES stepsby combination with other method1. further reduction of ri, with e.g. GMRES itself:GMRESR (VDV and Vuik ’94)2. precond of search direction: FGMRES (Saad ’93)GMRESR: more memory; robust, reductions addFGMRES: cheaper in memory; not robustGMRESRr0= b − Ax0;for i = 0, 1, 2, . . .zmapprox. sol. of Az = ri(nested GMRES)c = Azmfor k = 0, 1, . . . , i − 1zm= zm− αuk, c = c − αck(α = (ck, c))ci= c/kck2; ui= zm/kck2xi+1= xi+ (ci, ri)ui, ri+1= ri− (ci, ri)ci(FGMRES works with ui)Example of GMRESR; Navier StokesGMRES GMRES(50) GMRESR+GMRES(10)as ”solve”matvec 184 1220 198daxpy 17000 35000 1386ddot 17000 35000 1224memory 184 50 46cpu-time 7.2 17.0 1.2Choice of m in GMRESRGMRES: 177 it., GMRES(50): 1323 itGMRESR(m):m it m*it cpu-time memory4 45 180 1.15 948 23 184 0.89 5412 16 192 0.90 4416 12 192 1.05 4020 10 200 1.23 40Bi-CG and variantswith short recurrencies we can constructxisuch that ri⊥ Ki(AT; s0)• not optimal in Ki(A;
View Full Document
Unlocking...