在 Gilbert Strang 教授的 MIT 線性代數課程"Matrix Methods in Data Analysis, Signal Processing, and Machine Learning"的第一講中，我們遇到了一個基本而強大的矩陣分解方法：A=CR 分解。這個分解不僅是線性代數的基石，也是理解矩陣結構與向量空間的關鍵。本文將使用 Google Colab、NumPy 和 Matplotlib 來實現並視覺化這個概念，幫助我們從計算和幾何角度理解 A=CR 分解。
A=CR 分解的基本原理
A=CR 分解是線性代數中的第一個基本分解方法，它表示任何矩陣 A 都可以被分解為兩個特殊矩陣的乘積：一個包含獨立列的矩陣 C 和一個包含獨立行的矩陣 R。
具體來說：
- C 是一個 m×r 矩陣，包含 A 的 r 個線性獨立列，其中 r 是 A 的秩
- R 是一個 r×n 矩陣，包含 A 的簡約列梯式形式 (RREF) 的 r 個非零行
- R 的形式通常為[I F]，其中 I 是 r×r 單位矩陣，F 描述了 A 的依賴列如何由獨立列組合而成
這個分解揭示了線性代數的第一個重要定理：列秩等於行秩。也就是說，矩陣中線性獨立的列數等於線性獨立的行數。
在 Google Colab 上實現 A=CR 分解
首先，我們在 Google Colab 上設置環境並導入必要的庫：
import numpy as np
import matplotlib.pyplot as plt
from numpy.linalg import matrix_rank
from mpl_toolkits.mplot3d import Axes3D
# 設置中文顯示
!wget 'https://noto-website-2.storage.googleapis.com/pkgs/NotoSansCJKtc-hinted.zip'
!mkdir /tmp/fonts
!unzip -o NotoSansCJKtc-hinted.zip -d /tmp/fonts/
!mv /tmp/fonts/NotoSansMonoCJKtc-Regular.otf /usr/share/fonts/truetype/NotoSansMonoCJKtc-Regular.otf -f
!rm -rf /tmp/fonts
!rm NotoSansCJKtc-hinted.zip
import matplotlib.font_manager as font_manager
font_dirs = ['/usr/share/fonts/truetype/']
font_files = font_manager.findSystemFonts(fontpaths=font_dirs)
for font_file in font_files:
font_manager.fontManager.addfont(font_file)
plt.rcParams['font.family'] = "Noto Sans Mono CJK TC"
分解函數實現
我們首先實現一個計算 A=CR 分解的函數：
def cr_decomposition(A):
"""計算矩陣 A 的 CR 分解"""
A = np.array(A, dtype=float)
m, n = A.shape
r = matrix_rank(A)
# 創建用於跟蹤獨立列的列表
independent_cols = []
# 查找 r 個獨立列
temp_matrix = []
for j in range(n):
col = A[:, j]
temp_matrix.append(col)
temp_rank = matrix_rank(np.column_stack(temp_matrix))
if temp_rank > len(independent_cols):
independent_cols.append(j)
if len(independent_cols) == r:
break
# 構建 C 矩陣（包含獨立列）
C = A[:, independent_cols]
# 用最小二乘法解 R 矩陣：A = CR => R = (C^T C)^(-1) C^T A
R = np.linalg.lstsq(C, A, rcond=None)[0]
return C, R, independent_cols使用 Gilbert Strang 的例子
讓我們用 Gilbert Strang 在講座中使用的簡單整數例子來測試我們的分解：
# Gilbert Strang 講座中的例子
A = np.array([
[2, 1, 3],
[3, 1, 4],
[5, 7, 12]
])
print("原始矩陣 A:")
print(A)
# 計算 CR 分解
C, R, independent_cols = cr_decomposition(A)
print("\n矩陣 C（獨立列）:")
print(C)
print("\n矩陣 R:")
print(R)
print("\n獨立列索引：", independent_cols)
# 驗證 CR 乘積是否等於 A
CR = C @ R
print("\nC×R:")
print(CR)
print("\n原始矩陣 A 與 C×R 的差異：")
print(np.round(A - CR, 10)) # 四捨五入以處理浮點誤差
輸出應該顯示 C 包含 A 的前兩列（因為它們是線性獨立的），R 是一個 2×3 矩陣，其形式為[I F]，其中 F 表示第 3 列是前兩列的線性組合。我們可以看到第 3 列確實是第 1 列的 1 倍加上第 2 列，即[3][4][12] = 1*[2][3][5] + 1*[1][1][7]。
A=CR 分解在 3D 向量空間中的幾何意義
A=CR 分解的幾何意義可以通過矩陣的列空間來理解。矩陣 A 的列空間是所有可能的向量 Ax 的集合，其中 x 遍歷所有可能的向量。
在 Gilbert Strang 的例子中，A 是一個 3×3 矩陣，但其秩只有 2，這意味著它的列空間是 R³中的一個平面，而不是整個 R³空間。
讓我們通過可視化來理解這一點：
def visualize_column_space(A):
"""可視化矩陣 A 的列空間"""
# 創建 3D 圖形
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 繪製原點
ax.scatter([0], [0], [0], color='k', s=100, label='原點')
# 提取矩陣的列並繪製
for j in range(A.shape[1]):
col = A[:, j]
ax.quiver(0, 0, 0, col[0], col[1], col[2], color=f'C{j}',
label=f'列 {j+1}: [{col[0]}, {col[1]}, {col[2]}]')
# 隨機生成點在列空間中
num_points = 50
if matrix_rank(A) == 2: # 如果列空間是平面
# 找到兩個獨立列
independent_cols = []
temp_matrix = []
for j in range(A.shape[1]):
col = A[:, j]
temp_matrix.append(col)
temp_rank = matrix_rank(np.column_stack(temp_matrix))
if temp_rank > len(independent_cols):
independent_cols.append(j)
if len(independent_cols) == 2:
break
# 使用這兩個獨立列創建平面上的點
col1 = A[:, independent_cols[0]]
col2 = A[:, independent_cols[1]]
# 創建網格
t = np.linspace(-1, 1, 10)
s = np.linspace(-1, 1, 10)
T, S = np.meshgrid(t, s)
# 計算平面上的點
X = col1[0] * T + col2[0] * S
Y = col1[1] * T + col2[1] * S
Z = col1[2] * T + col2[2] * S
# 繪製平面
ax.plot_surface(X, Y, Z, alpha=0.3, color='cyan')
# 在平面上隨機生成點
coeffs1 = np.random.uniform(-1, 1, num_points)
coeffs2 = np.random.uniform(-1, 1, num_points)
points = np.outer(col1, coeffs1) + np.outer(col2, coeffs2)
ax.scatter(points[0, :], points[1, :], points[2, :], color='blue', alpha=0.6, s=20)
# 設置坐標軸標籤和標題
ax.set_xlabel('X 軸')
ax.set_ylabel('Y 軸')
ax.set_zlabel('Z 軸')
ax.set_title('矩陣的列空間可視化')
# 設置視角
ax.view_init(elev=20, azim=30)
# 添加圖例
ax.legend()
plt.tight_layout()
plt.show()
# 可視化 A 的列空間
visualize_column_space(A)
這個可視化展示了：
- 矩陣 A 的三個列向量
- 由前兩個獨立列向量生成的平面（列空間）
- 第三列向量是如何位於這個平面上，因為它是前兩列的線性組合
透過分解理解平面
A=CR 分解幫助我們理解矩陣 A 的列空間如何成為 R³中的平面：
-
矩陣 C 的列是 A 的列空間的基，在我們的例子中，是兩個線性獨立的向量，它們確定了一個平面
-
矩陣 R，尤其是其[I F]結構，告訴我們 A 的其他列如何由這些基向量表示
-
當我們計算 Ax 時，我們實際上是在計算 C(Rx)，意味著結果必須位於 C 的列空間內
讓我們通過一個簡單的實驗來確認這一點：
def verify_column_space(A, num_points=10):
"""驗證所有 Ax 向量都在列空間中"""
m, n = A.shape
C, R, _ = cr_decomposition(A)
# 生成隨機向量 x
X = np.random.uniform(-1, 1, (n, num_points))
# 計算 Ax 和 CR 分解乘積
AX = A @ X
CRX = C @ (R @ X)
# 計算差異
diff = np.linalg.norm(AX - CRX, axis=0)
print(f"隨機選擇{num_points}個向量 x，計算 Ax 和 C(Rx) 之間的差異：")
print(f"最大差異：{np.max(diff):.10f}")
print(f"平均差異：{np.mean(diff):.10f}")
# 3D 圖表可視化所有結果向量
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
# 繪製原點
ax.scatter([0], [0], [0], color='k', s=100, label='原點')
# 繪製列向量
for j in range(A.shape[1]):
col = A[:, j]
ax.quiver(0, 0, 0, col[0], col[1], col[2], color=f'C{j}',
label=f'列 {j+1}')
# 繪製結果向量 Ax
ax.scatter(AX[0, :], AX[1, :], AX[2, :], color='red', s=50, label='Ax 向量')
# 設置軸標籤和標題
ax.set_xlabel('X 軸')
ax.set_ylabel('Y 軸')
ax.set_zlabel('Z 軸')
ax.set_title('Ax 向量都在列空間（平面）中')
# 添加圖例
ax.legend()
plt.tight_layout()
plt.show()
# 驗證列空間
verify_column_space(A)
這個可視化確認了矩陣 A 的所有可能的輸出 Ax 都位於由其兩個獨立列所定義的平面上。
結論
A=CR 分解是理解矩陣和向量空間的強大工具。通過 Google Colab、NumPy 和 Matplotlib，我們能夠：
- 實現 A=CR 分解的計算
- 可視化矩陣的列空間是 R³中的一個平面
- 確認所有的 Ax 向量都位於這個平面上
這種分解揭示了線性代數的第一個基本定理：列秩等於行秩。它也幫助我們理解為什麼當矩陣的秩小於維度時，矩陣變換的結果會”壓縮”到更低維度的空間（在我們的例子中是平面）。
Gilbert Strang 教授的例子使用簡單整數有助於我們直觀地理解這些概念，而 Python 的實現則讓我們能夠通過計算和可視化來驗證這些理解。
通過本文的實踐，我們可以看到線性代數理論與計算實現的緊密結合，這對於數據分析、信號處理和機器學習等領域的應用至關重要。
