I need to multiply a series of matrices

Question

I need to perform a dot product on a large number of 2 x 2 matrices that I have defined to differ by even/odd subscripts (Ex: D-odd=x but D-even=y). The $\Pi$ function will not work because it does not perform dot products. Does anyone know how I can carry out this large series of matrix multiplications. Thanks

Answer 1

If you have a list of matrices, you can apply a dot multiplication to all the matrices by changing the Head of the list to Dot. To illustrate this, I'll first define a set of matrices symbolically:

matrixList = With[{numberOfMatrices = 2},
  Table[Array[Subsuperscript["M", Row[{#1, #2}], i] &, {2, 2}], {i, 
    numberOfMatrices}]
  ]

$\left\{\left( \begin{array}{cc} \text{M}_{11}^1 & \text{M}_{12}^1 \\ \text{M}_{21}^1 & \text{M}_{22}^1 \\ \end{array} \right),\left( \begin{array}{cc} \text{M}_{11}^2 & \text{M}_{12}^2 \\ \text{M}_{21}^2 & \text{M}_{22}^2 \\ \end{array} \right)\right\}$

This contains just two matrices for display reasons, but you can change numberOfMatrices to anything you like.

Now the answer to your question is just:

Apply[Dot, matrixList]

$\left( \begin{array}{cc} \text{M}_{11}^1 \text{M}_{11}^2+\text{M}_{12}^1 \text{M}_{21}^2 & \text{M}_{11}^1 \text{M}_{12}^2+\text{M}_{12}^1 \text{M}_{22}^2 \\ \text{M}_{11}^2 \text{M}_{21}^1+\text{M}_{21}^2 \text{M}_{22}^1 & \text{M}_{12}^2 \text{M}_{21}^1+\text{M}_{22}^1 \text{M}_{22}^2 \\ \end{array} \right)$

Edit

I noticed that there is a significant speedup over the above method, and also over the equivalent approach using Fold mentioned by J.M., if I partition the list of matrices beforehand and do the dot products in two steps:

With[{numberOfMatrices = 17},
  matrixList1 = RandomReal[{-100, 100}, {numberOfMatrices, 2, 2}]; 
  matrixList = 
   Table[Array[Subsuperscript["M", Row[{#1, #2}], i] &, {2, 2}], {i, 
     numberOfMatrices}]];


comparisonRule = Thread[Flatten[matrixList] -> Flatten[matrixList1]];

AbsoluteTiming[res1 = Apply[Dot, matrixList];]

(* ==> {0.060968, Null} *)

AbsoluteTiming[
 res2 = Apply[Dot, Dot @@@ Partition[matrixList, UpTo[4]]];]

(* ==> {0.00059, Null} *)

(res1 /. comparisonRule) == (res2 /. comparisonRule)

(* ==> True *)

The speed gain is more dramatic for larger lists of matrices, and one may have to play with the partition size. I used the command UpTo in Partition to allow the remainder partition to be of variable size.

If you want to test this for larger lists, I suggest omitting the test on the last line because it will be very slow (I already tried to speed it up by using random numbers in place of the symbolic matrix elements).

The original idea was to use ParallelMap on the partition, but the overhead due to parallelization eats up the speed gain with the tested example.

Answer 2

You can ParallelCombine as well. This has the added benefit of working on a list of more than the iteration limit of the session. My limit is 4096 so I would not be able to dot product more than 4096 matrices without using this or a similar method.

matrices = RandomReal[2, {10000, 2, 2}];

ParallelCombine[Dot[Sequence @@ ##] &, matrices, Dot]

Hope this helps.

Answer 3

For very long lists, it may be beneficial to iterate Jens' answer. This can be achieved using NestWhile

listmultiplier[list_, partitionwidth_: 5] := NestWhile[Dot @@@ Partition[#, partitionwidth, partitionwidth, 1, {}] &, list, Dimensions[#][[1]] > 1 &][[1]]

This code partitions the list into sublists containing at most partitionwidth matrices, applies Dot to each sublist and iterates until the final list contains a single matrix. As I am using version 10.2, the UpTo command is not available but can be replicated using the optional commands in Partition

For large lists, iterating can offer substantial savings:

matrixlist = RandomReal[{0, 1}, {4000, 20, 20}];

AbsoluteTiming[method1 = Apply[Dot, matrixlist];]

(* ==> {26.4274, Null} *)

AbsoluteTiming[method2 = Apply[Dot, Dot @@@ Partition[matrixlist, 5, 5, 1, {}]];]

(* ==> {5.50295, Null} *)

AbsoluteTiming[method3 = listmultiplier[matrixlist];]

(* ==> {0.204352, Null} *)

Chop[method2] == Chop[method1]

(* ==> True *)

Chop[method3] == Chop[method1]

(* ==> True *)

Answer 4

An alternative is to use Fold[] along with Dot[]. Using Jens's example, you can do Fold[Dot, IdentityMatrix[2], matrixList]. Note that your initial identity matrix must have dimensions that conform with the dimensions of your other matrices.