Matrix multiplication: exponentiation before summation

5 views (last 30 days)
Normal matrix multiplication would give C=A*B, such that C_ij = sum_k { A_ik * B_kj }
I want to compute a C* such that C*_ij = sum_k { exp( A_ik * B_kj ) }
That is, I want to exponentiate all the results of multiplication before they are summed into a single cell. How should I code this in matlab?

Accepted Answer

hyqneuron
hyqneuron on 13 Dec 2014
Thanks for your help John, Azzi, and Image Analyst.
John's answer is very thorough indeed. Here's something I wrote that's reasonably readable, equal in speed as C3, and avoids memory blowup completely.
C4=zeros(500,400);
for j=1:400
C4(:,j) = sum(exp(bsxfun(@times, A, B(:,j)')),2);
end
I'm only curious as to how to make it a bit faster. Exponentiating (500*400*1000) elements and doing a normal matrix multiplication (500,1000) with (1000, 400) only takes half as much time as C3 or C4. And exponentiation is certainly the most expensive part here.

More Answers (3)

Azzi Abdelmalek
Azzi Abdelmalek on 12 Dec 2014
Edited: Azzi Abdelmalek on 12 Dec 2014
EDIT
A=randi(5,2,2)
B=randi(5,2,3)
D=sum(exp(bsxfun(@times,permute(reshape(A',size(A',1),1,[]),[2 1 3]),B')),2)
out=D(:,:)'
%or with a for loop
A=randi(5,2,2)
B=randi(5,2,3)
n=size(A,1);
m=size(B,2);
out1=zeros(n,m);
for k=1:n
for p=1:m
out1(k,p)=sum(exp(A(k,:).*B(:,p)'));
end
end

John D'Errico
John D'Errico on 12 Dec 2014
Edited: John D'Errico on 12 Dec 2014
To be honest, I would not bother to answer this basic question if there was not something of some technical interest that can be a good teaching point.
You solve it WITH A LOOP. Actually, with loops. Given that you want to do something that is non-standard, you will need to write the code.
You will loop over each combination of i and j. Inside the double loop, take the element-wise products of row with columns, exponentiate, then sum. Do I really need to write this?
[r,n] = size(A);
[~,c] = size(B);
% note that I carefully preallocate C to the final size.
C = zeros(r,c);
for i = 1:r
for j = 1:c
C(i,j) = sum(exp(A(i,:).*B(:,j).'));
end
end
There is nothing wrong with simple loops when they are a logical option. The above is a basic matrix multiply, with an exponential inside, as you wanted. Effectively, I wrote it as a triply nested loop, with the outer loops explicit for loops. The innermost loop is an implicit one, where I create a vector of products inside the line, then exponentiate and sum them. So the above expression was mildly vectorized, since I could have written it all as an explicitly triply nested set of for loops.
Is there a fully "vectorized" alternative? (Why did I answer this question at all?) Yes, but you need to remember that vectorized code often trades off time or memory for "elegance" of the code, as well as replacing explicit loops with implicit loops. The trick is to recognize that a matrix multiply can be written using some alternatives. One way is to employ repmat calls, chosen to avoid the explicit loops. For example...
C2 = squeeze(sum(exp(repmat(A,[1 1 c]).*permute(repmat(B,[1 1 r]),[3 1 2])),2));
As you can see, the repmat vectorized code is far less readable than the simple loop. There are some who see elegance in the above line of code. But I also see an unreadable mess, that would be far less easy to debug in a crisis.
One problem with many codes that employ repmat is they tend to be less memory efficient, because we are forced to create additional matrices inline that could be quite large depending on the value of n. So a sometimes useful tool is bsxfun.
C3 = squeeze(sum(exp(bsxfun(@times,A,permute(B,[3 1 2]))),2));
Again, wholly a mess to read. No fun if I had to debug it. Did I get those parens in the right places? We can improve the readability somewhat by breaking that line down into several lines. Adding comments would help too.
For some simple test case matrices A and B...
A = rand(2,3);
B = rand(3,4);
C
C =
4.1459 5.3934 3.6706 4.1113
3.9077 4.7394 3.6274 3.7322
C2
C2 =
4.1459 5.3934 3.6706 4.1113
3.9077 4.7394 3.6274 3.7322
C3
C3 =
4.1459 5.3934 3.6706 4.1113
3.9077 4.7394 3.6274 3.7322
you can see the three methods all were the same in their end result. In fact, if we take out the exp call, this would have been a way to re-write a simple matrix multiply, A*B. Far less efficient of course, because MATLAB uses highly optimized code for that multiply.
Next, how did these variations compare in terms of time? I'll choose some medium large arrays to do a test.
A = rand(500,1000);
B = rand(1000,400);
tic
[r,n] = size(A);
[~,c] = size(B);
C = zeros(r,c);
for i = 1:r
for j = 1:c
C(i,j) = sum(exp(A(i,:).*B(:,j).'));
end
end
toc
Elapsed time is 5.247200 seconds.
tic
C2 = squeeze(sum(exp(repmat(A,[1 1 c]).*permute(repmat(B,[1 1 r]),[3 1 2])),2));
toc
Elapsed time is 7.145553 seconds.
tic
C3 = squeeze(sum(exp(bsxfun(@times,A,permute(B,[3 1 2]))),2));
toc
Elapsed time is 1.969375 seconds.
(Yes, I used tic and toc, even though normally I would recommend use of timeit. In this case, one of those alternatives used loops, so I would have had to encapsulate them into a function to use timeit. And since the time required was reasonably significant, tic and toc will give acceptable results.)
I had to be careful here. If I made A and B too large, then I might have run into memory problems. The intermediate arrays could be come HUGE. Here they are of size r*c*n elements, so
500*400*1000*8
ans =
1.6e+09
or roughly 1.6 gigabytes of RAM.
As you can see though, the bsxfun alternative was fastest, the repmat version slowest. Vectorized code need not ALWAYS be faster code! We can often gain time using vectorized code, but it does come at some costs. One cost is it makes your code unreadable. I'd strongly suggest that when you DO write a code fragment like those I have offered, you add some comments that clearly explain what was done in the line for future purposes during debugging.

Image Analyst
Image Analyst on 12 Dec 2014
You mean like this:
A = randi(9, 3, 5)
B = randi(9, 5, 2)
C = A * B
C_Star = sum(exp(C(:)))

Categories

Find more on Operating on Diagonal Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!