A way to speed up my custom interpolation algorithm?

I'm writing a simulation program which often has to read out data from small data tables.

The x and y values are provided as 1xn arrays:

x = [0, 2, 5, 9, 14, 67];
y = [0, 3, 8, 2, 45, 324];

The interpolation scheme i need is equivalent to

interp1(x,y,4,'previous','extrap');

But as interp1 is very slow i wrote my own function. According to Matlab Profiler it is 5x faster but still takes up 80% of my code runtime.

function y = table_lookup(x,table_y,table_x)    
    % Check if lookup tables are of the same length
    assert(length(table_y) == length(table_x))
    
    % initialize y
    y = NaN;
    
    % Iterate over lookup table
    i = 1;
    while i <= length(table_x)
        % If x value is exactly in the x_table assign corresponding y value
        if table_x(i) == x
            y = table_y(i);
            break;
    
        % If x value is grater than x value assign previous y value
        elseif table_x(i) > x && i ~= 1
            y = table_y(i-1);
            break;
        end
        i = i + 1;
    end
    
    % If the iteration completed over the whole table assign the last y value
    if i > length(table_x)
       y = table_y(end);
    end

end

Maybe someone has a suggestion for possible optimizations?

[Jonas]

can you try

yq=y(find(x<xq,1,'last'))

[valnim]

Thank you.

yq = y(find(x <= xq, 1, 'last')) 

is about 10% faster than my function.

Answer 1

[valnim]

I found

yq = y(find(x <= xq, 1, 'last'))

proposed by @Jonas works best for my usecase.

[Jon]

If @Jonas method worked the best for you, then you should accept his answer, rather than putting what should be a comment on your thread as an answer, and then accepting that as an answer. Note you can unaccept an answer and then pick a new one if you want to correct this.

Also, how did my approach using just logical indices compare to the find used by @Jonas. i had the impression that find could be a little slower if not needed but would be interested to if that were actually the case.

[Jon]

Oh I see that @Jonas never posted an answer only a comment, so you couldn't accept his answer. Sorry

Answer 2

[Jon]

You could try benchmarking this way too for comparison

xe = [x inf]

yq = y(xq >= xe(1:end-1) & xq <xe(2:end))

This returns empty if xq is less than the first element in x, not sure if you need to handle this edge case.

[Jon]

Did you benchmark this one. I'm curious how it performed relative to using the find proposed by @Jonas

[valnim]

The solution by @Jonas was on average three times faster than your solution.

For 1e5 iterations on 1x5 array's @Jonas's solution took 0.15 s and @Jon's solution took 0.45 s.

[Jon]

That's interesting. Thanks for checking it out. I had thought that perhaps the "find" function would be relatively expensive compared to just using logical indexing. Good to know that "find" is so efficient!