为何Numpy向量化实现蒙特卡洛求π比for循环更慢？

我实现了两个蒙特卡洛估算π的程序，原本预期纯Numpy向量化的method_1会比带for循环的method_2更快，但测试结果却相反，想请教背后的原因。

测试代码 #

import numpy as np
import time

def method_1():
    start = time.time()
    
    N = 100_000 #Number of points
    total = 1/N
    N_run = 100 #Number of runs
    pi = np.zeros(N_run)
    x = np.random.rand(N, N_run)
    y = np.random.rand(N, N_run)

    pi = np.sum(x*x + y*y < 1, axis = 0)*4*total
    
    stop = time.time()
    
    return stop-start


def method_2():
    start = time.time()

    N = 100_000 #Number of points
    total = 1/N
    N_run = 100 #Number of runs
    pi = np.zeros(N_run)
    x = np.random.rand(N)
    y = np.random.rand(N)
    
    for i in range(N_run):
        pi[i] = np.sum(x*x + y*y<1)*4*total
        x = np.random.rand(N)
        y = np.random.rand(N)
        
    stop = time.time()
    
    return stop-start



N = 10
m1=np.zeros(N)
m2=np.zeros(N)

for i in range(N):
    m1[i] = method_1() 
    m2[i] = method_2()

print(f'Method 1: {np.mean(m1)} sec, Method 2: {np.mean(m2)} sec')

测试输出 #

Method 1: 0.16888238191604615 sec, Method 2: 0.13313182592391967 sec

原因分析 #

• 内存缓存效率差异：method_1一次性生成两个100000×100的二维数组，总共有2000万浮点数，内存占用约1.6GB（每个float64占8字节），远远超出CPU缓存的容量（通常几MB到几十MB）。这会导致大量缓存失效，数据需要从内存中频繁读取，而内存访问速度比CPU缓存慢几个数量级，成为性能瓶颈。而method_2每次只生成两个10万元素的一维数组，内存仅约1.6MB，完全能被CPU缓存容纳，数据访问效率极高。
• 随机数生成的开销：Numpy生成大二维数组时，内存分配和初始化的开销比多次生成小一维数组更大。一次性生成大量随机数会触发更复杂的内存管理逻辑，而小批量生成的内存操作更轻量。
• 内存访问模式：method_1中计算x*x + y*y < 1并按列求和时，由于Numpy数组是行优先存储，按列访问会导致不连续的内存读取，进一步降低缓存命中率。而method_2的一维数组是连续内存访问，缓存命中率接近100%，计算效率更高。

内容的提问来源于stack exchange，提问作者TommyGiak