使用 Google Colab 理解線性代數中的 A=CR 分解：Gilbert Strang 講座解析

A=CR 分解的基本原理

Gilbert Strang 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

首先，我們在 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"

Google Colab Setup

分解函數實現

我們首先實現一個計算 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))  # 四捨五入以處理浮點誤差

Gilbert Strang Example

輸出應該顯示 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)