I am trying to solve a least squares problem where my independent variables are a handful of values at each of a few hundred control points, which are then interpolated to a huge number of other points via a RBFInterpolator (with neighbours <= 10), and then those interpolated values are used to compute the residuals – each point has some target values and the interpolated variables are used to generate predictions, from which I get residuals. These points are just coordinates in space, they are not on a regular grid or anything like that. The control points are a subset of the points.

The residual function used for the least squares optimizer is generated by a factory function and looks something like this (simplified):

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)
        <and so on ...>
        y_hat = model(x, x0, x1, <and so on...>)
        return (y_hat - y).flatten()

    return global_residual_function

In each iteration, what I am changing is only the d attribute of the RBFInterpolator. The points I interpolate onto are the same every time. The interpolation has to happen on the order of 100s to 1000s of times before the model will converge. It seems silly to me to invoke RBFInterpolator.call repeatedly when the the input points are unchanging. I’m sure there is a lot of computation being unnecessarily repeated to find the nearest X control points to each points, calculate the distances and feed into kernels, etc. In theory there should be a single (sparse) matrix that could be computed and then repeatedly used to compute the values of a scalar array on the points from the value of the scalar array on the control points.

My solution works – I’ve tested it on small test cases (~1000s of points with ~20 control points) and gotten correct results from the optimization. The problem is it is very inefficient and doesn’t scale well to larger problems.

Is there a more efficient solution for this in scipy? Perhaps some way to extract the (sparse) matrix I referred to above? Or a different interpolator class I should be using for this application?