Efficiently compute many inner products in matlab

Question

I am working on a project for which I need to compute a lot of inner products in high dimensions. I am aware that we should always try to vectorize operations in matlab, but I am not sure how to do this ...

Lets say we have two matrices, A and B of size N x d where N is the number of inner products to compute in d dimensions.

It is easy to implement using for-loops, but I suspect more efficient ways exist. A for-loop implementation might look like this:

innerprods=zeros(N,1);
for i=1:N
    innerprods(i)=A(i,:)*B(i,:)';
end

Does anyone have ideas how to vectorize this? I guess bsxfun should come into play at some point, but I cannot figure out how to use it for this ...

Thanks in advance!

Answer 1

Your suspicion regarding bsxfun is entirely justified! You can efficiently compute inner products using the following one-liner with bsxfun and the good old colon operator:

innerprods=sum(reshape(bsxfun(@times,A(:),B(:)),N,d),2);

For a more detailed explanation of what's going on here:

    bsxfun(@times,A(:),B(:)) --> element-wise product, returns Nd x 1 vector
    reshape( ^^^ , N, d)      --> reshape into N x d matrix
    sum( ^^^ , 2)             --> sum over columns to get N x 1 vector

EDIT: I did some timings to give some idea as to the performance difference. I used matlab R2010b (don't ask ...). The figure shows some empirical results using the loop in your question and the one-liner I suggest for N=500 and varying number of dimensions. Using loops is 2 to 8 times slower than the approach using bsxfun. The speed difference may be larger for newer versions of matlab due to better parallellization.

Answer 2

What about the simple:

sum(A.*B,2)

A simple test (R2013a WinVista 32bit Core duo):

N = 8e3;
A = rand(N);
B = rand(N);

tic
out = zeros(N,1);
for ii = 1:N
    out(ii) = A(ii,:)*B(ii,:)';
end
toc

tic
out2 = sum(A.*B,2);
toc

all(out-out2 < 1e5*eps) % note difference in precision

Times

loop     5.6 sec
multsum  0.8 sec

Additional tests on R2013a Win7 64 Xeon E5

Avg time loop:           2.00906 seconds
Avg time multSum:        0.18114 seconds
Avg time bsxfun:         0.18203 seconds
Avg time reshapeMultSum: 0.18088 seconds

The main take-away points:

• looping in this case is vastlly inefficient (expected);
• bsxfun() is totally redundant, although the overhead is not substantial (expected);
• Reshaping into a column and then back to a matrix is also expected to provide benefits due to the internal 'preference' of the MATLAB engine for columnwise operations (i.e. along rows first);
• for syntax clarity I still recommend multSum: sum(A.*B,2).

The testing suite (avg times on 100 trials, fixed matrix size 8e3, equality of results to 1e5*eps):

N = 8e3;
A = rand(N);
B = rand(N);

tic
for r = 1:100
    out = zeros(N,1);
    for ii = 1:N
        out(ii) = A(ii,:)*B(ii,:)';
    end
end
sprintf('Avg time loop: %.5f seconds', toc/100)

tic
for r = 1:100; out2 = sum(A.*B,2);                                  end
sprintf('Avg time multSum: %.5f seconds', toc/100)

tic
for r = 1:100; out3 = sum(reshape(bsxfun(@times,A(:),B(:)),N,N),2); end
sprintf('Avg time bsxfun: %.5f seconds', toc/100)

tic
for r = 1:100; out4 = sum(reshape(A(:).*B(:),N,N),2);               end
sprintf('Avg time reshapeMultSum: %.5f seconds', toc/100)