1d linear interpolation vectorization

Hi all, I have a vectorization problem I wonder if i could get some help with. I just can't seem to get it done.

I have a matrix A with size N*M, where N=200 and M=200000. The numbers in A are random, and cannot be sorted (for technical reasons elsewhere in my code). I have two vectors B and C with size 1*N. I want to use linear interpolation for all M, which can be done using something like

D=zeros(size(M))
Parfor m=1:200000
      D(:,)=interp1(M(:,m),A,B);
End

But needless to say this is super slow (it's all part of a large code run in several outer loops). The problem is that I have a new grid in each iteration, Griddedinterpolant would be faster if I had the same grid but new evaluation points since then I could define the interpolant outside of the loop. So I'd like to try to vectorize this code. I wrote a simple linear 1D interpolation routine (see eg loreen shures blog) but I just can't get it to work. Or even if I could loop over N instead of M that would be a huge improvement!

Any help would be amazing!!

Answer 1

[Matt J]

If B falls within the boundaries of all domain vectors M(:,m), you can vectorize it by pre-incrementing all the M(:,m) and stitching them together,

          Mmax=max(M);
          Mmin=min(M);
          z=cumsum([0,Mmax(1:end-1)-Mmin(2:end)+1]);

          Anew=repmat(A,1,size(M,2));
          Bnew=bsxfun(@plus,B,z);
          Mnew=bsxfun(@plus,M,z);

       D=interp1(Mnew(:),Anew(:),Bnew)

[Matt J]

Anders Österling commented:

Hi Matt, thanks a lot for your answer!! I tried it briefly and it is super fast! Unfortunately B doesn't always fall within the boundaries of M(:,m) so it doesn't work everywhere.

[Matt J]

You can easily detect those B that are out of bounds and post-process to your liking. By default, interp1 sets these to NaN and that's what I do below, but you can obviously post-process the out-of-bounds B in other ways if you wish.

outofbounds=bsxfun(@lt,Bnew, Mnew(1,:)) | bsxfun(@gt,Bnew, Mnew(end,:));

D(outofbounds)=NaN;

[Anders Österling]

Hi,

Excellent tip - it worked like a charm and is extremely fast!

Thanks a lot!! /A