Implementing binomial expansion of (non-commuting) matrices in Mathematica

Question

Given two matrices $A$ and $B$, their Binomial expansion is, in general, given by (see this page)

$(A+B)^n = \sum\limits_{k=0, k=n[2]}^n \Bigg( \sum\limits_{r=0}^k \binom{k}{r} A^r B^{k-r} \Bigg) \Bigg(\frac{-C}{2}\Bigg)^{\frac{n-k}{2}} \frac{n!}{k!(\frac{n-k}{2})!}$

Here, $C=AB-BA$, is the commutator of $A$ and $B$. For matM = {{m11, m12}, {m21, m22}}; matN = {{n11, n12}, {n21, n22}};, how can one implement the above formula in Mathematica, say for simple case with $n=2$?

Answer 1

commutator[a_, b_] := a.b - b.a
mtxbinomial[a_, b_, n_] := 
 Block[{c = commutator[a, b], zeros = ConstantArray[0, Dimensions[a]]},
  If[c == zeros, 
   Message[commutator::args, "Matrices A and B commute!"]; Abort[]
  ];
  Sum[
   If[Mod[k, 2] == Mod[n, 2],
    Sum[Binomial[k, r] MatrixPower[a, r].MatrixPower[b, k - r], {r, 0, k}]. 
     MatrixPower[-c/2, (n - k)/2]*n!/(k! ((n - k)/2)!),
    zeros],
   {k, 0, n}]]

Testing $(A+B)^2$ gives True for some random matrices:

a = RandomReal[{-1, 1}, {5, 5}];
b = RandomReal[{-1, 1}, {5, 5}];
expected = MatrixPower[a + b, 2];
mtxbinomial[a, b, 2] == expected