BUG: RBFInterpolator is much slower than Rbf for vector data

Describe your issue.

In the past ve been using the old RBF interpolator (Rbf) to interpolate between vectors and since it didn't support the interpolation of vectors, but only scalars, I've had to manually extract the matrix of distances from the Rbf object (see code below). Since the new RBFInterpolator now supports interpolation of vectors, I gave it a try and noticed that it is almost a factor of 10 slower in my application and testcase below ( 29.7 s vs 3.9s). While this is technically not a bug, but this is very surprising.

I see that the time is being spent mostly here

scipy/scipy/interpolate/_rbfinterp_pythran.py Line 175 15df208 inside hte pythran evaluation function. I can't immediately see anything obvious there that could be responsible for the regression though.

I also checked that if I set OMP_NUM_THREADS to 1, the performance difference is reduced but stays factor of 3.5 (10s vs 35s) (so the issue is not due to lack of parallelization of the pythran code)

Thanks

Reproducing Code Example

import numpy as np
import scipy.interpolate
import time
rstate = np.random.default_rng(44)
# initial point distribution
npt = 6000
xs, ys = rstate.normal(size=(2, npt))
# my data vectors are 4000dimensional
nx = 4000
xgrid = np.arange(nx)
data = (xgrid[None, :] * xs[:, None]**2 + ys[:, None]**2)
# number of points 
neval = 1000
xnew, ynew = rstate.normal(size=(2, neval))
t1 = time.time()
# use the old interface on 1d data to get coefficients to apply for the
# to the 4000d data vectors
RR = scipy.interpolate.Rbf(xs, ys, data[:, 0])
coeffs = scipy.linalg.solve(RR.A, data)
r = RR._call_norm(np.array([xnew, ynew]), RR.xi)
pred0 = np.dot(RR._function




    
(r), coeffs)
# use the new interface
t2 = time.time()
RR = scipy.interpolate.RBFInterpolator(np.array([xs, ys]).T, data)
pred1 = RR(np.array([xnew, ynew]).T)
t3 = time.time()
print(t2 - t1, t3 - t2)

Error message

SciPy/NumPy/Python version information

1.7.3 1.21.4 sys.version_info(major=3, minor=8, micro=10, releaselevel='final', serial=0)

@treverhines

@treverhines can you have a look at this issue?

@segasai

@segasai Thanks for pointing this out. I am getting similar times on my machine: 13 seconds for Rbf -based evaluation, and 35 seconds for RBFInterpolator with a single thread.

The issue is that, in this case, the matrix multiplication done to evaluate the interpolant is significantly faster with np.dot than the implementation in _evaluate . The difference in times seems to get larger with higher dimensional data. If you lower the data vector dimensions down to something like nx=100 , you should see that RBFInterpolator is either comparable or faster than the Rbf -based evaluation. So this issue is just for interpolating high-dimensional data.

If I replace this block:
scipy/scipy/interpolate/_rbfinterp_pythran.py Lines 220 to 222 15df208

with a call to np.dot , then BLAS will be used to perform a matrix-vector multiplication for each evaluation point, and that brings the time down from 35 seconds to about 20 seconds on my machine. So not as fast as your Rbf -based evaluation; however, it is more memory efficient since the Rbf -based evaluation requires holding a potentially large (neval, npts) matrix in memory.

If I were to make that change, _rbfinterp_pythran.py would need to be compiled with BLAS, and I am not sure if that is a deal-breaker for scipy because of added build complications. FWIW, I have been pip installing pythran-openblas ( https://github.com/serge-sans-paille/pythran-openblas ) and compiling the pythran code with --config compiler.blas=pythran-openblas . Perhaps we can add pythran-openblas as a build dependency? @rgommers

Can you try reversing loop order from j,k to k,j if it makes a difference ?

swapping the two for loops increased the time from 35 seconds to 154 seconds.

This change master...segasai:rbf_fix
allows me to reduce the time for the RBF evaluation from 27 seconds to 0.6 second. (basically I've moved np.dot out of the pythran file and only kept kernel/monomial calculations there)
Obviously there are potential memory implications here (we need to keep the matrix of coefficients in memory), but I'd say it's better to do chunking if we can get 50times better performance.

I was testing the performance of just this line in my code

pred1 = RR(np.array([xnew, ynew]).T)
and I was using OMP_NUM_THREADS=1

After playing around with this a bit, I think I prefer your suggested approach of building the kernel and polynomial matrix in chunks and dotting each chunk with coeffs. Although, I would suggest keeping the chunking and np.dot within _evaluate for slightly better performance (build issues aside). Are you interested in making a PR? Otherwise, I can probably work on it later this month.

I was planning to submit the PR based on what I have (i.e np.dot outside pythran). Personally I'd prefer not to delve into blas/lapack linking/building stuff/dependencies on pythran version, so if you think that np.dot has to be inside, I'll let you do it.

Navigation Menu

BUG: RBFInterpolator is much slower than Rbf for vector data #15150

Describe your issue.

Reproducing Code Example

Error message

SciPy/NumPy/Python version information

`Metadata`

`Assignees`

`Labels`

`Type`

`Projects`

`Milestone`

`Relationships`

`Development`

BUG: RBFInterpolator is much slower than Rbf for vector data #15150

Describe your issue.

Reproducing Code Example

Error message

SciPy/NumPy/Python version information

tupui commented on Dec 3, 2021

treverhines commented on Dec 4, 2021

segasai commented on Dec 5, 2021

treverhines commented on Dec 5, 2021

segasai commented on Dec 7, 2021

treverhines commented on Dec 8, 2021

segasai commented on Dec 8, 2021

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

`treverhines commented on Dec 8, 2021`

`segasai commented on Dec 8, 2021`

`Assignees`

`Labels`

`Type`

`Projects`

`Milestone`

`Relationships`

`Development`

`Issue actions`