---
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]
```