Memory efficient alternative for meshgrid?

16 次查看（过去 30 天）
Vogel
Vogel 2020-6-11
评论： Chuang 2024-9-17
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?

请先登录，再回答此问题。

采纳的回答

Fabio Freschi
Fabio Freschi 2020-6-11
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 个评论
Walter Roberson
Walter Roberson 2023-7-9
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(:))

请先登录，再进行评论。

更多回答（2 个）

Eran
Eran 2023-7-8
编辑：Eran 2023-7-8
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 个评论
Fabio Freschi
Fabio Freschi 2023-7-8
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
Bruno Luong 2023-7-9
编辑：Bruno Luong 2023-7-9
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;

类别

在 Help Center 和 File Exchange 中查找有关 Data Type Identification 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by