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

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 $$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).$$

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.

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"]},

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


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 SystemPrivateInternalSeries, so we teach this internal function about Dot and Inverse: Unprotect[SystemPrivateInternalSeries]; SystemPrivateInternalSeries[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)\$

• Wow, a very useful response. The step-by-step way to get the multi-dimensional Taylor expansion is fruitful.
– Faz
Commented Nov 16, 2018 at 17:48

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.