13  Lineare Algebra in Julia

using LinearAlgebra

Das LinearAlgebra-Paket liefert unter anderem:

13.1 Matrixtypen

A = SymTridiagonal(fill(1.0, 4), fill(-0.3, 3))
4×4 SymTridiagonal{Float64, Vector{Float64}}:
  1.0  -0.3    ⋅     ⋅ 
 -0.3   1.0  -0.3    ⋅ 
   ⋅   -0.3   1.0  -0.3
   ⋅     ⋅   -0.3   1.0
B = UpperTriangular(A)
4×4 UpperTriangular{Float64, SymTridiagonal{Float64, Vector{Float64}}}:
 1.0  -0.3   0.0   0.0
  ⋅    1.0  -0.3   0.0
  ⋅     ⋅    1.0  -0.3
  ⋅     ⋅     ⋅    1.0

Diese Typen werden platzsparend gespeichert. Die üblichen Rechenoperationen sind implementiert:

A + B
4×4 Matrix{Float64}:
  2.0  -0.6   0.0   0.0
 -0.3   2.0  -0.6   0.0
  0.0  -0.3   2.0  -0.6
  0.0   0.0  -0.3   2.0

Lesende Indexzugriffe sind möglich,

A[1,4]
0.0

schreibende nicht unbedingt:

A[1,3] = 17
ArgumentError: cannot set off-diagonal entry (1, 3)
Stacktrace:
 [1] setindex!(A::SymTridiagonal{Float64, Vector{Float64}}, x::Int64, i::Int64, j::Int64)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.3+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:466
 [2] top-level scope
   @ ~/Julia/23/Book-ansipatch/chapters/11_LinAlg.qmd:74

Die Umwandlung in eine ‘normale’ Matrix ist z.B. mit collect() möglich:

A2 = collect(A)
4×4 Matrix{Float64}:
  1.0  -0.3   0.0   0.0
 -0.3   1.0  -0.3   0.0
  0.0  -0.3   1.0  -0.3
  0.0   0.0  -0.3   1.0

13.1.1 Die Einheitsmatrix I

I bezeichnet eine Einheitsmatrix (quadratisch, Diagonalelemente = 1, alle anderen = 0) in der jeweils erforderlichen Größe

A + 4I
4×4 SymTridiagonal{Float64, Vector{Float64}}:
  5.0  -0.3    ⋅     ⋅ 
 -0.3   5.0  -0.3    ⋅ 
   ⋅   -0.3   5.0  -0.3
   ⋅     ⋅   -0.3   5.0

13.2 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|| \]

13.2.1 \(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.

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).

v = [3, 4]
w = [-1, 2, 33.2]

@show norm(v)  norm(v, 2)  norm(v, 1) norm(v, 4) norm(w, Inf);
norm(v) = 5.0
norm(v, 2) = 5.0
norm(v, 1) = 7.0
norm(v, 4) = 4.284572294953817
norm(w, Inf) = 33.2
  • 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.
A = [1 2 3
     4 5 6
     7 8 9]
norm(A)      # Frobenius norm
16.881943016134134

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).

Code 
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

Einheitskugeln im \(ℝ^2\) für verschiedene \(p\)-Normen: \(p\)=0.8; 1; 1.5; 2; 3.001 und 1000

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.

13.2.2 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.

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.

A = [ 0    1
      1.2  1.5 ]

@show opnorm(A, 1) opnorm(A, Inf) opnorm(A, 2) opnorm(A);
opnorm(A, 1) = 2.5
opnorm(A, Inf) = 2.7
opnorm(A, 2) = 2.0879899930905124
opnorm(A) = 2.0879899930905124

Das folgende Bild zeigt die Wirkung von \(A\) auf Einheitsvektoren. Vektoren gleicher Farbe werden aufeinander abgebildet. (Code durch anklicken sichtbar):

Code 
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

Bild der Einheitskugel unter \(v \mapsto Av\) mit \(||A||\approx 2.088\)

13.2.3 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 \]

@show cond(A, 1)  cond(A, 2) cond(A)  cond(A, Inf);
cond(A, 1) = 5.625
cond(A, 2) = 3.633085176038432
cond(A) = 3.633085176038432
cond(A, Inf) = 5.625000000000001

13.3 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\).

13.3.1 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):

A = [ 1   2  2
      3  -4  4
     -2   1  5]

L, U, _ = lu(A, NoPivot())
display(L)
display(U)
3×3 Matrix{Float64}:
  1.0   0.0  0.0
  3.0   1.0  0.0
 -2.0  -0.5  1.0
3×3 Matrix{Float64}:
 1.0    2.0   2.0
 0.0  -10.0  -2.0
 0.0    0.0   8.0

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.
F = lu(A)
typeof(F)
LU{Float64, Matrix{Float64}, Vector{Int64}}

Elemente des Objekts:

@show F.L F.U F.p;
F.L = [1.0 0.0 0.0; 0.3333333333333333 1.0 0.0; -0.6666666666666666 -0.5 1.0]
F.U = [3.0 -4.0 4.0; 0.0 3.333333333333333 0.6666666666666667; 0.0 0.0 8.0]
F.p = [2, 1, 3]

Man kann auch gleich auf der linken Seite ein entsprechendes Tupel verwenden:

L, U, p = lu(A);
p
3-element Vector{Int64}:
 2
 1
 3

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:

display(A)
display(A[p,:])
display(L*U)
3×3 Matrix{Int64}:
  1   2  2
  3  -4  4
 -2   1  5
3×3 Matrix{Int64}:
  3  -4  4
  1   2  2
 -2   1  5
3×3 Matrix{Float64}:
  3.0  -4.0  4.0
  1.0   2.0  2.0
 -2.0   1.0  5.0

Die Vorwärts/Rückwärtssubstition mit einem LU- erledigt der Operator \:

b = [1, 2, 3]
x = F \ b
3-element Vector{Float64}:
 -0.10000000000000002
 -0.012500000000000006
  0.5625

Probe:

A * x - b
3-element Vector{Float64}:
 0.0
 0.0
 0.0

In Julia verbirgt sich hinter dem \-Operator ein ziemlich universeller “matrix solver” und man kann ihn auch direkt anwenden:

A \ b
3-element Vector{Float64}:
 -0.10000000000000002
 -0.012500000000000006
  0.5625

Dabei wird implizit eine geeignete Faktorisierung durchgeführt, deren Ergebnis allerdings nicht abgespeichert.

13.3.2 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 abgespeichert. Umwandlung in eine “normale” Matrix ist bei Bedarf wie immer mit collect() möglich.

F = qr(A)
@show typeof(F) typeof(F.Q)
display(collect(F.Q))
display(F.R)
typeof(F) = LinearAlgebra.QRCompactWY{Float64, Matrix{Float64}, Matrix{Float64}}
typeof(F.Q) = LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}
3×3 Matrix{Float64}:
 -0.267261   0.872872  0.408248
 -0.801784  -0.436436  0.408248
  0.534522  -0.218218  0.816497
3×3 Matrix{Float64}:
 -3.74166  3.20713  -1.06904
  0.0      3.27327  -1.09109
  0.0      0.0       6.53197

13.3.3 Passende Faktorisierung

Die Funktion factorize() liefert eine dem Matrixtyp angepasste Form der Faktorisierung, siehe Dokumentation 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:

Af = factorize(A)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
  1.0        0.0  0.0
  0.333333   1.0  0.0
 -0.666667  -0.5  1.0
U factor:
3×3 Matrix{Float64}:
 3.0  -4.0      4.0
 0.0   3.33333  0.666667
 0.0   0.0      8.0
Af \ [1, 2, 3]
3-element Vector{Float64}:
 -0.10000000000000002
 -0.012500000000000006
  0.5625
Af \ [5, 7, 9]
3-element Vector{Float64}:
 0.4000000000000001
 0.42500000000000004
 1.875