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

Question

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.

Answer 1

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]

Your example is reproduced by:

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.

Answer 2

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]

Answer 3

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

Answer 4

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

Answer 5

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