Notebook 18
Practical iterative methods
Incomplete LU (ILU) preconditioning
Idea (Varga, Buleev,
LU factorization has large fill-in. For a preconditioner, just limit the fill-in to a fixed pattern.
Apply the standard LU factorization method, but calculate only a part of the entries, e.g. only those which are larger than a certain threshold value, or only those which correspond to certain predefined pattern.
Result: incomplete LU factors
, , remainder :What about zero pivots which prevent such an algoritm from being computable ?
Theorem (Saad, Th. 10.2): If
is an M-Matrix, then the algorithm to compute the incomplete LU factorization with a given pattern is stable, i.e. does not detriorate due to zero pivots (main diagonal elements) Moreover, where and is a regular splitting.
Discussion
Generally better convergence properties than Jacobi, Gauss-Seidel, though we mostly cannot apply the comparison theorem for regular splittings
Block variants are possible
ILU Variants:
ILUM: ("modified"): add ignored off-diagonal entries to main diagonal
ILUT: ("threshold"): zero pattern calculated dynamically based on drop tolerance
ILU0: Drop all fill-in
Incomplete Cholesky: symmetric variant of ILU
Dependence on ordering
Can be parallelized using graph coloring
Not much theory: experiment for particular systems and see if it works well
I recommend it as the default initial guess for a sensible preconditioner
Further approaches to preconditioning
These are based on ideas which are best explained and developed with multidimensional PDEs in mind.
Multigrid: gives indeed
optimal solver complexity in some situations. This is the holy grail method... I will try to discuss this later in the course.Domain decomposition - based on the idea the subdivision of the computational domain into a number of subdomains and subsequent repeated solution of the smaller subdomain problems
Iterative methods in Julia
Julia has some well maintained packages for iterative methods and preconditioning.
IterativeSolvers.jl: various Krylov subspace methods including conjugate gradients
IncompleteLU.jl: Incomplete LU factorizations
AlgebraicMultigrid.jl: Algebraic multigrid methods
Random sparse M-Matrices
We will test the methods with random sparse M matrices, so we define a function which gives us a random, strictly diagonally dominant M-Matrix which is not necessarily irreducible. For skew=0
it is also symmetric:
sprandm (generic function with 1 method)
Test the method a bit...
5
x1 | x2 | x3 | x4 | x5 | |
---|---|---|---|---|---|
Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | 2.6006 | -0.922413 | -0.165199 | -0.840878 | -0.462586 |
2 | -0.966885 | 1.1229 | 0.0 | -0.937833 | -0.0462712 |
3 | -0.566097 | 0.0 | 0.165822 | 0.0 | 0.0 |
4 | -0.703205 | 0.0 | 0.0 | 2.32171 | 0.0 |
5 | -0.363993 | -0.20023 | 0.0 | -0.542756 | 0.509066 |
Up to rounding errors, the inverse is nonnegative, as predicted by the theory. There are zero entries because it is not necessarily irreducible. Invertibility is guaranteed by strict diagonal dominance.
5×5 Matrix{Float64}:
308.942 308.846 307.78 308.839 308.807
363.27 364.062 361.903 363.534 363.193
1054.69 1054.36 1056.75 1054.34 1054.23
93.5729 93.544 93.221 93.9725 93.532
463.55 463.762 461.806 464.006 465.343
93.22100392714542
ρ_jacobi (generic function with 1 method)
0.9994562900651217