# Zassenhaus formula in Mathematica

I'm looking for Mathematica implementation of Zassenhaus formula -- given two matrices $$A,B$$, truncate the expansion of $$\exp(A+B)$$ from this paper:

$$e^{t(A+B)}= e^{tA}e^{tB}\prod_{n=2}^\infty e^{t^n Z_n(A,B)}$$

What's the most elegant way to implement this in Mathematica?

Below is an evaluation of "zero-th order" truncation of this expansion for $$A,B$$ appearing in my application (diagonal + rank1), trying to see how many terms I need for a reasonable approximation error.

d = 100;
h = Normalize[#^-1.1 & /@ Range[d], Tr];

a0 = 2 h (1 - h);
b0 = h;

ones = ConstantArray[1., d];
t = -d;
A = DiagonalMatrix[a0];
e[mat_] := MatrixExp[mat];
norm[mat_] := Total[mat, 2];
B = -{b0}\[Transpose] . {b0};
{{norm@e[t (A + B)], "exact"},
{norm[e[t A] . e[t B]], "order-0"}} // TableForm


Here is my old implementation based on M.WEYRAUCH, D.SCHOLZ, COMPUTER PHYSICS COMMUNICATIONS, 180, (2009), 1558-1565

Returns 'unfolded' or 'folded' (in terms of commutators):

ND$$ZASSENHAUS$$UNFOLDED[2, {"X", "Y"}]
ND$$ZASSENHAUS$$UNFOLDED[3, {"X", "Y"}]
ND$$ZASSENHAUS$$UNFOLDED[4, {"X", "Y"}]
ND$ZAS[4, X, Y, F] (* -((XY)/2)+(YX)/2 *) (* (XXY)/6-(XYX)/3-(XYY)/3+(YXX)/6+(2 YXY)/3-(YYX)/3 *) (* -((XXXY)/24)+(XXYX)/8+(XXYY)/8-(XYXX)/8-(XYXY)/4-(XYYY)/8+(YXXX)/24+(YXYX)/4+(3 YXYY)/8-(YYXX)/8-(3 YYXY)/8+(YYYX)/8 *) (* {X,Y,-(1/2) F[X,Y],-(1/6) F[F[X,Y],X]-1/3 F[F[X,Y],Y],-(1/24) F[F[F[X,Y],X],X]-1/8 F[F[F[X,Y],X],Y]-1/8 F[F[F[X,Y],Y],Y]} *)  Example for matrices: (* compute terms *) ClearAll[X, Y, F] ; order = 8 ; terms = ND$ZAS[order, X, Y, F] // Rest // Rest ;

(* set random input & compute exact output *)
SeedRandom[1] ;
t = 0.1 ;
a = RandomReal[{0.0, 1.0}, {2, 2}] ;
b = RandomReal[{0.0, 1.0}, {2, 2}] ;
c = MatrixExp[t * (a + b)]  ;

(* O(t^2) *)
result = MatrixExp[t*a] . MatrixExp[t * b] ;
Norm[result - c]
(* 0.001777447231126185 *)

(* higher order *)
result = Table[result = result . MatrixExp[terms[[i]] /. {X ->t*a, Y ->t* b, F -> Function[{x, y}, x . y -y . x]}] , {i, 1, Length[terms]}] ;
Map[Norm[# - c] &, result]
(* {0.00005544291090346927,2.773628721869079*^-6,7.421340761934133*^-8,1.375336002537211*^-9,3.593262198290658*^-11,1.2194115693175485*^-12,2.3818383250085455*^-14} *)


Main code:

(* ################################################################################ *)
(* BCH AND ZASSENHAUS (M.WEYRAUCH, D.SCHOLZ, COMPUTER PHYSICS COMMUNICATIONS, 180, (2009), 1558-1565) *)
(* ND$$BCH$$UNFOLDED[TERM,{LEFT,RIGHT}] *)
(* ND$$ZASSENHAUS$$UNFOLDED[TERM,{LEFT,RIGHT}] *)
(* ND$$CONVERT$$UNFOLDED[TYPE,TERM,{LEFT,RIGHT},HEAD] *)
(* ND$$BCH[TERM,LEFT,RIGTH,HEAD,MODE] *) (* ND$$ZAS[TERM,LEFT,RIGTH,HEAD] *)
(* ND$$TRUNCATED$$EXP[TERM,LEFT,RIGHT,HEAD,MODE] *)
(* ################################################################################ *)

(* ################################################################################ *)
(* AUXILIARY FUNCTION FOR BCH AND ZASSENHAUS *)
ClearAll[ND$$WORDS] ; ND$$WORDS[ND$$P_List, ND$$ALPH_List] := StringJoin[Apply[ConstantArray,Partition[Riffle[ND$$P,ND$$ALPH,List[1,Times[2,Length[ND$P]],2]],2],List[1]]] ; (* ################################################################################ *) (* ################################################################################ *) (* UNFOLDED BCH (GOLDBERG'S METHOD) *) ClearAll[ND$$G] ; ClearAll[ND$$C] ; ClearAll[ND$$BCH$$UNFOLDED] ; ND$$G[1] = 1 ; ND$$G[ND$$S_] := ND$$G[ND$$S] = Expand[1/ND$$S*D[ND$$T*(ND$$T-1)*ND$$G[ND$$S-1],ND$$T]]; ND$$C[ND$$W_] := ND$$C[ND$$W] = Block[{NDM = Length[NDW],NDP = Floor[Length[NDW]/2],NDQ = Floor[(Length[NDW]-1)/2],NDK},Integrate[ND$$T^ND$$P*(ND$$T-1)^ND$$Q*Product[ND$$G[ND$$W[[ND$$K]]],{ND$$K,ND$$M}],{ND$$T,0,1}]] ; ND$$BCH$$UNFOLDED::usage="ND$$BCH$$UNFOLDED[TER,{LEF,RIG}] " ; ND$$BCH$$UNFOLDED[ND$$N_,ND$$ALPH_] := Block[{NDP = Flatten[Map[Permutations,IntegerPartitions[NDN]],1]},Expand[Total[Map[ND$$C[Sort[#]](ND$$WORDS[#,ND$$ALPH]-(-1)^ND$$N ND$$WORDS[#,Reverse[ND$$ALPH]])&,ND$$P]]]] ; (* ################################################################################ *) (* ################################################################################ *) (* ZASSENHAUS, (COMPARISON METHOD) *) ClearAll[ND$$TPOWER] ; ClearAll[ND$$TPRODUCT] ; ClearAll[ND$$TSCALAR] ; ClearAll[ND$$TTRANSFORM] ; ClearAll[ND$$ZASSENHAUS$$UNFOLDED] ; ND$$TPOWER[ND$$X_,0] := List[List[1,""]] ; ND$$TPOWER[ND$$X_,1] := ND$$X ; ND$$TPOWER[ND$$X_,ND$$N_/;ND$$N>1] := Nest[Map[List[Times[Extract[#,{1,1}],Extract[#,{2,1}]],StringJoin[Extract[#,{1,2}],Extract[#,{2,2}]]]&,Flatten[Outer[List,#,ND$$X,1],1]]&,ND$$X,ND$$N-1] ; ND$$TPRODUCT[ND$$X_] := Fold[Map[List[Times[Extract[#,{1,1}],Extract[#,{2,1}]],StringJoin[Extract[#,{1,2}],Extract[#,{2,2}]]]&,Flatten[Outer[List,#1,#2,1],1]]&,First[ND$$X],Rest[ND$$X]] ; ND$$TSCALAR[ND$$C_,ND$$X_] := Map[List[ND$$C First[#],Last[#]]&,ND$$X] ; ND$$TTRANSFORM[ND$$X_List] := Total[Apply[Times,Flatten[ND$$X,1],List[1]]] ; ND$$TTRANSFORM[ND$$X_String] := List[List[1,ND$$X]] ; ND$$TTRANSFORM[ND$$X_Plus] := Apply[List,Apply[List,ND$$X],List[1]] ; ND$$ZASSENHAUS$$UNFOLDED::usage="ND$$ZASSENHAUS$$UNFOLDED[TER,{LEF,RIGH}] " ; ND$$ZASSENHAUS$$UNFOLDED[2,ND$$ALPH_List] := ND$$ZASSENHAUS$$UNFOLDED[2,ND$$ALPH] = - StringJoin[ND$$ALPH]/2+StringJoin[Reverse[ND$$ALPH]]/2 ; ND$$ZASSENHAUS$$UNFOLDED[ND$$N_/;ND$$N>2,ND$$ALPH_List] := ND$$ZASSENHAUS$$UNFOLDED[ND$$N,ND$$ALPH] = Block[ { NDA1,NDA2,NDA3,NDA4,NDA5,NDA6,NDA7,NDA8,NDA9 }, ND$$A5 = IntegerPartitions[ND$$N] ; ND$$A7 = Rest[Map[BinCounts[#,List[1,ND$$N,1]]&,ND$$A5]] ; ND$$A8 = Flatten[Map[(ND$$A3=ND$$A6=First[#];Table[Join[List[ND$$A4=ND$$A6--,ND$$A3-ND$$A4],Rest[#]],{NDA3+1}])&,ND$$A7],1] ; ND$$A2 = Map[ND$$TTRANSFORM,Join[ND$$ALPH,Table[ND$$ZASSENHAUS$$UNFOLDED[ND$$j,ND$$ALPH],{NDj,2,NDN-1}]]] ; ND$$A5 = Flatten[Map[Permutations,ND$$A5],1] ; ND$$A9 = Total[Join[Map[ND$$WORDS[#,ND$$ALPH]/Factorial[ND$$N]&,ND$$A5],Map[ND$$WORDS[#,Reverse[ND$$ALPH]]/Factorial[ND$$N]&,ND$$A5]]] ; ND$$A1 = MapThread[ND$$TSCALAR,{Map[Apply[Times,(1/Factorial[#])]&,ND$$A8],Map[ND$$TPRODUCT[MapThread[ND$$TPOWER,{NDA2,#}]]&,ND$$A8]}] ; ND$$A9 - ND$$TTRANSFORM[ND$A1]
] ;
(* ################################################################################ *)

(* ################################################################################ *)
(* POLYNOMIAL CONVERTER FOR BCH AND ZASSENHAUS, (I.M., 2014-2015) *)
ClearAll[ND$$CONVERT$$UNFOLDED] ;
ND$$CONVERT$$UNFOLDED::usage="ND$$CONVERT$$UNFOLDED[TYP,TER,{LEF,RIG},HEA] " ;
ND$$CONVERT$$UNFOLDED[
ND$$TYPE_, (* ND$$BCH$$UNFOLDED OR ND$$ZASSENHAUS$$UNFOLDED *) ND$$N_,        (* NUMBER OF LETTERS (INTEGER > 1) *)
ND$$LETTERS_, (* CHARACTERS FOR ALGEBRA ELEMENTS {,} *) ND$$HEAD_      (* COMMUTATOR HEAD *)
] := Block[
{
ND$$POLY, ND$$LISTA,
ND$$LISTB, ND$$A,
ND$$B, ND$$TMP,
ND$$SENTENCES, ND$$RESULT
},
ND$$HEAD[ND$$a_,ND$$b_,ND$$c__] := ND$$HEAD[ND$$HEAD[ND$$a,ND$$b],ND$$c] ; ND$$POLY = Apply[List,ND$$TYPE[ND$$N,ND$$LETTERS]] ; ND$$LISTA = Part[ND$$POLY,All,1] ; ND$$LISTB = Part[ND$$POLY,All,2] ; ND$$TMP = StringCases[ND$$LISTB,StringExpression[StartOfString,StringJoin[ND$$LETTERS],___]] ;
ND$$SENTENCES = Flatten[ND$$TMP] ;
ND$$A = Part[ND$$LISTA,Take[Flatten[Position[ND$$TMP,_String]],{1,-1,2}]] ; ND$$B = StringCount[ND$$SENTENCES,First[ND$$LETTERS]] ;
ND$$RESULT = Total[ND$$A/ND$$B Map[Apply[ND$$HEAD,ToExpression[StringCases[#, Repeated[_,1]]],List[0]]&,ND$$SENTENCES]] ; Clear[ND$$HEAD] ;
ND$$RESULT ] /; ND$$N > 1 ;
(* ################################################################################ *)

(* ################################################################################ *)
(* CONVERTED BCH *)
ClearAll[ND$$BCH] ; ND$$BCH::usage="ND$$BCH[TER,LEF,RIG,HEA,MOD] " ; ND$$BCH[
ND$$N_, (* -- NUMBER OF TERM *) ND$$A_,         (* -- LEFT *)
ND$$B_, (* -- RIGHT *) ND$$HEAD_,      (* -- COMMUTATOR HEAD *)
ND$$MODE_:True (* -- OUTPUT FORMAT *) (* True (DEFAULT) -- LIST OF TERMS UP TO N>  *)
(* False -- ONLY <ND$$N> TERM *) ] := Block[ { x, y, NDI1, NDI2 }, ND$$I$$1 = If[ND$$MODE,2,ND$$N] ; ND$$I$$2 = ND$$N ;
ReplaceAll[
If[
ND$$MODE, Join[List[ND$$A+ND$$B],#], First[#] ] & @ Map[ND$$CONVERT$$UNFOLDED[ND$$BCH$$UNFOLDED,#,{"x","y"},ND$$HEAD]&,Range[ND$$I$$1,ND$$I$$2]],
Thread[Rule[List[x,y],List[ND$$A,ND$$B]]]
]
] ;
(* ################################################################################ *)

(* ################################################################################ *)
(* CONVERTED ZAS *)
ClearAll[ND$$ZAS] ; ND$$ZAS::usage="ND$$ZAS[TER,LEF,RIG,HEA,MOD] " ; ND$$ZAS[
ND$$N_, (* -- TERM NUMBER *) ND$$A_,         (* -- LEFT *)
ND$$B_, (* -- RIGHT *) ND$$HEAD_,      (* -- COMMUTATOR HEAD *)
ND$$MODE_:True (* -- OUTPUT FORMAT *) (* True (DEFAULT) -- LIST OF TERMS UP TO N>  *)
(* False -- ONLY <ND$$N> TERM *) ] := Block[ { x, y, NDI1, NDI2 }, ND$$I$$1 = If[ND$$MODE,2,ND$$N] ; ND$$I$$2 = ND$$N ;
ReplaceAll[
If[
ND$$MODE, Join[List[ND$$A],List[ND$$B],#], First[#] ] & @ Map[ND$$CONVERT$$UNFOLDED[ND$$ZASSENHAUS$$UNFOLDED,#,{"x","y"},ND$$HEAD]&,Range[ND$$I$$1,ND$$I$$2]],
Thread[Rule[List[x,y],List[ND$$A,ND$$B]]]
]
] ;
(* ################################################################################ *)

• Thanks! background for this problem is given here Commented Apr 21, 2023 at 5:52