hellmund/JuliaKurs23
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Activity

LinearAlg extended

Browse Source
This commit is contained in:
Meik Hellmund 2024-04-22 14:51:37 +02:00
parent 81c688bc21
commit 4240d56455
1 changed files with 255 additions and 179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

442
chapters/11_LinAlg.qmd
View File

@ -1,4 +1,82 @@
# Lineare Algebra # 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 ## Normen
@ -15,6 +93,11 @@ $$
$$ $$
die die die euklidische Norm $p=2$ verallgemeinern. 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 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}}| $$ |x_{\text{max}}| \le ||\mathbf{x}||_p \le n^\frac{1}{p} |x_{\text{max}}|
$$ $$
@ -24,13 +107,12 @@ Damit folgt
$$ $$
\lim_{p \rightarrow \infty} ||\mathbf{x}||_p = |x_{\text{max}}| =: ||\mathbf{x}||_\infty. \lim_{p \rightarrow \infty} ||\mathbf{x}||_p = |x_{\text{max}}| =: ||\mathbf{x}||_\infty.
$$ $$
:::
In Julia definiert das `LinearAlgebra`-Paket eine Funktion `norm(v, p)`. In Julia definiert das `LinearAlgebra`-Paket eine Funktion `norm(v, p)`.
```{julia} ```{julia}
using LinearAlgebra
v = [3, 4] v = [3, 4]
w = [-1, 2, 33.2] w = [-1, 2, 33.2]
@ -52,10 +134,10 @@ norm(A) # Frobenius norm
Da Normen homogen unter Multiplikation mit Skalaren sind, 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. $||\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} ```{julia}
#| code-fold: false #| code-fold: true
#| fig-cap: "Einheitskugeln im $ℝ^2$ für verschiedene $p$-Normen: $p$=0.8; 1; 1.5; 2; 3.001 und 1000" #| fig-cap: "Einheitskugeln im $ℝ^2$ für verschiedene $p$-Normen: $p$=0.8; 1; 1.5; 2; 3.001 und 1000"
using Plots using Plots
@ -107,10 +189,10 @@ A = [ 0 1
@show opnorm(A, 1) opnorm(A, Inf) opnorm(A, 2) opnorm(A); @show opnorm(A, 1) opnorm(A, Inf) opnorm(A, 2) opnorm(A);
``` ```
Das folgende Bild zeigt die Wirkung von $A$ auf Einheitsvektoren (Code durch anklicken sichtbar): Das folgende Bild zeigt die Wirkung von $A$ auf Einheitsvektoren. Vektoren gleicher Farbe werden aufeinander abgebildet. (Code durch anklicken sichtbar):
```{julia} ```{julia}
#| code-fold: true #| code-fold: true
#| fig-cap: "Bild der Einheitskugel unter $v \\mapsto Av$" #| fig-cap: "Bild der Einheitskugel unter $v \\mapsto Av$ mit $||A||\\approx 2.088$"
using CairoMakie using CairoMakie
@ -123,7 +205,6 @@ tri = BezierPath([
A = [ 0 1 A = [ 0 1
1.2 1.5 ] 1.2 1.5 ]
colorlist = resample_cmap(:hsv, 12)
t = LinRange(0, 1, 100) t = LinRange(0, 1, 100)
xs = sin.(2π*t) xs = sin.(2π*t)
ys = cos.(2π*t) ys = cos.(2π*t)
@ -150,218 +231,213 @@ arrows!(fig2[1,3], x, y, Auv[1], Auv[2], arrowsize=10, arrowhead=tri, colormap=:
fig2 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 ## Matrixfaktorisierungen
Eine Basisaufgabe der numerischen linearen Algebra ist die Lösung von linearen Gleichungssystemen 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\mathbf{x} = \mathbf{b}. 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
$$ $$
Falls keine Lösung existiert, ist man oft an der bestmöglichen Annäherung an eine Lösung interessiert, d.h., gesucht wird der Vektor $\mathbf{x}$, für den
$$ $$
||A\mathbf{x} - \mathbf{b}|| \rightarrow \min. 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 _Linear Algebra_-Paket: eine kurze Auswahl 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]
- zusätzliche Matrix-Typen L, U, _ = lu(A, NoPivot())
display(L)
- Subtypen von `AbstractMatrix`: genauso verwendbar, wie andere Matrizen display(U)
- `Tridiagonal`
- `SymTridiagonal`
- `Symmetric`
- `UpperTriangular`
- ...
- Arithmetik: Matrixmultiplikation, `inv`, `det`, `exp`
- Lineare Gleichungssysteme: `\`
- Matrixfaktorisierungen
- `LU`
- `QR`
- `Cholesky`
- `SVD`
- ...
- Eigenwerte/-vektoren
- `eigen`, `eigvals`, `eigvecs`
- Zugriff auf BLAS/LAPACK-Funktionen
## Matrixtypen
```julia
using LinearAlgebra
``` ```
```julia
A = SymTridiagonal(fill(1.0, 4), fill(-0.3, 3)) #### 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}
```julia @show F.L F.U F.p;
B = UpperTriangular(A)
``` ```
Man kann auch gleich auf der linken Seite ein entsprechendes Tupel verwenden:
```{julia}
```julia L, U, p = lu(A);
A + B
```
### Einheitsmatrix `I`
`I` bezeichnet eine Einheitsmatrix (quadratisch, Diagonalelemente = 1, alle anderen = 0) in der jeweils erforderlichen Größe
```julia
A + 4I
```
## Faktorisierungen
### LU-Faktorisierung mit Zeilenpivoting
('Lower/Upper triangular matrix', im Deutschen auch oft "LR-Zerlegung" für "Linke/Rechte Dreiecksmatrix")
```julia
A = [0 22 1.
-1 2 3
77 18 19]
```
```julia
# Faktorisierungen geben eine spezielle Struktur zurück, die die Matrixfaktoren und weitere
# Informationen enthalten:
Af = lu(A);
@show Af.L Af.U Af.p;
```
```julia
# man kann auch gleich auf der linken Seite ein entsprechendes Tupel verwenden
l,u,p = lu(A)
l
```
```julia
u
```
```julia
p p
``` ```
Der Permutationsvektor `p` zeigt an, wie die Zeilen der Matrix permutiert wurden: 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)
```julia display(A[p,:])
A[p, :] # 3. Zeile, 1. Zeile, 2. Zeile von A display(L*U)
``` ```
Es gilt: $$ L\cdot U = PA$$
Die Vorwärts/Rückwärtssubstition mit einem `LU`- erledigt der Operator `\`:
```julia ```{julia}
l * u - A[p,:] 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 ### QR-Zerlegung
```julia 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
q, r = qr(A); [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}
q F = qr(A)
@show typeof(F) typeof(F.Q)
display(collect(F.Q))
display(F.R)
``` ```
### Passende Faktorisierung
```julia
r
```
### Singular value decomposition
```julia 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.
u, s, vt = svd(A); Wenn man Lösungen zu mehreren rechten Seiten $\mathbf{b_1}, \mathbf{b_2},...$ benötigt, sollte man die Faktorisierung nur einmal durchführen:
s
```
```julia ```{julia}
u
```
```julia
vt
```
```julia
evalues, evectors = eigen(A)
evalues
```
```julia
evectors
```
### Lineare Gleichungssysteme
Das lineare Gleichungssystem
$$Ax = b$$
kann man in Julia einfach lösen mit
```
x = A\b
```
```julia
b = [2,3,4]
A\b
```
Dabei wird eine geeignete Matrixfaktorisierung vorgenommen (meist LU). Wenn man Lösungen zu mehreren rechten Seiten $b_1, b_2,...$ benötigt, sollte man die Faktorisierung nur einmal durchführen:
```julia
Af = factorize(A) Af = factorize(A)
``` ```
```julia ```{julia}
Af\b Af \ [1, 2, 3]
``` ```
```julia ```{julia}
Af\[5,7,9] Af \ [5, 7, 9]
``` ```
```julia
```