Populating and storing a very large array in Mathematica

Question

I am using Mathematica 8 to generate a very large multidimensional array of floating-point numbers and then store it in a file. Typical array sizes are 10000 x 50 x 15 x 40, resulting in a 3-4 GB file. The entries of the array are generated recursively, after computing a few starting values. After computations finish (about 12 hours), Mathematica appears to hang while writing the array, for 8 or more hours. Sometimes the program doesn't finish within a 24 hour time limit (imposed by the systems I'm using), and no file is produced (not even a partial one). I am wondering: what are some strategies for dealing with these large arrays in Mathematica, and ensuring writing to a file happens successfully?

Some background: the array is allocated in memory, and then populated using some initial values and a recursion relation

prec = 50;
GenerateMyArray[iMax_,jMax_,kMax_,lMax_] := Module[{myArray},
  myArray = ConstantArray[SetPrecision[0.0, prec], {iMax, jMax, kMax, lMax}];

  (* do something to populate initial entries *)
  ...

  (* populate the rest with a recursion relation *)
  Do[
    myArray[[i,j,k,l]] = Evaluate[myRecursionRelation],
    {i,1,iMax},
    {j,1,jMax},
    {k,1,kMax},
    {l,1,lMax}
  ];
  myArray
];

(The calculations must be done at high precision, although the array is eventually stored at lower precision). Finally, I write the entries to a file

mp[x_] := NumberForm[SetPrecision[x, 16], 16, ExponentFunction -> (Null &)];    

myArray = GenerateMyArray[10000, 50, 15, 40];
w = OpenWrite[myFileName];
writeHeader[w];
Do[WriteString[w, mp[x], "\n"],
    {irow,myArray},
    {jrow,irow},
    {krow,jrow},
    {x,krow}];
Close[w];

My understanding of when things are being passed by value and when they're being passed by reference in Mathematica is pretty poor. Perhaps large amounts of data are being copied unnecessarily? I would appreciate help making this process more efficient and robust.

Additional note: The output will be read by a non-Mathematica program, so using .mx files is not an option.

Answer 1

The large data framework

Ok, here is a shameless plug for the large-data framework, which I developed some time ago for similar purposes. While it is not a proper and / or final place to store it, I put its code (straight from the linked post) to this gist, so that it is easier to copy and paste as a single piece of code (note that this is not yet a package).Here, I will illustrate how and why to use it for your case.

Why use it here?

Let us start with why. The short answer: because it was written to address exactly the (rather wide) class of use cases to which the case at hand seems to belong.

These are the advantages of using it here:

• Once you convert your array to the file-backed representation, you don't have to care any more about possible kernel crashes - you can always load things back, and load only those parts you need (so, you don't have to think much about the memory either)
• You don't have to convert your list all at once, neither do you have to wait until you have constructed all its elements. In fact, you can even build the list already in the file-backed form, by using the framework's API (it has two - higher-level for typical operations on lists such as appending a list, and lower-level one for accessing and possibly modifying individual parts of a list).
• The framework adds syntactic sugar, so that, when working with parts of your list, you use code which looks like it is a normal list. In other words: while the real data structure responsible for the file-backed list is more like a hash-table filled with class instances (with each class instance pointing at a specific file location where the appropriate part of the list is stored), it is hidden by an interface which overloads most list operations and implements the list abstraction in a rather complete way. So, you can work on a much higher, and familiar to most Mathematica users, level of abstraction, and not really be concerned with the details of how the parts of your list are actually backed by files, saved, loaded, etc (not that it is too complicated, but the whole point is in information hiding - these details are largely irrelevant from the point of view of the higher-level operations we perform on lists). For you, this will look just like a regular list, in most important respects.
• You can release memory used by the loaded parts of your list selectively, and right after you used this part - so the memory used will only be as much as many parts of your list you have to keep in memory at the same time. In the case of writing to a disk, this is just one part.
• You can set your writing-to-disk code up in such a way that a possible crash in the middle will only invalidate the last part being written, and you will be able to reload the kernel and start from where you stopped - and again, with minimal memory requirements.

How to use it here

First, run the code of the framework.

Then, run this initialization code, to switch to .mx files (which are very fast):

$fileNameFunction = mxFileName;
$importFunction = mxImport;
$exportFunction = mxExport;
$compressFunction = Identity;
$uncompressFunction = Identity;

Now, we create a test array:

myArray = RandomReal[10^6, {10000, 50, 15, 40}];

This took quite a bit of RAM, so you better have at least 5Gb of RAM available. Now, initialize the symbol used to represent your list in the framework:

ClearAll[test];
initList[test]

Now, convert your list to the file-backed list - this is the main operation (use the correct path on your machine. The directory must exist)

appendList[test,myArray,
   DestinationDirectory:>"C:\\Temp\\LargeData"
];//AbsoluteTiming

{327.4692831,Null}

Note that I have a very fast SSD (like 10x the usual HDD speed), so this part may be slower on your machine, if you have a regular HDD or SSD drive. In any case, several minutes for a file of a couple of Gb large is pretty decent I think.

Note also that, if you monitor the memory consumption, the additional memory used by this procedure is very little - below 50Mb according to what I saw.

Now, the last step: save the main symbol (again, use the correct path on your machine):

storeMainList[test,DestinationDirectory:>"C:\\Temp\\LargeData"]//AbsoluteTiming

{37.4501953,Null}

For the test list I generated, this file is about 7Mb. The write time is slow here because it is a normal .m file, not .mx file.

We can now quit the kernel:

Quit

On the fresh kernel, you have to execute the framework's code again. We can then reconstruct the top-level structure of our file-backed list:

retrieveMainList[test,
     DestinationDirectory:>"C:\\Temp\\LargeData"
]//AbsoluteTiming

{66.3115234,Null}

This does not take much memory - just a few extra Mb. You can now do a few tests:

Length[test]

10000

test[[1]]//Dimensions

{50,15,40}

Let us now write a few first elements to a disk, to illustrate how this is done. I will write the first 10 out of 10000 elements.

mp[x_] := NumberForm[SetPrecision[x, 16], 16, ExponentFunction -> (Null &)];

testLength = 10;
w=OpenWrite["C:\\Temp\\myLargeArray.dat"];
time = 
  AbsoluteTiming[
    Do[
      Do[WriteString[w,mp[x],"\n"],
        {jrow,test[[irow]]},{krow,jrow},{x,krow}
      ];
      releasePart[test,irow],
      {irow,1,testLength}
    ]
  ];
Close[w];
time

{7.4648438,Null}

Note that I was sloppy with stream closing (in the sense that I don't always guarantee that the stream will be close in case of exception or Abort in the middle) - but this is not my main point here.

Note also that I use the lower-level API provided by the framework, to release the memory used by a given part of the file-backed list, right after it was used. This makes the memory use minimal.

From the timing, it will take about 2 hours for you to write the entire file. What is important is that you can arrange your code to stop at any time and resume later, without the need to load the full huge array into the kernel each time.

Taking yet more advantage of the framework

Note also, that you can, instead of transferring the array to file-backed form only at the end, do so as you build an array. You can use the appendList function to add any number of elements to the file-backed representation of your array, so you can, for example, call it on every new element (which is, 3D array in your case), or every 10 new elements, etc.

One catch here is that, to be on the safe side, you will have to call the function storeMainList periodically as well, since those parts saved to disk can not be used by themselves, without this higher-level element. This takes time (for the final array it took about 40 sec. on my machine), but if your computation takes many hours, I think you can afford doing this every half an hour say. In this way, you will also protect your computation, so that, in the case of a crash, you will be able to resume from where you stopped, rather than starting from scratch.

Answer 2

As an alternative to all the mind-bendingly clever answers, here's a dumb trick that makes the writing 20% faster:

Do[WriteString[w, StringJoin @@ ({ToString[mp[#]], "\n"} & /@ krow)],
    {irow,myArray},
    {jrow,irow},
    {krow,jrow}];

Instead of writing each element individually, this constructs a string out of the deepest row then writes it all at once.