Memory efficient generation and selection of tuples

Question

The short of this question: I need to find a bunch of n-tuples of {0,1} which do not fail to satisfy a set of (non-linear polynomial) equations, without hogging up all of my memory by trying to generate every possible n-tuple before eliminating the bad ones. This set of 'not bad' n-tuples should be substantially smaller than every possible n-tuple for the problem I am working on.

The long of this question:

I have a system of non-linear (polynomial) equations I want to solve; indeed, multiple such systems, each one corresponding to a separate case of interest to a larger problem. Simply using Mathematica's built-in Solve function isn't very effective, as it quickly consumes memory and then run's at a snail's pace trying to work its way through the problem without consuming extra memory. Thankfully, the systems I'm working with have a couple of nice properties which give a more efficient workaround (indeed, my current method converts what could have taken days with Solve into mere minutes).

The one of relevance here is that, if I call the variables in the system d[i], then it turns out there are a fair number of equations of the form d[i]^2==d[i], which means all of those variables have to be either 0's or 1's. I am currently using Tuples to generate every possible substitution, and then proceeding to eliminate all of those which cannot satisfy the equations. However, this gets extremely memory intensive pretty fast. Even 'small examples' lead to insufficient memory errors at this stage, as they can have up to a few dozen of these squaring relations in each list of equations (and there can be half-dozen or more such lists in a 'small' example). There may be a number of ways to improve the efficiency of my code, as I'm still fairly new to advanced Mathematica programming, but mostly I am looking for ways to keep the memory demands down without drastically inflating computation time.

Before getting to this stage, I generate a nested list of equations I'm interested in, denoted CubeReduced2 (the equations end up being degree 2 multinomials; the reason for the 'cube' is not particularly germane here). The equations are all stored on the final level; the preceding levels effectively serve as markers of which case I'm looking at it, and do not store any other data. I am attempting to solve many possible systems without having to do them individually, especially when I have no advanced knowledge of how many systems there will be.

This is the code I am currently using to generate all possible combinations of 0/1 substitutions for these variables

SelfSquares = Flatten[Map[Reap[Do[If[MemberQ[#, d[i]^2 == d[i], Sow[i]], 
{i, 1, varnum}]][[2]]&, CubeReduced2, {2}], {{1}, {2}, {3, 4}}];
SquarePossibilities = 
 ParallelTable[
  Tuples[{0, 1}, Length[SelfSquares[[m]][[i]]]], {m, 1, 
   Length[CubeReduced2]}, {i, 1, Length[CubeReduced2[[m]]]}];
SquareReplaces = 
 ParallelTable[
  Map[Table[
     d[SelfSquares[[m]][[i]][[j]]]] -> #[[j]], {j, 1, 
      Length[SelfSquares[[m]][[i]]]}] &, 
   SquarePossibilities[[m]][[i]]], {m, 1, Length[SquarePossibilities]}, {i, 1, 
   Length[SquarePossibilities[[m]]]}];

The SquarePossibilities stage is where my system memory quickly vanishes. I Remove it after generating SquareReplaces, but when I can't even build the list that doesn't help any.

I would then use the following code to eliminate the tuples which definitely won't work:

SquareReplaces2 = 
  ParallelTable[
   Select[SquareReplaces[[m]][[
     i]], (And @@ 
       Map[! TrueQ[# == False] &, ((CubeReduced2[[m]][[
            i]]) /. (#))]) &], {m, 1, Length[CubeReduced2]}, {i, 1, 
    Length[CubeReduced2[[m]]]}];

There may be non-trivial equations left after apply the replacements, which will still need to be solved for. That's why I don't simply look for every equation to evaluate to True, but simply look for those tuples which don't produce a False. Hopefully this explanation is sufficient at this point to describe what I'm trying to accomplish.

I'd provide a sample list of equations I have applied this code to (succesfully), but right now my mathematica has the variables as $d_i$, which don't translate very well into an easy copy-paste. I'll try to fix that if a relevant example is desired.

Answer 1

Your explanations are very detailed and I have to admit it's too detailed for me to go through it, trying to understand without a minimal working example. Therefore, view this as some ideas for your first, short-version question. I will present two different approaches, where the first one takes very long, but consumes almost no memory and the second one is very fast, but uses much memory during the computation. Finally, I'll show you how you can combine them so that they work with your memory restrictions.

When I understood you correctly, then you have a function which tells you whether or not a combination {0,1,0,0,1,1,..} is relevant. I will call this function selector and what I use for the demonstration here is very simple: A combination is relevant if it contains exactly two 1's. A simple definition of this selector is therefore

selector[combination_] := Plus @@ combination === 2;

Whether or not my idea works for your problem depends partially on your selector. Here the question is for instance, can we include your selector in compiled code or how fast it is.

Another issue you'll see here is, that I don't make use of ParallelTable or related parallel constructs. The reason is simple: Since those functions distribute your problem to many subkernels they need to distribute parts of your data as well. If you are not absolutely sure that you know what you do there it can happen very easily that this consumes much memory while the speed up is often non-existent. In fact, I will make use of another kind of parallelization in my second approach.

General idea

The general idea bases on the fact that it is quite easy (and memory efficient) to create all tuples of {0,1} if you remember that all numbers are stored binary inside your computer. Therefore, to create all combinations of length 4 you only need to count from 0 to 15 and look at the binary representation. Mathematica has the IntegerDigits function to give you the representation of a number in a different base. Therefore an equivalent algorithm to Tuples[{0, 1}, 4] is the following

IntegerDigits[Range[0,15],2,4]
(*
{{0,0,0,0},{0,0,0,1},{0,0,1,0},{0,0,1,1},
 {0,1,0,0},{0,1,0,1},{0,1,1,0},{0,1,1,1},
 {1,0,0,0},{1,0,0,1},{1,0,1,0},{1,0,1,1},
 {1,1,0,0},{1,1,0,1},{1,1,1,0},{1,1,1,1}}
*)

The idea has two advantages. Firstly, you can just go through all combinations very easily by counting. Secondly, you don't need to store the {0,1,..} combinations explicitly because you can have them very memory efficient in a list of integers.

1. Approach: Slow but memory efficient

With the above idea the most direct algorithm to test all n-tuples of {0,1} is to go through all numbers 0 to 2^n-1 and test each binary representation whether for its relevance. If it is relevant then you store it. With this you consume (more or less) only the memory which your final result list would require anyway. Again, remember that you don't need to store a list of combinations as result because a list of integers representing the combinations is sufficient.

To try this we define a function iterativeSelect which implements the idea. It further sets a progess variable so that we have a feeling how long it will take. To store a relevant combination I Sow the integer representing the combination. All result combinations are collected using Reap at the outside. Reaping and sowing is a good way to collect results memory efficient and I strongly advice you not to use something like AppendTo.

iterativeSelect[n_, test_] := Last@Last@Reap[Do[
  If[Mod[i, 2^10] === 0, progress = i];
  If[test[IntegerDigits[i, 2, n]],
    Sow[i]
  ], {i, 0, 2^n - 1}
 ]
]

With selector defined as in the beginning of this answer, let's test this approach for n=25

n = 25;
ProgressIndicator[Dynamic[progress], {0, 2^n - 1}]
m1 = MemoryInUse[]/2^20.;
result = iterativeSelect[n, selector]; // AbsoluteTiming
m2 = MemoryInUse[]/2^20.;
m2 - m1

This takes about 156 seconds here and consumes almost no memory (~5kB). The equivalent run with Tuples

Select[Tuples[{0, 1}, 25], selector] // Length

occupies about 17GB of memory during the run (but frees it when it is finished) and takes maybe 1.5 minutes.

If we want to make an extrapolation how long this method would run for n=30 we could remember that

$$2^{30} = 32\cdot2^{25}$$

so I guess something around $156\cdot32$ seconds which is about 80 minutes.

2. Approach: Fast but memory consuming

The second approach I'm suggesting is to make a compiled function which takes a list of integers representing the combinations and tests the whole list in parallel. What you get is a list of True|False which you can use with Pick to select only relevant combinations.

This clearly is memory consuming because we need the integer list completely during the run. For n=30 this requires about 8GB

ByteCount[Range[2^30]]/2^20.

Remember, that we need a second list of True|False which is the result of the parallel selector. Furthermore, we don't now what Pick does internally.

Another restriction is that you have to be able to compile your selector function completely, otherwise it will not be of much use. Long story short, in the next code block the options to Compile are significant.

parallelSelector = Compile[{{n, _Integer}},
  Plus @@ IntegerDigits[n, 2] === 2,
  CompilationTarget -> "C",
  Parallelization -> True,
  RuntimeAttributes -> {Listable},
  RuntimeOptions -> "Speed"];

What you should always do is to use CompilePrint from the <<CompiledFunctionsTools` package the check whether your compiled function is free of any MainEvaluate call.

A speed test of the parallelSelector function only reveals how fast it is compared to the first approach

parallelSelector[Range[0, 2^25]]; // AbsoluteTiming
(* 3.094913 seconds *)

To create now all relevant combinations, we create the True|False list and select the True entries using Pick

result = Block[{data = Range[0, 2^25]},
  Pick[data, parallelSelector[data]]
]; // AbsoluteTiming

This takes about 5 seconds here and with 32GB of RAM I didn't see any noticeable memory consumption during the run, but a list of 2^25 integers needs only 256MB of memory. If I do the same calculation for n=28 on my machine it finishes very fast within 40 seconds but it needs a peak memory of about 11 GB during the run

Combining the two approaches

Combining the two approaches is now easy. Assuming that on your machine has enough memory to run parallelSelector for n=25 easily and you require to calculate n=30. Then you make a sequential loop that iterates to 2^30 in a stepsize of 2^25. In each step i (of the 32 steps) you calculate the interval [i,i+2^25-1] in parallel. This is kind of similar to what Sjoerd suggested in his first comment. We can even approximate the runtime quite exact. Since one parallel run for n=25 takes about 5 seconds on my machine and we need 32 steps to iterate to 2^30 I expect about 160 seconds. And indeed

sequencialParallelIterate[n_, step_, parallelTest_] :=
  Flatten@Last@Last@Reap@Do[
    Sow@Block[{data = Range[i, i + 2^step - 1]},
      Pick[data, parallelTest[data]]
    ], {i, 0, 2^n - 1, 2^step}
 ]

sequencialParallelIterate[30, 25, parallelSelector]; // AbsoluteTiming

needs 166 seconds here. If I'm allowed to refer to my system monitor once more:

First row shows the cpu usage and the second one again the memory in use. You see that the memory consumption is constantly low because it only needs what is used in a parallel run. When you correlate this with the cpu usage you see that the memory peaks are where the parallel selector runs. Additionally, you clearly see every single iteration because then all my core work very hard.

Conclusion

The benefit of the presented method is that you can tune it easily to your memory/runtime needs. If you are low on memory, you calculate smaller chunks in parallel and iterate more often. If you have enough memory you can make the chunks who are computed with a low-level parallelization larger and you gain speed.

The fast parallel method strongly depends that you don't do very complex things in your selector. If your selector depends on global variables they can most probably be included in the compiled function (using With or another strategy) but if you need something non-trivial like a call to Solve things get complicated if not impossible.

Answer 2

Here's the solution I came up with along the lines of running over tuples of a smaller set instead of the entire set, and building up the answer from there.

varnum is the maximum possible number of variables d[i], and CubeReduced2 is the set of equations under consideration.

I'm not sure how to best parallelize the code, so I'll use a sequential code. I've adjusted the code to use halirutan's idea of using binary representations of integers instead of actual tuples. It turns out, in the examples I've done so far, that it is extremely efficient to work on 5-tuples at a time (instead of n-tuples, with 2^n close to what your machine can handle). This is a feature of the exact equations I'm dealing with, rather than an innate feature of the code. It's convenient and somewhat unexpected. In any case, the size of the intermediate tuples can be adjusted up to suit the situation. Here's the code:

Clear[ValidTuples, SubBuild]
ValidTuples = {{}};
SubBuild[m_, n_] := d[m] -> n;
SetAttributes[SubBuild, Listable];

Module[{SelfSquares, SelfSquaresPartition, SquareReplaces, len}, 
 SelfSquares = 
  Flatten@Last@
    Reap@Do[If[
       MemberQ[CubeReduced2, d[i]^2 == d[i]] || 
        MemberQ[CubeReduced2, d[i]] == d[i]^2, Sow[i]], {i, 1, 
       varnum}];
 SelfSquaresPartition = Partition[SelfSquares, 5, 5, {1, 1}, {}];
 Do[
  len = Length[SelfSquaresPartition[[i]]];
  SquareReplaces = 
   Map[SubBuild[SelfSquaresPartition[[i]], 
      IntegerDigits[#, 2, len]] &, Range[0, 2^len - 1]];
  ValidTuples = 
   Flatten[Table[
     Map[Join[#, SquareReplaces[[x]]] &, ValidTuples], {x, 1, 
      Length[SquareReplaces]}], {{1, 2}, {3}}];
  ValidTuples = 
   Select[ValidTuples, 
    And @@ Map[! TrueQ[# == False] &, CubeReduced2 /. #] &], {i, 1, 
      Length[SelfSquaresPartition]}]]

The code builds up a list of all indices i such that d[i]^2==d[i], then Partitions them into smaller chunks. Very small chunks, as mentioned. It then successively builds up which of the 0/1 substitutions do not fail to satisfy the equations.

The final answer is contained in ValidTuples. Which are actually substitution rules, so that's a bit of a misnomer.

Since the code works on only 5-tuples at a time, it's possible using the binary trick is not worth it, but the algorithm is now robust enough to handle some larger examples, at which point a larger value than 5 may be desirable, and the binary digits much more productive.

This code (or rather, an altered version to deal with the various levels of the lists present in the actual data) has successfully worked it's way through an example for which my original code received insufficient memory errors. And it did so in a respectable 3.9 seconds.

Converting the variables from $d_i$ to d[i] also seems to have sped up the entire spreadsheet, as a bit of a side note. I'll definitely remember that.