Asked by Gustaf Lindberg
on 21 Feb 2013

I have a loop where I do many extrapolations using interp1 with the flags 'linear' and 'extrap' so I guess it's technically an extrapolation. Now I want to use the same code but do many pairs of extrapolations at the same time. I have two matrices, X and Y where each row is a matching set to be extrapolated for a value xi = 0 but I can't wrap my head around how to make it work. I've been told that it is possible.

Basically, what I want to do is a faster version of this:

for i = 1:n yi = interp1(X(i,:),Y(i,:),0,'linear','extrap') end

But without the loop using interp2. each row n in the matrices X and Y are a matching set to be evaluated independently.

*No products are associated with this question.*

Answer by Matt J
on 21 Feb 2013

Edited by Matt J
on 21 Feb 2013

[M,N]=size(X); yi=zeros(M,1); map=X>=0; [maxval,idx]=max(map,[],2);

a = idx==1 & maxval==1; %All X(i,:) are non-negative, extrapolate left b = idx==N | maxval==0; %All X(i,:) are nonpositive, extrapolate right c= ~(a|b); %Interpolation needed

%Extrapolate left e1=X(a,1); c1=Y(a,1); e2=X(a,2); c2=Y(a,2); slopes=(c2-c1)./(e2-e1); yi(a) = c1-slopes.*e1;

%Extrapolate right e1=X(b,end-1); c1=Y(b,end-1); e2=X(b,end); c2=Y(b,end); slopes=(c2-c1)./(e2-e1); yi(b) = c1-slopes.*e1;

%Interpolate (linearly) j = sub2ind([M,N],find(c),idx(c)); e1=X(j-M); c1=Y(j-M); e2=X(j); c2=Y(j); slopes=(c2-c1)./(e2-e1); yi(c) = c1+slopes.*abs(e1);

Matt J
on 21 Feb 2013

Had to make a few fixes, but it should work now. Here's some sample input/output data from my tests

X =

1.0000 2.0000 3.0000 4.0000 5.0000 -10.0000 -9.0000 -8.0000 -7.0000 -6.0000 -2.2500 -1.2500 -0.2500 0.7500 1.7500

Y =

3 4 5 6 7 -5 -4 -3 -2 -1 -2 -1 0 1 2

yi =

2.0000 5.0000 0.2500

Matt J
on 21 Feb 2013

Some more tests, this time a speed comparison with for-looped interp1. My version shows a factor of 300 speed-up, but of course I don't know what data sizes you're actually working with.

%Fake data n=3e4; Y=rand(n,300); X=bsxfun(@plus, sort(Y,2), 2*rand(n,1)+1 );

tic; % My code toc

Elapsed time is 0.017832 seconds.

yi2=yi; tic; for i = 1:n yi2(i) = interp1(X(i,:),Y(i,:),0,'linear','extrap'); end toc;

Elapsed time is 5.419498 seconds.

Answer by Teja Muppirala
on 21 Feb 2013

As others have also pointed out, this is **REALLY** not what INTERP2 is meant to do. In fact if anything, since there are 3 variables: X, Y, and i, you'd have to use INTERP3. But I think that would be rather inefficient.

Why do you think the FOR loop is a problem? How do you think doing something else will make it any more "stable"?

I can think of ways that might make it faster, but not really any more stable:

[n,m] = size(X) [XS,Xi] = sort(X,2);

YS = Y(bsxfun(@plus,(Xi(:,1:2)-1)*n,(1:n)')); y = YS(:,1) + diff(YS,1,2)./diff(XS(:,1:2),1,2).*(-XS(:,1));

Gustaf Lindberg
on 21 Feb 2013

José-Luis
on 21 Feb 2013

With *interp2()* you are interpolating to a plane, not a line.

I guess you could assign increasing y values for each row, call them *y =[ 1:number of rows]* and then interpolate for x and y, where *x = zeros(nRows,1)*, using *interp2()*. That could work.

Having said that, it is a very roundabout way of going about it and probably very inefficient. It takes more calculations to interpolate to a plane than to a line.

Opportunities for recent engineering grads.

## 4 Comments

## José-Luis (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/64379#comment_131259

Do you know in advance if it is always gonna be an extrapolation? Are the values in a row sorted? Is 0 always going to be inferior (or superior) to the values in X? A way to speed up a function is to get rid of the extra checks it performs. That can be dangerous, but you have to trim all the fat if you want speed.

## Gustaf Lindberg (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/64379#comment_131278

I know nothing in advance. There are many, many iterations done with this check and the values can vary widely and there are no real sorting. Since this is for a CFD solver, stability is more important than speed. I've been told that interp2 could do what I want but I can't figure out how. All examples just uses it as an look up table for a 2D surface but I want multiple, independent 1D vectors.

## José-Luis (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/64379#comment_131280

I don't think

interp2()is suited for what you want. If stability is more important than speed then I would stick tointerp1()and a loop. If I really want to speed up things and I know nothing about the data, then I would write my own routine, preferably as a mex file.You could also use

parfor()when evaluating in a loop to speed things along.## Gustaf Lindberg (view profile)

Direct link to this comment:http://www.mathworks.com/matlabcentral/answers/64379#comment_131282

I was more or less told to use interp2 by my supervisor so I don't have that much of a choice right now but your suggestions are appreciated.