Multiply slices of 3D matrices without duplicating in memory
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Ivan Bioli
el 21 de Oct. de 2023
Editada: Bruno Luong
el 21 de Oct. de 2023
Hi,
I have two 3d matrices: S of size [m, n, n] and W of size [n, n, n] and I need to compute the matrix X of size [m, n, n^2] such that
X(:, :, n * (i-1) + j) = squeeze(W(:, :, i)) * squeeze(S(:, :, j))
I would like to avoid doing a double for loop on i and j, in order to keep things as much efficient as possible. I was able to avoid one of the two for loops by doing something like
idx = 1:n;
for i = 1:n
X(:, :, n*(i-1)+idx) = pagemtimes(S, squeeze(W(:, :, i)));
end
However, profiling the code I noticed that most of the time is spent in memory management. Indeed, if rewrite the code as
idx = 1:n;
for i = 1:n
aux = pagemtimes(S, squeeze(W(:, :, i)));
X(:, :, n*(i-1)+idx) = aux;
end
25% of the time is effectively spent in computing the product, while 75% of the time is spent in assigning X(:, :, n*(i-1)+idx) = aux. In case it is useful, in my case m = 100, n = 10, but I get a similar ratio also if I use m = 10000, n = 10. I am not sure whether this can be avoided, but as a possible solution I was trying to assemble the matrix X all together. The only idea that came to my mind was to do something like
j = repmat(1:n, [1, n]);
i = repmat(1:n, [n, 1]); i = i(:)';
X = pagemtimes(S(:, :, j), W_hat(:, :, i));
which (I believe) copies the matrices S and W under the hood, and is not much more efficient for this reason. Is there a way to avoid execute this operation without duplicating memory and without the slow memory management I discussed above?
Thanks,
Ivan
0 comentarios
Respuesta aceptada
Bruno Luong
el 21 de Oct. de 2023
Editada: Bruno Luong
el 21 de Oct. de 2023
This would do
m = 100;
n = 10;
S=rand(m,n,n);
W=rand(n,n,n);
X=pagemtimes(S,reshape(W,[n n 1 n]));
X=reshape(X,[m n n^2]);
% Check againts the original code
idx = 1:n;
for i = 1:n
aux = pagemtimes(S, squeeze(W(:, :, i)));
X2(:, :, n*(i-1)+idx) = aux;
end
norm(X(:)-X2(:))/norm(X2(:))
0 comentarios
Más respuestas (0)
Ver también
Categorías
Más información sobre Operators and Elementary Operations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!