Vectorizing matrix multiplication inside a tensor

Question

I'm having some trouble vectorizing a part of my code. I have a (n,n,m) tensor, and I want to multiply each slice in m by a second (n by n) matrix (NOT element wise).

Here's what it looks like as a for-loop:

Tensor=zeros(2,2,3);
Matrix = [1,2; 3,4];

for j=1:n
    Matrices_Multiplied = Tensor(:,:,j)*Matrix;
    Recursive_Matrix=Recursive_Matrix + Tensor(:,:,j)/trace(Matrices_Multiplied);
end

How do I perform matrix multiplication for individual matrices inside a tensor in a vectorized manner? Is there a built-in function like tensor-dot that can handle this or is it more clever?

score 1 · Accepted Answer · edited May 23 '17 at 10:28

1

Bsxfunning and using efficient matrix-multiplication, we could do -

% Calculate trace values using matrix-multiplication
T = reshape(Matrix.',1,[])*reshape(Tensor,[],size(Tensor,3));

% Use broadcasting to perform elementwise division across all slices
out = sum(bsxfun(@rdivide,Tensor,reshape(T,1,1,[])),3);

Again, one can replace the last step with one more matrix-multiplication for possible further boost in performance. Thus, an all matrix-multiplication dedicated solution would be -

[m,n,r] = size(Tensor);
out = reshape(reshape(Tensor,[],size(Tensor,3))*(1./T.'),m,n)

Runtime test

Benchmarking code -

% Input arrays
n = 100; m = 100;
Tensor=rand(n,n,m);
Matrix=rand(n,n);
num_iter = 100; % Number of iterations to be run for

tic
disp('------------ Loopy woopy doops : ')
for iter = 1:num_iter
    Recursive_Matrix = zeros(n,n);
    for j=1:n
        Matrices_Multiplied = Tensor(:,:,j)*Matrix;
        Recursive_Matrix=Recursive_Matrix+Tensor(:,:,j)/trace(Matrices_Multiplied);
    end
end
toc, clear iter  Recursive_Matrix  Matrices_Multiplied

tic
disp('------------- Bsxfun matrix-mul not so dull : ')
for iter = 1:num_iter
    T = reshape(Matrix.',1,[])*reshape(Tensor,[],size(Tensor,3));
    out = sum(bsxfun(@rdivide,Tensor,reshape(T,1,1,[])),3);
end
toc, clear T out

tic
disp('-------------- All matrix-mul having a ball : ')
for iter = 1:num_iter
    T = reshape(Matrix.',1,[])*reshape(Tensor,[],size(Tensor,3));
    [m,n,r] = size(Tensor);
    out = reshape(reshape(Tensor,[],size(Tensor,3))*(1./T.'),m,n);
end
toc

Timings -

------------ Loopy woopy doops : 
Elapsed time is 3.339464 seconds.
------------- Bsxfun matrix-mul not so dull : 
Elapsed time is 1.354137 seconds.
-------------- All matrix-mul having a ball : 
Elapsed time is 0.373712 seconds.

edited May 23 '17 at 10:28

Community

1
1

answered Jul 17 '16 at 20:10

Divakar

218,885
19
262
358

Is there a way to vectorize out that second for loop (the one with "iter" as the index)? – Steven Sagona Jul 17 '16 at 23:58
@StevenSagona That loop : `for iter = 1:num_iter` is just for benchmarking so that we are running the code for enough time, so that the runtime measurements are as accurate as possible. If you are unfamiliar with timing MATLAB codes, take a look here for using timeit : http://stackoverflow.com/a/29719681/3293881 and here for another example on tic-toc with warm-up : http://stackoverflow.com/a/26137901/3293881 – Divakar Jul 18 '16 at 03:07
Sorry about that, I should've went through your code before commenting. I understand everything in the code except for this line: `(reshape(Tensor,[],size(Tensor,3))*(1./T.')` I believe this is doing something different than what the origional for-loop does. I ran this code when n, m were small and saw that the vectorized code doesn't produce the same result as the origional. I think it's the matrix multiplication at the end. I thought maybe switching to elementwise multiplication would work but it returns an error. – Steven Sagona Jul 18 '16 at 05:43
@StevenSagona Let me ask you - Aren't you initializing `Recursive_Matrix` as `zeros(n,n)` before going into the loop in your original code? – Divakar Jul 18 '16 at 07:26
Technically, before the for loop I initialize Recursive_Matrix as something different: Recursive_Matrix(1:length(Matrix),1:length(Matrix))=0; But I don't have to run my code to see a difference. If I run your code for n=2, I see that Recursive_Matrix and out produce different outputs for the same set of random numbers. – Steven Sagona Jul 18 '16 at 16:47
@StevenSagona Well at my end, results have matched for both approaches against the original loopy solution for sizes from small (n=2) to big (n=100). If you are seeing differences in values, even for `n = 2`, maybe post that data into the question as an edit section and show us the outputs? – Divakar Jul 18 '16 at 17:03
If m does not equal n, then they are different. I can fix it on my end if I make the for-loop end at m instead of n. But I think you're correct and I just made a typo with the for-loop in my question, because in my real code it goes to m. – Steven Sagona Jul 18 '16 at 21:55
@StevenSagona Well `n` being equal to or not to `m` shouldn't fail the listed approaches. – Divakar Jul 18 '16 at 22:53

Vectorizing matrix multiplication inside a tensor

1 Answers1