SLIDE 1
Summary of the course 2014h ed. Anton Evgrafov November 19, 2014 - - PowerPoint PPT Presentation
Summary of the course 2014h ed. Anton Evgrafov November 19, 2014 - - PowerPoint PPT Presentation
Summary of the course 2014h ed. Anton Evgrafov November 19, 2014 Linear algebra background: Gauss elimination, LU/Cholesky/LDL-factorization; Rank/null space/range/adjoint; Eigenvalues, -vectors, -spaces; algebraic and geometric
SLIDE 2
SLIDE 3
QR-factorization
A = QR, Q∗Q = I, Rij = 0, i ≥ j Used for:
◮ Solving least-squares problems (e.g. in
GMRES/DQGMRES/MINRES/SYMMLQ)
◮ QR-algorithm for eigenvalues
Computational algorithms:
◮ (modified) Gramm-Schmidt orthogonalization ◮ Hausholder reflections ◮ Givens rotations
SLIDE 4
Diagonalization methods
Rare instances, known:
- 1. Eigenvalues/vectors A = QΛQT
- 2. Fast algorithm (e.g. FFT) v ↔ Qv
Ax = b = ⇒ x = Q[Λ−1(QTb)]
SLIDE 5
Matrix-splitting methods
A = M − N xk+1 = M−1(Nxk + b) ek+1 = M−1Nek
Theorem
ek → 0 ∀x0 ⇐ ⇒ ρ(M−1N) < 1. Asymptotically ek+1/ek ∼ ρ(M−1N). Slow!
SLIDE 6
Matrix-splitting methods
Standard examples:
◮ Jacobi: M = diagA ◮ Gauss-Seidel: M = triuA ◮ Block variations
Use:
◮ Smoothers in multigrid ◮ Preconditioners: additive Schwarz ≈ block-Jacobi;
multiplicative Schwarz ≈ block-GS
SLIDE 7
Projection operators
P2 = P
◮ I − P also projector ◮ Null space/range: kerP ⊕ rangeP = Cn ◮ y = Px ⇐
⇒ y ∈ rangeP & x − y = (I − P)x ∈ kerP ⇐ ⇒ y ∈ rangeP & x − y = (I − P)x ⊥ [kerP]⊥
◮ Orthogonal projectors: PT = P
Used in this course: analysis of projection methods
SLIDE 8
Projection/Petrov-Galerkin framework
dim K = dim L << dim Cn Find: x ∈ x0 + K such that r = b − Ax ⊥ L
SLIDE 9
Projection/Petrov-Galerkin framework
dim K = dim L << dim Cn Find: x ∈ x0 + K such that r = b − Ax ⊥ L
◮ Error projection: A SPD, L = K
Theorem
x ∈ arg min
y∈K y − A−1bA
SLIDE 10
Projection/Petrov-Galerkin framework
dim K = dim L << dim Cn Find: x ∈ x0 + K such that r = b − Ax ⊥ L
◮ Error projection: A SPD, L = K
Theorem
x ∈ arg min
y∈K y − A−1bA ◮ Residual projection: L = AK
Theorem
x ∈ arg min
y∈K b − Ay2
SLIDE 11
1D variations:
◮ L = K = ek ⇐
⇒ Gauss-Seidel
◮ L = K = r = b − Ax ⇐
⇒ steepest descent
◮ . . .
SLIDE 12
Krylov subspace methods:
K = Km(A, r0) = spanr0, Aro, . . . , Am−1r0
◮ L = K: FOM, CG, SYMMLQ ◮ L = AK: GMRES, MINRES
SLIDE 13
Arnoldi:
Generating ON basis for Km(A, r0)! V T
m Vm = I, rangeVm = Km(A, r0)
AVm = Vm+1 ¯ Hm, [ ¯ Hm]ij = 0 if i > j + 1. Algorithms: (modified) Gramm-Schmidt, Hausholder Arnoldi used in: FOM, GMRES
SLIDE 14
Cannot store Arnoldi vectors:
◮ Restarts: FOM(k), GMRES(k) ◮ Incomplete (partial) orthogonalization: IOM/DIOM,
QGMRES/DQGMRES
SLIDE 15
Lanczos:
A = AT = ⇒ Hm = HT
m =
⇒ tri-diagonal! In Gramm-Schmidt, only need to orthogonolaze to two previous Arnoldi vectors (3-term recursion) = ⇒ computational savings O(m2n) → O(mn). Lanczos used in: D-LANCZOS (≈ CG), MINRES (≈ DQGMRES
- n symmetric A), SYMMLQ, Lanczos method for eigenvalues
(+Rayleigh-Ritz)
SLIDE 16
“Typical” convergence theorem:
Convergence, GMRES, diagonalizable case: A = XΛX −1
Theorem
rk2 ≤ κ(X) min
p max i
|p(λi)| p-polynomial degree k, p(0) = 1.
SLIDE 17
Preconditioning:
◮ Left: M−1Ax = M−1b ◮ Right: AM−1u = b, x = M−1u.
Goal: M−1A is easier for Krylov methods than A (e.g. has few clustered eigenvalues) M−1 ≈ A−1 but “simple” to apply (every Krylov iteration) ILU, SAI, multigrid, domain decomposition (DD). . .
SLIDE 18
2-grid algorithm:
- 1. (pre-)smooth error/residual: underrelaxed Jacobi, GS, . . .
- 2. restrict residual to coarse grid: rH = I H
h rh
- 3. appx. solve AHeH = −rH
- 4. interpolate onto fine grid: eh = I h
HeH
- 5. apply correction: xh = xh − eh
- 6. (post-)smooth error/residual after interpolations
- 7. repeat
Multigrid: recursion at step 3.
SLIDE 19
Domain decomposition:
◮ PDE-level: overlapping (Schwarz) vs non-overlapping (Schur).
Then discretize and solve! Use special block-structure of the discretized system when solving.
◮ Discrete level: additive/multiplicative Schwarz
SLIDE 20
SVD decomposition and applications
A = UΣV ∗
◮ U, V - unitary ◮ Σ - diagonal, positive
σ = [λ(AA∗)]1/2 = [λ(A∗A)]1/2 Uses:
◮ Least-squares ◮ Low rank approximation of A ◮ Range space ◮ Matrix norms
SLIDE 21
Eigenvalue algorithms for small problems:
Idea: eigenvalue revealing similarity transformation (Schur factorization): A = QRQT. Phase I: Hessenberg (non-symmetric)/tridiagonalization (symmetric) A = QHQT, QTQ = I, Hij = 0, i > j + 1 Phase I algorithms: Hausholder transformations
SLIDE 22
Eigenvalue algorithms for small problems:
Phase II algorithms:
◮ Rayleigh quotient: λ ≈ r(v) = [vTAv]/[vTv] ◮ Power iteration: vk = Akv0/ · ; eigenvector corr. largest in
magnitude eigenvalue; slow convergence
◮ Inverse iteration: vk = A−kv0/ · ; smallest in magn
eigenvalue
◮ Shifted inverse iteration: vk = (A − µI)−kv0/ · eigenvalue
closest to µ
◮ Rayleigh iteration: vk = (A − λkI)−kv0/ · , λk = r(vk).
Cubically convergent for symmetric matrices!
SLIDE 23
Eigenvalue algorithms for small problems:
Phase II algorithms:
◮ Simultaneous versions of power iterations; and their relation to ◮ QR-algorithm with shifts: QkRk = Ak − µkI;
Ak+1 = RkQk + µkI.
◮ Power/inverse power iteration for first/last column of Qk ◮ Wilkinson shift
SLIDE 24
Perturbation analysis for symmetric eigenvalue problems
◮ Eigenvalues (symm matrices):
|λ(A) − λ(A + E)| ≤ E2 |α − λ(A)| ≤ Av − αv2, ∀v2 = 1
◮ Eigenvectors (symm matrices):
1 2 sin(2 angle) ∼ perturbation2 eigenvalue gap
SLIDE 25
Eigenvalue algorithms for large problems:
◮ Lanczos for tridiagonalization: Tk = QT k AQk ◮ Rayleigh-Ritz idea: λ(A) ≈ λ(Tk); error estimates available ◮ Use Phase-II algorithms on Tk.
SLIDE 26