6.13. 線性代數¶

在相機上的線性代數屬於小型矩陣運算：將 IMU 取樣融合到世界座標系的 3x3 旋轉、校正鏡頭的標定矩陣、卡爾曼濾波器的狀態-共變異數更新、其正規方程式化為微型線性求解的多項式擬合。numpy.linalg 與 scipy.linalg 子模組正好涵蓋這個規模。

這些函式分布在兩個模組中。矩陣乘積運算位於 numpy 頂層；分解與矩陣反矩陣位於 numpy.linalg 之下；專用的線性系統求解器位於 scipy.linalg 之下。舉例來說，一個 2 乘 2 的矩陣相乘:

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. 可用的功能¶

dot() -- 矩陣或向量乘積。
cross() -- 三維向量外積。
trace() -- 對角線元素之和。
inv() -- 矩陣反矩陣。
det() -- 行列式。
cholesky() -- Cholesky 分解（輸入需為對稱正定矩陣）。
eig() -- 實對稱矩陣的特徵值與特徵向量。
norm() -- 向量或矩陣的 2-範數。
qr() -- QR 分解，可使用 mode='reduced'（預設）或 mode='complete'。
solve_triangular() -- 當 A 為三角矩陣時，求解 A @ x = b。
cho_solve() -- 給定 A 的 Cholesky 因子，求解 A @ x = b。

6.13.2. dot、cross、trace¶

dot() 是在相機上撰寫矩陣乘法的方式。桌面版 numpy 用於相同工作的 @ 運算子並未在 ndarray 上實作，因此每個矩陣-向量、矩陣-矩陣以及向量點積都透過 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

dot() 的結果其 dtype 永遠是 float。cross() 是兩個三維向量的外積，trace() 則是方陣主對角線之和。

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

反矩陣是透過高斯-喬丹消去法計算的，因此當矩陣為奇異矩陣時（消去過程中某個對角線元素變為零），inv() 會引發 ValueError。其 RAM 成本約為輸入大小的兩倍。

行列式重複使用相同的消去過程 -- 執行時間基本上與反矩陣相同。

當應用程式的目標是求解線性系統時，請勿先求反矩陣再相乘 -- 應優先使用下方的專用求解器。兩者都更快，且數值表現更佳。

6.13.4. cholesky¶

對於對稱正定矩陣 A，cholesky() 會傳回一個下三角矩陣 L，使得 A = L @ L.T:

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

若輸入並非正定或並非對稱，則會引發 ValueError。

Cholesky 因子的計算量只有 LU 分解的一半，對於任何已知矩陣為對稱正定的問題（共變異數更新、最小平方擬合的正規方程式），都是正確的起點。

6.13.5. eig¶

eig() 僅適用於實對稱矩陣。非對稱矩陣會引發 ValueError。它會傳回一個 2-元組 (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)

注意事項：

特徵值傳回時不具任何特定順序。當需要排序順序時，請套用 sort()（並透過 argsort() 對特徵向量套用相同的排列）。
特徵向量僅在差一個非零純量的範圍內是唯一的，因此個別特徵向量的正負號並非唯一確定。兩次正確的執行可能產生正負號相反的向量；這並無妨害。

6.13.6. norm¶

向量或矩陣的歐幾里得（Frobenius）範數:

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

選用的 axis= 關鍵字會沿著單一軸取範數，而非對整個陣列取範數。

6.13.7. qr¶

qr() 會將一個矩形矩陣 A（形狀為 (M, N)）分解為一個正交規範矩陣 Q 與一個上三角矩陣 R，使得 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)

此分解是透過連續的 Givens 旋轉實作的。對於矩陣非對稱的最小平方問題而言，這是正確的選擇。

6.13.8. 求解系統¶

ulab.scipy.linalg 之下的兩個專用求解器，都比 np.dot(np.linalg.inv(A), b) 更快且更精確：

solve_triangular(a, b, lower=False)() -- 假設 a 為三角矩陣，求解 a @ x = b:

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)() -- 給定 Cholesky 因子 L，求解 A @ x = b，其中 A = L @ L.T:
```
L = np.linalg.cholesky(A)
x = sp.linalg.cho_solve(L, b)
```

只要 A 的結構允許，就應使用這些求解器而非求反矩陣 -- 它們既能省下消去運算又能省下明確的反矩陣計算。

6.13.9. 求解小型線性系統¶

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

同樣的問題也可以透過分解 A 並呼叫適當的求解器來表達 -- 當 A 具有合適的結構時，會更快且更精確。

完整的引數層級參考資料，請參閱 numpy.linalg --- 線性代數常式與 scipy.linalg --- 線性代數常式。