Faster finding of 1D linear interpolation nodes and weights for each element in ND matrix

In a problem I'm working on now, I compute some values in a matrix x and I then for each element in x need to find the index of the closest element below in a monotonically increasing vector X as well as the relative proximity of the x elements to the first elements on their either side. (This is essentially linear interpolation without doing the actual interpolation.) I'm doing this maaaany times so I really super extra interested in it being as fast as possible.

I have written a function locate that I can call with some example data:

X = linspace(5, 300, 40);
x = randi(310, 5, 6, 7);

[ii, weights] = locate(x, X);

I have written two versions of locate. The first is for exposition and the second is my best attempt at speeding up the computations. Do you have any suggestions or alternative approaches for how I could accelerate performance further?

1. Exposition

function [ii, weights] = locate(x, X)
    % LOCATE Locate first node on grid below a given value.
    %
    %   [ii, weights] = locate(x, X) returns the first node in X that is below
    %   each element in x and the relative proximities to the two closest nodes.
    %
    %   X must be a monotonically increasing vector. x is a matrix (of any
    %   order).
    
    % Preallocate
    ii = ones(size(x));  % Indices of first node below (or 1 if no nodes below)
    weights = zeros([2, size(x)]);  % Relative proximity of the two closest nodes
    
    % Find indices and compute weights
    for ix = 1:numel(x)
        if x(ix) <= X(1)
            ii(ix) = 1;
            weights(:, ix) = [1; 0];
        elseif x(ix) >= X(end)
            ii(ix) = length(X) - 1;
            weights(:, ix) = [0; 1];
        else
            ii(ix) = find(X <= x(ix), 1, 'last');
            weights(:, ix) = ...
                [X(ii(ix) + 1) - x(ix); x(ix) - X(ii(ix))] / (X(ii(ix) + 1) - X(ii(ix)));
        end
    end
end

2. Best attempt

function [ii, weights] = locate(x, X)
    % LOCATE Locate first node on grid below a given value.
    %
    %   [ii, weights] = locate(x, X) returns the first node in X that is below
    %   each element in x and the relative proximities to the two closest nodes.
    %
    %   X must be a monotonically increasing vector. x is a matrix (of any
    %   order).
    
    % Preallocate
    ii = ones(size(x));  % Indices of first node below (or 1 if no nodes below)
    weights = zeros([2, size(x)]);  % Relative proximity of the two closest nodes
    
    % Find indices
    for iX = 1:length(X) - 1
        ii(X(iX) <= x) = iX;
    end
    
    % Find weights
    below = x <= X(1);
    weights(1, below) = 1;  % All mass on the first node
    weights(2, below) = 0;
    
    above = x >= X(end);
    weights(1, above) = 0;
    weights(2, above) = 1;  % All mass on the last node
    
    interior = ~below & ~above;
    xInterior = x(interior)';
    iiInterior = ii(interior);
    XBelow = X(iiInterior)';
    XAbove = X(iiInterior + 1)';
    weights(:, interior) = ...
        [XAbove - xInterior; xInterior - XBelow] ./ (XAbove - XBelow);
end

[Rik]

A few corrections to your array function:

%you need to add this to your array function:

    %force X to be a column vector
    X=reshape(X,1,[]);

%and replace the last block of code with this:
    interior = ~below & ~above;
    xInterior = x(interior)';
    iiInterior = ii(interior);
    XBelow = X(iiInterior);
    XAbove = X(iiInterior + 1);
    A=XAbove - xInterior;
    B=xInterior - XBelow;
    C=XAbove - XBelow;
    weights(1, interior) = A ./ C;
    weights(2, interior) = B ./ C;

Are these array sizes representative of your actual data? Using array operations can be slower for tiny loops, but much faster for longer loops. Unfortunately in this case the array version is as fast or slower than the looped version. In both cases the profiler indicates that the finding of the indices is the slow part. I tried to find a way around it, but I only found a slower version (which also requires lots of memory for larger arrays):

% Preallocate
ii = ones(size(x));

tic
% Find indices
for iX = 1:length(X) - 1
    ii(X(iX) <= x) = iX;
end
toc
tic
dims=[1 3:(ndims(x)+1) 2];
L=permute(X(1:(end-1)),dims)<=x;%implicit expansion
subs=cell(1,ndims(L));
[subs{:}]=ind2sub(size(L),find(L));
val=subs{end};subs(end)=[];subs=[subs{:}];
ii_=accumarray(subs,val,size(x),@max,1);
toc
max(abs(ii(:)-ii_(:)))%should return 0

[Ameer Hamza]

When I run you code with these inputs

X = 1:4;
x = 1:10;
[ii1, weights] = locate(x, X);

I get

ii1 =
     1     2     3     3     3     3     3     3     3     3

Shouldn't the output be

ii1 =
     1     2     3     4     4     4     4     4     4     4

[Fredrik P]

Thanks for giving it a go! And thanks for correcting my initial code! The array sizes are representative of my actual data. Or rather, in most function calls will x will be a 2D array (in the vicinity of 9x5) but for a sizeable share of the function calls x will be 3D (~9x5x5).

[Fredrik P]

@Ameer Hamza A very, very late reply to your comment. No, the results that you get are correct. As I need two nodes in X to compute the weights, the greatest index that I return is length(X) - 1.