6.13. Lineare Algebra¶

Lineare Algebra auf einer Kamera ist Arbeit mit kleinen Matrizen: eine 3x3-Rotation, die ein IMU-Sample in das Weltkoordinatensystem einbettet, eine Kalibriermatrix, die ein Objektiv entzerrt, das Zustands-Kovarianz-Update eines Kalman-Filters, eine Polynomanpassung, deren Normalgleichungen sich als winziges lineares Gleichungssystem ergeben. Die Untermodule numpy.linalg und scipy.linalg decken genau diese Größenordnung ab.

Die Funktionen liegen in zwei Modulen. Matrixprodukt-Operationen befinden sich auf der obersten Ebene von numpy; Zerlegungen und die Matrixinverse liegen unter numpy.linalg; die spezialisierten Löser für lineare Gleichungssysteme liegen unter scipy.linalg. Eine 2x2-Matrixmultiplikation zum Beispiel:

from ulab import numpy as np

A = np.array([[1, 2], [3, 4]], dtype=np.float)
B = np.array([[5, 6], [7, 8]], dtype=np.float)
np.dot(A, B)
# array([[19.0, 22.0],
#        [43.0, 50.0]])

6.13.1. Was verfügbar ist¶

dot() – Matrix- oder Vektorprodukt.
cross() – 3D-Vektorkreuzprodukt.
trace() – Summe der Diagonalen.
inv() – Matrixinverse.
det() – Determinante.
cholesky() – Cholesky-Zerlegung (symmetrisch positiv definite Eingabe).
eig() – Eigenwerte und Eigenvektoren einer reellen symmetrischen Matrix.
norm() – 2-Norm eines Vektors oder einer Matrix.
qr() – QR-Zerlegung mit mode='reduced' (Standard) oder mode='complete'.
solve_triangular() – löst A @ x = b, wenn A dreieckig ist.
cho_solve() – löst A @ x = b bei gegebenem Cholesky-Faktor von A.

6.13.2. dot, cross, trace¶

dot() ist die Art und Weise, wie Matrixmultiplikation auf der Kamera geschrieben wird. Der Operator @, den Desktop-numpy für dieselbe Aufgabe bereitstellt, ist auf ndarray nicht implementiert, sodass jedes Matrix-Vektor-, Matrix-Matrix- und Vektorskalarprodukt über den Aufruf dot() läuft:

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.dot(a, b)             # 32.0 (scalar product)
np.cross(a, b)           # array([-3.0, 6.0, -3.0])

m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
np.dot(m, a)             # matrix-vector product
np.dot(m, m)             # matrix-matrix product
np.trace(m)              # 1 + 5 + 9 = 15.0

Das Ergebnis von dot() hat immer den dtype float. cross() ist das Kreuzprodukt zweier 3er-Vektoren, trace() die Summe der Hauptdiagonalen einer quadratischen Matrix.

6.13.3. inv und det¶

m = np.array([[1, 2, 3, 4],
              [4, 5, 6, 4],
              [7, 9, 9, 4],
              [3, 4, 5, 6]])
print(np.linalg.inv(m))
print(np.linalg.det(m))

Die Inverse wird durch Gauß-Jordan-Elimination berechnet, daher löst inv() einen ValueError aus, wenn die Matrix singulär ist (ein Diagonalelement wird während der Elimination null). Der RAM-Bedarf beträgt etwa das Doppelte der Eingabegröße.

Die Determinante verwendet dieselbe Elimination wieder – die Laufzeit entspricht im Wesentlichen der der Inverse.

Wenn das Ziel der Anwendung darin besteht, ein lineares Gleichungssystem zu lösen, sollte man nicht invertieren und multiplizieren – bevorzugen Sie die spezialisierten Löser weiter unten. Beide sind schneller und numerisch gutartiger.

6.13.4. cholesky¶

Für eine symmetrisch positiv definite Matrix A liefert cholesky() ein unteres Dreieck L zurück, sodass A = L @ L.T:

a = np.array([[25, 15, -5],
              [15, 18,  0],
              [-5,  0, 11]])
L = np.linalg.cholesky(a)

Wenn die Eingabe nicht positiv definit oder nicht symmetrisch ist, wird ein ValueError ausgelöst.

Der Cholesky-Faktor erfordert die Hälfte des Aufwands einer LU-Zerlegung und ist der richtige Ausgangspunkt für jedes Problem, dessen Matrix bekanntermaßen symmetrisch positiv definit ist (Kovarianz-Updates, Normalgleichungen aus einer Kleinste-Quadrate-Anpassung).

6.13.5. eig¶

eig() funktioniert nur mit reellen symmetrischen Matrizen. Nicht-symmetrische Matrizen lösen einen ValueError aus. Es liefert ein 2-Tupel (eigenvalues, eigenvectors) zurück:

a = np.array([[1, 2, 1, 4],
              [2, 5, 3, 5],
              [1, 3, 6, 1],
              [4, 5, 1, 7]], dtype=np.uint8)
x, y = np.linalg.eig(a)

Hinweise:

Die Eigenwerte werden in keiner bestimmten Reihenfolge zurückgegeben. Wenden Sie sort() an (und über argsort() dieselbe Permutation auf die Eigenvektoren), wenn eine sortierte Reihenfolge wichtig ist.
Ein Eigenvektor ist nur bis auf einen von null verschiedenen Skalar eindeutig, daher ist das Vorzeichen einzelner Eigenvektoren nicht eindeutig festgelegt. Zwei korrekte Durchläufe können Vektoren mit entgegengesetztem Vorzeichen erzeugen; das ist unbedenklich.

6.13.6. norm¶

Die euklidische (Frobenius-)Norm eines Vektors oder einer Matrix:

v = np.array([1, 2, 3, 4, 5])
m = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

np.linalg.norm(v)        # 7.416...
np.linalg.norm(m)        # 16.881...

Das optionale Schlüsselwort axis= bildet die Norm entlang einer einzelnen Achse statt über das gesamte Array.

6.13.7. qr¶

qr() zerlegt eine rechteckige Matrix A (Form (M, N)) in ein orthonormales Q und ein oberes Dreieck R, sodass A == Q @ R:

A = np.arange(6).reshape((3, 2))

q, r = np.linalg.qr(A)
# mode='reduced' (default): q is (3, 2), r is (2, 2)

q, r = np.linalg.qr(A, mode='complete')
# q is (3, 3), r is (3, 2)

Die Zerlegung ist über aufeinanderfolgende Givens-Rotationen implementiert. Die richtige Wahl für Kleinste-Quadrate-Probleme, bei denen die Matrix nicht symmetrisch ist.

6.13.8. Gleichungssysteme lösen¶

Die beiden spezialisierten Löser unter ulab.scipy.linalg sind beide schneller und genauer als np.dot(np.linalg.inv(A), b):

solve_triangular(a, b, lower=False)() – löst a @ x = b unter der Annahme, dass a dreieckig ist:

A = np.array([[3, 0, 0, 0],
              [2, 1, 0, 0],
              [1, 0, 1, 0],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])
x = sp.linalg.solve_triangular(A, b, lower=True)

cho_solve(L, b)() – löst bei gegebenem Cholesky-Faktor L das System A @ x = b, wobei A = L @ L.T:
```
L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)
```

Greifen Sie zu diesen anstelle der Invertierung, wann immer die Struktur von A es zulässt – sie sparen sowohl den Eliminationsaufwand als auch die explizite Inverse.

6.13.9. Ein kleines lineares Gleichungssystem lösen¶

A = np.array([[3, 0, 1, 1],
              [0, 1, 0, 2],
              [1, 0, 1, 1],
              [1, 2, 1, 8]])
b = np.array([4, 2, 4, 2])

x = np.dot(np.linalg.inv(A), b)
print(x)
print(np.dot(A, x))         # should equal b

Dasselbe Problem lässt sich ausdrücken, indem man A zerlegt und den passenden Löser aufruft – schneller und genauer, wenn A die richtige Struktur hat.

Die vollständige Referenz auf Argumentebene finden Sie unter numpy.linalg — Routinen für lineare Algebra und scipy.linalg — Routinen der linearen Algebra.