# Memory efficient alternative for meshgrid?

27 views (last 30 days)
Vogel on 11 Jun 2020
Edited: Bruno Luong on 9 Jul 2023
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 on 11 Jun 2020

Fabio Freschi on 11 Jun 2020
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 on 9 Jul 2023
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(:))

Eran on 8 Jul 2023
Edited: Eran on 8 Jul 2023
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);
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 on 8 Jul 2023
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.

Bruno Luong on 9 Jul 2023
Edited: Bruno Luong on 9 Jul 2023
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;

### Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

R2020a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!