更改重塑顺序时Sobol向量的NumPy矩阵向量积性能问题及解决

问题背景 #

我在实现布朗桥路径构建时，希望通过NumPy矩阵乘法加速计算。具体场景：

• 3个标的资产，256个时间步长，生成2^16条路径
• 先抽取768维（3×256）的2^16个Sobol数，得到二维数组
• 将其重塑为形状(3,256,2^16)的三维数组，为了让布朗桥最大方差贡献对应首维度，使用order='F'参数
• 但发现用order='F'重塑后的数组执行矩阵积@时，速度比默认C顺序慢近100倍，完全抵消了矩阵乘法的性能优势

复现代码 #

from timeit import default_timer
from scipy.stats.qmc import Sobol
import numpy as np

np_gen = np.random.default_rng(seed=42)
matrix = np_gen.random((256, 256))
sobol_gen = Sobol(3*256, bits=64, seed=42)
sobol_nbrs = sobol_gen.random_base2(16).T

sobol_c_reshape = sobol_nbrs.reshape((3, 256, 2 ** 16), order='C')
sobol_f_reshape = sobol_nbrs.reshape((3, 256, 2 ** 16), order='F')

def perf_measure(input_nbrs):
    start = default_timer()
    res = matrix @ input_nbrs
    end = default_timer()
    print(end-start)

perf_measure(sobol_c_reshape)
perf_measure(sobol_f_reshape)

该问题并非Sobol序列特有，用随机数组也可复现：

sobol_nbrs = np.random.rand(65536, 768).T

性能差异的核心原因 #

NumPy数组的内存布局分为两种：

• C顺序（行优先）：内存中连续存储数组的每行元素，是NumPy默认的布局
• F顺序（列优先）：内存中连续存储数组的每列元素

矩阵乘法的底层实现（依赖BLAS等线性代数库）高度依赖内存连续性。当操作非连续内存的数组时，会频繁触发缓存未命中——CPU无法从高速缓存中读取所需数据，不得不从内存中加载，导致计算效率暴跌。

reshape(order='F')只是改变了数组的视图，并没有重新排列内存中的数据，因此sobol_f_reshape的内存是非连续的；而sobol_c_reshape是连续内存布局，符合BLAS库的优化预期，所以速度差距巨大。

解决方法 #

将重塑后的F顺序数组复制为C顺序的连续数组，让BLAS库可以高效利用缓存：

matrix @ sobol_f_reshape.copy(order='C')

复制操作会重新排列内存数据，生成连续的C顺序数组。虽然复制会带来少量额外开销，但相比矩阵乘法的性能提升可以忽略，最终计算速度仅比原生C顺序慢2倍左右。

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