Memory leak issues

Question

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.

Answer 1

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

Answer 2

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.