9.3 矩陣運算複雜度#

本節目標#

掌握各種矩陣運算的計算複雜度，知道什麼操作「貴」、什麼「便宜」
實測矩陣乘法的 \(O(n^3)\) 成長曲線（log-log 圖）
估算大矩陣的記憶體需求（含「1 億參數模型要幾 GB」）
實測你的電腦每秒能做多少 FLOPS

為什麼要關心複雜度？#

線性代數的數學性質沒有成本概念，但執行有。寫迴圈版本和向量化版本可能差 100 倍速度；選 inv 還是 solve 可能差 3 倍；矩陣從 \(1000 \times 1000\) 變 \(10000 \times 10000\) ，乘法時間不是慢 10 倍，而是慢 1000 倍。本章用實測取代背誦。

矩陣乘法是 \(O(n^3)\)#

\(C = AB\) ，其中每個 \(c_{ij} = \sum_k a_{ik} b_{kj}\) 。 \(C\) 有 \(n^2\) 個元素，每個要 \(n\) 次乘加，總計：

\[ \text{運算次數} \approx 2n^3 \quad (\text{n 次乘 + n 次加，共 } n^2 \text{ 個元素}) \]

理論講完，直接實測。規則： \(n\) 變 2 倍，時間變 8 倍：

import numpy as np
import timeit
import matplotlib.pyplot as plt

sizes = [128, 256, 512, 1024, 2048]
times = []

for n in sizes:
    A = np.random.rand(n, n)
    B = np.random.rand(n, n)
    t = timeit.timeit(lambda: A @ B, number=3) / 3
    times.append(t)
    print(f"n={n:5d}  時間={t:.4f}s  每倍增比={t/times[-2]:.1f}" if times[:-1] else f"n={n:5d}  時間={t:.4f}s")

# log-log 圖：O(n^p) 在 log-log 下是斜率 p 的直線
plt.loglog(sizes, times, "o-", label="matmul time")
plt.loglog(sizes, [times[0] * (n / sizes[0]) ** 3 for n in sizes],
           "--", label="slope = 3 reference")
plt.xlabel("n (log)")
plt.ylabel("seconds (log)")
plt.legend()
plt.title("矩陣乘法成長曲線：斜率應接近 3")
plt.show()

典型結果（桌機 CPU）：

n=  128  時間=0.0012s
n=  256  時間=0.0078s   (~6.5x)
n=  512  時間=0.0550s   (~7.0x)
n=1024  時間=0.4200s   (~7.6x)
n=2048  時間=3.1500s   (~7.5x)

倍增比接近 8 = \(2^3\) ，log-log 圖的斜率接近 3 —— 理論與實測吻合。（比值略低於 8 是因為大矩陣更能吃滿快取與 SIMD，BLAS 對此做了優化。）

各種運算的複雜度表#

運算	NumPy 呼叫	複雜度	備註
矩陣加法	`A + B`	\(O(n^2)\)	便宜
矩陣 × 向量	`A @ x`	\(O(n^2)\)	便宜，神經網路推論的主體
矩陣 × 矩陣	`A @ B`	\(O(n^3)\)	一切的重頭戲
解 \(A\mathbf{x}=\mathbf{b}\)	`np.linalg.solve`	\(\approx \tfrac{2}{3}n^3\)	LU 分解，比 `inv` 便宜 3 倍
求逆	`np.linalg.inv`	\(\approx 2n^3\)	又貴又（對條件數）不穩
行列式	`np.linalg.det`	\(\approx \tfrac{2}{3}n^3\)	基於 LU
SVD	`np.linalg.svd`	\(\approx O(mn^2)\) （ \(m \geq n\) ）	最貴但最穩
特徵值	`np.linalg.eig`	\(\approx O(n^3)\) 且常數大	對稱矩陣用 `eigvalsh` 更快

實測驗證「solve 比 inv 便宜」：

n = 2000
A = np.random.rand(n, n) + n * np.eye(n)   # 對角優勢，健康矩陣
b = np.random.rand(n)

t_inv = timeit.timeit(lambda: np.linalg.inv(A) @ b, number=3) / 3
t_solve = timeit.timeit(lambda: np.linalg.solve(A, b), number=3) / 3
print(f"inv @ b : {t_inv:.3f}s")
print(f"solve   : {t_solve:.3f}s   快 {t_inv/t_solve:.1f} 倍，而且誤差更小")

Strassen 演算法：一句話版本#

1969 年 Strassen 證明：把矩陣分成 \(2\times2\) 塊、用 7 次（而非 8 次）子矩陣乘法合併，複雜度降到 \(O(n^{\log_2 7}) \approx O(n^{2.81})\) ，且遞迴下去還能更低（理論最佳紀錄低於 \(n^{2.37}\) ）。但實務上 NumPy 不用它——常數因子大、快取不友善、數值穩定性較差，現代硬體上精調的 BLAS（ \(O(n^3)\) ）反而贏。知道它的存在，理解「理論複雜度 ≠ 實際速度」即可。

記憶體估算：1 億參數要多少 GB？#

時間之外，記憶體常是更早撞上的牆。float64 每個數 8 bytes、float32 4 bytes：

def gb(n, dtype=np.float64):
    return n * np.dtype(dtype).itemsize / 1024**3

# n x n 矩陣
for n in [1000, 5000, 10000, 20000]:
    print(f"{n}x{n} float64: {gb(n*n):6.1f} GB")

# 神經網路的重量級場景：1 億參數（1e8）
print(f"1e8 參數 float32: {gb(1e8, np.float32):.2f} GB")
print(f"1e8 參數 float64: {gb(1e8):.2f} GB")
print(f"1e8 參數 int1(bit): {1e8/8/1024**3:.4f} GB")

結果速查：

項目	記憶體
\(1000 \times 1000\) float64	0.007 GB（7.5 MB）
\(10000 \times 10000\) float64	0.75 GB
\(20000 \times 20000\) float64	3.0 GB
1 億參數模型（float32）	0.37 GB
1 億參數模型（float64）	0.75 GB

更狠的是：訓練時除了參數本身，還要存梯度和最佳化器狀態，實際占用常常是參數量的 3~4 倍——這就是大型模型需要多卡 GPU 的直接原因。用 np.zeros((20000, 20000)) 就能讓你自己的筆電立刻體會什麼叫記憶體壓力。

實用函式：A.nbytes 直接告訴你一個陣列占多少 bytes。

實驗：你的電腦每秒能做多少 FLOPS？#

FLOPS（每秒浮點運算次數）是硬體效能的通用量尺。矩陣乘法是標準基準： \(2n^3\) FLOPs / 實測秒數：

import numpy as np
import timeit

def measure_gflops(n=2048, repeats=3):
    A = np.random.rand(n, n)
    B = np.random.rand(n, n)
    flops = 2 * n**3
    t = timeit.timeit(lambda: A @ B, number=repeats) / repeats
    return flops / t / 1e9   # GFLOPS

print(f"float64 矩陣乘法: {measure_gflops(2048):.1f} GFLOPS")

# 對照：NumPy 的逐元素純 Python 端運算慢很多
x = np.random.rand(10_000_000)
t = timeit.timeit(lambda: np.sqrt(x), number=5) / 5
elem_flops = 1e7 / t   # sqrt 算 1 flop/element
print(f"逐元素 sqrt     : {elem_flops/1e9:.2f} GFLOPS")

# 再對照：純 Python 迴圈的慘況
import math
lst = x[:100_000].tolist()
t = timeit.timeit(lambda: [math.sqrt(v) for v in lst], number=3) / 3
print(f"Python 迴圈 sqrt: {100_000/t/1e6:.1f} MFLOPS（慢了約百倍）")

典型結果：矩陣乘法 30~200 GFLOPS（有 BLAS + 多核心 + SIMD加持），逐元素運算 1~5 GFLOPS，純 Python 迴圈只有 5~50 MFLOPS。 同樣一秒鐘，BLAS 矩陣乘法完成的運算是 Python 迴圈的數千倍—— 這就是下一節「向量化」要談的主題。

有 GPU 的讀者可以對照：消費級 GPU 的矩陣乘法動輒數十 TFLOPS，是 CPU 的百倍以上。深度學習硬體軍備競賽的根源就在這裡。

小結與練習#

矩陣乘法 \(O(n^3)\) ： \(n\) 倍增，時間 ×8；實測 log-log 斜率 ≈ 3
解方程 \(\tfrac{2}{3}n^3\) 比 inv 的 \(2n^3\) 便宜，且更穩定——優先用 solve / lstsq
SVD 與特徵值最貴；對稱問題用 eigvalsh 省時
記憶體：float64 每數 8 bytes； \(20000^2\) 矩陣 3 GB；1 億參數 float32 約 0.4 GB
你的 CPU 跑 BLAS 矩陣乘法可達數十 GFLOPS，純 Python 迴圈只有 MFLOPS 級

練習：

把 sizes 延伸到 4096，確認成長曲線斜率仍是 3；若時間突然暴增一個數量級，猜猜是不是矩陣大到超出快取/記憶體頻寬。
實測 A @ x（ \(O(n^2)\) ）在不同 \(n\) 的時間成長，驗證斜率 ≈ 2。
用 A.nbytes 計算你能開出最大的矩陣是多大（依你的 RAM），並實測它的乘法時間。
比較 np.linalg.eig 與 np.linalg.eigvalsh（對稱矩陣）在 \(n=1000\) 的時間差。