矩陣乘法演算法全解析：從朴素 O(n³) 到 Strassen O(n²·⁸¹) 與分治遞歸實作

前言

矩陣乘法是線性代數核心運算，廣泛應用於圖形學、機器學習、科學計算等領域。本文系統梳理三種經典演算法：一般矩陣乘法（O(n³））、傳統分治法（仍為 O(n³））、Strassen 演算法（O(n^log₂7) ≈ O(n²·⁸¹），並提供完整 C++ 實作碼與複雜度分析，助你掌握從理論到工程的完整鏈路。

一、一般矩陣乘法（朴素演算法）

定義與條件

設 A 為 m×n 矩陣，B 為 n×p 矩陣，乘積 C = AB 為 m×p 矩陣。元素計算公式：

核心代碼：

struct Mat { int n, m; double data[MAXN][MAXN]; };

int mul(Mat& c, const Mat& a, const Mat& b) {
    if (a.m != b.n) return 0;
    c.n = a.n; c.m = b.m;
    for (int i = 0; i < c.n; ++i)
        for (int j = 0; j < c.m; ++j)
            for (c.data[i][j] = 0, int k = 0; k < a.m; ++k)
                c.data[i][j] += a.data[i][k] * b.data[k][j];
    return 1;
}

時間複雜度： O(n³)

二、傳統分治法

分塊策略

將 n×n 方陣四等分為 n/2×n/2 子矩陣：

遞歸關係：

C₁₁ = A₁₁B₁₁ + A₁₂B₂₁
C₁₂ = A₁₁B₁₂ + A₁₂B₂₂
C₂₁ = A₂₁B₁₁ + A₂₂B₂₁
C₂₂ = A₂₁B₁₂ + A₂₂B₂₂

每層遞歸需 8 次子矩陣乘法 + 4 次矩陣加法，遞歸式：

T(2) = b
T(n) = 8T(n/2) + cn²  (n > 2)

解得 T(n) = O(n³)，未降低複雜度，但展示了分治思想：將大問題分解為相似子問題，再合併結果。

分治法適用問題四大特徵

1. 可縮性：規模縮小至可直接求解
2. 最優子結構：可分解為相同子問題
3. 可合性：子問題解可合併為原問題解
4. 獨立性：子問題間無公共子子問題

三、Strassen 演算法（突破 O(n³)）

核心創新：2×2 矩陣乘法僅需 7 次乘法

傳統 2×2 需 8 次乘法，Strassen 構造 7 個中間矩陣：

M₁ = A₁₁(B₁₂ - B₂₂)
M₂ = (A₁₁ + A₁₂)B₂₂
M₃ = (A₂₁ + A₂₂)B₁₁
M₄ = A₂₂(B₂₁ - B₁₁)
M₅ = (A₁₁ + A₂₂)(B₁₁ + B₂₂)
M₆ = (A₁₂ - A₂₂)(B₂₁ + B₂₂)
M₇ = (A₁₁ - A₂₁)(B₁₁ + B₁₂)

再以加減法組合得結果：

C₁₁ = M₅ + M₄ - M₂ + M₆
C₁₂ = M₁ + M₂
C₂₁ = M₃ + M₄
C₂₂ = M₅ + M₁ - M₃ - M₇

複雜度分析

遞歸式：

T(2) = b
T(n) = 7T(n/2) + cn²  (n > 2)

套用主定理：T(n) = O(n^log₂7) ≈ O(n²·⁸¹)，實現階級提升。

四、工程實作細節

資料結構選擇

• 二維陣列：簡單但需宏定義常數 N
• 二重指標：動態分配，擺脫固定大小限制

// 二重指標動態分配
int **c = (int**)malloc(n * sizeof(int*));
for (int i = 0; i < n; ++i)
    c[i] = (int*)malloc(n * sizeof(int));

核心輔助函數

// 矩陣加法
void add(int n, int A[][N], int B[][N], int R[][N]) {
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            R[i][j] = A[i][j] + B[i][j];
}

// 矩陣減法
void sub(int n, int A[][N], int B[][N], int R[][N]) {
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            R[i][j] = A[i][j] - B[i][j];
}

// 2×2 基礎乘法（Strassen 核心）
void multiply(int A[][N], int B[][N], int R[][N]) {
    int M1 = A[0][1] * (B[0][1] - B[1][1]);
    int M2 = (A[0][0] + A[0][1]) * B[1][1];
    int M3 = (A[1][0] + A[1][1]) * B[0][0];
    int M4 = A[1][1] * (B[1][0] - B[0][0]);
    int M5 = (A[0][0] + A[1][1]) * (B[0][0] + B[1][1]);
    int M6 = (A[0][1] - A[1][1]) * (B[1][0] + B[1][1]);
    int M7 = (A[0][0] - A[1][0]) * (B[0][0] + B[0][1]);
    R[0][0] = M5 + M4 - M2 + M6;
    R[0][1] = M1 + M2;
    R[1][0] = M3 + M4;
    R[1][1] = M5 + M1 - M3 - M7;
}

// 分塊函數
void divide(int A[][N], int n,
            int X11[][N], int X12[][N],
            int X21[][N], int X22[][N]) {
    int half = n / 2;
    for (int i = 0; i < half; ++i)
        for (int j = 0; j < half; ++j) {
            X11[i][j] = A[i][j];
            X12[i][j] = A[i][j + half];
            X21[i][j] = A[i + half][j];
            X22[i][j] = A[i + half][j + half];
        }
}

完整 Strassen 遞歸實作

void strassen(int n, int A[][N], int B[][N], int C[][N]) {
    if (n == 2) { multiply(A, B, C); return; }

    // 宣告子矩陣
    int A11[N][N], A12[N][N], A21[N][N], A22[N][N];
    int B11[N][N], B12[N][N], B21[N][N], B22[N][N];
    int C11[N][N], C12[N][N], C21[N][N], C22[N][N];
    int AA[N][N], BB[N][N];
    int M1[N][N], M2[N][N], M3[N][N], M4[N][N];
    int M5[N][N], M6[N][N], M7[N][N];

    int half = n / 2;
    // 分塊
    for (int i = 0; i < half; ++i)
        for (int j = 0; j < half; ++j) {
            A11[i][j] = A[i][j];           A12[i][j] = A[i][j + half];
            A21[i][j] = A[i + half][j];    A22[i][j] = A[i + half][j + half];
            B11[i][j] = B[i][j];           B12[i][j] = B[i][j + half];
            B21[i][j] = B[i + half][j];    B22[i][j] = B[i + half][j + half];
        }

    // 7 次遞歸乘法
    sub(half, B12, B22, BB);     strassen(half, A11, BB, M1);
    add(half, A11, A12, AA);     strassen(half, AA, B22, M2);
    add(half, A21, A22, AA);     strassen(half, AA, B11, M3);
    sub(half, B21, B11, BB);     strassen(half, A22, BB, M4);
    add(half, A11, A22, AA);
    add(half, B11, B22, BB);     strassen(half, AA, BB, M5);
    sub(half, A12, A22, AA);
    add(half, B21, B22, BB);     strassen(half, AA, BB, M6);
    sub(half, A11, A21, AA);
    add(half, B11, B12, BB);     strassen(half, AA, BB, M7);

    // 合併結果
    add(half, M5, M4, AA); sub(half, M6, M2, BB); add(half, AA, BB, C11);
    add(half, M1, M2, C12);
    add(half, M3, M4, C21);
    sub(half, M5, M3, AA); sub(half, M1, M7, BB); add(half, AA, BB, C22);

    // 寫回 C
    for (int i = 0; i < half; ++i)
        for (int j = 0; j < half; ++j) {
            C[i][j] = C11[i][j];
            C[i][j + half] = C12[i][j];
            C[i + half][j] = C21[i][j];
            C[i + half][j + half] = C22[i][j];
        }
}

五、手算矩陣乘法步驟範例

2×2 乘 2×3 → 2×3 結果：

|a b|   |e f g|   |ae+bh af+bi ag+bk|
|c d| × |h i k| = |ce+dh cf+di cg+dk|

規則：左矩陣第 i 列 × 右矩陣第 j 行對應元素相乘求和 → 結果第 i 列第 j 行元素。

六、總結與選型建議

現代高效能計算（如 OpenBLAS、Intel MKL）採用分塊矩陣乘法 + SIMD 向量化 + 多執行緒 + 緩存優化，實際效能遠超理論複雜度較低的 Strassen。理解演算法演進史，有助於在架構選型時做出權衡。