How to accelerate the speed of frequent interpolation？

For example,

I have a point set( A), where each point corresponds to a value of (Value_A). We can use "scatteredInterpolant" function to obtain the value (Value_B) on a point set B. The number of points in the point set B may be more or less than that in the point set A. However, this process is very time-consuming for a large number of points.

In my code, this interpolation process is required frequent operations, but the point sets A and B do not change, only the value (value_A) on the points will change.

In this case, can we extract the weight matrix A2B from A to B and directly obtain (Value_B) through Value_A * A2B. This way, the overall efficiency will be very high. I tried to check the underlying code of the scatteredInterpolant function, but it was not available.

Answer 1

[Balavignesh]

Hi XL,

Yes, you are correct. You can improve the efficiency of your interpolation process by precomputing a weight matrix that maps the values from point set A to point set B. This approach is better because it eliminates the need to perform the interpolation from scratch each time.

You can use Delaunay triangulation to triangulate the points in 'set A'. For each point in 'set B', find the corresponding simplex in the triangulation of 'set A'. Then, compute the barycentric coordinates of each point in 'set B' with respect to the vertices of the simplex it lies in. You can use these coordinates to form the weight matrix 'A2B'. After forming the weight matrix 'A2B', you can quickly compute the interpolated values for any set of values at the points in 'A' by multiplying the value vector with the weight matrix.

The following example code may help you understand this concept better:

% Dummy point sets A and B
A = [0, 0; 1, 0; 0, 1; 1, 1]; % Coordinates of points in set A
B = [0.5, 0.5; 0.75, 0.25; 0.25, 0.75]; % Coordinates of points in set B
Value_A = [10; 20; 30; 40]; % Values at points in set A

% Step 1: Triangulate point set A
tri = delaunayTriangulation(A);

% Step 2: Find the simplex containing each point in B
tess = pointLocation(tri, B);

% Step 3: Compute barycentric coordinates for points in B
numPointsB = size(B, 1);
numPointsA = size(A, 1);
A2B = zeros(numPointsB, numPointsA);

for i = 1:numPointsB
    if ~isnan(tess(i))
        % Vertices of the simplex containing point B(i)
        vertices = tri.ConnectivityList(tess(i), :);
        
        % Barycentric coordinates of B(i) with respect to these vertices
        lambda = cartesianToBarycentric(tri, tess(i), B(i, :));
        
        % Assign weights to the corresponding vertices
        A2B(i, vertices) = lambda;
    end
end

% Now, you can quickly compute Value_B for any Value_A
Value_B = A2B * Value_A;

% Display results
disp('Weight Matrix A2B:');

Weight Matrix A2B:

disp(A2B);

         0    0.5000    0.5000         0
         0    0.7500    0.2500         0
         0    0.2500    0.7500         0

disp('Interpolated Values at Points in B:');

Interpolated Values at Points in B:

disp(Value_B);

   25.0000
   22.5000
   27.5000

% Example usage for different Value_A
Value_A_new = [15; 25; 35; 45]; % New values at points in set A
Value_B_new = A2B * Value_A_new;

disp('Interpolated Values at Points in B for new Value_A:');

Interpolated Values at Points in B for new Value_A:

disp(Value_B_new);

   30.0000
   27.5000
   32.5000

This approach should significantly improve the efficiency of your interpolation process, especially when the point sets A and B remain constant but the values at the points in A change frequently.

Kindly have a look at the following documentation links to have more information on:

• Delaunay triangulation: https://www.mathworks.com/help/matlab/ref/delaunaytriangulation.html
• Point Loacation: https://www.mathworks.com/help/matlab/ref/delaunaytri.pointlocation.html

Hope that helps!

Balavignesh

[John D'Errico]

+1 of course, as this is correct. I might even go slightly further.

The weight matrix idea is a great one. It allows you to compute the interpolated values for new sets of outputs, requiring nothing more than a matrix multiply.

A problem is if there are many points to interpolate. Then the "weight matrix" will become a large one, and mostly composed of zeros. This will certainly be the case for a triangulated set of data. But a matrix multiply will then be making many multiplies by zeros.

You fix this problem by making that weight matrix a sparse one. When you multiply by a sparse matrix, MATLAB is smart, in that it knows to not bother performing the zero multiplies. And they do cost time, as it takes the same amount of time to perform x*0 as it does to compute x*y.

[xl z]

Thank you very much for your reply.

Answer 2

[Steven Lord]

How are you repeating your interpolation? Are you recreating the scatteredInterpolant object each time? If so, instead follow the "Replacement of Sample Values" example on that documentation page to use the same triangulation each time.

rng('default')
x = -2.5 + 5*rand([50 1]);
y = -2.5 + 5*rand([50 1]);
v = x.*exp(-x.^2-y.^2);
tic
F = scatteredInterpolant(x,y,v);
toc

Elapsed time is 0.013013 seconds.

xsample = -2.5:0.1:2.5;
ysample = -2.5:0.1:2.5;
tic
v1 = F(xsample, ysample);
toc

Elapsed time is 0.001524 seconds.

vnew = x.^2 + y.^2;
tic
F.Values = vnew;
v2 = F(xsample, ysample);
toc

Elapsed time is 0.003059 seconds.

Changing the values and interpolating does take longer than simply interpolating the original values, but it takes less time than creating the scatteredInterpolant object initially then interpolating. Granted in this case it's close because there are only 50 points, but try it with your own (I'm guessing larger) data set.

tic
F = scatteredInterpolant(x, y, vnew);
v3 = F(xsample, ysample);
toc

Elapsed time is 0.003083 seconds.

[xl z]

Thank you very much for your reply. There are many points to be interpolated, i.e., the size of the xsample and ysample matrices are large. This leads to a slower execution speed of ‘v1 = F(xsample, ysample)’ in each iteration. The xsample and ysample are variable. Namely, it requires frequent execution of 'v1 = F(xsample, ysample)' as the iterations.