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

• Wow, amazing. This is really helpful and complete. Commented Mar 22, 2016 at 10:22
• Can you extend your answer to solve this unanswered question: <mathematica.stackexchange.com/questions/109825/…> Commented Mar 22, 2016 at 10:23
• I can't get Mma to reproduce your analytic results there. Edit: actually, maybe I can. I'll take a look. Commented Mar 22, 2016 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. Commented Mar 22, 2016 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. Commented May 30, 2016 at 18:15

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[f, {b, 0, 4}]
(* a^-1 - a^-1 ** b ** a^-1 + a^-1 ** b ** a^-1 ** b ** a^-1 - a^-1 ** b ** a^-1 ** b ** a^-1 ** b ** a^-1 + a^-1 ** b ** a^-1 ** b ** a^-1 ** b ** a^-1 ** b ** a^-1 *)


Your other cases are calculated similarly:

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


A multivariate version can be computed in the same way.

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}} *)

• 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. Commented Mar 21, 2016 at 21:19
• @M.J.2 Well, you can. But it is nasty: mInvSeries[{{a, b}, {c, d}}, IdentityMatrix[2], 3] Commented Mar 21, 2016 at 21:29
• @M.J.2 Oh! you probably mean this mInvSeries[IdentityMatrix[2], {{a, b}, {c, d}}, 3] Commented Mar 21, 2016 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$. Commented Mar 21, 2016 at 21:47
• @M.J.2 Mathematica has no support for that notation. $I$ is the imaginary unit ... Commented Mar 21, 2016 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[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}}

• what "general" expression? If Ais the Identity, the result is a polynomial in B... there is no more "generality" than that AFAIK Commented Mar 21, 2016 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. Commented Mar 21, 2016 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! Commented Mar 21, 2016 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. Commented Mar 21, 2016 at 21:45
• @M.J.2 Please try running the following with my code mInvSeries[A, B, 3] Is that what you want? Commented Mar 21, 2016 at 21:45

This answer is an attempt to provide a solution for a hypothetical user that does not know the general formula beforehand and is unaware of any packages.

The answer below considers that the user is aware of the following:

• If $$f$$ is an expandable function that can be considered on the real line (for example $$f(x)=(1+x)^{-1}$$ or $$f(x)=\exp(x)$$), then, for a matrix variable $$A$$ and a real number valued variable $$a$$, $$f(A)$$=Replace[$$f(a),a \rightarrow A$$]. This point allows us to use Series for real valued variables.

• For $$A$$ and $$B$$ two matrices $$\left(A.B\right)^{-1}=B^{-1}.A^{-1}$$. This point is used to place the initial expression into the form $$f(x)$$ where $$f(x)=(1+x)^{-1}$$ and $$x=A^{-1}.B$$

Then if we define an identity tensor for the dot product:

Dot[s___, id, g___] ^:= Dot[s, g]


The formula may be obtained with the following code that attempts to do simple replacements sequentially (var below is a formal wrapper that formally replaces a matrix with a real valued variable) :

Fold[ReplaceAll,
Inverse[A + B] ,
{
Inverse[w_] :> Inverse[A . TensorExpand[Inverse@A . w]] ,
Inverse[a_ . b_] :> Inverse[b] . Inverse[a],
MatrixPower[_, 0] -> id,
Inverse[id + s_] :> ReplaceAll[

Normal@Series[(1 + epsilon*var[s])^(-1),
{epsilon, 0, 3}] /. epsilon -> 1

, {
var[a_]^m_:>Dot@@ConstantArray[a,m]
, var->Identity
, 1 -> id
}
]
}
]//TensorExpand


out: (* -MatrixPower[A, -1] . B . MatrixPower[A, -1] + MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] - MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] . B . MatrixPower[A, -1] + MatrixPower[A, -1] *)