Use interp2 for many simultanious 1D operations?
显示 更早的评论
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.
4 个评论
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
2013-2-21
I don't think interp2() is suited for what you want. If stability is more important than speed then I would stick to interp1() 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
2013-2-21
回答（2 个）
[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);
2 个评论
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
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.
Teja Muppirala
2013-2-21
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));
2 个评论
Gustaf Lindberg
2013-2-21
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Random Number Generation 的更多信息
2013-2-21
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
