Python中基于interp1d插值函数快速构建二维数组的优化方案咨询

• 核心思路是用NumPy的广播机制替代Python级嵌套循环，所有运算都在底层C层面执行，相比原循环速度可以提升几个数量级，不需要引入额外依赖。你之前使用numba优化失效的原因是numba无法编译scipy的interp1d函数实现，向量化方案不需要依赖numba即可获得极高的性能。

基础向量化实现代码 #

import numpy as np
from scipy.interpolate import interp1d

# 此处为模拟你已经生成的interp1d插值函数，实际使用时替换为你自己的f即可
x_sample = np.linspace(-50, 50, 1000)
y_sample = x_sample ** 3
f = interp1d(x_sample, y_sample)

# 你的原始y、z数组
y = np.arange(0, 100)
z = np.arange(0, 100)

# 广播计算所有i,j组合对应的x=(y[j]-z[i])/2
# y扩维为(1, N)，z扩维为(M, 1)，相除后自动广播为(M, N)的全量x数组
x_all = (y[np.newaxis, :] - z[:, np.newaxis]) / 2

# interp1d原生支持数组输入，直接传入全量x数组即可得到对应形状的结果
fnew = f(x_all)

• 若8000*8000的全量数组占用内存过高（单精度浮点数约占256MB，双精度约占512MB），可采用分批处理的方式降低内存占用，性能损失极小：

内存友好的分批实现 #

batch_size = 1000  # 可根据可用内存调整批次大小
fnew = np.zeros((len(z), len(y)))
for start in range(0, len(z), batch_size):
    end = min(start + batch_size, len(z))
    z_batch = z[start:end]
    x_batch = (y[np.newaxis, :] - z_batch[:, np.newaxis]) / 2
    fnew[start:end, :] = f(x_batch)

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