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

When $A$ and $B$ are matrices, we have the following Taylor expansion for the inversion function together with basic information on convergence: $$(A+B)^{-1}=A^{-1}-A^{-1}BA^{-1}+A^{-1}BA^{-1}BA^{-1}-A^{-1}BA^{-1}BA^{-1}BA^{-1}+...,$$ provided $Norm[A^{-1}B]<1$. So, my question is how we can write a code in Mathematica to write this approximation up to any desired term? For example, to obtain $(I+B)^{-1}$ up to six terms by keeping the fact that matrix multiplications are not commutative ($I$ is the identity matrix).

Can we also use three matrices inside, e.g. $(I+B+C)^{-1}$ and then compute the expansion using Mathematica?

Thank you in advance for any tips or tricks.

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[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[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]

• Wow, amazing. This is really helpful and complete. – M.J.2 Mar 22 '16 at 10:22
• Can you extend your answer to solve this unanswered question: <mathematica.stackexchange.com/questions/109825/…> – M.J.2 Mar 22 '16 at 10:23
• I can't get Mma to reproduce your analytic results there. Edit: actually, maybe I can. I'll take a look. – evanb Mar 22 '16 at 18:22
• Thank you. Please take a look. It is very similar to your above answer but with further matrix matrix multiplications. You have a better insight to this noncommutative products than me. – M.J.2 Mar 22 '16 at 22:02
• Evanb, have you found some time to take a look at this similar question "mathematica.stackexchange.com/questions/109825/…"? I hope you could extend your answer for that post as well. – M.J.2 May 30 '16 at 18:15
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, 150];

(a.ai ) // N

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

• Your answer is really nice in terms of numerical point of view. Please keep it but what I want is to obtain the expansion theoretically. I mean, we give A and B and obtain the matrix series up to 6 terms theoretically with consideration that matrix multiplication is not commutative. – M.J.2 Mar 21 '16 at 21:19
• @M.J.2 Well, you can. But it is nasty: mInvSeries[{{a, b}, {c, d}}, IdentityMatrix, 3] – Dr. belisarius Mar 21 '16 at 21:29
• @M.J.2 Oh! you probably mean this mInvSeries[IdentityMatrix, {{a, b}, {c, d}}, 3] – Dr. belisarius Mar 21 '16 at 21:30
• I wish to see the result of the code exactly as in the formulation of the question, I mean in terms of $A$ and $B$. Or when the inputs are $I$ and $A$, in terms of $I$ and $A$. – M.J.2 Mar 21 '16 at 21:47
• @M.J.2 Mathematica has no support for that notation. $I$ is the imaginary unit ... – Dr. belisarius Mar 21 '16 at 21:51

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.

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:=ReplaceAll[abSeries[a, ee b, 5].(a + ee b), {a -> PauliMatrix, b -> PauliMatrix}] // Simplify
Out:={{1 + ee^6, 0}, {0, 1 + ee^6}}

• what "general" expression? If Ais the Identity, the result is a polynomial in B... there is no more "generality" than that AFAIK – Dr. belisarius Mar 21 '16 at 21:37
• First of all notice in the question he commented the expression is only valid if $Norm[A^{-1}B]<1$. Second, I meant more general matrices like matrices of infinite dimensionality or for some generators of a Lie algebra etc. – Lior Blech Mar 21 '16 at 21:42
• In general $A$ is not always Identity matrix. I believe this answer is better than the above one but this is not the expansion of the question! Because it does not keep the non-commutative property of matrix multiplications. Please choose $n=5$ and see the differences! – M.J.2 Mar 21 '16 at 21:42
• I totally agree with the sentence "matrices of infinite dimensionality or for some generators of a Lie algebra etc." For example, I want to finally write the Taylor expansion of $(I+B+C)^{-1}$. So, the code should take $(B+C)$ as an input and then apply the Taylor expansion two times. – M.J.2 Mar 21 '16 at 21:45
• @M.J.2 Please try running the following with my code mInvSeries[A, B, 3] Is that what you want? – Dr. belisarius Mar 21 '16 at 21:45

Such formal power series can be computed in NCAlgebra by recursively calculating then evaluating directional derivatives. NCDirectionalD[f,{x,h}] calculates the directional derivative of f with respect to x in the direction h. The following will leverage that to produce a formal power series of f of degree n taken with respect to x around x0:

NCSeries[f_, {x_, x0_, n_}] := Block[{h},
SetNonCommutative[h];
Plus @@ (Table[1/i!, {i, 0, n}]*NestList[NCDirectionalD[#, {x, h}] &, f, n]) /. x -> x0 /. h -> x]


f = inv[a + b];

NCSeries[inv[1 + b], {b, 0, 6}]