Multiply slices of 3D matrices without duplicating in memory

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

Answer 1

[Bruno Luong]

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(:))

ans = 0