Random Sample from Very Large Table

Question

I'm currently trying to create some test data to check a numerical implementation against Mathematica's implementation of some functions.

In order to do this, I've been generating a broad range of arguments with a (Parallel)Table which works fine for a function with a few arguments, but I have one function I want to test which takes 4 integer arguments and 10 floating point arguments:

f[n1, n2, n3, n4, x0, x1, x2, x3, x4, x5, x6, x7, x8, x9]

The thing is creating a broad range of values use a table results in a simply enormous amount of data which cannot be held in memory. For example even:

realSample[n_] := Join[-10^Subdivide[-10, 10, n], 10^Subdivide[-10, 10, n], Subdivide[-10, 10, n]];
Table[{n1, n2, n3, n4, a, b, c, d, e, f, w, x, y, z},
  {n1, 0, 1}, {n2, 0, 1}, {n3, 0, 1}, {n4, 0, 1},
  {x0, realSample[1]}, ...
]

results in something that is just unmanageable even for n = 1 simply because of the sheer number of arguments.

So instead, I'd like to take a random sample of the whole table (say a million); however, RandomSample doesn't stream the table and the same problem as before still comes up.

Short of implementing a reservoir sampling, is there anything in Mathematica that allows me to take a random sample of the full table without loading the whole table into memory?

Answer 1

So one solution alluded to in the question is to use reservoir sampling which only holds in memory the desired final number of items while iterating through the table:

Attributes[ReservoirSample] = {HoldRest};
ReservoirSample[n_, arg_, iter__] := Block[{
    sample = {}, p = 0
  },

  Do[
    p += 1;
    If[Length[sample] < n,

      (* Fill the reservoir to begin with *)
      AppendTo[sample, arg]
      ,

      (* Otherwise we add the next argument with decreasing probability *)
      If[RandomInteger[{1, p}] <= n,
        sample[[RandomInteger[{1, n}]]] = arg;
      ];
    ];
    ,
    iter
  ];

  sample
];

Checking the output in a very simple case does seem to produce a flat distribution as desired:

Join @@ Table[ReservoirSample[
    10,
    x,
    {x, 10},
    {y, 10},
    {z, 10}
  ], 100] // Histogram

I have tried implementing a parallel reservoir sampling, but due to the way Mathematica subdivides the iteration between the kernels, I'm consistently getting samples which are very much non-uniform.

Answer 2

Following the suggestion of @HenrikSchumacher, here is another solution which creates the entries based on randomly selecting a value for each iteration variable at each step.

This is much faster than reservoir sampling in my other answer as it does not need to iterate through every possibility. The main drawback however is that this can produce duplicates, especially if the number of samples requested isn't much smaller than the total number of possibilities. The reservoir sampling is guaranteed to not do that.

Attributes[RandomAccessSample] = {HoldRest};
RandomAccessSample[n_, arg_, iter__] := Block[{
    iterArgs, iterLists
  },
  (* Separate the variables we're iterating over from the corresponding lists,
     and convert range specifications into lists for later. *)

  iterArgs = First /@ {iter};
  iterLists = Table[
    If[Head[i[[2]]] === List, i[[2]], Range @@ i[[2 ;;]]],
    {i, {iter}}
  ];

  Table[
    arg /. Thread[iterArgs -> RandomChoice /@ iterLists],
    n
  ]
];

The above can easily be parallelized by replace Table with ParallelTable and the results remain the same.

A quick check that it samples uniformly (though note here that it is repeating elements as it produces 100'000 samples from 1000 possibilities):

RandomAccessSample[100000,
  x,
  {x, 0, 10},
  {y, Range[10]},
  {z, 10}
] // Histogram