# Dot Product of Block Matrices

How to multiply two matrices, which elements themselves are matrices? The following works but is there a more elegant way?

(A = Array[Subscript[a, ##] &, {3, 3, 2, 2}]) // MatrixForm (B = Array[Subscript[b, ##] &, {3, 1, 2}]) // MatrixForm dotdot[A_, B_] := Table[Sum[A[[i, l]].B[[l, j]], {l, Dimensions[A][]}], {i, Dimensions[A][]}, {j, Dimensions[B][]}];
dotdot[A, B] // MatrixForm One way is to turn them into ordinary matrices, take the dot product, and then ArrayReshape them into the form you want.

Af = ArrayFlatten[A];
Bf = Flatten[B];

ArrayReshape[Af.Bf, {3, 1, 2}] // MatrixForm Just to check:

ArrayReshape[Af.Bf, {3, 1, 2}] == dotdot[A, B]

(* True *)


You could use TensorContract instead:

r1 = dotdot[A, B];
r2 = TensorTranspose[
TensorContract[TensorProduct[A, B], {{2, 5}, {4, 7}}],
{1, 3, 2}
];

r1 === r2


True

Comparison

One advantage of the TensorContract approach is that if A and B are sparse arrays, then the result is also a sparse array. This is not true for the other answers. For example:

A=SparseArray @ RandomChoice[{.9,.1}->{0,1}, {5, 6, 7, 8}];
B=SparseArray @ RandomChoice[{.9,.1}->{0,1}, {6, 9, 8}];

r1 = dotdot[A, B];
r2 = TensorTranspose[
TensorContract[TensorProduct[A, B], {{2, 5}, {4, 7}}],
{1, 3, 2}
];

Normal[r1] === Normal[r2]