interpolate over a multi-dimensional function

I want to use interpn to interpolate the values of function V:R^4->R^3. I.e. size(V) = [C C C C 3], but interpn requires that the last dimension be of size 1.

So instead of doing:

anew = interpn(C, M, Y, K, V, a(:,1), a(:,2), a(:,3), a(:,4));

I have to do this:

anew = zeros(size(a_cmyk));
for i=1:4   
    anew(:,i) = interpn(C, M, Y, K, V(:,:,:,:,i),a(:,1), a(:,2), a(:,3), a(:,4));
end

And it really slows me down. Am I doing something wrong? Is there a way around this?

Answer 1

[the cyclist]

I think you may be mistaken, about the last dimension of V needing to be 1. In the example in the help file, their V is length 6 in the final dimension, and it works fine. (They do use a loop over the final dimension in that example, but only for the visualization.)

[Lior Shapira]

Here's a simple test:

c = 1:100;

m = 1:100;

y = 1:100;

V = rand([100 100 100 4]);

[C M Y] = ndgrid(c,m,y);

vals = interpn(C,M,Y,V, 1:50,1:50,1:50);

It doesn't work, returns an error:

??? Error using ==> interpn at 155

Wrong number of input arguments or some dimension of V is less than 2.

Answer 2

[Alec]

It seems interpn only works with scalar functions. If you issue

open interpn

you can see where it parses the input it expects that you have as many X1,X2,... as dimensions in V. So V has to be a scalar function: R^n --> R. You'll have to interpolate each coordinate function independently.

In your original post, you write that V is a function from R^4 --> R^3, but then in the code V seems to be a function from R^4 --> R^4 (you're looping from 1 to 4 rather than 1 to 3). Was that just a typo?

In the example you posted in the comment. V is a function from R^3 to R^4. So you'll have to interpolate each of the 4 coordinate functions separately, perhaps like this:

c = 1:100; 
m = 1:100; 
y = 1:100; 
V = rand([100 100 100 4]); 
[C M Y] = ndgrid(c,m,y); 
vals = zeros([size(1:50) 4]);
for ii = 1:4
  vals(:,:,ii) = interpn(C,M,Y,V(:,:,:,ii), 1:50,1:50,1:50);
end

(Note: vals here has size(vals) = [ 1 50 4] because the 50 interpolated points are defined on a "grid" that only has one row )

As far as performance, the for loop isn't really a "bad" for loop in terms of matlab optimization. Presuming the number of coordinates is much much smaller than the number of data points or interpolated points.

[the cyclist]

If I am not mistaken, this is pretty much the same coding strategy that Lior included with her original question.

See my answer, posted (time-wise) just after yours, in which I think I have the correct vectorization of the last dimension. Performance wise, I am not sure if it will be faster or not, but it is worth a shot.

Answer 3

[the cyclist]

Lior,

Putting this as a new answer, because I did not have a perfect understanding of what you were trying to do before.

I am going to put this code in without explanation, because I think you are going to understand what it is doing. If you do not, add a comment, and I will try to explain a bit what is happening. The code that does what you want is the one commented "Four planes".

c = 1:20;
m = 1:20;
y = 1:20;
z = 1:4;
V = rand([20 20 20 4]); 
[C M Y Z] = ndgrid(c,m,y,z);

% One plane:
ci = 1:10;
mi = 1:10;
yi = 1:10;
zi = ones(1,10);

vals = interpn(C,M,Y,Z,V,ci,mi,yi,zi);

% Four planes:
ci = repmat(1:10,[1 4]);
mi = repmat(1:10,[1 4]);
yi = repmat(1:10,[1 4]);
zi = [ones(1,10), 2*ones(1,10), 3*ones(1,10), 4*ones(1,10)];

vals = interpn(C,M,Y,Z,V,ci,mi,yi,zi);

[Alec]

Hi the cyclist,

What is the function you are trying to interpolate? Is it some function V from R^3 to R^4?

Then if I understand, you separate the 4 R^3 -> R coordinate functions of V into different planes in R^4, so then you can redefine V to be a scalar function from R^4 to R. You can certainly interpolate this function, but the result might not be correct. The separate coordinate functions shouldn't influence each other.

I think you can get away with in practice though. As long as the planes are far enough away from each other, the answer should be the same as correctly interpolating each coordinate separately.

In terms of vectorizing the for loop this is a neat trick and I like it, though upon inspection it turns out to be slower than the for loop.

Think about this, to compute the interpolated values for the coordinate function i of some vector function V you should only have to worry about the data values of coordinate i. This is of course what happens in the for loop method.

In the four planes solution you're telling the interpolation scheme to consider all data coordinates even though you've separated them into planes in a way that they shouldn't influence each other.

So in the for loop method to interpolate at some new point, at most the interpolation scheme can consider n data values. Where n is the number of sample points. For the four planes method, at most the interpolation scheme will consider 4*n data values. It's just a multiplicative factor of the number of coordinates, but getting a speed up was the whole point to vectorizing this in the first place (I humbly assume anyway).

Here's a benchmark example:

c= (1:100);

m = (1:100);

y = (1:50);

% build a 3D grid, our data points

[C M Y] = ndgrid(c,m,y);

% the four coordinate functions of V

f1 = @(x,y,z) x;

f2 = @(x,y,z) y;

f3 = @(x,y,z) z;

f4 = @(x,y,z) zeros(size(x));

% define V over our 3D grid, so V is function from R^3 to R^4

V = cat(4,f1(C,M,Y),f2(C,M,Y),f3(C,M,Y),f4(C,M,Y));

% Build another 3D grid, our interpolated points

ci = (1:0.5:100);

mi = (1:0.5:100);

yi = (1:0.5:50);

[Ci Mi Yi] = ndgrid(ci,mi,yi);

tic;

% separated solution

vals_sep = zeros([size(Ci) 4]);

for ii = 1:4

vals_sep(:,:,:,ii) = interpn(C,M,Y,V(:,:,:,ii),Ci,Mi,Yi);

end

toc

% "four planes" solution

z = 1:4;

% 4D grid for data points of all 4 coordinate functions

[C M Y Z] = ndgrid(c,m,y,z);

% 4D grid for interpolated points of all 4 coordinate functions

zi = 1:4;

[Ci Mi Yi Zi] = ndgrid(ci,mi,yi,zi);

tic;

vals = interpn(C,M,Y,Z,V,Ci,Mi,Yi,Zi);

toc

if(all(vals(:) == vals_sep(:)))

fprintf('Solutions are the same\n');

else

fprintf('Solutions differ\n');

end

On my machine this prints:

Elapsed time is 10.378915 seconds.

Elapsed time is 54.006997 seconds.

Solutions are the same

I suspect the vectorized form does especially badly here because it needs to store all 4 coordinate functions (input and output) in memory during interpn. So maybe it's not so bad on a machine with more than 4GB ram.

I'd be curious if the vectorized version is faster if you have more/faster ram/processors.

[the cyclist]

Alec, interesting stuff. I have to admit that I have not thought it thoroughly through, but my gut does not agree with you that the planes need to be "far apart" to "get away with it". By construction, there is no variation in what I defined as the z-axis, so I think you will theoretically get strictly the same solution regardless of the (arbitrary) distance between the z-planes. Although there may be some computer-y round-off issues, I guess.

Regarding performance, I did a quick bit of testing, and on the largest arrays that did not bring my machine to a whimpering halt, the vectorized method was about twice as fast. This is definitely a "your mileage may vary" situation. Lior, I suggest you do some test runs with your method and mine, and see what is faster.

[Lior Shapira]

Thanks for the discussion guys. I'm running tests to see if this can help me. Also I'm writing the interpolation part as a MEX file, which will hopefully speed things up.