Maximally expand tensor series?

Question

I am trying to obtain a power series expansion in some real parameter, but all the terms are arbitrary products of tensors. I.e. I want to expand an expression containing sums/products/powers of this parameter (say a) and arbitrary sums/products/powers of tensors as much as possible. For example

a((MatrixPower[X.Y + Y.X, 2] + a X.Y.Y.X))

should give

a(X.Y.X.Y + X.Y.Y.X + Y.X.X.Y + Y.X.Y.X) + a^2 X.Y.Y.X

After trying combinations of TensorExpand, TensorReduce, Distribute, Collect, Refine, Replace, .. I'm lost. For example I get a term like

$Assumptions = a \[Element] Reals && A1 \[Element] Matrices[{4, 4}] && A2 \[Element] Matrices[{4, 4}] && A3 \[Element] Matrices[{4, 4}] && A4 \[Element] Matrices[{4, 4}];

a^3 (A1.A2/2 + A1.A3/2 - A2.A1/2 + A2.A3/2 - A3.A1/2 - A3.A2/2
+ 1/12 A3.(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2
- 2 A1.A2.A1 + A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)).
(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2 - 2 A1.A2.A1
+ A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)))

after Collect, which is obviously not only a^3, but also a^4, a^5, a^6. I cannot get this into a form like

a^3(...) + a^4(...) + a^5(...) + a^6(...)

TensorExpand / TensorReduce do not do what I want because they write MatrixProduct[.., ..] whenever possible which still contains factors of a inside and then cannot be properly used for Collect. Distribute also works only in the simplest cases.

Is it somehow possible to do what I want?

Answer 1

You can use my code to implement series expansions of tensor objects in the question How to force Series[] to compute expansions by considering non commutative multiplication?, repeated at the end of this question. Your expression:

expr = a^3 (A1.A2/2 + A1.A3/2 - A2.A1/2 + A2.A3/2 - A3.A1/2 - A3.A2/2
+ 1/12 A3.(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2
- 2 A1.A2.A1 + A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)).
(a A1 + a A2 + 1/2 a^2 (A1.A2 - A2.A1) + 1/12 a^3 (A1.A1.A2 - 2 A1.A2.A1
+ A2.A1.A1) + 1/12 a^3 (A1.A2.A2 - 2 A2.A1.A2 + A2.A2.A1)));

Using Series:

Series[expr, {a, 0, 6}] //TeXForm

$\frac{1}{2} a^3 (\operatorname{A1}.\operatorname{A2}-\operatorname{A2}.\operatorname{A1}+\operatorname{A1}. \operatorname{A3}-\operatorname{A3}.\operatorname{A1}+\operatorname{A2}.\operatorname{A3}-\operatorname{A3}.\operatorname{A2})+\frac{1}{12} a^5 (\operatorname{A3}.\operatorname{A1}.\operatorname{A2}+\operatorname{A3}.\operatorname{A2}. \operatorname{A1}+\operatorname{A3}.\operatorname{A1}.\operatorname{A1}+\operatorname{A3}.\operatorname{A2}.\operatorname{A2})+\frac{1}{24} a^6 (\operatorname{A3}.\operatorname{A1}.\operatorname{A1}.\operatorname{A2}+\operatorname{A3}. \operatorname{A1}.\operatorname{A2}.\operatorname{A2}-\operatorname{A3}.\operatorname{A2}.\operatorname{A1}.\operatorname{A1}-\operatorname{A3}.\operatorname{A2}.\operatorname{A2}.\operatorname{A1})+O\left(a^7\right)$

The code:

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

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
]