Partitioned matrix operations

I have a $(2 \times 2)$ matrix whose elements are themselves $(2 \times 2)$ matrices (i.e., a partitioned matrix), e.g: $$ M = \begin{pmatrix} A & B \\ C & D \end{pmatrix}, \quad \quad \text{with } A,B,C,D \in \text{Mat}_{2 \times 2}. $$ Problem: I want to compute some things with this matrix, in particular an expression for $M^{n}$, I have tried using MatrixPower[M, n] in Mathematica but this does not work, (the error message says that the argument at position 1 is not a non-empty square matrix). If I pretend that $A,B,C,D$ are scalars then Mathematica will automatically assume multiplicative commutativity when computing powers of $M$.

Additional info: Furthermore, matrix multiplication does not work, if I try $M.M$ then I expect to obtain a result which looks like: $$ \begin{pmatrix} A.A + B.C & A.B+B.D \\ C.A+D.C &C.B+D.D \end{pmatrix}, $$ however, instead, I obtain a result in which the entries of $M.M$ are $(2 \times 2)$ matrices of $(2 \times 2)$ matrices (a bit complicated to explain).

Many thanks if someone can tell me how to compute matrix powers of a matrix whose elements are themselves matrices.

Answer

We can use the symbolic representation of the block matrix and perform the matrix operations with Inner and NonCommutativeMultiply. At the end we replace the matrix blocks with numeric or symbolic matrices. Below is an example.

The package and articles at "Noncommutative Algebra Package and Systems" might be of interest if these kind of manipulations are going to be done often and in more general form.

Here is a block matrix in symbolic form:

blockMat = {{A11, A12}, {A21, A22}};
MatrixForm[blockMat]

Here is the matrix multiplication using Inner and NonCommutativeMultiply.

blockMat2 = 
  Outer[Inner[NonCommutativeMultiply, #1, #2, Plus] &, blockMat, 
   Transpose[blockMat], 1];

MatrixForm[blockMat2, 1]

Note that the use of Inner directly corresponds to the block-matrix row-by-column multiplication.

This replaces the symbols for the matrix blocks with numerical matrices:

nBlockMat2 = 
  blockMat2 /. {A11 -> RandomReal[{0, 1}, {2, 2}], 
    A12 -> RandomReal[{10, 11}, {2, 2}], 
    A21 -> RandomReal[{100, 101}, {2, 2}], 
    A22 -> RandomReal[{1000, 1001}, {2, 2}]};

This command makes the sub-matrix multiplications:

MatrixForm[nBlockMat2 /. NonCommutativeMultiply -> Dot]

The above result with ArrayFlatten applied to it:

MatrixForm[ArrayFlatten[nBlockMat2 /. NonCommutativeMultiply -> Dot]]

Here is an image of the above in a Mathematica session: