Does parallel programming use up large quantities of memory in Mathematica?

Question

When I start Mathematica with a fresh kernel and load a program I am developing to analyze biology data, the amount of used system memory is 1.9 GB (free memory is 14.5 GB). These memory values are reported by an application external to Mathematica.

After I run the code I have developed, the amount of used system memory is 13.2 GB (free memory is 3.2 GB). So, this means that running my code has consumed 13.2-1.9 = 11.3 GB of memory. After running the code, MaxMemoryUsed[] reports a value of 2.7 GB.

1. Even though Mathematica is reporting MaxMemoryUsed at 2.7Gb, my system is still reporting 13.2 GB of used memory. If this is true, why hasn't Mathematica cleared the 11.3 - 2.7 = 8.6 GB of memory it says it hasn't used and returned it to the system?

2. I have tried setting $HistoryLength = 0 and it has not made a difference in the reported memory usage.

3. My code (posted below) does use several ParallelTable statements. Do these statements send a copy of the memory content to each parallel kernel?

4. Do you have suggestions to solve this mystery of where is my memory?

Thank you for your thoughts. Todd

--------------- Mathematica Code ------------

(* load cel data *)
time0 = AbsoluteTime[];
SetDirectory[cellocation];
celfilenames = FileNames[]; (* names of data files to import *)
celvarnames = 
  Table[StringSplit[celfilenames[[i]], {"."}][[1]], {i, 1, 
    Length[celfilenames]}] ;(* names of datafiles, with extensions \
removed, to be used downstream *)
Table[microarray[celvarnames[[i]]] = Import[celfilenames[[i]]], {i, 1,
    Length[celfilenames]}];(* assign data to microarray variable \
using celvarnames as indexes to microarray; each chip's data is \
assigned to a different index to microarray *)

chipsize = 
 chipdimensions[
  microarray[celvarnames[[1]]]] (* determine chip dimensions *)

(* load affy cdf data *)
SetDirectory[affycdflocation];
cdffilenames = FileNames[];
cdffile = 
  Import[Flatten[StringCases[cdffilenames, ___ ~~ ".cdf" ~~ ___]][[
    1]]]; 

(* convert hybidization symbol names to string names *)
experimentchips = Map[ToString, experimentchips];
controlchips = Map[ToString, controlchips];
Print["Experimental condition chips: ", experimentchips]
Print[]
Print["Control condition chips: ", controlchips]

pmindexes = 
  Table[Select[cdffile[[i, 3, 2]], #[[5]] != #[[6]] &][[All, 4]], {i, 
    1, Length[
     cdffile]}]; (* perfect match probe indexes by probeset; same \
order as affy cdf file *)

mmapmindices = 
  Map[Transpose, 
   Thread[affyindextoMMAindices[pmindexes, 
     chipsize]]]; (* create Mathematica indices equivalent to Affy pm \
index positions by probeset *) 

pmtemp = ParallelTable[
   Extract[microarray[celvarnames[[i]]], mmapmindices[[j]]], {i, 1, 
    Length[celvarnames]}, {j, 1, 
    Length[mmapmindices]}]; (* get pm data from all the chips; data \
grouped by chip; within each chip,data is organized in same order as \
probes listed in cdf file *)

Table[pmsignal[celvarnames[[i]]] = pmtemp[[i]], {i, 1, 
   Length[celvarnames]}] ;(* put the pm data into the pmsignal \
variable, using celvarnames as indexes referring to the different \
chips *)

Clear[pmtemp, 
  microarray]; (* no need to retain duplicate information in memory, \
once data is assigned to specific chips *)

Table[pmstdev[celvarnames[[i]]] = 
   Map[StandardDeviation, pmsignal[celvarnames[[i]]]], {i, 1, 
   Length[celvarnames]}];
(* compute standard deviations by pm probeset,for each chip, and \
store it in pmstdev using celvarnames as indexes *)

subsetsize = 
  0.10*Length[
    pmstdev[celvarnames[[
      1]]]]; (* Use 10% of probesets to calculate clusters below *)

Table[pmstdevsubset[celvarnames[[i]]] = 
   BlockRandom[SeedRandom[1234]; 
    RandomChoice[pmstdev[celvarnames[[i]]], Floor[subsetsize]]] , {i, 
   1, Length[celvarnames]}];
(* select a random subset of pm probset stdev for establishing \
"normal" pm probeset standard deviations for each chip *)

pmtemp = ParallelTable[
   Sort[FindClusters[pmstdevsubset[celvarnames[[i]]], 3, 
      DistanceFunction -> EuclideanDistance, 
      Method -> "Agglomerate"]][[-1]], {i, 1, 
    Length[celvarnames]}]; (* Find the largest cluster for each pm \
standard deviation subset by chip; the mean of these clusters will \
set the "normal" pm variation for each chip; only find 3 clusters to \
make calculation quicker *)

Table[pmstdevcluster[celvarnames[[i]]] = pmtemp[[i]], {i, 1, 
   Length[pmtemp]}] ;(* assign pmtemp results to pmstdevcluster using \
celvarnames as keys *)

Clear[pmtemp]; (* no need to store same data twice *)

Table[pmpostocheck[celvarnames[[i]]] = 
   Flatten[Position[pmstdev[celvarnames[[i]]], 
     x_ /; x > Mean[pmstdevcluster[celvarnames[[i]]]]]], {i, 1, 
   Length[celvarnames]}];
 (* determine each chips pm probeset positions whose st devs are \
greater than the "normal" threshold, which is calculated by taking \
the mean of the largest clusters for each chip - stored in \
pmstdevcluster *)

pmexperimentemp = 
  ParallelTable[
   Flatten[Drop[
     errorcorrection[pmsignal, experimentchips[[i]], pmpostocheck], 
     1], 2], {i, 1, 
    Length[experimentchips]}];(*calculate improved pmsignal in \
experimental chip probesets that have unusually high signal variation \
using the same probesets in comparable experimental condition chips \
as a surrogate*)

pmcontroltemp = 
  ParallelTable[
   Flatten[Drop[
     errorcorrection[pmsignal, controlchips[[i]], pmpostocheck], 1], 
    2], {i, 1, 
    Length[controlchips]}];(*calculate improved pmsignal in control \
chip probesets that have unusually high signal variation using the \
same probesets in comparable control condition chips as a surrogate*)

Table[pmsignal[experimentchips[[i]]] = 
   ReplacePart[pmsignal[experimentchips[[i]]], 
    Thread[pmpostocheck[experimentchips[[i]]] -> 
      pmexperimentemp[[i]]]], {i, 1, 
   Length[experimentchips]}]; (* replace spurious pmsignal with \
replacement values in experimental condition chips *)

Table[pmsignal[controlchips[[i]]] = 
   ReplacePart[pmsignal[controlchips[[i]]], 
    Thread[pmpostocheck[controlchips[[i]]] -> 
      pmcontroltemp[[i]]]], {i, 1, 
   Length[controlchips]}]; (* replace spurious pmsignal with \
replacement values in control condition chips *)

Clear[pmexperimentemp, 
  pmcontroltemp]; (* don't store data longer than needed *)

Table[experimentchippmmean[experimentchips[[i]]] = 
   Mean[Flatten[pmsignal[experimentchips[[i]]]]], {i, 1, 
   Length[experimentchips]}];
 (* calculate each experimental chips pmsignal mean to use in \
standardization *)

Table[controlchippmmean[controlchips[[i]]] = 
   Mean[Flatten[pmsignal[controlchips[[i]]]]], {i, 1, 
   Length[controlchips]}]; (* calculate each control chips pmsignal \
mean to use in standardization *)

Table[experimentchippmstdev[experimentchips[[i]]] = 
   StandardDeviation[Flatten[pmsignal[experimentchips[[i]]]]], {i, 1, 
   Length[experimentchips]}]; (* calculate each experimental chips \
pmsignal standard deviation to use in standardization *)

Table[controlchippmstdev[controlchips[[i]]] = 
   StandardDeviation[Flatten[pmsignal[controlchips[[i]]]]], {i, 1, 
   Length[controlchips]}]; (* calculate each control chips pmsignal \
standard deviation to use in standardization *)

pmexpstandtemp = 
  ParallelTable[
   Map[(# - experimentchippmmean[experimentchips[[i]]])/
      experimentchippmstdev[experimentchips[[i]]] &, 
    pmsignal[experimentchips[[i]]]], {i, 1, 
    Length[experimentchips]}];(* standardize each experimental chips \
pmsignal *)

pmcontrolstandtemp = 
  ParallelTable[
   Map[(# - controlchippmmean[controlchips[[i]]])/
      controlchippmstdev[controlchips[[i]]] &, 
    pmsignal[controlchips[[i]]]], {i, 1, 
    Length[controlchips]}];(* standardize each control chips pmsignal \
*)

Table[experimentstandard[experimentchips[[i]]] = 
   pmexpstandtemp[[i]], {i, 1, Length[experimentchips]}];
 (* reassign data from pmexpstandtemp to proper indexed variable *)

Table[controlstandard[controlchips[[i]]] = 
   pmcontrolstandtemp[[i]], {i, 1, 
   Length[controlchips]}];(* reassign data from pmcontrolstandtemp to \
proper indexed variable *)

Clear[pmexpstandtemp, pmcontrolstandtemp]; (* free up memory *)

time1 = AbsoluteTime[] - time0

Answer 1

First off, you shouldn't be worried that Mathematica hasn't returned the memory to the system. The memory may as well be (and is being) held onto by the process that last had it, and if some other program needed this memory it would be handed back to the system without a problem.

MaxMemoryUsed reports the maximum amount of memory used, but it is only interacting to the current kernel. If you wanted to interact with the other kernels, you would have to use ParallelEvaluate. This paper contains a nice example of a table displaying the max memory used on each processor, along with the processor ID and some other information. (You will have to update some of the commands though, it was written in 2004.)

When I have used ParallelTable commands in the past I always was sending the bulk of the data in my current application for processing through ParallelTable, and so an entire copy of the data often was transferred to the memory of each parallel kernel. The answer by acl points out that the memory sent to each parallel kernel is dependent on what you send to it, and it may not necessarily need to send the entire memory state from the main kernel.

What OS are you running Mathematica on? In Linux, pulling up the System Monitor (or a command like top) should show the kernels loading in real-time, and if you sort by memory you should get a very clear picture of what's happening.

I don't have any tips to avoid this - it's a trade off, if you want faster processing of data and you have the RAM to spare then use the ParallelTable command. If you don't have the RAM to spare then stick with the regular Table command and brew some coffee.

As a last ditch attempt to scrape up some spare RAM, you can try throwing the Share[] command in every once in a while, but I haven't seen significant improvements from this myself either. As long as you're not actually running out of RAM, however, it's nothing to worry about. The fact that Mathematica hangs on to that extra RAM even after it's done processing the data is standard procedure, and you shouldn't concern yourself with it. If you want to prove that another application can access that RAM if needed, just run some RAM intensive process and watch in real time using the System Monitor / top.

---

I was reminded by Albert Retey that you can use distributed contexts to reduce the amount of information transferred to each kernel. The effectiveness of this approach will vary with your particular application, i.e. if you have a lot of information in memory before you perform ParallelTable on a small subset of that data.

Answer 2

I will only address a subquestion of yours, namely, "Does ParallelTable send a copy of the memory content to each parallel kernel?". The answer to that is a clear no: using ParallelTable does not mean that each kernel gets a copy of everything, but only of a subset of things (see the docs for DistributeDefinitions to see how the set of symbols to be distributed is determined).

Open some sort of system monitor utility (Activity Monitor or top on OS X, top on other Unices, the equivalent in windows) and watch for MathKernel. It should show a single kernel:

Try this:

Then launch 2 kernels using LaunchKernels[2]:

Then waste some memory:

tbl1 = Developer`FromPackedArray[ConstantArray[0, 5*10^7]];

As you can see, this increases the memory size of the master kernel but not the others.

Now send something to the kernels:

ParallelTable[$KernelID, {i, 1, 10}]
(*{2, 2, 2, 1, 1, 1, 2, 2, 1, 1}*)

Clearly this did get executed in the slave kernels. Their memory use has not increased:

So it looks like tbl1 did not get sent to the kernels; they don't get a complete copy of everything the master kernel has.

Answer 3

Others have given you some pointers about what's happening (and what not). I will focus on a simple self contained example which shows how you can set up a parallel computation which will not make copies of a large set of data but still make parts of that data accessible from the parallel kernels. Be aware that the approach I'm using here might have it's drawbacks and could not be very efficient, depending on the details of your problem. There might be better ways to achieve the same thing for your use case, but it's hard to say without a specific example that we can actually run. If you want better solutions, I would suggest you provide a self contained problem which is as simple as possible but still shows the main characteristics of your real problem (runtimes, datasize,...).

This is one way to make accessible parts of a large dataset to be used by kernels running in parallel without making a full copy on each parallel kernel:

LaunchKernels[2]

(*
==> {KernelObject[1, "local"], KernelObject[2, "local"]}
*)

{MaxMemoryUsed[], MemoryInUse[]}
ParallelEvaluate[{MaxMemoryUsed[], MemoryInUse[]}]

(*
==> {19909448, 19564448}
*)

(*
==> {{13627720, 13565704}, {13628376, 13566360}}
*)

data = Developer`FromPackedArray[RandomReal[{0, 1}, {5*10^6}]];

{MaxMemoryUsed[], MemoryInUse[]}
ParallelEvaluate[{MaxMemoryUsed[], MemoryInUse[]}]

(*
==> {220536232, 180507576}
*)

(*
==> {{13627720, 13565704}, {13628376, 13566360}}
*)

SetSharedFunction[getData];

getData[start_, end_] := data[[start ;; end]];

{MaxMemoryUsed[], MemoryInUse[]}
ParallelEvaluate[{MaxMemoryUsed[], MemoryInUse[]}]

(*
==> {220536232, 180519008}
*)

(*
==> {{13627720, 13571624}, {13628376, 13572280}}
*)

Total[
  ParallelTable[
   Total[getData[(i - 1)*5*10^5 + 1, i*5*10^5]], {i, 1, 10}, 
   DistributedContexts -> None,Method -> "FinestGrained"]] // AbsoluteTiming

(*
==> {1.2168021, 2.49971*10^6}
*)

{MaxMemoryUsed[], MemoryInUse[]}
ParallelEvaluate[{MaxMemoryUsed[], MemoryInUse[]}]

(*
==> {220536232, 180531696}
*)

(*
==> {{49573432, 29572336}, {49574088, 29572992}}
*)

Total[
  Table[Total[getData[(i - 1)*5*10^5 + 1, i*5*10^5]], {i, 1, 
    10}]] // AbsoluteTiming

(*
==> {0.9048016, 2.49971*10^6}
*)

Here are some explanations:

• I'm using SetSharedFunction to make the definitions of getData accessible by all kernels. Unlike using SetSharedVariable for the data itself, SetSharedFunction used this way makes no copies of the full data to the local kernels once accessed there but only those parts of the data that are actually needed. In fact this will even be true when you actually store the data in DownValues of the function that is shared. As stated in the documentation, the way I'm defining getData it will be evaluated on the master, so no copy will be made but only the result of getData will be transferred.

• Note that the efficiency of this approach will depend crucial on the partition of the data, only if these chunks have an appropriate size you can hope to see speedup. There is also some interdependency between the chunksize you use in the call to ParallelTable and it's Methods option, which I have set so that it's my partition of the data that is used.

• Using SetSharedFunction will lock the access to the data whenever one of the kernels tries to read or write. I'm not sure how smart the implementation of this locking actually is and whether it will recognize that the kernels actually only do read data, so no locking actually would be necessary. If not, it might well be that the overhead is much higher than it would need to be and it would pay off to find ways to avoid the locking.

• I have put some code in between the different parts so that you can check that the memory consumption is really what I am claiming (or see if that's not the case for you, which I can of course not exclude).

• For the example I show the serial code will be much faster -- the computational part is just too cheap to see any speedup, so that example is basically just measuring overhead. Only when the computational cost is high enough to compensate the overhead of sending around the data (and the additional time the kernels might be waiting due to the locking) you can hope to see any speedup (very first lesson from parallel coding trivia, but I think you either know that or could make it up yourself). On my computer the overhead for the above code seems to be in the order of half a second, so I'd expect speedup only when the computation for a single number takes at least a view milliseconds.

• In this example the result of each parallel computation is just a number. If the results would also be large datasets you could try using the same or a similar trick to get the output data back to the master. I have not actually tested whether that will also work as it does for the input case but would expect so.

The above is just a sketch on how you can proceed, I have no experience on how well this will work/perform with real use cases and it might well need some optimization. On the other hand I'm sure that distributing parts of a large dataset can be solved in an efficient way for every specific use case, it's probably just not easy or even possible to do it in a very general way. As explained above, I'm quite sure you'll get a lot more sophisticated answers once you post a self-contained example that people can play with...