Efficiently generating histograms with large numbers of bins

Question

I'm trying to generate a histogram for a large data set with the following lines of code:

MaxSideLength = Ceiling[Max[data]];
Export["Histogram.png",
 Image[ArrayPlot[HistogramList[data, {1, MaxSideLength, 1}][[2]],
   ColorRules -> {0 -> Black, 1 -> Red, 2 -> Yellow}], ImageSize -> MaxSideLength]];

Where ImageMatrix has something like $O(10^6)$ datapoints and I want to have maybe $50,000 \times 50,000$ bins (where each bin corresponds to a pixel in the output PNG image).

Unfortunately, my memory usage goes through the roof with the above script (exhausting over 200 GB of RAM). I can pull off having maybe $20{\rm k} \times 20{\rm k}$ bins, but anything above this fails.

Is there a more efficient way to proceed? If need by, I can properly round all of the values in data.

I know this comparison is unfair, but given that the typical output size for a $20,000 \times 20,000$ pixel output PNG is something like 40 Mb, it seems a little odd to me that the binning and image creation process would require over $2500x$ the memory in RAM.

Ok, let's run a little analysis script provided by ssch:

mmu = MaxMemoryUsed[];
n = 20*10^6;
{w, h} = {5000, 5000};
data = RandomVariate[NormalDistribution[], {n, 2}];
bins = BinCounts[data,
 {Min[data], Max[data], (Max[data] - Min[data])/w},
 {Min[data], Max[data], (Max[data] - Min[data])/h}];
i = Image[bins /. {0 -> {0, 0, 0}, 1 -> {1, 0, 0}, 2 -> {1, 1, 0},
                   3 -> {0, 0, 1}, _?IntegerQ -> {1, 1, 1}}];
mmu2 = MaxMemoryUsed[];
mmu2 - mmu

The result for my value of n is $\approx 600*10^6$. Increasing n by an order of magnitude to $20*10^7$ yields a memory usage of $\approx 12*10^9$, however, if we run the script again, overwriting the "data" data structure, the memory usage reported is again only $\approx 600*10^6$.

However, if we increase the bin sizes to {w, h} = {10000, 10000}, the memory usage jumps to $\approx 3.5*10^9$ (for the first run), and then oddly $\approx 2.5*10^9$ for subsequent runs (since I suppose we're overwriting previous data structures). My conclusion is that an increase in the number of bins is responsible for the blowup in RAM usage.

Answer 1

This seems to work -- well, on the small example. I don't have 200GB to play with.

Th OP's setup (without the image):

Quit[]

SeedRandom[2];
n = 20*10^6;
{w, h} = {5000, 5000};
data = RandomVariate[NormalDistribution[], {n, 2}];

mmu = MaxMemoryUsed[];
bins = BinCounts[
    data, {Min[data], 
     Max[data], (Max[data] - Min[data])/w}, {Min[data], 
     Max[data], (Max[data] - Min[data])/h}]; // Timing

mmu2 = MaxMemoryUsed[];
mmu2 - mmu

(* {12.107534, Null} *)
(* 960687112 *)

I adapted the OP's color function to create a packed array for the image data:

cf = Compile[{{bins, _Integer, 2}},
    Map[Piecewise[{
        {{0., 0., 0.}, # == 0},
        {{1., 0., 0.}, # == 1},
        {{1., 1., 0.}, # == 2},
        {{0., 0., 1.}, # == 3}},
       {1., 1., 1.}] &,
     bins, {2}],
    RuntimeOptions -> "Speed", CompilationTarget -> "C"
   ];

One might convert the image to "Byte" to save space, considering the colors of the color function.

(img = Image[Image @ cf @ bins, "Byte"]) // Timing

It does save space:

ByteCount @ img
(* 75000400 *)

ByteCount @ Image @ cf @ bins
(* 600000384 *)

Export finishes in a reasonable time (less than half the time it took for the larger Image@cf@bins).

With[{t0 = AbsoluteTime[]},
 Export["1Example.png", img];
 AbsoluteTime[] - t0]
(* 3.000048 *)

Possible alternative

There is another alternative, based on the answers by kguler and ruebenko to these questions respectively:

1. Optimising 2D binning code

2. https://stackoverflow.com/questions/8178714/mathematica-fast-2d-binning-algorithm

It uses SparseArray and it can be fast. Indeed when the OP's n was originally 2 * 10^6 it seemed promising. But when it changed to 20 * 10^6, the performance was considerably worse. Perhaps, if in the actual application, n is much smaller than w * h, it may still be worth considering. It seemed to use a lot a memory in constructing the SparseArray, but the final result can be quite a bit smaller than bins.

See the linked questions about the option "TreatRepeatedEntries". I present three variations depending on how one wants to treat Max[data]. The OP's code tends not to count the maximum, which lands just at the beginning of the bin after the last one, unless due to round-off error it just happens to land in the last one.

This is a direct adaptation of the answers to the linked questions. I assume h = w as they happen to be equal in the OP's question. The original answers fudge a little with epsilon, but the variations below avoid this.

epsilon = 10^-12;
indexes1 = 1 + Floor[(1 - epsilon) w Rescale[data]];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
bins1 = SparseArray[indexes1 -> 1, {w, w}];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];

Here is one that does what the OP's BinCounts code does. This one lops off the maximums with ArrayPad.

indexes2 = 1 + Floor[w Rescale[data]];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
bins2 = ArrayPad[SparseArray[indexes2 -> 1, 1 + {w, w}], {0, -1}];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];

One can put the maximum in the last bin with this:

indexes2 = # + Unitize[#-w] &@ Floor[w Rescale[data]];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
bins2 = SparseArray[indexes2 -> 1, {w, w}]];
SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];

Note if w != h1, then these variations would need to be adapted.