6.13. Lineární algebra¶

Lineární algebra na kameře je práce s malými maticemi: rotace 3x3, která fúzuje vzorek z IMU do souřadnic světa, kalibrační matice, která rektifikuje objektiv, aktualizace stavové kovariance Kalmanova filtru, polynomiální proložení, jehož normální rovnice vyjdou jako drobné lineární řešení. Submoduly numpy.linalg a scipy.linalg pokrývají přesně tento rozsah.

Funkce sídlí ve dvou modulech. Operace maticového součinu jsou na nejvyšší úrovni numpy; dekompozice a inverze matice jsou pod numpy.linalg; vyhrazené řešiče lineárních soustav jsou pod scipy.linalg. Například násobení matic 2x2:

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. Co je k dispozici¶

dot() – maticový nebo vektorový součin.
cross() – vektorový součin 3D vektorů.
trace() – součet diagonály.
inv() – inverze matice.
det() – determinant.
cholesky() – Choleského rozklad (symetrický pozitivně definitní vstup).
eig() – vlastní čísla a vlastní vektory reálné symetrické matice.
norm() – 2-norma vektoru nebo matice.
qr() – QR rozklad s mode='reduced' (výchozí) nebo mode='complete'.
solve_triangular() – řeší A @ x = b, když je A trojúhelníková.
cho_solve() – řeší A @ x = b při zadaném Choleského faktoru matice A.

6.13.2. dot, cross, trace¶

dot() je způsob, jakým se na kameře zapisuje násobení matic. Operátor @, který desktopový numpy poskytuje pro stejnou úlohu, není na ndarray implementován, takže každý maticově-vektorový, maticově-maticový a vektorový skalární součin prochází voláním dot()

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

Výsledek dot() je vždy typu float. cross() je vektorový součin dvou 3-vektorů, trace() součet hlavní diagonály čtvercové matice.

6.13.3. inv a 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))

Inverze se počítá Gaussovou-Jordanovou eliminací, takže inv() vyvolá ValueError, když je matice singulární (diagonální prvek se během eliminace stane nulovým). Náklady na RAM jsou zhruba dvojnásobkem velikosti vstupu.

Determinant znovu využívá tutéž eliminaci – doba běhu je v podstatě stejná jako u inverze.

Když je cílem aplikace vyřešit lineární soustavu, neinvertujte a nenásobte – dejte přednost vyhrazeným řešičům níže. Oba jsou rychlejší a numericky se chovají lépe.

6.13.4. cholesky¶

Pro symetrickou pozitivně definitní matici A vrací cholesky() dolní trojúhelníkovou matici L takovou, že A = L @ L.T

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

Pokud vstup není pozitivně definitní nebo není symetrický, vyvolá se ValueError.

Choleského faktor představuje polovinu práce LU faktorizace a je správným výchozím bodem pro každý problém, jehož matice je známa jako symetrická pozitivně definitní (aktualizace kovariance, normální rovnice z proložení metodou nejmenších čtverců).

6.13.5. eig¶

eig() funguje pouze na reálných symetrických maticích. Nesymetrické matice vyvolají ValueError. Vrací dvojici (eigenvalues, eigenvectors)

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)

Poznámky:

Vlastní čísla se vracejí bez konkrétního pořadí. Když na seřazeném pořadí záleží, použijte sort() (a stejnou permutaci aplikujte na vlastní vektory pomocí argsort()).
Vlastní vektor je jednoznačný pouze až na nenulový skalár, takže znaménko jednotlivých vlastních vektorů není jednoznačně určeno. Dva správné běhy mohou vyprodukovat vektory opačného znaménka; to je neškodné.

6.13.6. norm¶

Eukleidovská (Frobeniova) norma vektoru nebo matice:

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

Volitelné klíčové slovo axis= počítá normu podél jediné osy namísto přes celé pole.

6.13.7. qr¶

qr() rozkládá obdélníkovou matici A (tvar (M, N)) na ortonormální Q a horní trojúhelníkovou R takovou, že 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)

Rozklad je implementován pomocí postupných Givensových rotací. Správná volba pro problémy nejmenších čtverců, kde matice není symetrická.

6.13.8. Řešení soustav¶

Dva vyhrazené řešiče pod ulab.scipy.linalg jsou rychlejší i přesnější než np.dot(np.linalg.inv(A), b):

solve_triangular(a, b, lower=False)() – řeší a @ x = b za předpokladu, že a je trojúhelníková:

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)() – při zadaném Choleského faktoru L řeší A @ x = b, kde A = L @ L.T
```
L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)
```

Sáhněte po nich namísto invertování, kdykoli to struktura A dovolí – ušetří práci s eliminací i explicitní inverzi.

6.13.9. Řešení malé lineární soustavy¶

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

Tentýž problém lze vyjádřit faktorizací A a voláním vhodného řešiče – rychlejší a přesnější, když má A správnou strukturu.

Úplnou referenci na úrovni argumentů najdete v numpy.linalg — Rutiny lineární algebry a scipy.linalg — Rutiny lineární algebry.