1969-Strassen矩陣乘法#
案件摘要#
1969 年,德國數學家 Volker Strassen 發表了一篇標題即宣言的論文:〈Gaussian elimination is not optimal〉。兩百年來,矩陣乘法 \(O(n^3)\) 被視為顯然的最優——三個巢狀迴圈,一個不多。 Strassen 卻證明:只要把 \(2 \times 2\) 矩陣乘法所需的 8 次標量乘法巧妙地降到 7 次,分治遞迴就會把整體複雜度降到 \(O(n^{\log_2 7}) \approx O(n^{2.807})\) 。 這樁案子的破案時刻不只是一個更快的演算法,而是一個哲學宣言:「顯然最優」也可以被打破。
前因 -- 為什麼會有這個案子#
- 矩陣乘法的 \(O(n^3)\) 定義自高斯消去法的傳統: \(n \times n\) 矩陣相乘需要 \(n^3\) 次標量乘法,自 19 世紀以來無人質疑。
- Strassen 的動機恰恰是想證明「高斯消去法不是最優」——他在研究消除法的複雜性下界時,意外發現上界反而可以更低。
- 分治法在 1960 年代已經成熟:von Neumann 的合併排序(1945)展示 \(O(n \log n)\) ;Karatsuba 的整數乘法(1960)展示 \(O(n^{1.585})\) ,證明「乘法次數」也可以削減,為 Strassen 提供了方法論模板。
線索與推理 -- 數學式、程式、理論#
線索一:Karatsuba 的先聲#
兩個 \(n\) 位數相乘,樸素法需要 \(n^2\) 次單位乘法。Karatsuba(1960)把每數切成兩半 \(x = x_1 2^m + x_0\) 、 \(y = y_1 2^m + y_0\) ,樸素法需要 4 次半長乘法( \(x_1 y_1, x_1 y_0, x_0 y_1, x_0 y_0\) ),Karatsuba 只用 3 次:
遞迴關係:
「少一次乘法、多幾次加法」在漸進意義下大賺——這正是 Strassen 的靈感來源。
線索二:Strassen 的 7 次乘法#
\(2 \times 2\) 矩陣乘法 \(C = AB\) ,樸素分治需要 8 次子矩陣乘法。Strassen 把 \(A\) 、 \(B\) 各分成四塊 \(A = \begin{pmatrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{pmatrix}\) , \(B\) 同理,然後只做 7 次乘法:
結果由加減法組合(無需再乘): \(C_{11} = M_1 + M_4 - M_5 + M_7\) , \(C_{12} = M_3 + M_5\) , \(C_{21} = M_2 + M_4\) , \(C_{22} = M_1 - M_2 + M_3 + M_6\) (注意 \(M_3\) 的符號)。
線索三:主定理的應用#
遞迴關係與主定理(master theorem)的應用: \(a = 7\) , \(b = 2\) , \(f(n) = O(n^2)\) ,且 \(n^{\log_2 7} \approx n^{2.807}\) 支配 \(n^2\) ,故
每省一次標量乘法,指數就下降一截:8 次 ⇒ \(n^3\) ,7 次 ⇒ \(n^{2.807}\) 。
線索四:後續戰場與理論極限#
- Pan(1978)降到 \(O(n^{2.795})\) ,Coppersmith–Winograd(1990)降到 \(O(n^{2.376})\) ;2010 年代的組合改進(Vassilevska Williams、Le Gall 等)壓到 \(O(n^{2.373})\) 左右;理論下界是 \(n^2\) (輸出本身就有 \(n^2\) 個元素),但至今無人達到。
- AlphaTensor(DeepMind,2022)用強化學習搜尋張量分解,為 \(4 \times 4\) 矩陣找到 47 次乘法的方案。
線索五:Python 實作與實測#
以下程式實作 Strassen 矩陣乘法(7 次乘法遞迴),統計乘法次數並對照樸素法的 \(n^3\) :
mul_count = {"strassen": 0}
def mat_mul(A, B):
n = len(A)
return [[sum(A[i][k] * B[k][j] for k in range(n)) for j in range(n)]
for i in range(n)]
def add(A, B):
return [[a + b for a, b in zip(ra, rb)] for ra, rb in zip(A, B)]
def sub(A, B):
return [[a - b for a, b in zip(ra, rb)] for ra, rb in zip(A, B)]
def strassen(A, B):
n = len(A)
if n == 1:
mul_count["strassen"] += 1
return [[A[0][0] * B[0][0]]]
m = n // 2
half = lambda M, i, j: [row[j*m:(j+1)*m] for row in M[i*m:(i+1)*m]]
A11, A12, A21, A22 = (half(A, i, j) for i in (0, 1) for j in (0, 1))
B11, B12, B21, B22 = (half(B, i, j) for i in (0, 1) for j in (0, 1))
M1 = strassen(add(A11, A22), add(B11, B22))
M2 = strassen(add(A21, A22), B11)
M3 = strassen(A11, sub(B12, B22))
M4 = strassen(A22, sub(B21, B11))
M5 = strassen(add(A11, A12), B22)
M6 = strassen(sub(A21, A11), add(B11, B12))
M7 = strassen(sub(A12, A22), add(B21, B22))
C11 = add(sub(add(M1, M4), M5), M7)
C12 = add(M3, M5)
C21 = add(M2, M4)
C22 = add(add(sub(M1, M2), M3), M6)
return [r1 + r2 for r1, r2 in zip(C11, C12)] + \
[r1 + r2 for r1, r2 in zip(C21, C22)]
n = 16
A = [[i * n + j + 1 for j in range(n)] for i in range(n)]
B = [[(i + j) % 7 + 1 for j in range(n)] for i in range(n)]
print("results equal:", mat_mul(A, B) == strassen(A, B))
print(f"naive: {n**3:>6} multiplications (theory)")
print(f"Strassen: {mul_count['strassen']:>6} multiplications (counted)")
理論預測: \(n = 16\) 時樸素法做 \(16^3 = 4096\) 次乘法,Strassen 做 \(7^{\log_2 16} = 7^4 = 2401\) 次。執行程式驗證結果相等且次數吻合。
結案 -- 後果與影響#
- 第一次證明「顯然最優」可以打破: \(O(n^3)\) 不是矩陣乘法的宿命,也動搖了其他「顯然」下界的信心。
- 分治遞迴的教科書案例: \(T(n) = 7T(n/2) + O(n^2)\) 成為主定理教學的標準例題。
- 矩陣乘法指數成為獨立的研究戰場:從 Strassen 的 2.807 到 Coppersmith–Winograd 的 2.376,再到 2024 年理論極限 \(n^2\) 仍未達成;實務上 BLAS 等函式庫也採用了它的變體( \(n > 10^3\) 時可勝過樸素法)。
- Strassen 之後轉向機率演算法:Solovay–Strassen 質數測試(1977)是 RSA 時代機率質數測試的先驅。
關鍵人物與文獻#
- Volker Strassen:德國數學家(波昂大學、蘇黎世大學),矩陣乘法與機率質數測試的先驅。
- Anatoly Karatsuba:蘇聯數學家,1960 年發明分治整數乘法;Don Coppersmith 與 Shmuel Winograd(IBM)於 1990 年將指數壓到 2.376。
主要文獻:
- Strassen, V. (1969). "Gaussian elimination is not optimal". Numerische Mathematik, 13(4), 354–356.
- Karatsuba, A., & Ofman, Y. (1962). "Multiplication of Many-Digital Numbers by Automatic Computers". Doklady Akademii Nauk SSSR, 145, 293–294.
- Coppersmith, D., & Winograd, S. (1990). "Matrix Multiplication via Arithmetic Progressions". Journal of Symbolic Computation, 9(3), 251–280.
- Fawzi, A., et al. (2022). "Discovering Faster Matrix Multiplication Algorithms with Reinforcement Learning". Nature, 610, 47–53.(AlphaTensor)
- Solovay, R., & Strassen, V. (1977). "A Fast Monte-Carlo Test for Primality". SIAM Journal on Computing, 6(1), 84–85.