466 lines
14 KiB
Plaintext
466 lines
14 KiB
Plaintext
---
|
||
engine: julia
|
||
---
|
||
|
||
```{julia}
|
||
#| error: false
|
||
#| echo: false
|
||
#| output: false
|
||
using InteractiveUtils
|
||
```
|
||
|
||
# Lineare Algebra in Julia
|
||
|
||
```{julia}
|
||
using LinearAlgebra
|
||
```
|
||
|
||
Das `LinearAlgebra`-Paket liefert unter anderem:
|
||
|
||
- zusätzliche Subtypen von `AbstractMatrix`: genauso verwendbar, wie andere Matrizen, z.B.
|
||
|
||
- `Tridiagonal`
|
||
- `SymTridiagonal`
|
||
- `Symmetric`
|
||
- `UpperTriangular`
|
||
|
||
|
||
- zusätzliche/erweiterte Funktionen: `norm`, `opnorm`, `cond`, `inv`, `det`, `exp`, `tr`, `dot`, `cross`, ...
|
||
|
||
- einen universellen Solver für lineare Gleichungssysteme: `\`
|
||
- `x = A \ b` löst $A \mathbf{x}=\mathbf{b}$ durch geeignete Matrixfaktorisierung und Vorwärts/Rückwärtssubstition
|
||
|
||
- [Matrixfaktorisierungen](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#man-linalg-factorizations)
|
||
- `LU`
|
||
- `QR`
|
||
- `Cholesky`
|
||
- `SVD`
|
||
- ...
|
||
|
||
- Berechnung von Eigenwerte/-vektoren
|
||
|
||
- `eigen`, `eigvals`, `eigvecs`
|
||
|
||
|
||
- Zugriff auf BLAS/LAPACK-Funktionen
|
||
|
||
## Matrixtypen
|
||
|
||
|
||
|
||
```{julia}
|
||
A = SymTridiagonal(fill(1.0, 4), fill(-0.3, 3))
|
||
```
|
||
|
||
|
||
```{julia}
|
||
B = UpperTriangular(A)
|
||
```
|
||
|
||
Diese Typen werden platzsparend gespeichert. Die üblichen Rechenoperationen sind implementiert:
|
||
|
||
```{julia}
|
||
A + B
|
||
```
|
||
|
||
Lesende Indexzugriffe sind möglich,
|
||
```{julia}
|
||
A[1,4]
|
||
```
|
||
|
||
schreibende nicht unbedingt:
|
||
```{julia}
|
||
#| error: true
|
||
A[1,3] = 17
|
||
```
|
||
|
||
Die Umwandlung in eine 'normale' Matrix ist z.B. mit `collect()` möglich:
|
||
```{julia}
|
||
A2 = collect(A)
|
||
```
|
||
|
||
### Die Einheitsmatrix `I`
|
||
|
||
`I` bezeichnet eine Einheitsmatrix (quadratisch, Diagonalelemente = 1, alle anderen = 0) in der jeweils erforderlichen Größe
|
||
|
||
|
||
```{julia}
|
||
A + 4I
|
||
```
|
||
|
||
|
||
## Normen
|
||
|
||
Um Fragen wie Kondition oder Konvergenz eines Algorithmus studieren zu können, brauchen wir eine Metrik. Für lineare Räume ist es zweckmäßig, die Metrik über eine Norm zu definieren:
|
||
$$
|
||
d(x,y) := ||x-y||
|
||
$$
|
||
|
||
### $p$-Normen
|
||
|
||
Eine einfache Klasse von Normen im $ℝ^n$ sind die $p$-Normen
|
||
$$
|
||
||\mathbf{x}||_p = \left(\sum |x_i|^p\right)^\frac{1}{p},
|
||
$$
|
||
die die die euklidische Norm $p=2$ verallgemeinern.
|
||
|
||
|
||
:::{.callout-note}
|
||
|
||
## Die Max-Norm $p=\infty$
|
||
|
||
Sei $x_{\text{max}}$ die _betragsmäßig_ größte Komponente von $\mathbf{x}\in ℝ^n$. Dann gilt stets
|
||
$$ |x_{\text{max}}| \le ||\mathbf{x}||_p \le n^\frac{1}{p} |x_{\text{max}}|
|
||
$$
|
||
(Man betrachte einen Vektor, dessen Komponenten alle gleich $x_{\text{max}}$ sind bzw. einen Vektor, dessen Komponenten außer $x_{\text{max}}$ alle gleich Null sind.)
|
||
|
||
Damit folgt
|
||
$$
|
||
\lim_{p \rightarrow \infty} ||\mathbf{x}||_p = |x_{\text{max}}| =: ||\mathbf{x}||_\infty.
|
||
$$
|
||
:::
|
||
|
||
|
||
In Julia definiert das `LinearAlgebra`-Paket eine Funktion `norm(v, p)`.
|
||
|
||
```{julia}
|
||
v = [3, 4]
|
||
w = [-1, 2, 33.2]
|
||
|
||
@show norm(v) norm(v, 2) norm(v, 1) norm(v, 4) norm(w, Inf);
|
||
```
|
||
|
||
- Wenn das 2. Argument `p` fehlt, wird `p=2` gesetzt.
|
||
- Das 2. Argument kann auch `Inf` (also $+\infty$) sein.
|
||
- Das 1. Argument kann ein beliebiger Container voller Zahlen sein. Die Summe $\sum |x_i|^p$ erstreckt sich über *alle* Elemente des Containers.
|
||
- Damit ist für eine Matrix `norm(A)` gleich der _Frobenius-Norm_ der Matrix `A`.
|
||
|
||
```{julia}
|
||
A = [1 2 3
|
||
4 5 6
|
||
7 8 9]
|
||
norm(A) # Frobenius norm
|
||
```
|
||
|
||
|
||
|
||
Da Normen homogen unter Multiplikation mit Skalaren sind,
|
||
$||\lambda \mathbf{x}|| = |\lambda|\cdot||\mathbf{x}||$, sind sie durch die Angabe der Einheitskugel vollständig bestimmt. Subadditivität der Norm (Dreiecksungleichung) ist äquivalent zur Konvexität der Einheitskugel
|
||
(Code durch anklicken sichtbar).
|
||
```{julia}
|
||
#| code-fold: true
|
||
#| fig-cap: "Einheitskugeln im $ℝ^2$ für verschiedene $p$-Normen: $p$=0.8; 1; 1.5; 2; 3.001 und 1000"
|
||
using Plots
|
||
|
||
colors=[:purple, :green, :red, :blue,:aqua, :black]
|
||
x=LinRange(-1, 1, 1000)
|
||
y=LinRange(-1, 1, 1000)
|
||
|
||
fig1=plot()
|
||
for p ∈ (0.8, 1, 1.5, 2, 3.001, 1000)
|
||
contour!(x,y, (x,y) -> p * norm([x, y], p), levels=[p], aspect_ratio=1,
|
||
cbar=false, color=[pop!(colors)], contour_labels=true, ylim=[-1.1, 1.1])
|
||
end
|
||
fig1
|
||
```
|
||
Wie man sieht, muß $p\ge 1$ sein, damit die Einheitskugel konvex und $||.||_p$ eine Norm ist.
|
||
|
||
Die Julia-Funktion `norm(v, p)` liefert allerdings für beliebige Parameter `p` ein Ergebnis.
|
||
|
||
|
||
### Induzierte Normen (Operatornormen)
|
||
|
||
Matrizen $A$ repräsentieren lineare Abbildungen $\mathbf{v}\mapsto A\mathbf{v}$. Die von einer Vektornorm Induzierte Matrixnorm beantwortet die Frage:
|
||
|
||
> _„Um welchen Faktor kann ein Vektor durch die Transformation $A$ maximal gestreckt werden?“_
|
||
|
||
Auf Grund der Homogenität der Norm unter Multiplikation mit Skalaren reicht es aus, das Bild der Einheitskugel unter der Transformation $A$ zu betrachten.
|
||
|
||
::: {.callout-tip}
|
||
## Definition
|
||
Sei $V$ ein Vektorraum mit einer Dimension $0<n<\infty$ und
|
||
$A$ eine $n\times n$-Matrix. Dann ist
|
||
$$
|
||
||A||_p = \max_{||\mathbf{v}||_p=1} ||A\mathbf{v}||_p
|
||
$$
|
||
:::
|
||
|
||
Induzierte Normen lassen sich für allgemeines $p$ nur schwer berechnen. Ausnahmen sind die Fälle
|
||
|
||
- $p=1$: Spaltensummennorm
|
||
- $p=2$: Spektralnorm und
|
||
- $p=\infty$: Zeilensummennorm
|
||
|
||
Diese 3 Fälle sind in Julia in der Funktion `opnorm(A, p)` aus dem `LinearAlgebra`-Paket implementiert, wobei wieder `opnorm(A) = opnorm(A, 2)` gilt.
|
||
|
||
```{julia}
|
||
A = [ 0 1
|
||
1.2 1.5 ]
|
||
|
||
@show opnorm(A, 1) opnorm(A, Inf) opnorm(A, 2) opnorm(A);
|
||
```
|
||
|
||
Das folgende Bild zeigt die Wirkung von $A$ auf Einheitsvektoren. Vektoren gleicher Farbe werden aufeinander abgebildet. (Code durch anklicken sichtbar):
|
||
|
||
```{julia}
|
||
#| error: false
|
||
#| echo: false
|
||
#| output: false
|
||
|
||
using CairoMakie
|
||
CairoMakie.activate!(type = "svg")
|
||
# set_theme!(size=(1600, 360))
|
||
```
|
||
|
||
```{julia}
|
||
#| code-fold: true
|
||
#| fig-cap: "Bild der Einheitskugel unter $v \\mapsto Av$ mit $||A||\\approx 2.088$"
|
||
|
||
using CairoMakie
|
||
|
||
# Makie bug https://github.com/MakieOrg/Makie.jl/issues/3255
|
||
# Würgaround https://github.com/MakieOrg/Makie.jl/issues/2607#issuecomment-1385816645
|
||
tri = BezierPath([
|
||
MoveTo(Point2f(-0.5, -1)), LineTo(0, 0), LineTo(0.5, -1), ClosePath()
|
||
])
|
||
|
||
A = [ 0 1
|
||
1.2 1.5 ]
|
||
|
||
t = LinRange(0, 1, 100)
|
||
xs = sin.(2π*t)
|
||
ys = cos.(2π*t)
|
||
Axys = A * [xs, ys]
|
||
|
||
u = [sin(n*π/6) for n=0:11]
|
||
v = [cos(n*π/6) for n=0:11]
|
||
x = y = zeros(12)
|
||
|
||
Auv = A * [u,v]
|
||
|
||
fig2 = Figure(size=(800, 400))
|
||
lines(fig2[1, 1], xs, ys, color=t, linewidth=5, colormap=:hsv, axis=(; aspect = 1, limits=(-2,2, -2,2),
|
||
title=L"$\mathbf{v}$", titlesize=30))
|
||
arrows!(fig2[1,1], x, y, u, v, arrowsize=10, arrowhead=tri, colormap=:hsv, linecolor=range(0,11), linewidth=3)
|
||
|
||
Legend(fig2[1,2], MarkerElement[], String[], L"⟹", width=40, height=30, titlesize=30, framevisible=false)
|
||
|
||
lines(fig2[1,3], Axys[1], Axys[2], color=t, linewidth=5, colormap=:hsv, axis=(; aspect=1, limits=(-2,2, -2,2),
|
||
title=L"$A\mathbf{v}$", titlesize=30))
|
||
arrows!(fig2[1,3], x, y, Auv[1], Auv[2], arrowsize=10, arrowhead=tri, colormap=:hsv, linecolor=range(0,11),
|
||
linewidth=3)
|
||
|
||
fig2
|
||
```
|
||
|
||
### Konditionszahl
|
||
|
||
Für p = 1, p = 2 (default) oder p = Inf liefert `cond(A,p)` die Konditionszahl in der $p$-Norm
|
||
$$
|
||
\text{cond}_p(A) = ||A||_p \cdot ||A^{-1}||_p
|
||
$$
|
||
|
||
```{julia}
|
||
@show cond(A, 1) cond(A, 2) cond(A) cond(A, Inf);
|
||
```
|
||
|
||
## Matrixfaktorisierungen
|
||
|
||
Basisaufgaben der numerischen linearen Algebra:
|
||
|
||
- Löse ein lineares Gleichungssystem $A\mathbf{x} = \mathbf{b}$.
|
||
- Falls keine Lösung existiert, finde die beste Annäherung, d.h., den Vektor $\mathbf{x}$, der $||A\mathbf{x} - \mathbf{b}||$ minimiert.
|
||
- Finde Eigenwerte und Eigenvektoren $A\mathbf{x} = \lambda \mathbf{x}$ von $A$.
|
||
|
||
Diese Aufgaben sind mit Matrixfaktorisierungen lösbar. Einige grundlegende Matrixfaktorisierungen:
|
||
|
||
- **LU-Zerlegung** $A=L\cdot U$
|
||
- faktorisiert eine Matrix als Produkt einer _lower_ und einer _upper_ Dreiecksmatrix
|
||
- im Deutschen auch LR-Zerlegung (aber die Julia-Funktion heisst `lu()`)
|
||
- geht (eventuell nach Zeilenvertauschung - Pivoting) immer
|
||
- **Cholesky-Zerlegung** $A=L\cdot L^*$
|
||
- die obere Dreiecksmatrix ist die konjugierte der unteren,
|
||
- halber Aufwand im Vergleich zu LU
|
||
- geht nur, wenn $A$ hermitesch und positiv definit ist
|
||
- **QR-Zerlegung** $A=Q\cdot R$
|
||
- zerlegt $A$ als Produkt einer orthogonalen (bzw. unitären im komplexen Fall) Matrix und einer oberen Dreiecksmatrix
|
||
- $Q$ ist längenerhaltend (Drehungen und/oder Spiegelungen); die Stauchungen/Streckungen werden durch $R$ beschrieben
|
||
- geht immer
|
||
- **SVD** _(Singular value decomposition)_: $A = U\cdot D \cdot V^*$
|
||
- $U$ und $V$ sind orthogonal (bzw. unitär), $D$ ist eine Diagonalmatrix mit Einträgen in der Diagonale $σ_i\ge 0$, den sogenannten _Singulärwerten_ von $A$.
|
||
- Jede lineare Transformation $\mathbf{v} \mapsto A\mathbf{v}$ läßt sich somit darstellen als eine Drehung (und/oder Spiegelung) $V^*$, gefolgt von einer reinen Skalierung $v_i \mapsto \sigma_i v_i$ und einer weitere Drehung $U$.
|
||
|
||
|
||
### LU-Faktorisierung
|
||
|
||
LU-Faktorisierung ist Gauß-Elimination. Das Resultat der Gauß-Elimination ist die obere Dreiecksmatrix $U$. Die untere Dreiecksmatrix $L$ enthält Einsen auf der Diagonale und die nichtdiagonalen Einträge $l_{ij}$ sind gleich minus den Koeffizienten, mit denen im Gauß-Algorithmus Zeile $Z_j$ multipliziert und zu Zeile $Z_i$ addiert wird.
|
||
Ein Beispiel:
|
||
$$
|
||
A=\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
3 &-4& 4 \\
|
||
-2 & 1 & 5
|
||
\end{array}\right]
|
||
~ \begin{array}{c}
|
||
~\\
|
||
Z_2 \mapsto Z_2 \mathbin{\color{red}-}\textcolor{red}{3} Z_1\\
|
||
Z_3 \mapsto Z_3 + \textcolor{red}{2} Z_1
|
||
\end{array} \quad \Longrightarrow\
|
||
\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
&-10& -2 \\
|
||
& 5 & 9
|
||
\end{array}\right]
|
||
~ \begin{array}{c}
|
||
~\\
|
||
~\\
|
||
Z_3 \mapsto Z_3 + \textcolor{red}{\frac{1}{2}} Z_2
|
||
\end{array} \quad \Longrightarrow\
|
||
\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
&-10& -2 \\
|
||
& & 8
|
||
\end{array}\right] = U
|
||
$$
|
||
$$
|
||
A = \left[
|
||
\begin{array}{ccc}
|
||
1 && \\
|
||
\mathbin{\color{red}+}\textcolor{red}{3} &1 & \\
|
||
\mathbin{\color{red}-}\textcolor{red}{2} & \mathbin{\color{red}-}\textcolor{red}{\frac{1}{2}}& 1
|
||
\end{array}
|
||
\right] \cdot
|
||
\left[
|
||
\begin{array}{ccc}
|
||
1 &2 &2 \\
|
||
&-10& -2 \\
|
||
& & 8
|
||
\end{array}\right]
|
||
$$
|
||
|
||
|
||
- Häufig in der Praxis: $A\mathbf{x}=\mathbf{b}$ muss für ein $A$ und viele rechte Seiten $\mathbf{b}$ gelöst werden.
|
||
- Die Faktorisierung, deren Aufwand kubisch $\sim n^3$ mit der Matrixgröße $n$ wächst, muss nur einmal gemacht werden.
|
||
- Der anschliessende Aufwand der Vorwärts/Rückwärtssubstition für jedes $\mathbf{b}$ ist nur noch quadratisch $\sim n^2$.
|
||
|
||
Das `LinearAlgebra`-Paket von Julia enthält zur Berechnung einer LU-Zerlegung die Funktion `lu(A, options)`:
|
||
```{julia}
|
||
A = [ 1 2 2
|
||
3 -4 4
|
||
-2 1 5]
|
||
|
||
L, U, _ = lu(A, NoPivot())
|
||
display(L)
|
||
display(U)
|
||
```
|
||
|
||
|
||
|
||
#### Pivoting
|
||
|
||
Sehen wir uns einen Schritt der Gauß-Elimination an:
|
||
$$
|
||
\left[
|
||
\begin{array}{cccccc}
|
||
* &\cdots & * & * & \cdots & * \\
|
||
&\ddots & \vdots &\vdots && \vdots \\
|
||
&& * & * &\cdots & * \\
|
||
&& & \textcolor{red}{a_{ij}}&\cdots & a_{in}\\
|
||
&& & \textcolor{blue}{a_{i+1,j}}&\cdots & a_{i+1,n}\\
|
||
&&& \textcolor{blue}{\vdots} &&\vdots \\
|
||
&& & \textcolor{blue}{a_{mj}}&\cdots & a_{mn}
|
||
\end{array}
|
||
\right]
|
||
$$
|
||
Ziel ist es, als nächstes den Eintrag $a_{i+1,j}$ zum Verschwinden zu bringen, indem zur Zeile $Z_{i+1}$ ein geeignetes Vielfaches von Zeile $Z_i$ addiert wird. Das geht nur, wenn das _Pivotelement_ $\textcolor{red}{a_{ij}}$ nicht Null ist. Falls $\textcolor{red}{a_{ij}}=0$, müssen wir Zeilen vertauschen um dies zu beheben.
|
||
|
||
Darüber hinaus ist die Kondition des Algorithmus am besten, wenn man bei jedem Schritt die Matrix so anordnet, dass das Pivotelement das betragsmäßig größte
|
||
in der entsprechenden Spalte der noch zu bearbeitenden Umtermatrix ist. Beim (Zeilen-)Pivoting wird bei jedem Schritt durch Zeilenvertauschung sichergestellt, dass
|
||
gilt:
|
||
$$
|
||
|\textcolor{red}{a_{ij}}|=\max_{k=i,...,m} |\textcolor{blue}{a_{kj}}|
|
||
$$
|
||
|
||
#### LU in Julia
|
||
|
||
- Die Faktorisierungen in Julia geben ein spezielles Objekt zurück, das die Matrixfaktoren und weitere
|
||
Informationen enthält.
|
||
- Die Julia-Funktion `lu(A)` führt eine LU-Faktorisierung mit Pivoting durch.
|
||
|
||
|
||
```{julia}
|
||
F = lu(A)
|
||
typeof(F)
|
||
```
|
||
Elemente des Objekts:
|
||
```{julia}
|
||
@show F.L F.U F.p;
|
||
```
|
||
Man kann auch gleich auf der linken Seite ein entsprechendes Tupel verwenden:
|
||
```{julia}
|
||
L, U, p = lu(A);
|
||
p
|
||
```
|
||
|
||
Der Permutationsvektor zeigt an, wie die Zeilen der Matrix permutiert wurden. Es gilt: $$ L\cdot U = PA$$. Die Syntax der indirekten Indizierung erlaubt es, die Zeilenpermutation durch die Schreibweise `A[p,:]` anzuwenden:
|
||
```{julia}
|
||
display(A)
|
||
display(A[p,:])
|
||
display(L*U)
|
||
```
|
||
|
||
|
||
Die Vorwärts/Rückwärtssubstition mit einem `LU`- erledigt der Operator `\`:
|
||
```{julia}
|
||
b = [1, 2, 3]
|
||
x = F \ b
|
||
```
|
||
Probe:
|
||
```{julia}
|
||
A * x - b
|
||
```
|
||
|
||
In Julia verbirgt sich hinter dem `\`-Operator ein ziemlich universeller "matrix solver" und man kann ihn auch direkt anwenden:
|
||
```{julia}
|
||
A \ b
|
||
```
|
||
Dabei wird implizit eine geeignete Faktorisierung durchgeführt, deren Ergebnis allerdings nicht abgespeichert.
|
||
|
||
|
||
### QR-Zerlegung
|
||
|
||
|
||
Die Funktion `qr()` liefert ein epezielles QR-Objekt zurück, das die Komponenten $Q$ und $R$ enthält. Die orthogonale (bzw. unitäre) Matrix $Q$ ist
|
||
[in einer optimierten Form](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#man-linalg-abstractq) abgespeichert. Umwandlung in eine "normale" Matrix ist bei Bedarf wie immer mit `collect()` möglich.
|
||
```{julia}
|
||
F = qr(A)
|
||
@show typeof(F) typeof(F.Q)
|
||
display(collect(F.Q))
|
||
display(F.R)
|
||
```
|
||
|
||
### Passende Faktorisierung
|
||
|
||
|
||
Die Funktion `factorize()` liefert eine dem Matrixtyp angepasste Form der Faktorisierung, siehe [Dokumentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.factorize) für Details.
|
||
Wenn man Lösungen zu mehreren rechten Seiten $\mathbf{b_1}, \mathbf{b_2},...$ benötigt, sollte man die Faktorisierung nur einmal durchführen:
|
||
|
||
|
||
```{julia}
|
||
Af = factorize(A)
|
||
```
|
||
|
||
|
||
```{julia}
|
||
Af \ [1, 2, 3]
|
||
```
|
||
|
||
|
||
```{julia}
|
||
Af \ [5, 7, 9]
|
||
```
|
||
|