How to write the Taylor expansion of the sum of several matrices?

dotk[ab_, k_] := Dot @@ ConstantArray[ab, k]
dotk[ab_, 0] := IdentityMatrix[Length@ab];
mInvSeries[a_, b_, n_] := Sum[(-1)^k Inverse[a].dotk[(b.Inverse[a]), k], {k, 0, n}]

Usage

a = 10 {{6, 3}, {2, 8}};
ai = mInvSeries[a, IdentityMatrix[2], 150];

(a.ai ) // N

(* {{0.981341, 0.00691085}, {0.00460723, 0.985948}} *)

How about using NonCommutativeMultiply as matrix multiplication?

Then, we can let -1 escape the noncommutative multiply

Unprotect[NonCommutativeMultiply];
NonCommutativeMultiply[H___, Times[-1, M_], T___] := - NonCommutativeMultiply[H, M, T]

And do the Taylor series about the big matrix by the small matrix by hand:

inv[n_][big_, small_] := Total@FoldList[
     #1 ** #2 &, 
     Inverse[big], 
     -(small ** Inverse[big]) & /@ Range[n]
]

Then, inv[3][A, B] gives

Inverse[A] - Inverse[A] ** B ** Inverse[A] + 
Inverse[A] ** B ** Inverse[A] ** B ** Inverse[A] - 
Inverse[A] ** B ** Inverse[A] ** B ** Inverse[A] ** B ** Inverse[A]

and inv[3][A, B + C] gives

Inverse[A] - Inverse[A] ** (B + C) ** Inverse[A] + 
Inverse[A] ** (B + C) ** Inverse[A] ** (B + C) ** Inverse[A] - 
Inverse[A] ** (B + C) ** Inverse[A] ** (B + C) ** 
Inverse[A] ** (B + C) ** Inverse[A]

In the same spirit to Dr. belisarius' answer:

abSeries[A_, B_, a_, n_] := Series[Inverse[A + a B], {a, 0, n}]

After obtaining the series you can get the matrix without the a by using

ReplaceAll[Normal[abSeries[A, B, a, n]],a->1]

There is definitely a way to generate the expression you quoted for "general" matrices in mathematica, but since I don't know the math in detail I don't know how to teach mathematica to derive it either. Of course programming a known expression is easy, similarly to what Dr. belisarius did.

Addition

This might be more akin to what was originally intended:

abSeries[A_, B_, n_] := Plus @@ NestList[-Dot[Inverse[A].B, #] &, Inverse[A], n]

It gives an expression for general expressions A and B meaning you do not have to input a literal matrix for the expressions for this to work. Of course to get the final literal result you still need to put in the matrices using ReplaceAll after the calculation.

If you want to use three matrices you can simply put in the addition of two matrices as A or B. i.e. B=C+D. You can then use Distribute and Factor to put the resulting expression into different forms.

Example:

In[1]:=ReplaceAll[abSeries[a, ee b, 5].(a + ee b), {a -> PauliMatrix[1], b -> PauliMatrix[2]}] // Simplify
Out[1]:={{1 + ee^6, 0}, {0, 1 + ee^6}}