How to force Series[] to compute expansions by considering non commutative multiplication?

Question

I wish to compute the Taylor series expansion of the following iteration method $x_{k+1}=x_k-f'(x_k)^{-1}f(x_k)$ up to four terms of error ($e_k=x_k-\alpha$). When this is a scalar iteration, I very simply write the following

ClearAll["Global`*"]
f[e_] := df*(e^1 + c2 e^2 + c3 e^3 + c4 e^4 + c5 e^5);
fe = f[e];
f1e = f'[e];
Series[f1e^-1, {e, 0, 4}]*df // Simplify
u = e - Series[fe/f1e, {e, 0, 4}] // Simplify

and obtain correct results $$f(x_k)^{-1}=[1-2 \text{c2} e_k+ \left(4 \text{c2}^2-3 \text{c3}\right)e_k^2$$ $$-4 (2 c2^3 - 3 c2 c3 + c4) e_k^3 + (16 c2^4 - 36 c2^2 c3 + 9 c3^2 + 16 c2 c4 - 5 c5) e_k^4]f'({\alpha})^{-1}.$$ But, how to do this for the multi-dimensional case. That is, when the coefficients $c2,c3,c4$, and even $df=f'({\alpha})^{-1}$ are all matrices (note that e.g. $c2=\frac{1}{2!}f'(\alpha)^{-1}f^{(2)}(\alpha)$). By hand, I obtain the following correct results: $$ f(x_k)^{-1}=\left(I-2c2e_k+(4c2^2-3c3)e_k^2+(6c3c2+6c2c3-8c2^3-4c4)e_k^3\\ + (8c4c2+9c3^2+8c2c4-5c5-12c3c2^2-12c2c3c2-12c2^2c3+16c2^4)e_k^4\right)f'(\alpha)^{-1} $$ and \begin{equation} e_{k+1}=x_{k+1}-\alpha=-c2e_k^2+(2c2^2-2c3)e_k^3+(-4c2^3+4c2c3+3c3c2-3c4)e_k^4+\mathcal{O}(e_k^5). \end{equation}

So, my question is how to handle this non-commutative multiplication inside Series[] (e.g. $c2c3$ is not equal to $c3c2$ in this case)? I also tried to apply $**$, but I failed. Any suggestions or tricks to solve this problem will be appreciated fully.

I also saw the following posts in MathematicaStackExchange, but they were not useful for this problem, How to make noncommutative multiplication agree with commutative multiplication and Noncommutative multiply- expand expression.

Answer 1

One idea is to use the function MatrixD from my answer to question det simplification that takes scalar derivatives of functions of matrices:

MatrixD[expr_, x__] := With[
    {old = OptionValue[SystemOptions[], "DifferentiationOptions"->"ExcludedFunctions"]},

    Internal`WithLocalSettings[
        SetSystemOptions["DifferentiationOptions"->"ExcludedFunctions"->Join[old, {Det, Inverse, Tr}]];

        Unprotect[D];
        (* handle list derivatives *)
        D[h:((Det|Tr|Inverse)[m_]), {z_, n_Integer}] := Nest[D[#, Replace[z, _List :> {z}]]&, h, n];
        D[h:((Det|Tr|Inverse)[m_]), {z_List}] := D[h, #]& /@ z;
        D[h:((Det|Tr|Inverse)[m_]), z_, y___] := D[D[h, z], y];

        (* define derivatives for Det, Tr, and Inverse *)
        D[Det[m_], z:Except[_List]] := Det[m] Tr[Inverse[m] . D[m,z]];
        D[Tr[m_], z:Except[_List]] := Tr[D[m,z]];
        D[Inverse[m_], z:Except[_List]] := -Inverse[m] . D[m, z] . Inverse[m],

        D[expr, x],

        SetSystemOptions["DifferentiationOptions"->"ExcludedFunctions"->old];
        Clear[D];
        Protect[D]
    ]
]

Here are your expressions:

f[e_] := df*(IdentityMatrix[d] e^1 + c2 e^2 + c3 e^3 + c4 e^4 + c5 e^5);
fe = f[e];
f1e = f'[e];

I added the IdentityMatrix and assumed that df is just a scalar. Then:

MatrixD[Inverse[f1e], e] /. e->0

-Inverse[df IdentityMatrix[d]].(2 c2 df).Inverse[df IdentityMatrix[d]]

In order to make progress, we need to include assumptions about the coefficients:

$Assumptions = (c2 | c3 | c4 | c5) ∈ Matrices[{d, d}] && df ∈ Complexes;

The first few derivatives of Inverse[f1e] tensor expanded using the above assumptions:

TensorExpand[MatrixD[df Inverse[f1e], e] /. e->0]
TensorExpand[MatrixD[df Inverse[f1e], e, e] /. e->0]
TensorExpand[MatrixD[df Inverse[f1e], e, e, e] /. e->0]

-2 c2

-6 c3 + 8 MatrixPower[c2, 2]

-24 c4 + 36 c2.c3 + 36 c3.c2 - 48 MatrixPower[c2, 3]

So, the series expansion is:

TensorExpand[
    (df Inverse[f1e] /. e->0) +
    (MatrixD[df Inverse[f1e], e] /. e->0) e +
    (MatrixD[df Inverse[f1e], e, e] /. e->0) e^2/2! + 
    (MatrixD[df Inverse[f1e], e, e, e] /. e->0) e^3/3! +
    (MatrixD[df Inverse[f1e], e, e, e, e] /. e->0) e^4/4!
] + O[e]^5 //TeXForm

$\operatorname{IdentityMatrix}[d]-2 \operatorname{c2} e+e^2 (4 \operatorname{MatrixPower}[\operatorname{c2},2]-3 \operatorname{c3})+e^3 (-8 \operatorname{MatrixPower}[\operatorname{c2},3]+6 \operatorname{c2}.\operatorname{c3}+6 \operatorname{c3}.\operatorname{c2}-4 \operatorname{c4})+e^4 (-12 \operatorname{c3}.\operatorname{MatrixPower}[\operatorname{c2},2]-12 \operatorname{MatrixPower}[\operatorname{c2},2].\operatorname{c3}+16 \operatorname{MatrixPower}[\operatorname{c2},4]+9 \operatorname{MatrixPower}[\operatorname{c3},2]-12 \operatorname{c2}.\operatorname{c3}.\operatorname{c2}+8 \operatorname{c2}.\operatorname{c4}+8 \operatorname{c4}.\operatorname{c2}-5 \operatorname{c5})+O\left(e^5\right)$

which is in basic agreement with your inverse computation. A similar computation can be done for Inverse[f1e] . fe. However, I will instead present an approach teaching Series to use MatrixD for matrix functions. The internal function used by Series is System`Private`InternalSeries, so we teach this internal function about Dot and Inverse:

Unprotect[System`Private`InternalSeries];
System`Private`InternalSeries[a_Inverse|a_Dot, {e_,e0_,n_}] := SeriesData[
    e,
    e0,
    TensorReduce[TensorExpand[NestList[MatrixD[#,e]&,a,n]/.e->e0]]/Range[0,n]!,
    0,
    n+1,
    1
]

Now, we can use Series directly:

Series[df Inverse[f1e], {e, 0, 4}] //TeXForm

$\operatorname{IdentityMatrix}[d]-2 \operatorname{c2} e+e^2 (4 \operatorname{MatrixPower}[\operatorname{c2},2]-3 \operatorname{c3})+e^3 (-8 \operatorname{MatrixPower}[\operatorname{c2},3]+6 \operatorname{c2}.\operatorname{c3}+6 \operatorname{c3}.\operatorname{c2}-4 \operatorname{c4})+e^4 (-12 \operatorname{c3}.\operatorname{MatrixPower}[\operatorname{c2},2]-12 \operatorname{MatrixPower}[\operatorname{c2},2].\operatorname{c3}+16 \operatorname{MatrixPower}[\operatorname{c2},4]+9 \operatorname{MatrixPower}[\operatorname{c3},2]-12 \operatorname{c2}.\operatorname{c3}.\operatorname{c2}+8 \operatorname{c2}.\operatorname{c4}+8 \operatorname{c4}.\operatorname{c2}-5 \operatorname{c5})+O\left(e^5\right)$

in agreement with the previous result. Finally:

Series[Inverse[f1e] . fe, {e, 0, 4}] - e IdentityMatrix[d] //TeXForm

$-\operatorname{c2} e^2-2 e^3 (\operatorname{c3}-\operatorname{MatrixPower}[\operatorname{c2},2])+e^4 (-4 \operatorname{MatrixPower}[\operatorname{c2},3]+4 \operatorname{c2}.\operatorname{c3}+3 \operatorname{c3}.\operatorname{c2}-3 \operatorname{c4})+O\left(e^5\right)$

in agreement with your result.

Answer 2

This can be calculated by using NCAlgebra and this answer another question. The function

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]

will calculate a formal power series of rational noncommutative functions. Your example is mixed, in the sense that e seems to commute but the coefficients do not. I propose that you lift e into a noncommutative variable ee

SNC[c1, c2, c3, c4, c5, ee] 
f[e_] := df*(e^1 + c2 ** e^2 + c3 ** e^3 + c4 ** e^4 + c5 ** e^5);
fe = f[ee]
f1e = f'[ee]   (* this works as expected in this simple case only! *)

use NCSeries to calculate the formal power series

series = NCSeries[ee - inv[f1e] ** fe, {ee, 0, 4}]

which you can then project back by replacing ee by the commutative e:

SetCommutative[e];
Collect[NCExpand[series /. ee -> e], e]

which will give the result you're looking for.