如何在Scipy中优化大量点的重复RBF插值以提升最小二乘效率？

我正在求解一个最小二乘问题：自变量为数百个控制点上的若干数值，通过RBFInterpolator（邻居数≤10）插值至约15万个其他点；插值结果用于计算残差——每个点存在目标值，插值变量生成预测值后得到残差。这些点为空间坐标，非规则网格，控制点是点集的子集。

残差函数由工厂函数生成，简化代码如下：

def create_global_residual_function(
    model: Callable[[np.ndarray, float, ...], np.ndarray],
    control_points: np.ndarray,
    points: np.ndarray,
    y: np.ndarray,
    x: np.ndarray,
    thin_plate_spline_smoothing: float,
    thin_plate_spline_neighbors: int,
    thin_plate_spline_degree: int,
    residual_boost_factor: float = 2.0,
) -> Callable[[Tuple[float, ...]], np.ndarray]:
    
    num_control_points = control_points.shape[0]
    num_points = points.shape[0]

    x = x.reshape(1, len(x))

    tps = RBFInterpolator(
        control_points,
        np.zeros((num_control_points, 5), dtype=float),
        smoothing = thin_plate_spline_smoothing,
        neighbors = thin_plate_spline_neighbors,
        kernel = "thin_plate_spline",
        degree = thin_plate_spline_degree
    )

    def global_residual_function(params: Tuple[float, ...]) -> np.ndarray:
        tps.d = np.array(params).reshape(num_control_points, <number of parameters>)
        model_parameters = tps(points)
        x0 = model_parameters[:, 0].reshape(num_points, 1)
        x1 = model_parameters[:, 1].reshape(num_points, 1)
        # ... 其他参数提取逻辑
        y_hat = model(x, x0, x1, <and so on...>)
        return (y_hat - y).flatten()

    return global_residual_function

每次迭代仅修改RBFInterpolator的d属性，插值目标点固定，但需执行数百至数千次插值才会收敛。重复调用RBFInterpolator.__call__会重复计算最近邻、距离、核函数等冗余操作，理论上可预计算稀疏矩阵复用，实现从控制点值到目标点值的高效计算。

我的方案在小测试用例（约数千点+20个控制点）中验证有效，但效率极低，无法扩展至大规模问题。

请问在Scipy中是否存在更高效的解决方案？例如提取上述稀疏矩阵，或使用其他适配的插值类？

针对大规模迭代插值的效率问题，Scipy生态下有以下几种直接的优化方式：

1. 预计算插值权重稀疏矩阵 #

RBFInterpolator的核心是控制点值到目标点值的线性变换，当控制点、目标点、邻居数、核函数及平滑参数固定时，这个变换的权重矩阵是恒定的，可提前计算并复用：

• 用scipy.spatial.cKDTree预计算每个目标点的thin_plate_spline_neighbors个最近控制点，得到邻居索引数组；
• 对每个目标点的邻居，预计算薄板样条核函数值（r² * log(r)）和多项式基函数值；
• 结合平滑参数，为每个目标点构建局部线性方程组，求解得到该点对邻居控制点的权重；
• 将所有目标点的权重整理成scipy.sparse.csr_matrix（形状为(num_points, num_control_points)），后续迭代仅需执行model_parameters = W @ params_reshaped即可完成插值，彻底避免重复计算最近邻、距离和核函数。

2. 复用RBFInterpolator的内部计算结果 #

查看scipy.interpolate.RBFInterpolator的源码可知，它初始化时会构建KDTree，但每次调用__call__仍会重复查询最近邻和计算核矩阵。你可以手动提取并缓存这些中间结果：

• 初始化tps后，调用tps._query_neighbors(points)获取邻居索引和距离；
• 预计算核函数矩阵和多项式基矩阵，结合平滑参数求解局部权重；
• 将权重存储为稀疏矩阵，后续直接用矩阵乘法替代tps(points)调用。

3. 切换到线性插值类（若场景允许） #

如果你的问题可以接受线性插值的平滑特性，scipy.interpolate.LinearNDInterpolator在固定点集时的效率更高：它会预计算Delaunay三角剖分，后续迭代仅需查询三角单元并计算线性组合，避免了RBF的核函数计算开销。但需注意这会改变插值结果的平滑度，需根据需求权衡。

4. 分块并行计算 #

若预计算全量权重矩阵的内存占用过高（15万点×数百控制点的稀疏矩阵通常内存友好，但极端场景可考虑），可将目标点分块处理：每次迭代仅对一块点执行插值，用scipy.sparse的并行矩阵乘法加速，平衡内存占用和计算速度。

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