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

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1.