diff --git a/chapters/11_LinAlg.qmd b/chapters/11_LinAlg.qmd index b160e59..5b604e1 100644 --- a/chapters/11_LinAlg.qmd +++ b/chapters/11_LinAlg.qmd @@ -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 @@ -15,6 +93,11 @@ $$ $$ 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}}| $$ @@ -24,13 +107,12 @@ 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} -using LinearAlgebra - v = [3, 4] w = [-1, 2, 33.2] @@ -52,10 +134,10 @@ 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. - +$||\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: 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" using Plots @@ -107,10 +189,10 @@ A = [ 0 1 @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} #| 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 @@ -123,7 +205,6 @@ tri = BezierPath([ A = [ 0 1 1.2 1.5 ] -colorlist = resample_cmap(:hsv, 12) t = LinRange(0, 1, 100) xs = sin.(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 ``` -## Matrixfaktorisierungen +### Konditionszahl -Eine Basisaufgabe der numerischen linearen Algebra ist die Lösung von linearen Gleichungssystemen +Für p = 1, p = 2 (default) oder p = Inf liefert `cond(A,p)` die Konditionszahl in der $p$-Norm $$ -A\mathbf{x} = \mathbf{b}. -$$ -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. +\text{cond}_p(A) = ||A||_p \cdot ||A^{-1}||_p $$ +```{julia} +@show cond(A, 1) cond(A, 2) cond(A) cond(A, Inf); +``` +## Matrixfaktorisierungen -# Das _Linear Algebra_-Paket: eine kurze Auswahl +Basisaufgaben der numerischen linearen Algebra: -- zusätzliche Matrix-Typen +- 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$. - - Subtypen von `AbstractMatrix`: genauso verwendbar, wie andere Matrizen - - `Tridiagonal` - - `SymTridiagonal` - - `Symmetric` - - `UpperTriangular` - - ... +Diese Aufgaben sind mit Matrixfaktorisierungen lösbar. Einige grundlegende Matrixfaktorisierungen: - -- Arithmetik: Matrixmultiplikation, `inv`, `det`, `exp` - -- Lineare Gleichungssysteme: `\` - -- Matrixfaktorisierungen - - - `LU` - - `QR` - - `Cholesky` - - `SVD` - - ... +- **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$. -- Eigenwerte/-vektoren - - - `eigen`, `eigvals`, `eigvecs` - - -- Zugriff auf BLAS/LAPACK-Funktionen -## Matrixtypen +### 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$. -```julia -using LinearAlgebra +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) ``` -```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) ``` - - -```julia -B = UpperTriangular(A) +Elemente des Objekts: +```{julia} +@show F.L F.U F.p; ``` - - -```julia -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 +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) ``` -Der Permutationsvektor `p` zeigt an, wie die Zeilen der Matrix permutiert wurden: - -```julia -A[p, :] # 3. Zeile, 1. Zeile, 2. Zeile von A +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 ``` -Es gilt: $$ L\cdot U = PA$$ - - -```julia -l * u - A[p,:] +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 -```julia -q, r = qr(A); - -q +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) ``` - -```julia -r -``` - -### Singular value decomposition +### Passende Faktorisierung -```julia -u, s, vt = svd(A); - -s -``` +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 -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 +```{julia} Af = factorize(A) ``` -```julia -Af\b +```{julia} +Af \ [1, 2, 3] ``` -```julia -Af\[5,7,9] +```{julia} +Af \ [5, 7, 9] ``` - -```julia - -```