Grundlagen der linearen Algebra in

Über Grundlagen der linearen Algebra

Was ist „Lineare Algebra"?

Es gibt jede Menge technischer Definitionen, auf Wikipedia, in Lehrbüchern und auf hervorragenden Websites wie 3Blue1Brown, aber für unsere Zwecke können wir es etwas lockerer angehen.

Die Maintainer von Julia (und Python) lieben solche Dinge. Die Verantwortlichen bei Exercism eher weniger.

Dieses Konzept beschränkt sich auf die einfacheren Aspekte eines riesigen Themas, aber sei gewarnt: Das hier ist zwangsläufig ein ziemlich mathematisches Konzept.

Was ist ein Vektor?

Das hängt davon ab, wen du fragst und wie du etwas ziemlich Abstraktes veranschaulichen willst.

• In Julia ist Vector{T} nur ein Alias für Array{T, 1}, und der eltype T kann alles Mögliche sein.
• In (weiten Teilen) der Physik ist ein Vektor ein Pfeil mit length und direction im N-dimensionalen Raum, aber ohne feste Position.
• In der linearen Algebra ist der Pfeil der Physiker im Raum fixiert, mit seinem Fuß im Ursprung. Dadurch wird der Schaft des Pfeils überflüssig, und wir können den Vektor als point im N-dimensionalen Raum darstellen, wobei die N Elemente den Abstand vom Ursprung entlang jeder der N Achsen angeben.

Für unsere Zwecke ignorieren wir Vektoren von Strings oder Zeichen. In diesem Konzept arbeiten wir mit numerischen Typen: Int, Float oder Complex.

Was ist eine Matrix?

Auch hier wirst du eine Vielzahl von Antworten finden.

• Ein rechteckiges 2-D-Array von Zahlen (ohne ausgefranste Ränder). Eine quadratische Matrix ist ein häufiger Sonderfall.
• Eine linear combination von Spaltenvektoren, nebeneinander gestapelt.
• Eine transformation, die auf einen Vektor angewendet werden kann, so wie eine function auf andere Eingabetypen angewendet wird.

Einige Arten von quadratischen Matrizen sind so verbreitet, dass sie eigene Namen haben.

Diagonalmatrix: Alle von Null verschiedenen Einträge liegen auf der main diagonal (oben links nach unten rechts).

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

Einheitsmatrix: Eine Diagonalmatrix, bei der nur Einsen auf der Diagonalen stehen (aus Gründen, die später klarer werden). Oft als I abgekürzt.

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

Obere Dreiecksmatrix: Von Null verschiedene Werte auf und über der Diagonalen, Nullen darunter. Untere Dreiecksmatrix kannst du dir denken.

Eine Matrix transponieren

Wenn du wirklich Zeilen mit Spalten vertauschen willst, erledigt das die Funktion permutedims() und liefert dir eine neue Matrix.

Das erfordert eine Kopie, was bei großen Matrizen langsam ist und viel Speicher braucht.

Für Zwecke der linearen Algebra ist die Funktion transpose() nützlicher, denn sie erzeugt schnell einen Lazy-Wrapper um die ursprüngliche Matrix.

Noch nützlicher ist die Funktion adjoint(), die zusätzlich das Vorzeichen des Imaginärteils in komplexen Zahlen umkehrt (wenn du dich fragst, warum das so wichtig ist, lautet eine schnelle Antwort: Quantenmechanik). Diese Operation ist so verbreitet, dass wir einfach einen Apostroph ' an den Variablennamen hängen können, um die Adjungierte zu erzeugen.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

Multiplikation

Der Rest dieses Dokuments nutzt Funktionalität aus dem Modul LinearAlgebra. Eine Zeile using LinearAlgebra holt es in den Namensraum, aber wir werden das in den Beispielen nicht ständig wiederholen (zu viel visueller Ballast).

Vektoren, elementweises Produkt

Das wurde im Konzept Vektoroperationen besprochen. Der Operator ist .*, der paarweise auf die Eingabevektoren wirkt und eine Ausgabe derselben Größe und desselben Typs wie die Eingaben liefert.

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

Vektoren, Skalarprodukt

Diese extrem verbreitete Operation wird in Lehrbüchern als u ⋅ v geschrieben und „Skalarprodukt" genannt.

Ein Skalarprodukt entspricht der Summe des elementweisen Produkts.

Zwei Vektoren lassen sich mit dem üblichen *-Operator multiplizieren, aber nur, wenn der linke Vektor eine adjoint ist: praktischerweise geschrieben als u' * v. Das dürfte im späteren Abschnitt über Matrixmultiplikation klarer werden.

Der hochgestellte (zentrierte) Punkt ist in Julia verfügbar (eingegeben mit \cdot und Tab) als syntaktischer Zucker für die Funktion dot(). Die Adjungierte musst du nicht angeben, dieses Detail wird automatisch behandelt.

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

Bei komplexwertigen Vektoren muss der linke Vektor das Konjugierte sein (Vorzeichen des Imaginärteils umgekehrt). Die Funktion dot macht das automatisch, die Syntax u' macht es explizit, aber sum(u .* v) schlägt fehl, wenn u und v komplex sind.

Vektoren, Kreuzprodukt

Ein bisschen Nostalgie für alle, die früher mal eine Vorlesung über Elektrizität und Magnetismus gehört haben! Und für Ingenieure, die mit der Berechnung von Drehmoment- oder Drehimpulsvektoren vertraut sind.

Während das Skalarprodukt zwei Vektoren in einen scalar umwandelt, macht das Kreuzprodukt aus zwei 3-Vektoren einen dritten 3-Vektor.

In der geometrischen Darstellung des Vektorraums steht der neue Vektor senkrecht auf der Ebene, die die beiden Eingabevektoren enthält. Wenn die beiden Eingaben parallel sind (bis aufs Vorzeichen), definieren sie keine Ebene, also ist die Ausgabe [0, 0, 0].

Die Reihenfolge ist wichtig: u × v == -(v × u). Das ist die berühmte Rechte-Hand-Regel, bei der schon viele von uns auf Daumen und zwei Finger gestarrt haben, während sie sie im Raum verdrehten (oft mit ratlosem Gesichtsausdruck).

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

Beachte, dass diese Operation auf Vektoren der Länge 3 beschränkt ist (äquivalent zum euklidischen Raum mit orthogonalen x, y, z-Achsen). Mathematiker definieren die Operation gerne auch für andere Dimensionen, mit einem unverständlichen Wortsalat (Tensoren höherer Ordnung, Keilprodukt, äußeres Produkt, Multivektoren ...) als Ergebnis. Die meisten von uns rennen an dieser Stelle weg oder wechseln schnell das Thema!

Normen

Wie „groß" ist ein Vektor?

Die norm ist der Versuch, das einzufangen, indem der Vektor auf einen passenden Skalar reduziert wird.

Es ist eine ganze Familie von Normen definiert, aber bei Weitem am gebräuchlichsten ist die 2-Norm, die √(v ⋅ v) ist.

Diese Wurzel-Mittelwert-Operation liefert den pythagoreischen Abstand vom Ursprung (im N-dimensionalen Raum). Wenn wir den Vektor als Pfeil mit seinem Fuß im Ursprung veranschaulichen, ist die 2-Norm die Länge des Pfeils.

Jede p-Norm lässt sich berechnen, indem du p als zweites Argument angibst. Die 1-Norm ist manchmal nützlich: Sie ist einfach die Summe der Beträge der Einträge (also sehr schnell und einfach zu berechnen).

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

Matrixmultiplikation

Manchmal scheint es, als hätte jeder angewandte Mathematiker auf der Welt einen Großteil der letzten 80 Jahre damit verbracht, jede Berechnung in eine Reihe von Matrixmultiplikationen zu verwandeln.

Das ist eine Operation, in der Computer sehr gut sind:

• Sie ist stark repetitiv.
• Sie lässt sich effizient parallelisieren.
• Es wurde viel spezialisierte Hardware entwickelt, um sie schneller zu machen, vom 80 Millionen Dollar teuren Cray-1 in den 1970ern bis zur GPU, die wahrscheinlich in deinem Laptop steckt.

Die Details sind ziemlich einfach, wenn auch auf den ersten Blick nicht sehr intuitiv.

Betrachte eine Matrix A, die mit einem Vektor v multipliziert wird und eine Ausgabe w ergibt (nach Konvention der linearen Algebra verwenden wir Großbuchstaben für Matrizen, Kleinbuchstaben für Vektoren).

Die oberste Zeile von A wird mit v skalar multipliziert, um das erste Element von w zu erhalten, die zweite Zeile ergibt das zweite Element, und so weiter nach unten.

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

Eine Matrix-mal-Matrix-Multiplikation erstreckt dies auf die Spalten der Matrix auf der rechten Seite.

Für C = A * B können wir uns vorstellen, dass A jede Spalte von B multipliziert und die entsprechende Spalte von C ergibt: eine Reihe von Matrix-Vektor-Multiplikationen.

Gleichwertig könnten wir sagen, dass die erste Zeile von A mit jeder Spalte von B skalar multipliziert wird und die oberste Zeile von C ergibt, die zweite Zeile ergibt die zweite Zeile, und so weiter nach unten.

Es gibt mehrere solcher Denkmodelle, und interessierte Lernende können sich eine ganze MIT-Vorlesung ansehen, die sie bespricht.

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53

julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

Wie im obigen Beispiel gezeigt, ist die Matrixmultiplikation nicht kommutativ: Es gibt keine einfache Beziehung zwischen A*B und B*A.

Matrixmultiplikation ist wahrscheinlich für jeden, der neu darin ist, schwer zu veranschaulichen, wenn man nur die Worte liest. Auf YouTube gibt es jede Menge Videos, die sie grafisch demonstrieren. Suche also nach „Matrixmultiplikation" und wähle eines mit deinem bevorzugten Stil, Detaillierungsgrad und deiner Sprache.

Dimensionen

Das Skalarprodukt zweier Vektoren setzt voraus, dass sie gleich lang sind.

Analog dazu muss bei der Matrixmultiplikation die Anzahl der Spalten der linken Matrix mit der Anzahl der Zeilen der rechten Matrix übereinstimmen.

Drücken wir die Größen als Tupel (nrows, ncols) aus, wie sie size(A) ausgibt, dann gilt (a, b) * (b, c) -> (a, c). Die „inneren" Dimensionen, hier b und b, sind mit der Bildung von Skalarprodukten kompatibel. Die „äußeren" Dimensionen, hier a und c, bestimmen die Dimensionen der Ausgabe.

Ein Beispiel mit rechteckigen Matrizen:

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

Vektorpaare auf Matrix-Art zu multiplizieren hat zwei Möglichkeiten.

Konventionell würden wir u' * v als Äquivalent zum Skalarprodukt verwenden. Die Dimensionen sind (1, 3) * (3, 1), und Julia vereinfacht die Ausgabe (1, 1) zu einem Skalar (anders als zum Beispiel R).

Alternativ könnten wir u * v' mit den Dimensionen (3, 1) * (1, 3) -> (3, 3) verwenden, um alle möglichen Elementpaare zu multiplizieren und die Vektoren zu einer Matrix zu erweitern.

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * v wird manchmal das innere Produkt genannt.

Das (seltenere) u * v' ist entsprechend das äußere Produkt. Das hängt mit dem Tensorprodukt zusammen, das aber weit über unseren Rahmen hinausgeht.

Rotationen

Eine besonders häufige Art von Matrixmultiplikation betrifft Rotationsmatrizen.

In 2D gibt es eine relativ einfache Matrix, um einen Vektor um θ Radiant gegen den Uhrzeigersinn zu drehen.

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0

Rotationen in 3D brauchen eine komplexere Matrix mit zwei Winkeln. Das lässt sich in einem Markdown-Dokument schwer darstellen, ohne viel LaTeX einzubetten, also schau für die Formel bei Wikipedia nach.

In der 3-D-Grafik, etwa bei OpenGL und seinen Nachfolgern, arbeitet man konventionell mit homogenen Koordinaten.

Jeder Scheitelpunkt wird durch einen 4-Vektor dargestellt (oder gleichwertig eine Matrixspalte): [x, y, z, 1.0] für einen Punkt bei (x, y, z).

Das ermöglicht aufwendigere Transformationsmatrizen. Die Rotation steckt weiterhin in A[1:3, 1:3], Translationen sind [Δx, Δy, Δz] in A[1:3, 4], die Skalierung liegt auf der Diagonalen, und für Scherung, Perspektive usw. gibt es weitere Möglichkeiten.

Performance

An Informatik-Fakultäten auf der ganzen Welt gibt es viele Doktorarbeiten mit einem oder mehreren Kapiteln über kleine, schrittweise Verbesserungen von Algorithmen für die Matrixmultiplikation. Das ist wirklich wichtig, und verschiedene Organisationen sind bereit, die Forschung zu ihrem eigenen Nutzen zu finanzieren.

Performance ist ein großes Thema, in das wir nicht tief einsteigen können, aber zur Veranschaulichung können wir versuchen, Zufallsmatrizen verschiedener Größen zu multiplizieren.

Der folgende Code ist eine schnelle, grobe Schätzung (nutze BenchmarkTools.jl für einen besseren Ansatz).

Das verwendete System war ein kleiner PC für 430 US-Dollar: Ryzen-9-Prozessor, 32 GB RAM, Linux Mint 22.1, Julia 1.11.6.

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

Grob gesagt brauchten zwei Matrizen mit einer Million Elementen ein paar Millisekunden, Matrizen mit 100 Millionen Elementen ein paar Sekunden. Ein paar weitere Größenordnungen draufzulegen erfordert bessere Hardware ...