How to vectorize a product of a tensor and a vector.

How could I vectorize the following loop?

A=randn(10,10,40);
b=randn(10,1);
a=zeros(40,1);
for i=1:40
    a(i)=b'*A(:,:,i)*b;
end

[madhan ravi]

Have you considered using James Tursa's FEX?:

https://www.mathworks.com/matlabcentral/fileexchange/25977-mtimesx-fast-matrix-multiply-with-multi-dimensional-support

[madhan ravi]

Another possiblity:

A=randn(10,10,40);
[m,n,p]=size(A);
b=randn(10,1);
AA=reshape(A,m,[],1)';
B=mat2cell(AA,repmat(n,1,size(AA,1)/n));
R=cellfun(@(x)b'*x*b,B);

Test result:

Elapsed time is 0.001497 seconds. % using a simple loop
Elapsed time is 0.062664 seconds. % using cellfun

From the test it's quite clear that the loop you have is at it's best shape and is by nature the fastest in this case , so why bother?

[Moslem Zamani]

Thanks so much for your favor.

Answer 1

[David Goodmanson]

Hi Moslem

Here is one way, and it does test out faster than the for loop. How much faster depends on the sizes of what I called m and n, which in your example are 10 and 40. But 2 - 80 times faster, using the method of total tic and toc elapsed time on a lot of repetitions of the same code, and around 15 times faster in the 10 and 40 case. (I know some people on this site are down on tic and toc, but total time for a lot of repetitions seems to me to be a valid real world test). There is more of an advantage the larger n is, which is the length of the for loop.

The code uses reshape, which can be pretty time consuming. [ I thought. See Guillaume's comment ]. There are any number of situations where reshapes and transposes take significantly more time than just using a for loop, but this seems to be a case where matrix manipulation wins out. (There is a transpose at the end, but it's only on a vector).

I tried in on complex to make sure it checked out vs the original.

m = 10;
n = 40;
% A=rand(m,m,n);
% b=rand(m,1);
A = 2*(rand(m,m,n) + i*rand(m,m,n)) - (1+i);   
b = 2*(rand(m,1) + i*rand(m,1)) - (1+i);

B = reshape(A,m,m*n);
bB = reshape(b'*B,m,n);
anew = b.'*bB;
anew = anew(:);
   

[Guillaume]

The code uses reshape, which can be pretty time consuming

reshape should always be near instantaneous as it doesn't affect the matrix elements. All it does is rewrite the matrix header from this matrix has m dimensions of size a,b,... to this matrix has n dimensions of size x,y,...

transpose and permute on the other hand are expensive operations, they need to reorder all the elements of the matrix (as well as rewrite the header).

[Moslem Zamani]

Hi David,

Thanks so much for your favor.

[David Goodmanson]

Hi Moslem,

You're welcome. And thanks to Guillaume for pointing out the reason behind why reshape is quick.