Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Advanced Topics from Scientific Computing
TU Berlin Winter 2022/23
Notebook 08
Creative Commons Lizenzvertrag Jürgen Fuhrmann

html"""
<b> Advanced Topics from Scientific Computing<br>
TU Berlin Winter 2022/23<br>
Notebook 08<br>
<a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/"><img alt="Creative Commons Lizenzvertrag" style="border-width:0" src="https://i.creativecommons.org/l/by-sa/4.0/80x15.png" /></a>
Jürgen Fuhrmann
</b> <br>
"""
3.5 ms
begin
using Graphs
using ExtendableSparse
using GraphPlot
using Colors
using LinearAlgebra
using SparseArrays
using IterativeSolvers
end
4.8 s

Iterative methods for linear systems

md"""
# Iterative methods for linear systems
"""

115 ms

Let $V=\mathbb R^n$ be equipped with the inner product $(\cdot,\cdot)$. Let $A$ be an $n\times n$ nonsingular matrix.

Solve $Au=b$ iteratively. For this purpose, two components are needed:

  • Preconditioner: a matrix $M\approx A$ "approximating" the matrix $A$ but with the property that the system $Mv=f$ is easy to solve

  • Iteration scheme: algorithmic sequence using $M$ and $A$ which updates the solution step by step

md"""
Let $V=\mathbb R^n$ be equipped with the inner product $(\cdot,\cdot)$.
Let $A$ be an $n\times n$ nonsingular matrix.
Solve $Au=b$ iteratively. For this purpose, two components are
needed:
- __Preconditioner__: a matrix $M\approx A$ "approximating" the matrix
$A$ but with the property that the system $Mv=f$ is easy to solve
- __Iteration scheme__: algorithmic sequence using $M$ and $A$ which
updates the solution step by step

"""

517 ms

Simple iteration scheme

md"""
## Simple iteration scheme
"""

6.6 ms

Assume we know the exact solution $\hat u$: $A\hat u=b$.

Then it must fulfill the identity

$$ \begin{aligned} \hat u= \hat u - M^{-1}(A \hat u -b) \end{aligned}$$

$\Rightarrow$ iterative scheme: put the "old" value on the right hand side and the ""new" value on the left hand side:

$$ \begin{aligned} u_{k+1}= u_k - M^{-1}(A u_k -b) \quad (k=0,1\dots) \end{aligned}$$

Obviously, if $u_k=\hat u$, the process would be stationary.

Otherwise it leads to a sequence of approximations

$$u_0, u_1,\dots, u_k, u_{k+1},\dots$$

md"""
Assume we know the exact solution $\hat u$: $A\hat u=b$.

Then it must fulfill the identity
```math
\begin{aligned}
\hat u= \hat u - M^{-1}(A \hat u -b)
\end{aligned}
```
$\Rightarrow$ iterative scheme: put the "old" value on the right hand
side and the ""new" value on the left hand side:
```math
\begin{aligned}
u_{k+1}= u_k - M^{-1}(A u_k -b) \quad (k=0,1\dots)
\end{aligned}
```
Obviously, if $u_k=\hat u$, the process would be stationary.

Otherwise it leads to a sequence of approximations

```math
u_0, u_1,\dots, u_k, u_{k+1},\dots
```
"""

370 μs

Implementation: solve $Au=b$ with tolerance $\varepsilon$:

  1. Choose initial value $u_0$, set $k=0$

  2. Calculate residuum $r_k=Au_k -b$

  3. Test convergence: if $||r_k||<\varepsilon$ set $u=u_k$, finish

  4. Calculate update: solve $M v_k = r_k$

  5. Update solution: $u_{k+1}=u_k - v_k$, set $k=k+1$, repeat with step 2.

md"""
Implementation: solve $Au=b$ with tolerance $\varepsilon$:
1. Choose initial value $u_0$, set $k=0$
2. Calculate _residuum_ $r_k=Au_k -b$
3. Test convergence: if $||r_k||<\varepsilon$ set $u=u_k$, finish
4. Calculate _update_: solve $M v_k = r_k$
5. Update solution: $u_{k+1}=u_k - v_k$, set $k=k+1$, repeat with step 2.

"""

173 ms

General convergence theorem

md"""
### General convergence theorem
"""

5.6 ms

Let $\hat u$ be the solution of $Au=b$.

Let $e_k= u_k - \hat u$ be the error of the $k$-th iteration step. Then:

$$ \begin{aligned} u_{k+1}&= u_k - M^{-1}(A u_k -b) \\ & = (I-M^{-1}A) u_k + M^{-1}b\\ u_{k+1}-\hat u &= u_k -\hat u - M^{-1}(A u_k - A \hat u)\\ &=(I-M^{-1}A)( u_k - \hat u)\\ &=(I-M^{-1}A)^k( u_0 - \hat u) \end{aligned}$$

resulting in

$$ \begin{aligned} e_{k+1}= (I-M^{-1}A)^k e_0 \end{aligned}$$

  • So when does $(I-M^{-1}A)^k$ converge to zero for $k\to\infty$ ?

  • Denote $B=I-M^{-1}A$

md"""
Let $\hat u$ be the solution of $Au=b$.

Let $e_k= u_k - \hat u$ be the error of the $k$-th iteration step. Then:
```math
\begin{aligned}
u_{k+1}&= u_k - M^{-1}(A u_k -b) \\
& = (I-M^{-1}A) u_k + M^{-1}b\\
u_{k+1}-\hat u &= u_k -\hat u - M^{-1}(A u_k - A \hat u)\\
&=(I-M^{-1}A)( u_k - \hat u)\\
&=(I-M^{-1}A)^k( u_0 - \hat u)
\end{aligned}
```
resulting in
```math
\begin{aligned}
e_{k+1}= (I-M^{-1}A)^k e_0
\end{aligned}
```
- So when does $(I-M^{-1}A)^k$ converge to zero for $k\to\infty$ ?
- Denote $B=I-M^{-1}A$
"""

3.1 ms

Definition The spectral radius $\rho(B)$ is the largest absolute value of any eigenvalue of $B$: $\rho(B)= \max_{\lambda\in\sigma(B)} |\lambda|$.

definition"""
__Definition__
The spectral radius $\rho(B)$ is the largest absolute
value of any eigenvalue of $B$:
$\rho(B)= \max_{\lambda\in\sigma(B)} |\lambda|$.
"""

10.9 ms

Sufficient condition for iterative method convergence:

$$\rho(I-M^{-1}A)<1$$

important"""
__Sufficient condition for iterative method convergence:__
```math
\rho(I-M^{-1}A)<1
```
"""

73.7 μs

Asymptotic convergence factor $\rho_{it}$ can be estimated via the spectral radius:

$$ \begin{aligned} \rho_{it}&= \lim_{k\to\infty}\left(\max_{u_0}\frac{|| (I-M^{-1}A)^k(u_0-\hat u)||}{||u_0-\hat u||}\right)^{\frac1k}\\ &= \lim_{k\to\infty}||(I-M^{-1}A)^k||^{\frac1k}\\ &=\rho(I-M^{-1}A) \end{aligned}$$

Depending on $u_0$, the rate may be faster, though

important"""
__Asymptotic convergence factor__ ``\rho_{it}`` can be estimated via the spectral
radius:
```math
\begin{aligned}
\rho_{it}&= \lim_{k\to\infty}\left(\max_{u_0}\frac{|| (I-M^{-1}A)^k(u_0-\hat u)||}{||u_0-\hat u||}\right)^{\frac1k}\\
&= \lim_{k\to\infty}||(I-M^{-1}A)^k||^{\frac1k}\\
&=\rho(I-M^{-1}A)
\end{aligned}
```
Depending on ``u_0``, the rate may be faster, though

"""
2.4 ms

Convergence estimate for symmetric positive definite A,M

md"""
### Convergence estimate for symmetric positive definite A,M
"""
152 μs

Matrix preconditioned Richardson iteration: $M$, $A$ spd.

Scaled Richardson iteration with preconditioner $M$

$$ u_{k+1}= u_k- \alpha M^{-1}(Au_k -b)$$

Spectral equivalence estimate

$$ 0<\gamma_{min}(Mu,u) \leq (Au,u) \leq \gamma_{max}(Mu,u)$$

$\Rightarrow$ $\gamma_{min}\leq \lambda_i \leq \gamma_{max}$

$\Rightarrow$ optimal parameter $\alpha=\frac{2}{\gamma_{max}+\gamma_{min}}$

$\Rightarrow$ convergence rate with optimal parameter: $\rho_{opt}\leq\frac{\kappa(M^{-1}A)-1}{\kappa(M^{-1}A)+1}$ where $\kappa(M^{-1}A)\leq \frac{\gamma_{max}}{\gamma_{min}}$

statement"""
__Matrix preconditioned Richardson iteration:__ ``M``, ``A`` spd.

Scaled Richardson iteration with preconditioner ``M``
```math
u_{k+1}= u_k- \alpha M^{-1}(Au_k -b)
```

Spectral equivalence estimate
```math
0<\gamma_{min}(Mu,u) \leq (Au,u) \leq \gamma_{max}(Mu,u)
```

``\Rightarrow`` ``\gamma_{min}\leq \lambda_i \leq \gamma_{max}``

``\Rightarrow`` optimal parameter ``\alpha=\frac{2}{\gamma_{max}+\gamma_{min}}``

``\Rightarrow`` convergence rate with optimal parameter: ``\rho_{opt}\leq\frac{\kappa(M^{-1}A)-1}{\kappa(M^{-1}A)+1}`` where ``\kappa(M^{-1}A)\leq \frac{\gamma_{max}}{\gamma_{min}}``


"""

169 μs

Regular splittings

md"""
### Regular splittings
"""
143 μs

Definiton

  • $A=M-N$ is a regular splitting if

    • $M$ is nonsingular

    • $M^{-1}\geq 0$, $N\geq 0$ are element-wise nonnegative

definition"""
__Definiton__
-
``A=M-N`` is a regular splitting if

-
``M`` is nonsingular
-
``M^{-1}\geq 0``, ``N\geq 0`` are element-wise nonnegative
"""

133 μs

Just remark that in this case $M^{-1}N=I-M^{-1}A$, and that we don't assume symmetry.

md"""
Just remark that in this case ``M^{-1}N=I-M^{-1}A``, and that we don't assume symmetry.
"""
94.1 ms

Theorem: Assume $A$ is nonsingular, $A^{-1}\geq 0$, and $A=M-N$ is a regular splitting. Then $\rho(M^{-1}N)<1$.

statement"""
__Theorem__: Assume ``A`` is nonsingular, ``A^{-1}\geq 0``, and
``A=M-N`` is a regular splitting. Then ``\rho(M^{-1}N)<1``.

"""

85.0 μs

With this theory we cannot say much about the value of the convergence rate, but we have a comparison theorem:

md"""
With this theory we cannot say much about the value of the convergence rate, but we have a comparison theorem:
"""
169 μs

Theorem: Let $A^{-1}\geq 0$, $A=M_1-N_1$ and $A=M_2-N_2$ be regular splittings.

If $N_2\geq N_1$, then $1 >\rho(M_2^{-1}N_2) \geq \rho(M_1^{-1}N_1)$.

statement"""
__Theorem:__ Let ``A^{-1}\geq 0``, ``A=M_1-N_1`` and ``A=M_2-N_2`` be
regular splittings.


If ``N_2\geq N_1``, then
``1 >\rho(M_2^{-1}N_2) \geq \rho(M_1^{-1}N_1)``.

"""
87.4 μs

What can we say about inverse nonnegative matrices ?

md"""
What can we say about inverse nonnegative matrices ?
"""
144 μs
Loading...