为何RegularGridInterpolator无法返回R^d维度的函数输出值？

最小可复现示例（含可行与不可行输出） #

import numpy as np
from scipy.interpolate import RegularGridInterpolator, griddata

def f(x1, x2, x3):
    return x1 + 2*x2 + 3*x3, x1**2, x2

# 定义输入点序列
xi = [np.linspace(0, 1, 5), np.linspace(0, 1, 5), np.linspace(0, 1, 5)]

# 生成规则网格
x1, x2, x3 = np.meshgrid(*xi, indexing='ij')

# 计算多维度输出
y = f(x1, x2, x3)
assert (y[0][1][1][3] == (0.25 + 2*0.25 + 3*0.75))
assert (y[1][1][1][3] == (0.25**2))
assert (y[2][1][1][3] == 0.25)

#### 可行：仅插值单个输出维度
# 在(0.3, 0.3, 0.4)处插值第一个输出维度
interp = RegularGridInterpolator(xi, y[0])
print(interp([0.3, 0.3, 0.4]))  # 预期输出2.1，实际符合预期


#### 不可行：尝试同时插值所有输出维度
# 在(0.3, 0.3, 0.4)处插值全部输出维度
interp = RegularGridInterpolator(xi, y)
print(interp([0.3, 0.3, 0.4]))  # 未输出预期的array([2.1, 0.1, 0.3])

griddata的对比实现（支持多输出插值） #

# griddata的多输出插值实现
grid_for_griddata = np.array([x1.flatten(), x2.flatten(), x3.flatten()]).T
assert (grid_for_griddata.shape == (125, 3))
y_for_griddata = np.array([y[0].flatten(), y[1].flatten(), y[2].flatten()]).T
assert (y_for_griddata.shape == (125, 3))
# 在(0.3, 0.3, 0.4)处插值
result = griddata(grid_for_griddata, y_for_griddata, [0.3, 0.3, 0.4], method='linear')[0]
print(result)  # 如预期输出array([2.1, 0.1, 0.3])

疑问 #

• 是否错误使用了RegularGridInterpolator？
• 已知griddata可实现需求，但由于数据是规则网格，理论上RegularGridInterpolator速度更快，应该优先使用它，对吗？

性能对比 #

内容的提问来源于stack exchange，提问作者FluidMechanics Potential Flows