11
$\begingroup$

I am working with noncommuting objects and basically I am using Mathematica to sort large expressions by normal ordering and simplifying them. To this end, I first create an object called Boson (These satisfy $[a,a^\dagger]=1$):

Clear[Boson, BosonC, BosonA]
Boson /: MakeBoxes[Boson[cr : (True | False), sym_], fmt_] := 
With[{sbox = If[StringQ[sym], sym, ToBoxes[sym]]}, 
With[{abox = 
 If[cr, SuperscriptBox[#, FormBox["\[Dagger]", Bold]], #] &@sbox},
InterpretationBox[abox, Boson[cr, sym]]]]
BosonA[sym_: String "a"] := Boson[False, sym]
BosonC[sym_: String "a"] := Boson[True, sym]

Next I alias the noncommutative product with CenterDot as follows:

Unprotect[NonCommutativeMultiply];
Clear[NonCommutativeMultiply, CenterDot];
CenterDot[a__] := NonCommutativeMultiply[a];
NonCommutativeMultiply /: 
MakeBoxes[NonCommutativeMultiply[a__], fmt_] := 
With[{cbox = ToBoxes[HoldForm[CenterDot[a]]]}, 
InterpretationBox[cbox, NonCommutativeMultiply[a]]]
Protect[NonCommutativeMultiply];
Clear[CRule]
CRule = {NonCommutativeMultiply[a_] :> a};

Then I define the Function clean which does the normal ordering. It is defined in terms of expand. The main definition for expand is

ClearAll@expand
SetAttributes[expand, HoldAll]

Unevaluated[
  expand[expr_] :=
   Block[
    {NonCommutativeMultiply (*or times*)},

    expr //. (*ReplaceRepeated instead of ReplaceAll*)
     {times[left___, cnum_ /; FreeQ[cnum, (_Boson)], 
        right___] :> cnum*times[left, right], 
      times[left___, 
        cnum_ /; (! FreeQ[cnum, Times[n___?NumericQ, ___Boson]]), 
        right___] :> 
       Times @@ Apply[Power, Drop[FactorList[cnum], -1], 2]*
        times[left, First[Last[FactorList[cnum]]], right], 
      times[left___, Boson[False, s_], Boson[True, s_], right___] :> 
       times[left, right] + 
        times[left, Boson[True, s], Boson[False, s], right], 
      times[left___, fst : Boson[_, s_], sec : Boson[_, t_], 
        right___] :> 
       times[left, sec, fst, right] /; 
        FreeQ[Ordering[{s, t}], {1, 2}], times[b_Boson] :> b,
      times[] -> 1
      }
    ]
  ] /. {HoldPattern[times] -> NonCommutativeMultiply }

More definitions for expand are

expand[Alternatives[NonCommutativeMultiply, CenterDot][
a1_, (a2_ + a3_)]] := expand[a1 ** a2] + expand[a1 ** a3]
expand[Alternatives[NonCommutativeMultiply, CenterDot][(a1_ + a2_), 
a3_]] := expand[a1 ** a3] + expand[a2 ** a3]

Finally, the definition of clean is

Clear[clean]
clean = Simplify[FixedPoint[expand, Distribute //@ #] //. CRule] &;

This has been working very well for computing quartic terms but now the problem is that I am computing higher order terms which can involve products of up to 16 oscillators. When I use the clean[] function on those terms, the computation goes okay for a while but then the kernel starts using pretty much all the memory and it just stops doing anything else and I have quit the kernel. I am assuming that the Module[] function is not deleting the temporary variables and is leaking memory. I tried setting $HistoryLength=0 but that didn't help either. What else can I do prevent this?

Update: I think the problem is in the implementation of FixedPoint in clean. If I manually apply expand and Distribute and then clear the cache and keep doing it till the computation doesn't take long (I take that as an indication that not much is changing), then I can finish the computation. This is a "dirty" fix but it would be nice to know if there is a way I can implement clearing the cache in clean along with FixedPoint.

Update: The data in one of the computations became so large that I kept running out of memory for everything and I had to keep dividing the expression into smaller and smaller parts to the point where it became hard to manage. Could be just that expressions are too big and I just don't have enough memory? Also edited the last block of code to remove a stray piece of code.

Update: Here is the full working notebook in pastebin.

$\endgroup$
11
  • $\begingroup$ Would there be the same problem if Module to be changed to Block? $\endgroup$
    – Andrew
    Dec 25, 2013 at 11:30
  • $\begingroup$ Yes I already tried that. I am still not very sure of the exact differences but when I use Module, the memory seems to be around 300-400 MB for a minute and then quickly shoots up to 2.5GB or so and stalls. With Block it stays around 1.3GB and then shoots up and stalls. $\endgroup$
    – Bilentor
    Dec 25, 2013 at 14:15
  • $\begingroup$ The memory use behavior sounds very similar to what I stumbled on in a completely different context. Support suggested I try this support.wolfram.com/kb/11466 and that solved my problem. Can you try that and see whether it makes a difference? It doesn't require code changes and it is easy to undo the changes after the experiment. $\endgroup$ Jan 6, 2014 at 4:15
  • $\begingroup$ @BillSimpson: I tried that and it didn't really help. I have seen this kind of problem on SE in some other questions without any clear solutions. $\endgroup$
    – Bilentor
    Jan 6, 2014 at 11:23
  • $\begingroup$ Am I correct in thinking that you replace NonCommutativeMultiply by times in order to avoid issues with the attribute Flat? To avoid this, you can simply use expand[expr_] := Block[{NonCommutativeMultiply}, replacementCode], which temporarily disregards any attributes (as any other definitions) of NonCommutativeMultiply. $\endgroup$ Jan 6, 2014 at 13:42

2 Answers 2

9
+150
$\begingroup$

This is a very confusing mess of code. It seems to conmingle formatting with functionality, and that certainly does not make it straightforward to see what is wanted, let alone provide code. So what I show below could well be off the mark.

I basically copied and modified code for working with commutators avaliable here

(Note to any colleagues of mine who might be reading this: Searching the site at library.wolfram.com remains an unfriendly task, many years after such things have become a "solved problem" in the "real world". Instead of navigating to the cite and trying my bad luck with a search, I should have just used Google. Unhappy that I'm calling you search engineer folks out on this? Feel free to drop by my office and punch me before I leave this evening.)

(Note to other readers: It's -10 F in Champaign today, or -23 C, and the place is empty. So I'm actually pretty safe. Also my office was move last August and nobody knows where to find me anymore. Or maybe nobody cares. Safe either way.)

So back to the code. I'm not sure all these rules are needed and some might be amiss or missing.

canonicalizeRules = {NonCommutativeMultiply[x___, Boson[False, y_], 
     Boson[True, y_], z___] :> 
    NonCommutativeMultiply[x, Boson[True, y], Boson[False, y], z] + 
     NonCommutativeMultiply[x, z],
   NonCommutativeMultiply[w___, Boson[tf1_, x_], Boson[tf2_, y_], 
      z___] /; (tf1 === tf2 || x =!= y) && Not[OrderedQ[{x, y}]] :>
      NonCommutativeMultiply[w, Boson[tf2, y], Boson[tf1, x], z],
   NonCommutativeMultiply[x___, y_, z___] /; 
     FreeQ[y, NonCommutativeMultiply | Boson] :> y*NonCommutativeMultiply[x, z],
   NonCommutativeMultiply[x___, y_NonCommutativeMultiply, z___] :> 
    NonCommutativeMultiply[x, Apply[Sequence, y], z],
   NonCommutativeMultiply[x___, y_*y2__, z___] /; 
     FreeQ[y, NonCommutativeMultiply | Boson] :> 
    y*NonCommutativeMultiply[x, y2, z],
   NonCommutativeMultiply[a_] /; Head[a] =!= Times :> a,
   NonCommutativeMultiply[x___, y_Plus, z___] :> 
    Map[NonCommutativeMultiply[x, #, z] &, y],
   NonCommutativeMultiply[] -> 1,
   NonCommutativeMultiply[a_, b_] /; 
     FreeQ[{a, b}, 
       NonCommutativeMultiply] && (FreeQ[a, Boson] || 
        FreeQ[b, Boson]) :> a*b,
   NonCommutativeMultiply[x___, y___, z___] /; FreeQ[{x, z}, Boson] :>
     x*z*NonCommutativeMultiply[y]};

I wrote it for NonCommutativeMultiply, then realized we're really using CenterDot. I guess. Anyway, I just changed the rules using a rule rule, and things started to evaluate as expected.

canonicalizeRules2 = 
  canonicalizeRules /. NonCommutativeMultiply -> CenterDot;

With these rules, expand and clean can be done as below.

expand[expr_] := Expand[expr //. canonicalizeRules2]
clean[expr_] := FixedPoint[expand, expr]

CommOp[x_, y_] := clean[(x\[CenterDot]y - y\[CenterDot]x)]; 
ACommOp[x_, y_] := clean[(x\[CenterDot]y + y\[CenterDot]x)];

I made another change in how various fellas were defined, mostly to use List rather than DownValues. From the interface view, the difference is in the double vs single bracket indexing. Under the hood of course is a different matter but that's outside the scope here.

GG = Table[
   Sum[\[CapitalSigma][i, j][[m, 
       n]] \[Chi][[1, m]]\[CenterDot]\[Psi][[n, 1]], {m, 1, 4}, {n, 1,
      4}], {i, 1, Dim}, {j, 1, Dim}];

MM = Table[GG[[i, j]], {i, 1, 4}, {j, 1, 4}];
KK = Table[GG[[i, 6]] - GG[[i, 5]], {i, 1, 4}];
\[CapitalPi] = Table[GG[[i, 6]] + GG[[i, 5]], {i, 1, 4}];
\[CapitalDelta] = -GG[[5, 6]];

One last thing to observe is that since g is a diagonal matrix, many sums can be shortened considerably by contracting indices. I show such a computation with Cas2b.

Cas2b = clean[
   Sum[g[[ind1, ind1]]\[CenterDot]GG[[ind1, ind2]]\[CenterDot]g[[ind2,
        ind2]]\[CenterDot]GG[[ind2, ind1]], {ind1, 1, 6}, {ind2, 1, 
     6}]];

The result is failry short.

6 - (3 Boson[True, a]\[CenterDot]Boson[False, a])/2 - (
 3 Boson[True, b]\[CenterDot]Boson[False, b])/2 - (
 3 Boson[True, x]\[CenterDot]Boson[False, x])/2 - (
 3 Boson[True, y]\[CenterDot]Boson[False, y])/2 - 
 3 Boson[True, a]\[CenterDot]Boson[False, a]\[CenterDot]Boson[
   True, b]\[CenterDot]Boson[False, b] + 
 3 Boson[True, a]\[CenterDot]Boson[False, a]\[CenterDot]Boson[
   True, x]\[CenterDot]Boson[False, x] + 
 3 Boson[True, a]\[CenterDot]Boson[False, a]\[CenterDot]Boson[
   True, y]\[CenterDot]Boson[False, y] - 
 3/2 Boson[True, a]\[CenterDot]Boson[True, a]\[CenterDot]Boson[
   False, a]\[CenterDot]Boson[False, a] + 
 3 Boson[True, b]\[CenterDot]Boson[False, b]\[CenterDot]Boson[
   True, x]\[CenterDot]Boson[False, x] + 
 3 Boson[True, b]\[CenterDot]Boson[False, b]\[CenterDot]Boson[
   True, y]\[CenterDot]Boson[False, y] - 
 3/2 Boson[True, b]\[CenterDot]Boson[True, b]\[CenterDot]Boson[
   False, b]\[CenterDot]Boson[False, b] - 
 3 Boson[True, x]\[CenterDot]Boson[False, x]\[CenterDot]Boson[
   True, y]\[CenterDot]Boson[False, y] - 
 3/2 Boson[True, x]\[CenterDot]Boson[True, x]\[CenterDot]Boson[
   False, x]\[CenterDot]Boson[False, x] - 
 3/2 Boson[True, y]\[CenterDot]Boson[True, y]\[CenterDot]Boson[
   False, y]\[CenterDot]Boson[False, y]

I do not know if it is correct but at least it is not in any obvious (to myself) way wrong.

--- edit ---

I also mention that it could be safer/faster to use a table, then total it, with clean applied to the summands. This way we know that no individual term is overly large, hence one can hope for relatively straightforward, fast, pattern matching.

Here is an example. I also use the contracting of indices.

Timing[cas4c = 
   Table[clean[
     g[[ind1, ind1]] g[[ind2, ind2]] g[[ind3, ind3]] g[[ind4, 
        ind4]] GG[[ind1, ind2]]\[CenterDot]GG[[ind2, 
         ind3]]\[CenterDot]GG[[ind3, ind4]]\[CenterDot]GG[[ind4, 
         ind1]]], {ind1, 1, 6}, {ind2, 1, 6}, {ind3, 1, 6}, {ind4, 1, 
     6}];]

(* {423.500000, Null} *)

It's big...

cas4c // LeafCount

(* Out[362]= 2916979 *)

... but it collapses nicely.

cas4d = Total[Flatten[cas4c]];
cas4d // LeafCount

(* Out[364]= 1689 *)

There are ways to get further contraction by grouping identical consecutive NonCommutativeMultiply factors in powers, say.

Also I will repeat that the rules I used might be overkill and might contain subtle (or not so subtle) bottlenecks. So there is likely to be room for improvement in terms of speed as well as size of result.

--- end edit ---

$\endgroup$
6
  • $\begingroup$ Next time I'm in Champaign, I'll have to hunt you down. :D $\endgroup$
    – rcollyer
    Jan 6, 2014 at 22:04
  • 2
    $\begingroup$ @rcollyer We've relocated to the Siberian Steppes, to hear the weather folks describe it. My kid says it's an ordinary day for Ft Mac and we whine too much. Both of which, sad to say, are true. $\endgroup$ Jan 6, 2014 at 22:17
  • $\begingroup$ @DanielLichtblau I think there is some extra stuff in the canonicalizeRules2 definition. I tried to use your modifications to see the results but the product is not behaving distributively. Its strange because the answer that you computed for Cas2b is correct and same as what I had gotten previously. I exactly copied your rules but still I don't know what's going wrong in my notebook. $\endgroup$
    – Bilentor
    Jan 7, 2014 at 2:43
  • $\begingroup$ I'll fix the cut/paste error. If that doesn't do the job then I'll recheck tomorrow. I went through several debug cycles and it's possible the code I posted for canonicalizeRules was not the final version. $\endgroup$ Jan 7, 2014 at 3:29
  • $\begingroup$ @Daniel The distributive rule didn't work till I put it at the top of the list in 'canonicalizeRules'. Could you tell me what caused that? $\endgroup$
    – Bilentor
    Jan 7, 2014 at 14:50
1
$\begingroup$

This is really an extended comment, not an answer to the question. But who knows, maybe this solves the problem anyway. Also I guess this solves a more important question than the one asked by the OP as I guess the code was broken before.

I would rewrite your last code block as

ClearAll@expand
(*the attribute is new*)
SetAttributes[expand, HoldAll]

Unevaluated[
  expand[expr_] :=
   Block[
    {NonCommutativeMultiply (*or times*)}
    ,
    expr /.
     {
      times[left___, cnum_ /; FreeQ[cnum, (_Boson)], right___] :> 
       cnum*times[left, right]
      ,
      times[left___, 
        cnum_ /; (! 
           FreeQ[cnum, Unevaluated@Times[n___?NumericQ, ___Boson]]), 
        right___] :>
       Times @@ Apply[Power, Drop[FactorList[cnum], -1], 2]*
        times[left, FactorList[cnum][[-1,1]], right]
      ,
      times[left___, Boson[False, s_], Boson[True, s_], right___] :> 
       times[left, right] + 
        times[left, Boson[True, s], Boson[False, s], right]
      ,
      times[left___, fst : Boson[_, s_], sec : Boson[_, t_], 
        right___] :> 
       times[left, sec, fst, right] /; 
        FreeQ[Ordering[{s, t}], {1, 2}]
      ,
      times[b_Boson] :> b
      ,
      times[] -> 1
      }
    ]
  ] /. {HoldPattern[times] -> NonCommutativeMultiply }

And

(*Not equivalent to your code, this is a fix*)
expand[Alternatives[NonCommutativeMultiply, CenterDot][
   a1_, (a2_ + a3_)]] := expand[a1 ** a2] + expand[a1 ** a3]
expand[Alternatives[NonCommutativeMultiply, CenterDot][(a1_ + a2_), 
   a3_]] := expand[a1 ** a3] + expand[a2 ** a3]

This fixes only minor issues with memory, but it makes the last definitions for expand useful. Note that your definition already had the nice property that expr was evaluated inside the Module (and now inside Block), so that any aliases Centerdot will be replaced by NonCommutativeMultiply before the rules in terms of times are applied.

If you saw that I was hinting at using With, and you are curious about how stuff stuff work ReplaceAll and Unevaluated works, my explanation is that this basically an alternative way to write With. The reason I do this is that With is a scoping construct, and it tries to resolve conflicts between variables, which was undesirable in this case.

$\endgroup$
4
  • 1
    $\begingroup$ If you agree that this is nice, please replace the code in your question and I will delete the answer, as well as any comments that hinted at this, in order to clean things up and attack the real problem. $\endgroup$ Jan 6, 2014 at 14:38
  • $\begingroup$ Jacob, I agree this is neater looking and more efficient. I will replace this in my code on pastebin. $\endgroup$
    – Bilentor
    Jan 6, 2014 at 14:43
  • $\begingroup$ Jacob, I just noticed the changes and was in the process of changing things. It should be fine now. $\endgroup$
    – Bilentor
    Jan 6, 2014 at 15:06
  • $\begingroup$ @Karan I will keep the answer for now. Also note that I have added an Unevaluated around Times inside FreeQ $\endgroup$ Jan 6, 2014 at 17:11

Your Answer

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

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