## How to efficiently multiply many pairs (k-tuples) of matrices?

Asked by Tamar Friedlander

### Tamar Friedlander (view profile)

on 30 Jun 2013
Accepted Answer by the cyclist

### the cyclist (view profile)

I have a code in which I have to multiply many pairs of matrices. This part repeats many times so it is important to make it as efficient as possible.

Now I'm storing this as 2 3D matrices: for example:

```A_Pop = rand(5,5,500);
B_Pop = rand(5,5,500);
```

Now I'd like to multiply A(:,:,i)*B(:,:,i).

` tic; for ii=1:500, C_Pop(:,:,ii) = A_Pop(:,:,ii) * B_Pop(:,:,ii); end; toc;`
` Elapsed time is 0.001915 seconds.`

-----------------------------------------------------------------------------

Using cellfun or arrayfun was even slower than this loop.

``` multmat = @(x,y) x*y;
tic; C = cellfun(multmat, A, B,'UniformOutput', false);toc;```
` Elapsed time is 0.003202 seconds.`

-----------------------------------------------------------------------------

The most efficient way I found so far was to store A_Pop and B_Pop as blocks in a block-diagonal matrices, (and then store these huge matrices as sparse) and simply multiply the 2 matrices by each other:

```A_block = []; B_block = [];
for ii=1:500, A=rand(5,5); B=rand(5,5); A_block = blkdiag(A_block,A);
B_block = blkdiag(B_block,B);end;
A_bl = sparse(A_block); B_bl = sparse(B_block);
```

And then multiply:

` tic; C_bl=A_bl*B_bl; toc`
` Elapsed time is 0.000425 seconds.`

-----------------------------------------------------------------------------

Does anyone have an idea how to implement this more efficiently than the block-diagonal version?

That will be very helpful. Thanks, --Tamar.

Matt J

### Matt J (view profile)

on 30 Jun 2013

This is incidental to your question, but looping over

`    A_block = blkdiag(A_block,A);`

is an unnecessarily slow and painful way to build a multi-block block diagonal sparse matrix. Better options are

(1)

```    Acell=num2cell(A,[1,2]);
Acell = cellfun(@sparse,Acell,'uni',0);
A_block=blkdiag(Acell{:});```

(2)

```   template=kron(speye(N), ones(5));
idx=template>0;
A_block=template;
A_block(idx)=A(:);```

## Products

### the cyclist (view profile)

Answer by the cyclist

### the cyclist (view profile)

on 30 Jun 2013

#### 1 Comment

Tamar Friedlander

### Tamar Friedlander (view profile)

on 2 Jul 2013

Thanks, looks like this is what I need.

#### Join the 15-year community celebration.

Play games and win prizes!

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi