Lineare Gleichungen lösen in

Über Lineare Gleichungen lösen

Erinnerst du dich noch daran, wie du zum ersten Mal in der Schule Algebra gelernt hast?

Normalerweise bekamen wir ein Gleichungssystem vorgesetzt und sollten es nach x und y auflösen.

4x - 3y = -2
2x + 7y = 16

Zwei Gleichungen, zwei Unbekannte: ein wenig Einsetzen, und du siehst bald, dass x = 1, y = 2 gilt.

Schau dir noch einmal die linke Seite des vorherigen Beispiels an. Sie sieht aus wie eine Matrix [4 -3; 2 7], die einen unbekannten Vektor [x, y] multipliziert und den Vektor [-2, 16] ergibt.

In der Notation der linearen Algebra: A x = b, nach der üblichen Konvention, Großbuchstaben für Matrizen und Kleinbuchstaben für Vektoren zu verwenden.

A x = 0 und der Nullraum

Ein wichtiger Spezialfall ist, wenn alle unsere Gleichungen auf der rechten Seite null haben.

4x - 3y = 0
2x + 7y = 0

Im obigen Beispiel ist die einzige Lösung x = y = 0, was normalerweise nicht interessant ist.

Nichttriviale Lösungen für A x = 0 gibt es nur, wenn A eine Determinante von null hat.

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

Betrachte stattdessen diese Gleichungen:

4x - 3y = 0
8x + 6y = 0

Die zweite Gleichung ist nur das Doppelte der ersten und fügt keine neue Information hinzu. Jegliche Werte, bei denen x = 0.75y gilt, sind eine Lösung.

Diese Werte (x, y) liegen auf einer Geraden im zweidimensionalen Raum. Es gibt unendlich viele Lösungen, und sie bilden den Nullraum der Matrix.

Die Zeilen der Matrix sind in diesem Fall nicht linear unabhängig, und der Rang der Matrix ist kleiner als die Anzahl der Zeilen oder Spalten.

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

Die Funktion nullspace() gibt einen einzelnen Punkt zurück, der so gewählt ist, dass er zu einem Einheitsvektor normiert ist, aber jedes skalare Vielfache dieses Vektors liegt ebenfalls im Nullraum.

A x = b lösen

Nimm nun an, du hast 1000 Gleichungen mit 1000 Unbekannten? Wenn das albern klingt, bedenke, dass drahtlose digitale Messaufnehmer heute billig und vielseitig sind (sie messen Dehnung, Windgeschwindigkeit, dreiachsige Beschleunigung, was auch immer …). 1000 davon an einem modernen Bauwerk wie einer Hängebrücke zu haben, ist absolut vernünftig. Es muss ein Steuerungssystem geben, das den Datenstrom auswertet und bereit ist, Alarm zu schlagen, wenn die Dinge beunruhigend werden.

Wieder gilt A x = b, wobei wir A (aus dem technischen Entwurf) und b (die gemessenen Werte der Messaufnehmer) haben, aber x finden müssen.

Ein erster Gedanke ist für alle, die damit nicht vertraut sind, irgendwie „durch A zu teilen“, um es auf die rechte Seite zu bringen.

Tatsächlich hat A eine Inverse (die meisten, aber nicht alle quadratischen Matrizen haben eine: die Determinante muss ungleich null sein), und die Berechnung funktioniert (langsam! ).

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

Leider ist die Berechnung der Inversen bei kleinen Matrizen langsam und bei großen im Schneckentempo. Nicht ideal, wenn deine Brücke während der Berechnung einstürzt!

Glücklicherweise sind andere Algorithmen dramatisch schneller (in diesem Fall die Gauß-Elimination).

Julia (in Anlehnung an Matlab) verwendet einfach einen Backslash für den Löser (technisch gesehen die „Linksdivision“).

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

Das ist ein triviales Beispiel, aber es gibt ein paar Details, auf die du achten solltest.

• Die Koeffizienten müssen in Spalten ausgerichtet sein, also füge bei fehlenden Termen Nullen ein.
• Die Gleichungen (und damit die Zeilen der Matrix) müssen linear unabhängig sein. Wenn Zeile 3 nur die Summe der Zeilen 1 und 2 ist, fügt sie keine neue Information hinzu und das Problem ist unterbestimmt.

Der Befehl rank() gibt die Anzahl der linear unabhängigen Zeilen/Spalten an, also strebe an, dass diese gleich der Anzahl der unbekannten Variablen ist.

julia> rank(A)
2

Rechteckige Matrizen

Irgendwann als Kind hast du wahrscheinlich gelernt, dass du ein System von N Gleichungen brauchst, um N Unbekannte zu lösen.

Wie bei den meisten Dingen in der Mathematik ist die Realität etwas differenzierter (nicht nur das Detail der linearen Unabhängigkeit, das im vorherigen Abschnitt erwähnt wurde).

„Zu wenige“ Zeilen

Mit N-1 Gleichungen (Zeilen in der Matrix A) ist das Problem unterbestimmt. Geben wir verzweifelt auf?

Es kommt darauf an!

Die N-1 Gleichungen enthalten immer noch eine Menge Informationen. Geometrisch betrachtet brauchen wir N, um die Lösung als Punkt im N-dimensionalen Raum zu bestimmen, aber N-1 sagt uns, dass die Lösung irgendwo auf einer Geraden liegen muss.

Damit bleiben uns unendlich viele Lösungen, aber ebenso unendlich viel Parameterraum, in dem es keine Lösungen gibt.

Verläuft deine Lösungsgerade durch einen Bereich des Parameterraums, über den ein seriöser Ingenieur sich Sorgen machen sollte? Oder der rechtfertigen könnte, das Bauwerk für die Öffentlichkeit zu sperren (Hängebrücken können bei starkem Seitenwind lebhaft werden)? Vielleicht lohnt es sich, dafür noch ein paar weitere Berechnungen anzustellen!

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

Der Backslash-Löser liefert uns eine Lösung! Inoffiziell ist dies offenbar die Lösung mit der kleinsten Norm (geometrisch am nächsten am Ursprung).

Um die anderen Lösungen zu erhalten, brauchen wir den nullspace von A2.

Addiert man ein beliebiges skalares Vielfaches des Nullraums zur ursprünglichen Lösung, erhält man eine weitere gültige Lösung.

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

Ähnlich schränken N-2 Gleichungen die (unendlichen) Lösungen auf eine Ebene ein, und dieselben Prinzipien gelten.

„Zu viele“ Zeilen

Die entgegengesetzte Situation tritt auf, wenn es mehr Gleichungen als Unbekannte gibt, sie aber trotzdem linear unabhängig sind.

In der realen Ingenieurspraxis ist das völlig normal und wird als etwas Gutes angesehen!

Messaufnehmer haben eine begrenzte Genauigkeit, der b-Vektor hat eine begrenzte Genauigkeit, und in der Berechnung steckt Rauschen. Jetzt ist die Lösung eine Anpassung nach der Methode der kleinsten Quadrate an die verrauschten Daten.

Die einfachste Technik verwendet die pseudoinverse der Matrix, die Julia als Funktion pinv() implementiert.

Diese lässt sich ähnlich wie die Inverse einer (nichtsingulären) quadratischen Matrix verwenden und liefert eine Schätzung der Variablen nach der Methode der kleinsten Quadrate statt einer exakten Lösung.

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

Eigenwerte und Eigenvektoren

Der vorherige Abschnitt beschrieb Lösungen für A x = b, was der vertrauten Algebra entspricht.

Dieser Abschnitt handelt von Lösungen für A x = λ x, wobei λ ein Skalar ist.

Das ist ein typisch lineare-Algebra-Konzept, aber wir können es geometrisch interpretieren:

Für geeignete Werte von x und λ streckt die quadratische Matrix A den Vektor x ungleich null um den Faktor λ in der Länge, ohne seine Richtung zu ändern (außer dass ein negatives λ sie umkehrt).

Das klingt nach einem Nischenthema, aber es stellt sich als lächerlich nützlich heraus!

Terminologie: Gültige Werte von λ sind die Eigenwerte von A, und die entsprechenden Werte von x sind die Eigenvektoren von A. Leider müssen wir eben mit Wörtern leben, die mitten im Wort die Sprache wechseln.

Studierenden wird normalerweise beigebracht, wie man Eigenwerte/-vektoren für 2×2-Matrizen von Hand berechnet, aber einen Computer zu benutzen ist viel einfacher (wenn auch immer noch ziemlich langsam und im allgemeinen Fall mit O(n^3) für eine n×n-Matrix skalierend).

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

Im Allgemeinen hat eine n×n-Matrix n Eigenwerte, allerdings nicht immer verschiedene Werte. Betrachte sie als die Wurzeln eines Polynoms n-ten Grades (das charakteristische Polynom), die mehrfach vorkommen können und selbst bei einer reellwertigen Matrix oft komplex sind.

Jeder Eigenvektor repräsentiert eine Richtung, und jedes skalare Vielfache ist ebenfalls ein gültiger Eigenvektor. Zur Vereinfachung weiterer Berechnungen gibt Julia Einheitsvektoren mit einer Norm von 1 zurück.

Anwendungen

Zu erklären, warum Eigenvektoren wichtig sind, ist eigentlich eine Aufgabe für ein 500-seitiges Lehrbuch und nicht für ein paar kurze Absätze. Dieses Thema durchdringt sehr viel der modernen angewandten Mathematik.

Auf oberster Ebene repräsentieren Eigenvektoren die „wichtigsten“ Achsen (Richtungen) in einem Datensatz. Eigenwerte geben die „relative Wichtigkeit“ jeder Achse an (unter bestimmten Annahmen über die Normierung der Eingaben).

Was das tatsächlich bedeutet, hängt von der Anwendung ab.

Hauptkomponentenanalyse

In einem typischen Data-Science-Problem haben wir vielleicht 100 oder mehr „Merkmale“, die als Datenspalten gespeichert sind. Fast zwangsläufig gibt es Rauschen, Redundanz und unerwünschte Korrelationen.

Wir müssen die Dimension reduzieren, und die PCA ist eine Möglichkeit, Ordnung ins Chaos zu bringen:

• Berechne die Kovarianzmatrix der Daten.
• Jetzt haben wir eine quadratische Matrix, also berechne als Nächstes die Eigenwerte und Eigenvektoren.
• Sortiere die Eigenvektoren absteigend nach ihren Eigenwerten (genauer gesagt nach dem Absolutwert |λ|).
• Die obersten k Eigenvektoren sind nun deine principal components, wobei k deutlich kleiner ist als die ursprüngliche Anzahl der Merkmale.
• Projiziere den vollständigen Datensatz auf diese k Achsen und beginne, nach interessanten Mustern zu suchen.

In dieser Form wurde die PCA von so unterschiedlichen Gruppen genutzt wie Medizinwissenschaftlern, die molekulare Eigenschaften im Wirkstoffdesign verbessern wollen, und politischen Wahlkämpfern, die aus Umfragedaten die Präferenzen der Wähler verstehen wollen. Wahrscheinlich auch von Marketinggruppen, die dir etwas verkaufen wollen, aber jede Technologie kann zum Guten oder zum Schlechten eingesetzt werden.

Die PCA ist nicht Teil einer minimalen Julia-Installation, aber (außerhalb von Exercism) enthält das Paket MultivariateStats, was du brauchst.

Bildverarbeitung

Wenn man die PCA-Idee einen Schritt weiterführt, sind digitale Bilder einfach Matrizen von Pixelwerten, und wir können Hauptkomponenten dafür berechnen.

Das wird seit Jahrzehnten in der Bildkompression verwendet, wo die PCA ein Leitfaden dafür ist, was beim Verkleinern der Datei am wichtigsten zu erhalten ist.

Zunehmend ist die PCA ein wesentlicher Bestandteil von Bildklassifikatoren wie der Gesichtserkennung. Wenn du das nächste Mal durch einen Flughafen gehst, schaut Big Brother nicht nur zu, er verwendet auch lineare Algebra, um zu verstehen, was er sieht!

Maschinenbau

Weniger kontrovers: Die Hauptachsen eines mechanischen Bauteils sind die Eigenvektoren seines Trägheitstensors I (einer Matrix mit einem etwas anderen Namen).

Jedes Rad an deinem Auto hat wahrscheinlich ein oder mehrere Auswuchtgewichte, die so justiert sind, dass die Nichtdiagonalelemente in diesem Tensor I zu null werden. Der Mechaniker in der Werkstatt rechnet das sicher nicht durch (im Gegensatz zu Flugzeugkonstrukteuren und Raketeningenieuren), aber „Wobbeln“ ist nur ein Alltagswort für die Kreuzterme im Tensor, die in diesem Fall deine Fahrt weniger angenehm machen und den mechanischen Verschleiß der Lager erhöhen würden.