# How can we multiply nested matrices?

If we have nested matrices, such as:

$$\begin{bmatrix} \begin{bmatrix} a_{1,1} & a_{1,2}\\ a_{2,1} & a_{2,2}\\ \end{bmatrix} \begin{bmatrix} b_{1,1} & b_{1,2}\\ b_{2,1} & b_{2,2}\\ \end{bmatrix} \\ \begin{bmatrix} c_{1,1} & c_{1,2}\\ c_{2,1} & c_{2,2}\\ \end{bmatrix} \begin{bmatrix} d_{1,1} & d_{1,2}\\ d_{2,1} & d_{2,2}\\ \end{bmatrix} \\ \end{bmatrix}$$

...and we wish to multiply two of these matrices together, how can we accomplish this?

Here I assume that the conventional rules of matric multiplication apply. For example, we treat each "sub-matrix" (or block matrix, if you prefer), for lack of a better term, as an entry in the main matrix. Then we recursively apply the rules of matrix multiplication.

AN ATTEMPT TO EXPLAIN

So, for example, the matrix from above could be represented as:

$$\begin{bmatrix} a & b\\ c & d\\ \end{bmatrix}$$

...and multiplied by another matrix to get:

$$\begin{bmatrix} a & b\\ c & d\\ \end{bmatrix} * \begin{bmatrix} e & f\\ g & h\\ \end{bmatrix} = \begin{bmatrix} (ae+bg) & (af+bh)\\ (ce+dg) & (cf+dh)\\ \end{bmatrix}$$

Where each entry above is actually another matrix. For example, the entry $(ae+bg)$ is actually:

$$(ae+bg) = \begin{bmatrix} a_{1,1} & a_{1,2}\\ a_{2,1} & a_{2,2}\\ \end{bmatrix} * \begin{bmatrix} e_{1,1} & e_{1,2}\\ e_{2,1} & e_{2,2}\\ \end{bmatrix} + \begin{bmatrix} b_{1,1} & b_{1,2}\\ b_{2,1} & b_{2,2}\\ \end{bmatrix} * \begin{bmatrix} g_{1,1} & g_{1,2}\\ g_{2,1} & g_{2,2}\\ \end{bmatrix}$$

• How do you have this stored? – Andrew Jaffe Nov 22 '13 at 16:53
• Mathematica works with general tensors of arbitrary dimensions, not just (2-dimensional) matrices. Dot can operate on any pair of tensors (for as long as their sizes match). It will contract the last index of the first tensor with the first index of the last tensor. TensorContract allows contracting arbitrary indices. – Szabolcs Nov 22 '13 at 17:01
• @AndrewJaffe: I'm open to any storage method. My idea is to store the smallest matrices as variables, and work from there. But again, any method that works is a big help. I haven't tried very much yet, so I'm open to ideas. – Matt Groff Nov 22 '13 at 17:02
• @Szabolcs: I don't know much about tensors. Could you incorporate them into an answer? – Matt Groff Nov 22 '13 at 17:03
• @MattGroff Not at this moment for lack of time. But this is not a difficult topic, and I think you should really read up on it. Everything will become clear and easy once you start using the Einstein convention to write out these contractions. (No need to learn about covariant and contravariant for your application, avoid books with those for now.) – Szabolcs Nov 22 '13 at 17:06

Update note: All but matrixDot2 work on general matrices with matrix entries. On the other hand, all but matrixDot2 unpack packed arrays. Only matrixDot and matrixDot4 return packed arrays. So perhaps there's still yet a method that that is more efficient.

matrixDot4[mat1_, mat2_] :=
Transpose[TensorContract[mat1.Transpose[mat2, {2, 3, 1, 4}], {{2, 4}}], {1, 3, 2, 4}];

matrixDot3[mat1_, mat2_] := Block[{Identity},
Inner[Dot, Map[Identity, mat1, {2}], Map[Identity, mat2, {2}], Plus]
];

matrixDot2[{{a_, b_}, {c_, d_}}, {{e_, f_}, {g_, h_}}] := (* similar to Jason B's *)
Evaluate[{{a, b}, {c, d}} . {{e, f}, {g, h}} /. Times -> Dot];

matrixDot[mat1_, mat2_] :=
First@Transpose[
ListCorrelate[Transpose@mat1, Transpose /@ mat2, {1, -1}, 0, Dot, Plus, 1],
{1, 2, 4, 3, 5}]


Example matrices:

Clear[f1, f2, a, n, c, d, e, f, g];
f1[1, 1, x__] := a[x]; f1[1, 2, x__] := b[x]; f1[2, 1, x__] := c[x]; f1[2, 2, x__] := d[x];
f2[1, 1, x__] := e[x]; f2[1, 2, x__] := f[x]; f2[2, 1, x__] := g[x]; f2[2, 2, x__] := h[x];

mat1 = Array[f1, {2, 2, 2, 2}];
mat2 = Array[f2, {2, 2, 2, 2}];

(* {{{{a[1, 1], a[1, 2]}, {a[2, 1], a[2, 2]}},
{{b[1, 1], b[1, 2]}, {b[2, 1], b[2, 2]}}},
{{{c[1, 1], c[1, 2]}, {c[2, 1], c[2, 2]}},
{{d[1, 1], d[1, 2]}, {d[2, 1], d[2, 2]}}}}  *)

(* {{{{e[1, 1], e[1, 2]}, {e[2, 1], e[2, 2]}},
{{f[1, 1], f[1, 2]}, {f[2, 1], f[2, 2]}}},
{{{g[1, 1], g[1, 2]}, {g[2, 1], g[2, 2]}},
{{h[1, 1], h[1, 2]}, {h[2, 1], h[2, 2]}}}}  *)


The product:

matrixDot[mat1, mat2]

(* {{
{{a[1, 1] e[1, 1] + a[1, 2] e[2, 1] + b[1, 1] g[1, 1] + b[1, 2] g[2, 1],
a[1, 1] e[1, 2] + a[1, 2] e[2, 2] + b[1, 1] g[1, 2] + b[1, 2] g[2, 2]},
{a[2, 1] e[1, 1] + a[2, 2] e[2, 1] + b[2, 1] g[1, 1] + b[2, 2] g[2, 1],
a[2, 1] e[1, 2] + a[2, 2] e[2, 2] + b[2, 1] g[1, 2] + b[2, 2] g[2, 2]}},

{{a[1, 1] f[1, 1] + a[1, 2] f[2, 1] + b[1, 1] h[1, 1] + b[1, 2] h[2, 1],
a[1, 1] f[1, 2] + a[1, 2] f[2, 2] + b[1, 1] h[1, 2] + b[1, 2] h[2, 2]},
{a[2, 1] f[1, 1] + a[2, 2] f[2, 1] + b[2, 1] h[1, 1] + b[2, 2] h[2, 1],
a[2, 1] f[1, 2] + a[2, 2] f[2, 2] + b[2, 1] h[1, 2] + b[2, 2] h[2, 2]}}},
{
{{c[1, 1] e[1, 1] + c[1, 2] e[2, 1] + d[1, 1] g[1, 1] + d[1, 2] g[2, 1],
c[1, 1] e[1, 2] + c[1, 2] e[2, 2] + d[1, 1] g[1, 2] + d[1, 2] g[2, 2]},
{c[2, 1] e[1, 1] + c[2, 2] e[2, 1] + d[2, 1] g[1, 1] + d[2, 2] g[2, 1],
c[2, 1] e[1, 2] + c[2, 2] e[2, 2] + d[2, 1] g[1, 2] + d[2, 2] g[2, 2]}},

{{c[1, 1] f[1, 1] + c[1, 2] f[2, 1] + d[1, 1] h[1, 1] + d[1, 2] h[2, 1],
c[1, 1] f[1, 2] + c[1, 2] f[2, 2] + d[1, 1] h[1, 2] + d[1, 2] h[2, 2]},
{c[2, 1] f[1, 1] + c[2, 2] f[2, 1] + d[2, 1] h[1, 1] + d[2, 2] h[2, 1],
c[2, 1] f[1, 2] + c[2, 2] f[2, 2] + d[2, 1] h[1, 2] + d[2, 2] h[2, 2]}}
}}    *)


matrixDot[mat1, mat2] == matrixDot2[mat1, mat2] ==
matrixDot3[mat1, mat2] == matrixDot4[mat1, mat2]
(* True *)

• I don't know if that fits what he wants exactly. In your result, the [[1,1]] element isn't purely a.e + b.g, but a mixture of a.e + b.g and a.f + b.h. – Jason B. Nov 22 '13 at 20:37
• These are ok. Although I ended up with 4-tuples after these were evaluated. Anyways, thanks for the ideas! – Matt Groff Nov 23 '13 at 3:43
• @MattGroff 4-tuples? I get an array with Dimensions of {2, 2, 2, 2}. – Michael E2 Nov 23 '13 at 3:45
• What I mean is that I get something like f1[1,1,1,1] or f2[2,1,2,2]. For some reason, it won't modify these to be a[1,1] or g[2,2], for example. – Matt Groff Nov 23 '13 at 3:50
• @MattGroff That's odd. Check what you get when you execute ?f1. You should see the definition of f1. If not, you didn't execute the definition. – Michael E2 Nov 23 '13 at 3:52

I spent some time looking at the documentation for TensorProduct and TensorContract but it was all greek to me (I wish I was more mathematically minded). But some trial and error gave me what I think you are looking for.

For the innerdot function, I have to wrap it in Quiet since it gives error messages, but still outputs what I want.

Edit: I changed it to work with square matrices that are an arbitrary size, but still have TensorRank of 4.

mat[lab_] :=
Table[ToExpression[lab <> IntegerString[n] <> IntegerString[m]], {n,
2}, {m, 2}];

m1 = {{mat["a"], mat["b"]}, {mat["c"], mat["d"]}};
m2 = {{mat["e"], mat["f"]}, {mat["g"], mat["h"]}};

MatrixDot[mat1_, mat2_] := Module[{dim, d1, d2, dummyfunc},
If[Dimensions[mat1] != Dimensions[mat2] ||
Dimensions[mat1][] != Dimensions[mat1][],
Print["Matrices not square or not the same size"];,
dim = Length[mat1];
dummyfunc[d1_, d2_] =
Quiet[Inner[Dot, Table[d1[[n, m]], {n, dim}, {m, dim}],
Table[d2[[n, m]], {n, dim}, {m, dim}], Plus]];
dummyfunc[mat1, mat2]]];

MatrixDot[m1, m2]

(*  {{{{a11 e11+a12 e21+b11 g11+b12 g21,a11 e12+a12 e22+b11 g12+b12 g22},
{a21 e11+a22 e21+b21 g11+b22 g21,a21 e12+a22 e22+b21 g12+b22 g22}},
{{a11 f11+a12 f21+b11 h11+b12 h21,a11 f12+a12 f22+b11 h12+b12 h22},
{a21 f11+a22 f21+b21 h11+b22 h21,a21 f12+a22 f22+b21 h12+b22 h22}}},
{{{c11 e11+c12 e21+d11 g11+d12 g21,c11 e12+c12 e22+d11 g12+d12 g22},
{c21 e11+c22 e21+d21 g11+d22 g21,c21 e12+c22 e22+d21 g12+d22 g22}},
{{c11 f11+c12 f21+d11 h11+d12 h21,c11 f12+c12 f22+d11 h12+d12 h22},
{c21 f11+c22 f21+d21 h11+d22 h21,c21 f12+c22 f22+d21 h12+d22 h22}}}}  *)


Or, in MatrixForm It seems as though this should work:

Inner[Dot, m1, m2, Plus]


but that returns a rank 6 tensor.