Memory efficient alternative for meshgrid?

I am generating a meshgrid to be able to calculate my result fast:

% x, y, z are some large vectors
[a,b,c] = meshgrid(x,y,z);

% s, t are constants, M some matrix
result = (((c*s - b*t).^2)./(a.^2 + b.^2 + c.^2)).*M;

This is actually working quite nicely. Unfortunately, for very large x,y,z, the meshgrid function is running out of memory.

How do I rewrite the meshgrid function to be memory efficient?

I had thought of three loops like this:

result = zeros(length(x), length(y), length(z));
for i = 1:lenght(x)-1
    for j = y = 1:lenght(y)-1
        for k = z = 1:lenght(z)-1
            b = ??
            c = ??
            result(i,j,k) = (((c*s - b*t).^2)./(x(i)^2 + y(j)^2 + z(k).^2));
        end
    end
end
result = result.*M;

What are the values for b and c?

How can I turn the outer for into a parfor?

[Stephen23]

https://www.mathworks.com/matlabcentral/fileexchange/74956-ndgridvecs

Answer 1

[Fabio Freschi]

This is how to make the three-loop version analogous to the meshgrid version

% some dummy values
N = 300;
x = linspace(1,10,N);
y = linspace(1,10,N);
z = linspace(1,10,N);
s = 1;
t = 1;
M = rand(N,N,N);

%% meshgrid
tic
% x, y, z are some large vectors
[a,b,c] = meshgrid(x,y,z);
% s, t are constants, M some matrix
result = (((c*s - b*t).^2)./(a.^2 + b.^2 + c.^2)).*M;
toc


%% three-loops
tic
% preallocation
result2 = zeros(length(x), length(y), length(z));

% note the order of the for-loop indices to mimic meshgrid
for iz = 1:length(x)
    for jx = 1:length(y)
        for ky = 1:length(z)
            result2(iz,jx,ky) = M(iz,jx,ky)*(((z(iz)*s - y(ky)*t).^2)./(x(jx)^2 + y(ky)^2 + z(iz).^2));
        end
    end
end
toc

% check results
norm(result(:)-result2(:))./norm(result(:))

However I don't see how you can avoid running out of memory: meshgrid is creating a N*N*N (with my notation) matrix, if it runs out of memory, also the preallocation of the result matrix will

result2 = zeros(length(x), length(y), length(z));

It is however true that in the second version you only have 2 matrices with dimensions N*N*N (M and result2) whereas in the first case you have 5 (a, b, c, M, result).

Note that according to my tests, the meshgrid version with vectorization is always faster than the version with three loops

[Walter Roberson]

Slightly more efficiently:

% some dummy values
N = 300;
x = linspace(1,10,N);
y = linspace(1,10,N);
z = linspace(1,10,N);
s = 1;
t = 1;
M = rand(N,N,N);

%% meshgrid
tic
% x, y, z are some large vectors
[a,b,c] = meshgrid(x,y,z);
% s, t are constants, M some matrix
result = (((c*s - b*t).^2)./(a.^2 + b.^2 + c.^2)).*M;
toc


%% three-loops
tic
% preallocation
result2 = zeros(length(x), length(y), length(z));

% note the order of the for-loop indices to mimic meshgrid
for iz = 1:length(x)
    X = x(iz);
    for jx = 1:length(y)
        Y = y(jx);
        for ky = 1:length(z)
            Z = z(ky);
            result2(iz,jx,ky) = M(iz,jx,ky)*(((Z*s - Y*t).^2)./(X^2 + Y^2 + Z.^2));
        end
    end
end
toc

% check results
norm(result(:)-result2(:))./norm(result(:))

Answer 2

[Eran]

You can do it faster and without the meshgrid memory allocation:

Common code

M=1; t=1; s=1;
x=(1:200);
y=(1:200);
z=(1:200);

your code:

tic; [a,b,c] = meshgrid(x,y,z); result = (((c*s - b*t).^2)./(a.^2 + b.^2 + c.^2)).*M;toc                
Elapsed time is 0.047975 seconds.

Save memory and time using:

tic; x=x(:); y=y(:).'; z=permute(z,[3 1 2]); result1 = (((c.*s - b.*t).^2)./(a.^2 + b.^2 + c.^2)).*M;toc
Elapsed time is 0.015030 seconds.

Verify that you get the same results:

all(result==result1,'all')

[Fabio Freschi]

Note that you are using a,b,c from the previous computation. Your method should read

tic; x=x(:); y=y(:).'; z=permute(z,[3 1 2]); result1 = (((z.*s - y.*t).^2)./(x.^2 + y.^2 + z.^2)).*M;toc

Still faster, anyway

Note also that the OP said M is a matrix.

Answer 3

[Bruno Luong]

Compute with auto-expansion capability after reshapeing vectors in appropriate dimensions rather than meshgrid/ndgrid

% [a,b,c] = meshgrid(x,y,z);
a = reshape(x, 1, [], 1);
b = reshape(y, [], 1, 1);
c = reshape(z, 1, 1, []);
result = (((c*s - b*t).^2)./(a.^2 + b.^2 + c.^2)).*M;