How do I efficiently replace a function with a lookup?

Question

I am trying to increase the speed of code that operates on large datasets. I need to perform the function out = sinc(x), where x is a 2048-by-37499 matrix of doubles. This is very expensive and is the bottleneck of my program (even when computed on the GPU).

I am looking for any solution which improves the speed of this operation. I expect that this might be achieved by pre-computing a vector LookUp = sinc(y) where y is the vector y = min(min(x)):dy:max(max(x)), i.e. a vector spanning the whole range of expected x elements.

How can I efficiently generate an approximation of sinc(x) from this LookUp vector?

I need to avoid generating a three dimensional array, since this would consume more memory than I have available.

Here is a test for the interp1 solution:

a = -15;
b = 15;
rands = (b-a).*rand(1024,37499) + a;

sincx = -15:0.000005:15;
sincy = sinc(sincx);

tic
res1 = interp1(sincx,sincy,rands);
toc

tic
res2 = sinc(rands);
toc'

sincx = gpuArray(sincx);
sincy = gpuArray(sincy);
r = gpuArray(rands);

tic
r = interp1(sincx,sincy,r);
toc

r = gpuArray(rands);

tic
r = sinc(r);
toc

Elapsed time is 0.426091 seconds.
Elapsed time is 0.472551 seconds.
Elapsed time is 0.004311 seconds.
Elapsed time is 0.130904 seconds.

Corresponding to CPU interp1, CPU sinc, GPU interp1, GPU sinc respectively

Answer 1

Not sure I understood completely your problem. But once you have LookUp = sinc(y) you can use the Matlab function interp1

out = interp1(y,LookUp,x)

where x can be a matrix of any size

Answer 2

I came to the conclusion, that your code can not be improved significantly. The fastest possible lookup table is based on simple indexing. For a performance test, lets just perform the test based on random data:

%test data:
x=rand(2048,37499);
%relevant code:
out = sinc(x);

Now the lookup based on integer indices:

a=min(x(:));
b=max(x(:));
n=1000;
x2=round((x-a)/(b-a)*(n-1)+1);
lookup=sinc(1:n);
out2=lookup(x2);

Regardless of the size of the lookup table or the input data, the last lines in both code blocks take roughly the same time. Having sinc evaluate roughly as fast as a indexing operation, I can only assume that it is already implemented using a lookup table.

Answer 3

I found a faster way (if you have a NVIDIA GPU on your PC) , however this will return NaN for x=0, but if, for any reason, you can deal with having NaN or you know it will never be zero then:

if you define r = gpuArray(rands); and actually evaluate the sinc function by yourself in the GPU as:

tic
r=rdivide(sin(pi*r),pi*r);
toc

This generally is giving me about 3.2x the speed than the interp1 version in the GPU, and its more accurate (tested using your code above, iterating 100 times with different random data, having both methods similar std).

This works because sin and elementwise division rdivide are also GPU implemented (while for some reason sinc isn't) . See: http://uk.mathworks.com/help/distcomp/run-built-in-functions-on-a-gpu.html

Answer 4

m = min(x(:));
y = m:dy:max(x(:));
LookUp = sinc(y);

now sinc(n) should equal

LookUp((n-m)/dy + 1)

assuming n is an integer multiple of dy and lies within the range m and max(x(:)). To get to the LookUp index (i.e. an integer between 1 and numel(y), we first shift n but the minimum m, then scale it by dy and finally add 1 because MATLAB indexes from 1 instead of 0.

I don't know what that wll do for you efficiency though but give it a try.

Also you can put this into an anonymous function to help readability:

sinc_lookup = @(n)(LookUp((n-m)/dy + 1))

and now you can just call

sinc_lookup(n)