block matrix multiplication with matlab cells

Question

I want to do block matrix multiplication of A and B using cells in Matlab. More specifically suppose

a=

 1     1     2     2
 1     1     2     2
 3     3     4     4
 3     3     4     4

b=

 2     2     4     4
 2     2     4     4
 6     6     8     8
 6     6     8     8

We can now convert a and b to cell arrays containing their blocks.

A = mat2cell(a,[2,2],[2,2])

ans = 

[2x2 double]    [2x2 double]
[2x2 double]    [2x2 double]

B = mat2cell(b,[2,2],[2,2])

ans = 

[2x2 double]    [2x2 double]
[2x2 double]    [2x2 double]

I need a function C=foo(A,B) that will return a cell array C such that the blocks of C are the blocks of matrix product A*B, for example in this case:

C{1,1} = A{1,1}*B{1,1} + A{1,2}*B{2,1}
C{1,2} = A{1,1}*B{1,2} + A{1,2}*B{2,2}
...

and cell2mat(C) should return:

ans =

28    28    40    40
28    28    40    40
60    60    88    88
60    60    88    88

• The reason I can't just do cell2mat(A)*cell2mat(B) is because in my application most blocks are zero and that would be inefficient.

• Even though most blocks are zero, I also can't do sparse(cell2mat(A))*sparse(cell2mat(B)) because the blocks that are non-zero are fully dense so that would also be inefficient.

Is there a better way to do this problem without using clunky and slow for loops? Thanks!

Edit: I wrote a small code that illustrates what I want to do. It is slow, however, and I'm wondering if there's a better way.

function C = celltimes(A,B,nn,blocksize)

C = cell(nn);
[C{:}] = deal(sparse(blocksize,blocksize));

for ii = 1:nn
   for jj = 1:nn
      row = A(ii,:);
      col = B(:,jj);
      for kk = 1:nn
         if ( nnz(row{kk}~=0) && nnz(col{kk}~=0) )
            C{ii,jj} = C{ii,jj}+row{kk}*col{kk};
         end
      end
   end
end

and a test code:

%test

nn = 3; %number of blocks
blocksize = 3; %block size

a = randi([0,10],nn*blocksize)
b = randi([0,10],nn*blocksize)

A = mat2cell(a,repmat(blocksize,[1,nn]),repmat(blocksize,[1,nn]));
B = mat2cell(b,repmat(blocksize,[1,nn]),repmat(blocksize,[1,nn]));

C = celltimes(A,B,nn,blocksize);

%verify result
c = a*b;
max(max(abs(cell2mat(C)-c)))

Answer 1

Your C from a, b will be obtained using:

C = a*b ;

a=[ 1     1     2     2
 1     1     2     2
 3     3     4     4
 3     3     4     4] ;

b=[ 2     2     4     4
 2     2     4     4
 6     6     8     8
 6     6     8     8] ;

C = a*b