# I need to multiply a series of matrices

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

-
Could you please show your starting code? It's hard to have nothing to start from. – Vitaliy Kaurov Dec 2 '12 at 6:01
Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2)Read the FAQs! 3) When you see good Q&A, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. ALSO, remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign – chris Dec 2 '12 at 11:09

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.

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}]];

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.

-
Yes, you could also take a List and change its head to Times if you want to get the normal multiplication of all elements in a list. But just to clarify: Dot is really very different from Times because conventional matrix multiplication is not commutative, whereas Times doesn't care about the order of the factors. BTW - strange that @NasserM.Abbasi appears here as your name. That's not you, right? – Jens Dec 2 '12 at 7:18
@NasserM.Abbasi lol, it's too late at night here... – Jens Dec 2 '12 at 7:51
is there a way to parallelize this? I have a similar problem where I need to multiply a lot of 2x2 matrices in a list. I used your method and there doesn't seem to be a difference between that and using a Do loop to multiply them. This method is more elegant but not faster. – lucian Jan 14 at 12:20
If your matrices are numerical, you could try to speed things up by using machine-precision arithmetic, e.g. by wrapping the list in N before doing the product. That would be much more of a gain than trying to parallelize the product. – Jens Jan 14 at 17:18
@lucian You can use the new UpTo to fix that problem. I'll add something to my answer shortly. – Jens Jan 14 at 19:28

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.

-
Interesting observation about the limit of factors (+1). I didn't dare to try it for 10000 symbolic matrices, but I chose 100 to compare to my method. With that, I can still make my (non-parallel) method run about twice as fast as ParallelCombine by choosing UpTo[10] in Partition. The more matrices, the better ParallelCombine seems to get. – Jens Jan 15 at 4:30

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.

-
However, it should be mentioned that the Fold approach is (on my machine, for matrixList= = RandomReal[2, {10000, 2, 2}]) an order of magnitude slower than Jens' approach using Partition. – Lukas Mar 31 at 11:42

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 *)

-