6.13. Đại số tuyến tính

Đại số tuyến tính trên camera là công việc với ma trận nhỏ: một phép quay 3x3 hợp nhất một mẫu IMU vào khung thế giới, một ma trận hiệu chỉnh chỉnh lại một thấu kính, cập nhật hiệp phương sai trạng thái của bộ lọc Kalman, một bài khớp đa thức mà các phương trình chuẩn tắc cho ra một hệ tuyến tính nhỏ. Các mô-đun con numpy.linalg và scipy.linalg xử lý đúng quy mô đó.

Các hàm nằm trong hai mô-đun. Các phép toán tích ma trận nằm ở cấp trên của numpy; các phép phân tích và nghịch đảo ma trận nằm dưới numpy.linalg; các bộ giải hệ tuyến tính chuyên dụng nằm dưới scipy.linalg. Ví dụ, nhân hai ma trận 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. Những gì khả dụng

• dot() -- tích ma trận hoặc tích vectơ.

• cross() -- tích có hướng vectơ 3 chiều.

• trace() -- tổng đường chéo chính.

• inv() -- nghịch đảo ma trận.

• det() -- định thức.

• cholesky() -- phân tích Cholesky (đầu vào xác định dương đối xứng).

• eig() -- trị riêng và vectơ riêng của ma trận đối xứng thực.

• norm() -- chuẩn 2 của vectơ hoặc ma trận.

• qr() -- phân tích QR với mode='reduced' (mặc định) hoặc mode='complete'.

• solve_triangular() -- giải A @ x = b khi A là ma trận tam giác.

• cho_solve() -- giải A @ x = b từ nhân tử Cholesky của A.

6.13.2. dot, cross, trace

dot() là cách viết phép nhân ma trận trên camera. Toán tử @ mà numpy trên máy tính để bàn cung cấp cho cùng công việc đó không được triển khai trên ndarray, vì vậy mọi tích vectơ-ma trận, ma trận-ma trận, và tích vô hướng đều qua lệnh gọi 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

Kết quả của dot() luôn có dtype float. cross() là tích có hướng của hai vectơ 3 chiều, trace() là tổng đường chéo chính của ma trận vuông.

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

Nghịch đảo được tính bằng phương pháp khử Gauss-Jordan, vì vậy inv() ném ValueError khi ma trận suy biến (một phần tử đường chéo trở thành 0 trong quá trình khử). Chi phí RAM xấp xỉ gấp đôi kích thước đầu vào.

Định thức tái sử dụng cùng phép khử đó -- thời gian chạy về cơ bản bằng với nghịch đảo.

Khi mục tiêu của ứng dụng là giải một hệ tuyến tính, đừng nghịch đảo rồi nhân -- hãy ưu tiên các bộ giải chuyên dụng bên dưới. Cả hai đều nhanh hơn và ổn định số hơn.

6.13.4. cholesky

Với ma trận xác định dương đối xứng A, cholesky() trả về ma trận tam giác dưới L sao cho A = L @ L.T

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

Nếu đầu vào không xác định dương hoặc không đối xứng, ValueError sẽ được ném ra.

Nhân tử Cholesky chỉ cần nửa công việc của một phân tích LU và là điểm khởi đầu đúng đắn cho bất kỳ bài toán nào mà ma trận được biết là xác định dương đối xứng (cập nhật hiệp phương sai, phương trình chuẩn tắc từ bài khớp bình phương nhỏ nhất).

6.13.5. eig

eig() chỉ hoạt động trên ma trận đối xứng thực. Ma trận không đối xứng ném ValueError. Nó trả về một bộ 2 phần tử (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)

Lưu ý:

• Các trị riêng trả về không theo thứ tự cụ thể nào. Áp dụng sort() (và cùng hoán vị cho các vectơ riêng qua argsort()) khi cần thứ tự đã sắp xếp.

• Một vectơ riêng chỉ duy nhất đến một số vô hướng khác không, vì vậy dấu của các vectơ riêng riêng lẻ không được xác định duy nhất. Hai lần chạy đúng có thể tạo ra các vectơ có dấu ngược nhau; điều này vô hại.

6.13.6. norm

Chuẩn Euclid (Frobenius) của vectơ hoặc ma trận:

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

Từ khóa tùy chọn axis= lấy chuẩn theo một trục duy nhất thay vì trên toàn bộ mảng.

6.13.7. qr

qr() phân tích ma trận chữ nhật A (hình dạng (M, N)) thành Q trực chuẩn và R tam giác trên sao cho 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)

Phân tích được thực hiện qua các phép quay Givens liên tiếp. Lựa chọn phù hợp cho các bài toán bình phương nhỏ nhất khi ma trận không đối xứng.

6.13.8. Giải hệ phương trình

Hai bộ giải chuyên dụng trong ulab.scipy.linalg đều nhanh hơn và chính xác hơn np.dot(np.linalg.inv(A), b):

• solve_triangular(a, b, lower=False)() -- giải a @ x = b giả sử a là ma trận tam giác:

  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)() -- từ nhân tử Cholesky L, giải A @ x = b với A = L @ L.T

  L = np.linalg.cholesky(A)
  x = sp.linalg.cho_solve(L, b)
  

Hãy dùng những cái này thay vì nghịch đảo bất cứ khi nào cấu trúc của A cho phép -- chúng tiết kiệm công việc khử và tránh phải tính nghịch đảo tường minh.

6.13.9. Giải một hệ tuyến tính nhỏ

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

Cùng bài toán đó có thể được biểu diễn bằng cách phân tích A và gọi bộ giải phù hợp -- nhanh hơn và chính xác hơn khi A có cấu trúc phù hợp.

Để tham khảo đầy đủ cấp đối số, xem numpy.linalg --- Các hàm đại số tuyến tính và scipy.linalg --- Các thủ tục đại số tuyến tính.