9.4 向量化效能#

本節目標#

理解 NumPy 為什麼快：C 底層、SIMD、記憶體連續性、BLAS
實測 Python 迴圈 vs 向量化的速度差距，畫出速度比
學會用 broadcasting 消除迴圈（距離矩陣重構實例）
建立撰寫高效 NumPy 的準則清單，並重構自己的舊程式碼

NumPy 快的四個理由#

C 底層：np.sum(a) 的迴圈在編譯好的 C 程式碼裡跑，沒有 Python 每圈都要「建立整數物件、查型別、分派方法」的開銷。
SIMD：現代 CPU 有一條指令處理多筆資料（AVX2 一次算 4 個 float64）， NumPy/BLAS 會自動使用。
記憶體連續：NumPy 陣列是一整塊連續記憶體，CPU 快取預取器可以一口氣搬； Python list 是「指標的陣列」，資料散落各處。
BLAS：矩陣乘法等線性代數操作呼叫底層 BLAS 函式庫（OpenBLAS、MKL、Accelerate），數十年最佳化 + 多核心並行。

記住一個量級感：Python 迴圈是 MFLOPS 級，NumPy 向量化是 GFLOPS 級，BLAS 矩陣乘法可達數十 GFLOPS（上節實測過）。差距是 100~1000 倍。

實驗一：迴圈 vs np.sum#

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

sizes = [10**3, 10**4, 10**5, 10**6, 10**7]
t_loop, t_vec = [], []

for n in sizes:
    a = np.random.rand(n)
    t_loop.append(timeit.timeit(lambda: sum(a.tolist()), number=3) / 3)
    t_vec.append(timeit.timeit(lambda: a.sum(), number=3) / 3)
    print(f"n={n:>8}  迴圈={t_loop[-1]*1e3:8.2f}ms  向量={t_vec[-1]*1e3:8.3f}ms  快 {t_loop[-1]/t_vec[-1]:5.0f} 倍")

plt.loglog(sizes, t_loop, "o-", label="Python loop")
plt.loglog(sizes, t_vec, "s-", label="np.sum")
plt.xlabel("n"); plt.ylabel("seconds (log)")
plt.legend(); plt.title("迴圈 vs 向量化：兩條平行線，差距約 100 倍")
plt.show()

典型輸出：

n=    1000  迴圈=    0.05ms  向量=    0.006ms  快     9 倍
n=   10000  迴圈=    0.45ms  向量=    0.008ms  快    56 倍
n=  100000  迴圈=    4.5 ms  向量=    0.05 ms  快    90 倍
n= 1000000  迴圈=   45   ms  向量=    0.5  ms  快    90 倍
n=10000000  迴圈=  450   ms  向量=    5    ms  快    90 倍

log-log 圖上兩條線平行（都是線性掃描，斜率 1），只是相差約兩個數量級。注意 sum(a.tolist()) 甚至比 sum(a) 更公平（先把 list 化，避免對 NumPy 陣列逐項迭代反而更慢的干擾）。

實驗二：broadcasting 消除迴圈——距離矩陣#

經典任務：給 \(n\) 個二維點，算所有點對之間的距離矩陣 \(D_{ij} = \lVert \mathbf{p}_i - \mathbf{p}_j \rVert\) 。

初學者版（雙重迴圈）：

def dist_loops(P):
    n = len(P)
    D = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            D[i, j] = np.sqrt(np.sum((P[i] - P[j])**2))
    return D

broadcasting 版：

def dist_vec(P):
    diff = P[:, None, :] - P[None, :, :]   # (n,1,2) - (1,n,2) -> (n,n,2)
    return np.sqrt((diff**2).sum(axis=-1))

# 驗證正確性
rng = np.random.default_rng(0)
P = rng.standard_normal((50, 2))
print(np.allclose(dist_loops(P), dist_vec(P)))   # True

# 效能對決（n=300 時雙重迴圈要跑 9 萬次 Python 迴圈）
P = rng.standard_normal((300, 2))
import timeit
print(timeit.timeit(lambda: dist_loops(P), number=3) / 3)   # ~1.5 s
print(timeit.timeit(lambda: dist_vec(P),    number=3) / 3)  # ~0.01 s，快百倍

broadcasting 的規則：比較兩個形狀時，由尾端對齊，各維度要嘛相等、要嘛其中一邊是 1（可複製擴展）、要嘛其中一邊不存在。 (n,1,2) 與 (1,n,2) 對齊後廣播成 (n,n,2)，一次做完所有點對的減法。 外層沒有任何 Python 迴圈，全部落到 C。

還有更省記憶體的進階版（不配置 \(n \times n \times 2\) 的中間陣列）：

def dist_fast(P):
    sq = (P**2).sum(axis=1)
    d2 = sq[:, None] + sq[None, :] - 2 * P @ P.T   # 展開 |a-b|^2 = |a|^2+|b|^2-2a·b
    np.maximum(d2, 0, out=d2)                      # 浮點誤差可能產生微小負值！
    return np.sqrt(d2)

這裡順便複習 9.1： \(d^2\) 理論上非負，但浮點誤差會讓它出現 \(-10^{-17}\) ， sqrt 會給 nan —— 所以先 np.maximum(d2, 0) 夾住。效能優化與浮點誠實經常一起出現。

np.dot 與 BLAS#

矩陣操作盡量走 BLAS 通道：

A = np.random.rand(1000, 1000)
x = np.random.rand(1000)

# 這些都會呼叫 BLAS：
y = A @ x            # gemv：矩陣乘向量
C = A @ A            # gemm：矩陣乘矩陣
s = x @ x            # dot：內積

# 這些不是 BLAS，是逐元素廣播，慢於 gemm 但仍遠快於迴圈：
row_dot = (A * x).sum(axis=1)

形狀小知識：np.dot(x, y) 對一維是內積；x @ y 同義且可讀性更好。 A @ B 對二維就是矩陣乘法，背後是 dgemm。

%timeit 用法#

在 Jupyter 裡量測效能的標配：

# In [1]:
%timeit np.sum(np.random.rand(1_000_000))
# 輸出範例：1.23 ms ± 45.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# In [2]: 只跑一次的大操作用 %%time
%%time
A = np.random.rand(3000, 3000)
B = A @ A

%timeit 會自動重複執行多次取平均與標準差，比手刻 time.time() 可靠得多。腳本裡則用 timeit.timeit(lambda: f(), number=...)。

高效 NumPy 準則清單#

禁止 Python 層的 for/while 遍歷陣列元素——先想 broadcasting 能不能做。
盡量讓資料以「整塊陣列運算」流動：a.sum()、a.max(axis=1)、a ** 2。
矩陣乘法交給 @，不要手寫內積迴圈。
中間陣列別太濫用：a - b 再 **2 再 sum 會配置多個臨時陣列；熱點程式碼可用 np.einsum 或 out= 參數省配置。
用對 dtype：float32 讓記憶體減半、SIMD 吞吐加倍（精度夠用時）。
保持記憶體連續：轉置後馬上做大運算前考慮 np.ascontiguousarray。
先測量、再優化：%timeit 找熱點，別憑感覺優化沒人跑的程式碼。
組合技：向量化讓單核快 100 倍，還能再搭 numba / Cython / GPU 上萬倍。

實驗三：重構自己的舊程式碼#

拿一個「你一定寫過」的典型初學者函式——用迴圈算每個資料點到中心的距離並分類：

import numpy as np
import timeit

rng = np.random.default_rng(1)
X = rng.standard_normal((100_000, 3))          # 十萬筆三維資料
center = np.array([0.5, -0.2, 0.1])

# ── 重構前：Python 迴圈思維 ──
def assign_loops(X, center):
    labels = np.empty(len(X), dtype=int)
    for i in range(len(X)):
        d = 0.0
        for k in range(3):
            d += (X[i, k] - center[k])**2
        labels[i] = 0 if d < 1.0 else 1
    return labels

# ── 重構後：broadcasting 思維 ──
def assign_vec(X, center):
    d2 = ((X - center) ** 2).sum(axis=1)
    return (d2 >= 1.0).astype(int)

assert np.array_equal(assign_loops(X, center), assign_vec(X, center))
print(timeit.timeit(lambda: assign_loops(X, center), number=3) / 3)  # ~0.5 s
print(timeit.timeit(lambda: assign_vec(X, center),    number=3) / 3) # ~0.001 s

# 進階挑戰：把它重構成 k 個中心的版本（k-means 的一小步）
centers = rng.standard_normal((4, 3))
# 距離矩陣 (100000, 4)：X[:,None,:] - centers[None,:,:] 會占 100000*4*3*8B ≈ 9.6 MB，可接受
d2 = ((X[:, None, :] - centers[None, :, :]) ** 2).sum(axis=2)
labels = d2.argmin(axis=1)
print(np.bincount(labels))    # 各群數量 —— 這就是最近中心分類

重構的心法只有一句話：把「對每個元素做 X」翻譯成「對整個陣列做 X，讓形狀說話」。翻譯訣竅是先寫下輸出形狀，再倒推需要什麼 broadcasting。

小結與練習#

NumPy 快的原因：C 迴圈、SIMD、連續記憶體、BLAS；量級差距 100~1000 倍
迴圈 vs np.sum：log-log 圖兩條平行線，恆定約兩個數量級差距
broadcasting 消除迴圈：(n,1,d) - (1,n,d) 一次算出所有點對差值
效能優化常伴隨浮點陷阱：d2 可能出現微小負值，記得夾非負
準則：先向量化、矩陣乘交給 @、用 %timeit 測量後再優化

練習：

把 dist_fast 與 dist_vec 在 \(n = 100, 500, 2000\) 比較時間與記憶體用量（用 A.nbytes 估算中間陣列大小）。
實測 float32 版的 A @ A 比 float64 快多少？快的原因是什麼？
重構：給一個用雙重迴圈計算兩矩陣乘積的函式，改成 @ 並量測 \(n = 200\) 的加速比。
用 %timeit 找出 np.einsum('ij,ij->i', X, X) 與 (X**2).sum(axis=1) 誰快，並解釋為什麼 einsum 可能省下中間陣列配置。