Skip to main content

programming - Memory leak issues


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



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


Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...