Absolute value preconditioning for symmetric linear systems and - - PowerPoint PPT Presentation

Absolute value preconditioning for symmetric linear systems and eigenvalue problems

Eugene Vecharynski (joint work with Andrew Knyazev)

Georgia State University Department of Mathematics and Statistics

Scientific Computing Group Meeting Emory University, Atlanta, GA October 17, 2012

1 / 46

Two problems

1

Symmetric indefinite linear system Ax = b, A = A∗

Mixed FE discretizations of PDEs in fluid and solid mechanics Discretizations of PDEs in acoustics Inner steps of interior point methods in optimization Inner steps of symmetric eigensolvers (e.g., Jacobi-Davidson)

2

Symmetric interior eigenvalue problem Kv = λMv, K = K ∗, M = M∗ > 0

Vibration modes analysis Electronic structure calculations

2 / 46

Ax = b and Kv = λMv

General setting Large size, ill-conditioned Sparse matrices or matrix-free environment Direct methods are inapplicable. Iterate! Use preconditioners to improve convergence

3 / 46

Ax = b and Kv = λMv

General setting Large size, ill-conditioned Sparse matrices or matrix-free environment Direct methods are inapplicable. Iterate! Use preconditioners to improve convergence ⇒ This talk concerns the construction of symmetric positive definite (SPD) preconditioners for both problems

4 / 46

Preconditioning

Traditionally, For Ax = b, define T ≈ A−1 For Kv = λMv, define T ≈ (K − c2M)−1 ⇒ T is not SPD!

5 / 46

Outline

In this talk ... Advocate SPD preconditioning Motivate and introduce absolute value (AV) preconditioning Construct an AV preconditioner for a model linear system and eigenvalue problem

6 / 46

AV preconditioning

Symmetric indefinite linear systems

7 / 46

How to solve Ax = b, A = A∗?

Preconditioned system TAx = Tb

8 / 46

How to solve Ax = b, A = A∗?

Preconditioned system TAx = Tb T is symmetric indefinite or nonsymmetric (e.g., T ≈ A−1) ⇒ TA is nonsymmetric in any inner product.

The symmetry is lost. Methods: GMRES, BiCG, BiCGstab, QMR, etc.

9 / 46

How to solve Ax = b, A = A∗?

Preconditioned system TAx = Tb T is symmetric indefinite or nonsymmetric (e.g., T ≈ A−1) ⇒ TA is nonsymmetric in any inner product.

The symmetry is lost. Methods: GMRES, BiCG, BiCGstab, QMR, etc.

T is SPD. ⇒ TA is symmetric in the T −1-based inner product (·, ·)T −1 = (·, T −1·)

The symmetry is preserved. Method: PMINRES (optimal, short-term recurrent).

10 / 46

How to solve Ax = b, A = A∗?

Preconditioned system TAx = Tb T is symmetric indefinite or nonsymmetric (e.g., T ≈ A−1) ⇒ TA is nonsymmetric in any inner product.

The symmetry is lost. Methods: GMRES, BiCG, BiCGstab, QMR, etc.

T is SPD. ⇒ TA is symmetric in the T −1-based inner product (·, ·)T −1 = (·, T −1·)

The symmetry is preserved. Method: PMINRES (optimal, short-term recurrent).

⇒ How to correctly define an SPD preconditioner T?

11 / 46

Matrix absolute value

The eigenvalue decomposition A = VΛV ∗ The matrix absolute value of A is defined as |A| = V |Λ| V ∗, |Λ| = diag{

• λj
• }.

The matrix sign of A is defined as sign(A) = Vsign(Λ)V ∗, sign(Λ) = diag

• sign(λj)
• .

The polar decomposition A = |A| sign(A) = sign(A) |A| .

12 / 46

An ideal SPD preconditioner

Let T = |A|−1. Then TAx = Tb is sign(A)x = |A|−1 b. The matrix TA = sign(A) has two distinct eigenvalues: −1 and 1. ⇒ PMINRES converges in at most two steps

13 / 46

An ideal SPD preconditioner

Let T = |A|−1. Then TAx = Tb is sign(A)x = |A|−1 b. The matrix TA = sign(A) has two distinct eigenvalues: −1 and 1. ⇒ PMINRES converges in at most two steps T = |A|−1 is an ideal SPD preconditioner Construction of the exact |A|−1 is prohibitively expensive AV preconditioning: Construct an SPD preconditioner T as an approximation to |A|−1, i.e., T ≈ |A|−1 .

14 / 46

Bounds on eigenvalues of TA

Theorem Given a symmetric indefinite A, an SPD T, and constants δ1 ≥ δ0 > 0, such that δ0(v, T −1v) ≤ (v, |A| v) ≤ δ1(v, T −1v), ∀v. Then Λ(TA) ⊂ [−δ1, −δ0]

• [δ0, δ1] .

15 / 46

Challenges

For practical applications: Can we construct T ≈ |A|−1 using A (matvecs with A)? Can we construct T ≈ |A|−1 at an optimal cost O(nnz(A))?

16 / 46

AV preconditioning

Eigenvalue problems

17 / 46

Kv = λMv

18 / 46

Kv = λMv

19 / 46

Kv = λMv

LOBPCG: Λ(k) = (V (k))∗KV (k), (V (k))∗MV (k) = I V (k+1) = V (k)C(k)

V

+ T

• KV (k) − MV (k)Λ(k)

C(k)

R

+ P(k)C(k)

P

P(k+1) = V (k+1) − V (k)C(k)

V

20 / 46

Kv = λMv

LOBPCG: Λ(k) = (V (k))∗KV (k), (V (k))∗MV (k) = I V (k+1) = V (k)C(k)

V

+ T

• KV (k) − MV (k)Λ(k)

C(k)

R

+ P(k)C(k)

P

P(k+1) = V (k+1) − V (k)C(k)

V

Choice of the preconditioner T ≈ K −1 ⇒ emphasizes the extreme (smallest) eigenpairs T ≈

• K − c2M

−1 ⇒ indefinite preconditioning Our approach: T ≈

• K − c2M
• −1

21 / 46

The same challenges!

For practical applications: Can we construct T ≈ |K − c2M|−1 using K and M? Can we construct T ≈ |K − c2M|−1 at an optimal cost O(nnz(K)+nnz(M))?

22 / 46

Model problems

The “shifted Laplacian” equation −∆u(x, y) − c2u(x, y) = f(x, y), (x, y) ∈ Ω = (0, 1) × (0, 1) u|Γ = 0. The eigenvalue problem (find eigenpairs near c2) −∆u(x, y) = λu(x, y), (x, y) ∈ Ω = (0, 1) × (0, 1) u|Γ = 0.

23 / 46

Model problems

After discretization (FD): Symmetric indefinite linear system (L − c2I)x = b ⇒ T ≈ |L − c2I|−1 Symmetric eigenvalue problem (find eigenpairs near c2). Lv = λv ⇒ T ≈ |L − c2I|−1 ⇒ Apply T by approximately solving |L − c2I|w = r. Use MG!

24 / 46

MG

25 / 46

Two-grid scheme for |L − c2I|w = r

Algorithm (Two-grid scheme for |L − c2I|w = r) Input r. Output w = Tr.

1

• Presmoothing. Apply ν smoothing steps,

w(i+1) = w(i) + τ(r − |L − c2I|w(i)), i = 0, . . . , ν − 1, w(0) = 0. Set wpre = w(ν), ν ≥ 1.

2

Coarse grid correction. Restrict residual to the coarse grid (R), apply

• LH − c2IH
• −1, prolongate to the fine grid (P), and add to wpre,

wcgc = wpre + P|LH − c2IH|−1R

• r − |L − c2I|wpre

.

3

• Postsmoothing. Apply ν smoothing steps,

w(i+1) = w(i) + τ(r − |L − c2I|w(i)), i = 0, . . . , ν − 1, w(0) = wcgc. Return w = w(ν).

26 / 46

Two-grid scheme for |L − c2I|w = r

Algorithm (Two-grid scheme for |L − c2I|w = r) Input r. Output w = Tr.

1

• Presmoothing. Apply ν smoothing steps,

w(i+1) = w(i) + τ(r − |L − c2I|w(i)), i = 0, . . . , ν − 1, w(0) = 0. Set wpre = w(ν), ν ≥ 1.

2

Coarse grid correction. Restrict residual to the coarse grid (R), apply

• LH − c2IH
• −1, prolongate to the fine grid (P), and add to wpre,

wcgc = wpre + P|LH − c2IH|−1R

• r − |L − c2I|wpre

.

3

• Postsmoothing. Apply ν smoothing steps,

w(i+1) = w(i) + τ(r − |L − c2I|w(i)), i = 0, . . . , ν − 1, w(0) = wcgc. Return w = w(ν).

Idea: Replace |L − c2I| ≈ B, where y = Bu is easy to compute

27 / 46

Two-grid AV preconditioner

Algorithm (Two-grid AV preconditioner) Input r, B ≈ |L − c2I|. Output w = Tr.

1

• Presmoothing. Apply ν smoothing steps,

w(i+1) = w(i) + τ(r − Bw(i)), i = 0, . . . , ν − 1, w(0) = 0. Set wpre = w(ν), ν ≥ 1.

2

Coarse grid correction. Restrict residual to the coarse grid (R), apply

• LH − c2IH
• −1, prolongate to the fine grid (P), and add to wpre,

wcgc = wpre + P|LH − c2IH|−1R

• r − Bwpre

.

3

• Postsmoothing. Apply ν smoothing steps,

w(i+1) = w(i) + τ(r − Bw(i)), i = 0, . . . , ν − 1, w(0) = wcgc. Return w = w(ν).

28 / 46

MG AV preconditioner (V-cycle)

Algorithm (AV-MG(rl): the MG AV preconditioner) Input rl, output wl.

1

Set Bl ≈ |Ll − c2Il|.

2

• Presmoothing. Apply ν smoothing steps,

w(i+1)

l

= w(i)

l

+ τl(rl − Blw(i)

l ), i = 0, . . . , ν − 1, w(0) l

= 0. Set wpre

l

= w(ν)

l

, ν ≥ 1.

3

Coarse grid correction. Restrict residual to coarser grid (Rl−1), apply AV-MG on this level, prolongate to the finer grid (Pl), add to wpre

l

, wcgc

l

= wpre

l

+ Pl AV-MG(Rl−1

• rl − Blwpre

l

• ), if l > 1;

wcgc

1

= wpre

1

+ P1|L0 − c2I0|−1R0

• r1 − B1wpre

1

• , if l = 1.

4

• Postsmoothing. Apply ν smoothing steps,

w(i+1)

l

= w(i)

l

+ τl(rl − Blw(i)

l ), i = 0, . . . , ν − 1, w(0) l

= wcgc

l

. Return wl = w(ν)

l

.

29 / 46

How to choose Bl ≈ |Ll − c2Il|?

Observations: On sufficiently fine grids (chl << 1), |Ll − c2Il| ≈ Ll. Since |Ll − c2Il| =

• 2h(Ll − c2Il) − I

Ll − c2Il

• ,

|Ll − c2Il| ≈

• 2qm−1(Ll − c2Il) − I

Ll − c2Il

• pm(Ll−c2Il)

, where qm−1 is a least-squares polynomial approximation of the Heaviside step function h (polynomial filtering). On coarser grids, evaluation of y = pm(Ll − c2Il)u is inexpensive.

30 / 46

How to choose Bl ≈ |Ll − c2Il|?

Observations: On sufficiently fine grids (chl << 1), |Ll − c2Il| ≈ Ll. Since |Ll − c2Il| =

• 2h(Ll − c2Il) − I

Ll − c2Il

• ,

|Ll − c2Il| ≈

• 2qm−1(Ll − c2Il) − I

Ll − c2Il

• pm(Ll−c2Il)

, where qm−1 is a least-squares polynomial approximation of the Heaviside step function h (polynomial filtering). On coarser grids, evaluation of y = pm(Ll − c2Il)u is inexpensive. ⇒ Given a “switching parameter” δ ∈ (0, 1), define Bl = Ll, chl < δ; pm(Ll − c2Il),

• therwise.

31 / 46

Summary of theoretical analysis

The preconditioner is SPD under mild assumptions on the smoother, restriction and prolongation. The size of the coarsest grid h0 should be ∼ 1/c. Two-grid preconditioner reduces the error for |L − c2I|w = r in the directions of almost all Laplacian eigencomponents. For more details, see

Eugene Vecharynski and Andrew Knyazev: Absolute value preconditioning for symmetric indefinite linear systems. Available at http://arxiv.org/abs/1104.4530, submitted

32 / 46

Numerical experiments

In our tests Standard coarsening Full weighting for restriction and piecewise multilinear interpolation for prolongation Choose δ = 1/3 and m = 10, so that Bl = Ll, chl < 1/3; p10(Ll − c2Il),

• therwise.

Computations of p10 are based on filtering with Chebyshev polynomials. Smoother for Bl = Ll: 1 Richardson’s iteration, τl = h2

l /5

Smoother for Bl = pm(Ll − c2Il): 5 Richardson’s iterations, τl = h2

l /(5 − c2h2 l )

33 / 46

Numerical experiments

AV preconditioning for a model linear systems

34 / 46

Spectra of preconditioned matrices

Eigenvalues of T|L − c2 I| and T(L − c2 I) for different shifts 10 10

1

300 400 1500 3000 Λ(T(L − c2 I))>0 −10

0 −10 −1

300 400 1500 3000 Λ(T(L − c2 I))<0 10

010 1

300 400 1500 3000 Shift value c2 Λ(T|L − c2 I|)

Figure : Spectrum of T

• L − c2I
• (left), negative eigenvalues of T
• L − c2I
• (center), and positive eigenvalues of T
• L − c2I
• (right); n = 16129

35 / 46

Available preconditioning for (L − c2I)x = b

Compare AV-MG with SPD preconditioning: T = L−1 (Bayliss, Goldstein, Turkel 83) Indefinite preconditioning: T = MG V-cycle for (L − c2I)x = b Many related works: e.g., Elman, Ernst, O’Leary 01; Bramble, Leyk, Pasciak 93; van Gijzen, Erlangga, Vuik 07, etc.

36 / 46

Comparisons with available preconditioned schemes

10 20 30 40 50 60 10

−10

10

−5

10 10

5

Shift value c2 = 300 (19 negative eigenvalues) Iteration number Euclidean norm of error

Bi−CGSTAB−MG GMRES(5)−MG MINRES−AV−MG MINRES−Laplace 20 40 60 80 10

−10

10

−5

10

Shift value c2 = 400 (26 negative eigenvalues) Iteration number Euclidean norm of error

Bi−CGSTAB−MG GMRES(10)−MG MINRES−AV−MG MINRES−Laplace

Figure : Coarsest problems of size 225; n = 65025

37 / 46

Comparisons with available preconditioned schemes

100 200 300 400 500 10

−10

10

−5

10 10

5

Shift value c2 = 1500 (108 negative eigenvalues) Iteration number Euclidean norm of error

Bi−CGSTAB−MG GMRES(20)−MG MINRES−AV−MG MINRES−Laplace 100 200 300 400 500 600 700 10

−5

10

Shift value c2 = 3000 (222 negative eigenvalues) Iteration number Euclidean norm of error

Bi−CGSTAB−MG GMRES(150)−MG MINRES−AV−MG MINRES−Laplace

Figure : Coarsest problems of size 961; n = 65025

38 / 46

Mesh-independent convergence

Number of steps performed to decrease error norm by 10−8.

h = 2−7 h = 2−8 h = 2−9 h = 2−10 h = 2−11 c2 = 300 31 30 30 30 30 c2 = 400 38 37 37 37 37 c2 = 1500 97 89 88 89 90 c2 = 3000 222 279 256 257 256 Table : Mesh-independent convergence of MINRES-AV-MG.

39 / 46

Numerical experiments

AV preconditioning for a model eigenvalue problem

40 / 46

The eigenvalue problem Lv = λv

41 / 46

The eigenvalue problem Lv = λv

42 / 46

Convergence to eigenpairs (λ25, v25) through (λ30, v30)

2 4 6 8 10 12 14 10

−2

10 10

2

10

4

LOBPCG with different preconditioners Euclidean norm of residuals Iteration number

— “Standard” SPD preconditioning: T ≈ L−1 (MG V-cycle for L)

43 / 46

Convergence to eigenpairs (λ25, v25) through (λ30, v30)

2 4 6 8 10 12 14 10

−4

10

−2

10 10

2

10

4

LOBPCG with different preconditioners Euclidean norm of residuals Iteration number

— “Standard” SPD preconditioning: T ≈ L−1 (MG V-cycle for L) — Indefinite preconditioning: T ≈ (L − c2I)−1 (MG V-cycle for L − c2I)

44 / 46

Convergence to eigenpairs (λ25, v25) through (λ30, v30)

2 4 6 8 10 12 14 10

−6

10

−4

10

−2

10 10

2

10

4

LOBPCG with different preconditioners Euclidean norm of residuals Iteration number

— “Standard” SPD preconditioning: T ≈ L−1 (MG V-cycle for L) — Indefinite preconditioning: T ≈ (L − c2I)−1 (MG V-cycle for L − c2I) — AV preconditioning: T ≈ |L − c2I|−1 (AV-MG)

45 / 46

Current and future work

Domain decomposition AV preconditioning Algebraic MG type AV preconditioning AV preconditioning for saddle point problems Construction of novel AV preconditioned iterative schemes for interior eigenvalue problems Application to electronic structure calculations HPC software development Thank you!

46 / 46