4
$\begingroup$

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
Improve this question
$\endgroup$
0

1 Answer 1

Reset to default
9
$\begingroup$

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[{ND$M = Length[ND$W],ND$P = Floor[Length[ND$W]/2],ND$Q = Floor[(Length[ND$W]-1)/2],ND$K},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[{ND$P = Flatten[Map[Permutations,IntegerPartitions[ND$N]],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[
    {
        ND$A1,ND$A2,ND$A3,ND$A4,ND$A5,ND$A6,ND$A7,ND$A8,ND$A9
    },
    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[#]],{ND$A3+1}])&,ND$A7],1] ;
    ND$A2 = Map[ND$TTRANSFORM,Join[ND$ALPH,Table[ND$ZASSENHAUS$UNFOLDED[ND$j,ND$ALPH],{ND$j,2,ND$N-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,{ND$A2,#}]]&,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 {<CHAR1>,<CHAR2>} *)
    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 <ND$N>  *)
    (* False -- ONLY <ND$N> TERM *)
] := Block[
    {
        x,
        y,
        ND$I$1,
        ND$I$2
    },
        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 <ND$N>  *)
    (* False -- ONLY <ND$N> TERM *)
] := Block[
    {
        x,
        y,
        ND$I$1,
        ND$I$2
    },
        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]]]
        ]
] ;
(* ################################################################################ *)
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ Thanks! background for this problem is given here $\endgroup$
    – Yaroslav Bulatov
    Commented Apr 21, 2023 at 5:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.