Phân rã ma trận¶

Các phép phân rã ma trận chia nhỏ các ma trận phức tạp thành những nhân tử đơn giản hơn để giải hệ phương trình, tính nghịch đảo, và nén dữ liệu. File này bao quát khử Gauss, LU, QR, Cholesky, phân rã trị riêng, và SVD — những thuật toán đứng sau PCA, hệ thống gợi ý, và độ ổn định số trong ML.

Một phép phân rã (hay phân tích nhân tử) ma trận chia ma trận thành những mảnh đơn giản hơn dễ làm việc hơn. Hãy tưởng tượng nó giống như phân tích một số ra thừa số: \(12 = 3 \times 4\) dễ lập luận hơn so với số 12 một mình.
Chúng ta phân rã ma trận để giải các hệ phương trình nhanh hơn, tính toán nghịch đảo ổn định, tìm trị riêng, nén dữ liệu, và hiểu hình học của các phép biến đổi.
Kỹ thuật cơ bản nhất là khử Gauss (Gaussian elimination) (rút gọn hàng). Ý tưởng rất đơn giản: cho một hệ \(A\mathbf{x} = \mathbf{b}\), dùng ba phép toán được phép để đơn giản hóa \(A\) cho đến khi câu trả lời hiển nhiên.
Các phép toán là: hoán đổi hai hàng, nhân một hàng với một vô hướng khác không, hoặc cộng một bội số của hàng này vào hàng khác.
Ví dụ, để khử cột đầu tiên bên dưới trụ, trừ các bội số của hàng 1 từ các hàng bên dưới:

\[ \begin{bmatrix} 2 & 1 & 5 \\ 4 & 3 & 7 \\ 6 & 5 & 9 \end{bmatrix} \xrightarrow{R_2 - 2R_1} \begin{bmatrix} 2 & 1 & 5 \\ 0 & 1 & -3 \\ 6 & 5 & 9 \end{bmatrix} \xrightarrow{R_3 - 3R_1} \begin{bmatrix} 2 & 1 & 5 \\ 0 & 1 & -3 \\ 0 & 2 & -6 \end{bmatrix} \]

Mục tiêu là dạng bậc thang hàng (REF): các số không dưới mỗi trụ (phần tử khác không đầu tiên trong mỗi hàng), với mỗi trụ ở bên phải trụ của hàng phía trên nó. Ma trận trở thành hình dạng cầu thang.

Khử Gauss: các phép toán hàng tạo ra dạng tam giác, sau đó giải từ dưới lên

Đi xa hơn đến dạng bậc thang hàng rút gọn (RREF) làm cho mỗi trụ bằng 1 và là phần tử khác không duy nhất trong cột của nó. Mọi ma trận có một RREF duy nhất.
Một khi ở dạng tam giác, ta giải bằng thế ngược: hàng dưới cùng cho biến cuối trực tiếp, rồi làm ngược lên.
Đây là nền tảng mà tất cả các phép phân rã khác xây dựng lên, mục tiêu của phân rã là đưa ma trận về dạng tam giác, để ta có thể thế ngược và giải cho các biến.
Phân rã LU hình thức hóa khử Gauss bằng cách phân tích một ma trận vuông thành \(A = LU\) (hoặc \(A = PLU\) với hoán đổi hàng), với \(L\) là tam giác dưới và \(U\) là tam giác trên.

Phân rã LU: chia một ma trận khó thành hai ma trận tam giác dễ

Để giải \(A\mathbf{x} = \mathbf{b}\): trước tiên giải \(L\mathbf{y} = \mathbf{b}\) bằng thế xuôi (từ trên xuống dưới), sau đó giải \(U\mathbf{x} = \mathbf{y}\) bằng thế ngược (từ dưới lên trên). Hai phép giải tam giác dễ thay vì một phép giải tổng quát khó.
Lợi thế so với khử Gauss thô là khả năng tái sử dụng. Một khi bạn có \(L\) và \(U\), bạn có thể giải cho nhiều vector \(\mathbf{b}\) khác nhau mà không cần làm lại phép phân rã.
Nếu bạn cần giải cùng một hệ với 1000 vế phải khác nhau (phổ biến trong mô phỏng), bạn phân rã một lần và tái sử dụng.
Khi một ma trận đối xứng và xác định dương (như ma trận hiệp phương sai), ta có thể làm tốt hơn nữa.
Phân rã Cholesky phân tích nó thành \(A = LL^T\), với \(L\) là tam giác dưới. Ví dụ:

\[ \begin{bmatrix} 4 & 2 \\ 2 & 5 \end{bmatrix} = \begin{bmatrix} 2 & 0 \\ 1 & 2 \end{bmatrix} \begin{bmatrix} 2 & 1 \\ 0 & 2 \end{bmatrix} \]

Cách này nhanh gấp đôi LU và được đảm bảo ổn định về số. Hãy nghĩ về nó như "căn bậc hai" của ma trận.
Nếu phép phân rã thất bại (một giá trị âm dưới căn bậc hai), ma trận không phải xác định dương. Cholesky do đó cũng đóng vai trò một phép kiểm tra tính xác định dương.
Vector riêng (eigenvectors) của một ma trận vuông \(A\) là những hướng đặc biệt mà phép biến đổi chỉ co giãn hoặc thu nhỏ, không xoay. Trị riêng (eigenvalue) là hệ số co giãn:

\[A\mathbf{x} = \lambda\mathbf{x}\]

Vector riêng nằm yên trên cùng một đường thẳng (chỉ bị co giãn), vector thông thường bị xoay

Hầu hết các vector thay đổi hướng khi nhân với một ma trận. Nhưng vector riêng thì đặc biệt: đầu ra chỉ cùng hướng với đầu vào, chỉ được co giãn bởi \(\lambda\). Nếu \(\lambda = 2\), vector riêng tăng gấp đôi độ dài. Nếu \(\lambda = -1\), nó lật hướng. Nếu \(\lambda = 0\), nó bị ép xuống không.
Ví dụ, với:

\[ A = \begin{bmatrix} 3 & 1 \\ 0 & 2 \end{bmatrix} \]

vector \([1, 0]^T\) là một vector riêng với \(\lambda = 3\) vì \(A[1, 0]^T = [3, 0]^T = 3[1, 0]^T\).
Để tìm trị riêng, giải đa thức đặc trưng \(\det(A - \lambda I) = 0\). Các nghiệm là các trị riêng. Sau đó thay từng \(\lambda\) trở lại vào \((A - \lambda I)\mathbf{x} = \mathbf{0}\) để tìm các vector riêng tương ứng.
Các tính chất chính:
- Trace của \(A\) bằng tổng các trị riêng của nó.
- Định thức của \(A\) bằng tích các trị riêng của nó.
- Ma trận đối xứng có các vector riêng vuông góc nhau và trị riêng thực.
- Ma trận xác định dương có tất cả trị riêng dương.
- Ma trận hiệp phương sai (mà ta sẽ gặp trong thống kê) luôn xác định bán dương.
Tính toán trị riêng qua đa thức đặc trưng là không thực tế cho các ma trận lớn. Thay vào đó, các phương pháp lặp được sử dụng:
- Lặp lũy thừa (power iteration): nhân liên tục với \(A\) và chuẩn hóa. Hội tụ về vector riêng trội (trị riêng lớn nhất). Đơn giản nhưng chỉ tìm được một cặp trị riêng.
- Thuật toán QR: phương pháp chủ lực. Phân rã và tổ hợp lại nhiều lần dùng phân rã QR cho đến khi ma trận hội tụ về dạng tam giác, tiết lộ tất cả trị riêng trên đường chéo.
- Lặp ngược (inverse iteration): tìm vector riêng gần nhất với một giá trị mục tiêu cho trước. Hữu ích khi bạn biết gần đúng trị riêng bạn muốn.
- Với ma trận thưa lớn, lặp Arnoldi và Lanczos khai thác tính thưa để đạt hiệu quả.
Nếu một ma trận vuông có một tập hợp đầy đủ các vector riêng độc lập tuyến tính, nó có thể được chéo hóa (diagonalised): \(A = PDP^{-1}\), với \(D\) là ma trận đường chéo các trị riêng và các cột của \(P\) là các vector riêng.
Tại sao điều này hữu ích? Ma trận đường chéo rất dễ làm việc. Cần \(A^{100}\)? Thay vì nhân \(A\) với chính nó 100 lần, hãy tính \(PD^{100}P^{-1}\), và nâng một ma trận đường chéo lên lũy thừa chỉ là nâng từng phần tử trên đường chéo. Điều này biến một phép toán đắt đỏ thành một phép toán rẻ.
Một cơ sở trị riêng (eigenbasis) là một cơ sở được tạo hoàn toàn từ các vector riêng. Trong cơ sở này, ma trận trở thành đường chéo và phép biến đổi chỉ là co giãn độc lập dọc theo mỗi hướng vector riêng. Điều này giống như việc tìm ra hệ tọa độ tự nhiên cho phép biến đổi.
Phân rã QR phân tích bất kỳ ma trận nào \(A\) thành \(A = QR\), với \(Q\) là trực giao (các cột của nó trực chuẩn) và \(R\) là tam giác trên. Hãy nghĩ về nó như việc tách thông tin "hướng" (\(Q\)) khỏi thông tin "co giãn và pha trộn" (\(R\)).
Quy trình Gram-Schmidt xây dựng \(Q\) từng cột một. Lấy cột đầu tiên của \(A\) và chuẩn hóa nó. Lấy cột thứ hai, trừ đi phép chiếu của nó lên cột thứ nhất (để làm nó vuông góc), và chuẩn hóa. Lặp lại cho mỗi cột. Kết quả là một tập hợp các vector trực chuẩn.
Phân rã QR là động cơ đằng sau thuật toán QR cho trị riêng. Nó cũng được dùng trực tiếp để giải các bài toán bình phương tối thiểu: khi \(A\mathbf{x} = \mathbf{b}\) không có nghiệm chính xác (nhiều phương trình hơn ẩn), QR tìm câu trả lời xấp xỉ tốt nhất.
SVD (Singular Value Decomposition) là phép phân rã tổng quát nhất và được tranh luận là quan trọng nhất. Mọi ma trận (mọi hình dạng, mọi hạng) đều có SVD: \(A = U\Sigma V^T\)
- \(V^T\) (\(n \times n\), trực giao): xoay đầu vào
- \(\Sigma\) (\(m \times n\), đường chéo): co giãn dọc theo các trục trực giao (các giá trị singular, không âm, theo thứ tự giảm dần)
- \(U\) (\(m \times m\), trực giao): xoay đầu ra

SVD: mọi phép biến đổi = xoay, rồi co giãn, rồi xoay lần nữa

Về mặt hình học, SVD nói rằng mọi phép biến đổi tuyến tính, dù phức tạp đến đâu, chỉ là một phép xoay, tiếp theo là co giãn dọc theo các trục, rồi một phép xoay khác. Một hình tròn trở thành một hình elip.
Các giá trị singular (\(\sigma_1 \geq \sigma_2 \geq \ldots\)) tiết lộ "tầm quan trọng" của mỗi hướng. Các giá trị singular lớn tương ứng với các hướng quan trọng nhất. Hạng của \(A\) bằng số lượng giá trị singular khác không.
Xấp xỉ hạng thấp: bằng cách chỉ giữ lại \(k\) giá trị singular lớn nhất và đặt phần còn lại về không, bạn có được xấp xỉ hạng \(k\) tốt nhất có thể của \(A\). Đây là cách nén ảnh hoạt động: một ảnh \(1000 \times 1000\) có thể chỉ cần \(k = 50\) giá trị singular để trông gần như giống hệt, nén nó xuống 20 lần.
SVD cũng cung cấp giả nghịch đảo: \(A^+ = V\Sigma^+U^T\), với \(\Sigma^+\) nghịch đảo các giá trị singular khác không.
Trong khi phân rã trị riêng chỉ hoạt động với ma trận vuông, SVD hoạt động với bất kỳ ma trận nào. Đây là lợi thế chính của nó.
PCA (Phân tích thành phần chính) dùng phân rã trị riêng (hoặc SVD) để giảm số chiều.
Hãy tưởng tượng một tập dữ liệu với 100 đặc trưng mỗi mẫu (vector 100 chiều được xếp chồng thành một ma trận). Nhiều đặc trưng trong số đó có tương quan và dư thừa.
PCA tìm các hướng mà dữ liệu thực sự biến thiên, cho phép bạn chỉ giữ lại những gì quan trọng.

PCA tìm các hướng có phương sai lớn nhất trong dữ liệu

Thành phần chính đầu tiên (PC1) là hướng có phương sai lớn nhất.
Thành phần chính thứ hai (PC2) nắm bắt nhiều phương sai nhất của phần còn lại, và vuông góc với thành phần thứ nhất.
Nếu hầu hết phương sai chỉ nằm dọc theo một vài hướng, bạn có thể chiếu dữ liệu xuống các chiều đó và bỏ phần còn lại với tổn thất tối thiểu.
Các bước:
- Chuẩn hóa dữ liệu (trừ giá trị trung bình, chia cho độ lệch chuẩn) để tất cả đặc trưng đóng góp như nhau
- Tính ma trận hiệp phương sai
- Tìm các trị riêng và vector riêng của nó
- Chọn \(k\) vector riêng với trị riêng lớn nhất (đây là các thành phần chính)
- Chiếu dữ liệu lên các thành phần này
Chuẩn hóa rất quan trọng: nếu không có nó, một đặc trưng được đo bằng kilomet sẽ chi phối một đặc trưng được đo bằng centimet, bất kể tầm quan trọng thực tế.
Trong thực tế, PCA được dùng để trực quan hóa (chiếu dữ liệu nhiều chiều xuống 2D hoặc 3D), giảm nhiễu (loại bỏ các hướng phương sai thấp chủ yếu là nhiễu), và tăng tốc các mô hình ML bằng cách giảm số lượng đặc trưng đầu vào.
Nhân PCA (Kernel PCA) mở rộng PCA cho các mối quan hệ phi tuyến. Nó ánh xạ dữ liệu qua một hàm nhân vào một không gian nhiều chiều hơn nơi cấu trúc trở thành tuyến tính, sau đó áp dụng PCA chuẩn và chiếu ngược lại.
Phân rã Schur phân tích một ma trận vuông thành \(A = QTQ^\ast\), với \(Q\) là unita và \(T\) là tam giác trên. Mọi ma trận vuông đều có phân rã Schur, ngay cả khi nó không thể chéo hóa được.
Phân tích ma trận không âm (NMF) phân rã một ma trận thành hai ma trận không âm: \(A \approx WH\), với tất cả các phần tử trong \(W\) và \(H\) đều \(\geq 0\). Không giống như SVD có thể tạo ra các phần tử âm, NMF chỉ cộng, không bao giờ trừ. Điều này làm cho các phần có thể giải thích được: trong mô hình hóa chủ đề, \(W\) cho trọng số chủ đề mỗi tài liệu và \(H\) cho trọng số từ mỗi chủ đề, tất cả đều không âm, phù hợp với cách chúng ta nghĩ về "mỗi tài liệu chứa bao nhiêu phần của mỗi chủ đề."
Định lý phổ (spectral theorem) phát biểu rằng các ma trận đối xứng (hoặc Hermite) luôn có thể được chéo hóa bằng một ma trận trực giao (hoặc unita). Trị riêng của chúng luôn thực và các vector riêng của chúng luôn trực giao. Đây là nền tảng lý thuyết đằng sau PCA.

Bài tập lập trình (dùng CoLab hoặc notebook)¶

Tính trị riêng và vector riêng của một ma trận đối xứng. Xác minh các vector riêng vuông góc và tái tạo ma trận từ phân rã trị riêng của nó.

import jax.numpy as jnp

A = jnp.array([[4.0, 2.0],
               [2.0, 3.0]])

eigenvalues, eigenvectors = jnp.linalg.eigh(A)
print(f"Eigenvalues: {eigenvalues}")
print(f"Eigenvectors orthogonal: {jnp.dot(eigenvectors[:,0], eigenvectors[:,1]):.6f}")

# Reconstruct: A = P D P^T
D = jnp.diag(eigenvalues)
A_reconstructed = eigenvectors @ D @ eigenvectors.T
print(f"Reconstruction matches: {jnp.allclose(A, A_reconstructed)}")

Triển khai lặp lũy thừa để tìm trị riêng lớn nhất, và lặp ngược để tìm trị riêng nhỏ nhất. So sánh với jnp.linalg.eigh. Sau đó thử triển khai thuật toán QR.

import jax.numpy as jnp

A = jnp.array([[4.0, 2.0],
               [2.0, 3.0]])

# Power iteration: finds the LARGEST eigenvalue
v = jnp.array([1.0, 0.0])
for _ in range(20):
    v = A @ v
    v = v / jnp.linalg.norm(v)
print(f"Largest eigenvalue:  {v @ A @ v:.4f}")

# Inverse iteration: multiply by A^{-1} instead of A, finds the SMALLEST eigenvalue
v = jnp.array([1.0, 0.0])
for _ in range(20):
    v = jnp.linalg.solve(A, v)
    v = v / jnp.linalg.norm(v)
print(f"Smallest eigenvalue: {1.0 / (v @ jnp.linalg.solve(A, v)):.4f}")

print(f"jnp.linalg.eigh:    {jnp.linalg.eigh(A)[0]}")

Tính SVD của một ma trận, sau đó tái tạo nó chỉ dùng top-k giá trị singular và quan sát chất lượng xấp xỉ thay đổi như thế nào với k.

import jax.numpy as jnp

A = jnp.array([[1.0, 2.0, 3.0],
               [4.0, 5.0, 6.0],
               [7.0, 8.0, 9.0]])

U, S, Vt = jnp.linalg.svd(A)

for k in [1, 2, 3]:
    approx = U[:, :k] @ jnp.diag(S[:k]) @ Vt[:k, :]
    error = jnp.linalg.norm(A - approx)
    print(f"k={k}, reconstruction error: {error:.4f}")