Preconditioning

Preconditioning

Recall that preconditioning is the act of creating an “affordable” operation “$P^{-1}$” such that $P^{-1} A$ (or $A P^{-1}$) is is well-conditoned or otherwise has a “nice” spectrum. We then solve the system

$$ P^{-1} A x = P^{-1} b \quad \text{or}\quad A P^{-1} \underbrace{(P x)}_y = b $$

in which case the convergence rate depends on the spectrum of the iteration matrix $$ I - \omega P^{-1} A . $$

• The preconditioner must be applied on each iteration.
• It is not merely about finding a good initial guess.

Classical methods

We have discussed the Jacobi preconditioner

$$ P_{\text{Jacobi}}^{-1} = D^{-1} $$

where $D$ is the diagonal of $A$. Gauss-Seidel is

$$ P_{GS}^{-1} = (L+D)^{-1} $$

where $L$ is the (strictly) lower triangular part of $A$. The upper triangular part may be used instead, or a symmetric form

$$ P_{SGS}^{-1} = (D+U)^{-1} D (L+D)^{-1} . $$

Domain decomposition

Given a linear operator $A : V \to V$, suppose we have a collection of prolongation operators $P_i : V_i \to V$. The columns of $P_i$ are “basis functions” for the subspace $V_i$. The Galerkin operator $A_i = P_i^T A P_i$ is the action of the original operator $A$ in the subspace.

Define the subspace projection

$$ S_i = P_i A_i^{-1} P_i^T A . $$

• $S_i$ is a projection: $S_i^2 = S_i$
• If $A$ is SPD, $S_i$ is SPD with respect to the $A$ inner product $x^T A y$
• $I - S_i$ is $A$-orthogonal to the range of $P_i$

These projections may be applied additively

$$ I - \sum_{i=0}^n S_i, $$

multiplicatively

$$ \prod_{i=0}^n (I - S_i), $$

or in some hybrid manner, such as

$$ (I - S0) (I - \sum{i=1}^n S_i) . $$ In each case above, the action is expressed in terms of the error iteration operator.

Examples

• Jacobi corresponds to the additive preconditioner with $P_i$ as the $i$th column of the identity
• Gauss-Seidel is the multiplicate preconditioner with $P_i$ as the $i$th column of the identity
• Block Jacobi corresponds to labeling “subdomains” and $P_i$ as the columns of the identity corresponding to non-overlapping subdomains
• Overlapping Schwarz corresponds to overlapping subdomains
• $P_i$ are eigenvectors of $A$
• A domain is partitioned into interior $VI$ and interface $V\Gamma$ degrees of freedom. $PI$ is embedding of the interior degrees of freedom while $P\Gamma$ is “harmonic extension” of the interface degrees of freedom. Consider the multiplicative combination $(I - S_\Gamma)(I - S_I)$.

Convergence theory

The formal convergence is beyond the scope of this course, but the following estimates are useful. We let $h$ be the element diameter, $H$ be the subdomain diameter, and $\delta$ be the overlap, each normalized such that the global domain diameter is 1. We express the convergence in terms of the condition number $\kappa$ for the preconditioned operator.

• (Block) Jacobi: $\delta=0$, $\kappa \sim H^{-2} H/h = (Hh)^{-1}$
• Overlapping Schwarz: $\kappa \sim H^{-2} H/\delta = (H \delta)^{-1}$
• 2-level overlapping Schwarz: $\kappa \sim H/\delta$

Hands-on with PETSc: demonstrate these estimates

• Symmetric example: src/snes/examples/tutorials/ex5.c
• Nonsymmetric example: src/snes/examples/tutorials/ex19.c
• Compare preconditioned versus unpreconditioned norms.
• Compare BiCG versus GMRES
• Compare domain decomposition and multigrid preconditioning
  • -pc_type asm (Additive Schwarz)
  • -pc_asm_type basic (symmetric, versus restrict)
  • -pc_asm_overlap 2 (increase overlap)
  • Effect of direct subdomain solver: -sub_pc_type lu
  • -pc_type mg (Geometric Multigrid)
• Use monitors:
  • -ksp_monitor_true_residual
  • -ksp_monitor_singular_value
  • -ksp_converged_reason
• Explain methods: -snes_view
• Performance info: -log_view

Examples

mpiexec -n 4 ./ex19 -lidvelocity 2 -snes_monitor -da_refine 5 -ksp_monitor -pc_type asm -sub_pc_type lu