Memory efficient alternative for meshgrid?

33 views (last 30 days)
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?

Accepted Answer

Fabio Freschi
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
  1 Comment
Walter Roberson
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(:))

Sign in to comment.

More Answers (2)

Eran
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);
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')
  1 Comment
Fabio Freschi
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.

Sign in to comment.


Bruno Luong
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

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!