Scan through (partial) tuples

I will show how similar tasks can be done using (undocumented) Streaming framework. The usual caveat applies: since this is undocumented functionality, there is no guarantee that it will exist in the future versions in the same exact form, or at all. But I thought it may be a nice application to illustrate some of the ideas behind Streaming, as well as to have a somewhat deeper look into its internals and possible use cases.

Preparation

Download and install a patch for Streaming, if don't yet have it

To start with, one would have to follow this Q/A to patch Streaming, which is necessary to run it on modern versions of Mathematica. Specifically, you need to execute:

Import[StringJoin[
  "https://raw.githubusercontent.com/", 
  "lshifr/StreamingPatch/master/StreamingPatchBootstrap.m"]
]
DownloadAndInstallStreamingPatch[]

where I used StringJoin because the bug in SE code editor would not allow me to enter long strings. This step has to be done only once on a given machine.

Streaming and LazyTuples

Since my post has become very long, I have collected essential parts into a gist, which can be easily imported. For lazy tuples, I will use the function TuplesFunction from this excellent answer by Carl Woll, with a few minor changes.

The initialization code is then:

$HistoryLength = 0;
Get["StreamingPatch`"];
Import @ StringJoin[
  "https://gist.githubusercontent.com/lshifr/",
  "7d5d3ef064b2eb1a1e9b997726923331/raw/LazyTuples.m"
]

Example data

We start by defining sample lists for tuple elements:

lists = Table[RandomSample[Range[30], 20], 5]


(* 
  {
   {4, 5, 7, 16, 20, 17, 14, 6, 22, 25, 12, 15, 19, 1, 2, 13, 29, 26, 10, 3}, 
   {7, 4, 21, 28, 29, 3, 10, 9, 15, 12, 17, 14, 1, 11, 16, 2, 23, 30, 6, 5}, 
   {1, 7, 4, 18, 17, 21, 16, 30, 26, 11, 3, 22, 6, 14, 24, 19, 12, 23, 9, 28}, 
   {28, 23, 30, 4, 11, 16, 27, 20, 6, 14, 24, 18, 12, 3, 7, 21, 2, 25, 26, 1}, 
   {11, 17, 5, 1, 18, 25, 4, 12, 16, 24, 2, 9, 10, 8, 3, 23, 28, 26, 19, 15}
  }
*)

This defines tuples with the length 20^5 == 3200000 and total ByteCount about 128Mb:

Times @@ Length /@ lists
ByteCount[Tuples[lists]]

(* 3200000 *)
(* 128000208 *)

Main example

LazyList brief intro

Evaluating LazyTuples[lists] will create a LazyList object, which is a lazy representation of a list of tuples, with a default chunk size 10000.

Sections below have more detailed explanations of the internals, but for practical purposes one can in many ways use LazyList in place of a normal List. Many operations on LazyList, such as e.g. Select and Take, are lazy - they produce a new LazyList with very little computation. The real work is done when Normal is applied to a LazyList object, which converts it to an equivalent normal List.

Illustration

This creates a LazyList object (delayed assignment is deliberate, for the purposes of benchmarking we want a new instance to be created every time):

lt := LazyTuples[lists]

One can convert it to a normal list

Normal[lt] // Length // AbsoluteTiming

(* {1.3135, 3200000} *)

and verify that it produces the result identical to standard Tuples:

Normal[lt] == Tuples[lists]

(* True *)

In the following, we will search for all tuples which sum exactly to 100, and measure performance and memory footprint:

Normal @ Select[lt,  Total[#]==100&]//Short//AbsoluteTiming 
Normal @ Select[lt,  Total[#]==100&]//MaxMemoryUsed
Normal @ Take[Select[lt,  Total[#]==100&], 1]//AbsoluteTiming 

(*
   {5.91346,{{4,21,17,30,28},{4,21,21,28,26},{4,21,21,30,24},{4,21,21,26,28},<<23635>>,{3,30,28,14,25},{3,30,28,24,15},{3,30,28,21,18}}}

  7595088

  {0.108224,{{4,21,17,30,28}}}
*)

This shows that it takes about 6 seconds to process entire list of tuples in this case, and about 0.1 sec to get just the first match, while memory use does not exceed 8Mb.

We now compare that to in-memory computations using Tuples:

Select[Tuples[lists], Total[#]==100&]//Short//AbsoluteTiming 
Select[Tuples[lists], Total[#]==100&]//Short//MaxMemoryUsed 
Select[Tuples[lists], Total[#] == 100 &, 1]//Short//AbsoluteTiming 

(* 
   {4.55724,{{4,21,17,30,28},{4,21,21,28,26},{4,21,21,30,24},{4,21,21,26,28},<<23635>>,{3,30,28,14,25},{3,30,28,24,15},{3,30,28,21,18}}}

  692115784

  {1.09674,{{4,21,17,30,28}}}
*)

which shows that in-memory computation is about two orders of magnitude more memory-hungry, about 1.5x - 2x faster for all matches, and significantly slower for a single match search.

Note that run-time and memory efficiency of lazy computations in general do depend on the chunk size. One can specify it using "ChunkSize" option to LazyTuples. For example, this would create a LazyList with larger chunks:

LazyTuples[lists, "ChunkSize" ->  50000]

Previous editions of this post contain more in-depth explanations of the internals of LazyTuples.

Speeding things up

Using compiled predicate

We can compile the predicate in Select:

selFun = Compile[{{ints, _Integer, 1}}, Total[ints] == 100, "CompilationTarget" -> "C"]

which would give about 1.5x - 2x speedup (for both lazy and in-memory computations):

Normal @ Select[lt,  selFun];//AbsoluteTiming 
Normal @ Select[lt,  selFun]//MaxMemoryUsed
Normal @ Take[Select[lt,  selFun], 1]//AbsoluteTiming 

(*
   {3.28482,Null}

   4941080

   {0.075051,{{4,21,17,30,28}}}
*)

Compiling entire Select

This section will show how Streaming can actually beat in-memory computations not only by memory efficiency, but also by speed.

Let's start by asking, whether one can further speed things up. One obvious improvement is to try compiling entire Select, rather than just the predicate.

The following function is a generator of compiled Select functions, which take (preferably compiled) predicate and return compiled Select:

ClearAll[createTupleSelect]
createTupleSelect[criteria_] := createTupleSelect[criteria, All]
createTupleSelect[criteria_, numberOfResults_] :=
  Replace[
    If[
        numberOfResults === All, 
        Hold[Select[tuples, criteria]], 
        Hold[Select[tuples, criteria, numberOfResults]]
    ],
    Hold[select_] :> Compile @@ Hold[
        {{tuples, _Integer, 2}}, 
        select, 
        CompilationTarget->"C", CompilationOptions -> {"InlineCompiledFunctions"-> True}
    ]
]   

Let us use it to create two compiled functions, one that selects all matching elements, and the other that selects just the first one:

sel = createTupleSelect[selFun]
selFst = createTupleSelect[selFun, 1]

---

We can test these easily by using them for the in-memory computation:

sel[Tuples[lists]]//Short//AbsoluteTiming
selFst[Tuples[lists]]//Short//AbsoluteTiming

(*
  {0.413269,{{4,21,17,30,28},{4,21,21,28,26},{4,21,21,30,24},<<23636>>,{3,30,28,14,25},{3,30,28,24,15},{3,30,28,21,18}}}

  {0.283024,{{4,21,17,30,28}}}
*)

We see that this version is about 5-6 times faster than the one where only the predicate has been compiled.

But what is remarkable here, is that the time to get a single result is almost the same as the time to get all results. This is because, for the in-memory computation, one first has to generate all tuples, and in the compiled setup this operation takes the majority of time. This is where lazy computations potentially have an advantage. Let's see whether we actually can have it with Streaming.

---

To test Streaming in this regime, we will switch to larger chunks - that would reduce the top-level overhead of Streaming chunk manipulation:

ltLrg := LazyTuples[lists, "ChunkSize" ->  50000]

Unfortunately, the implementation of Streaming Select does not have an option to provide an entire select function for a chunk, rather than a predicate. I have created a gist where I added such capability. You can also look at it if you are interested in how Select is implemented in Streaming currently, since it is mostly exact same code.

We will import it now:

Import[StringJoin[
  "https://gist.githubusercontent.com/lshifr/",
  "2fa8ff7c950595d2ca32adc41b76bd7d/raw/LazyListSelect.m"
]]

which will make a new LazyListSelect[list, pred, entireSelect] function available (again, the only reason for StringJoin is the M SE editor bug).

We can now use it:

Normal @ LazyListSelect[ltLrg,  None, sel]//Short//AbsoluteTiming 
Normal @ LazyListSelect[ltLrg,  None, sel]//MaxMemoryUsed
Normal @ Take[LazyListSelect[ltLrg,  None, selFst],1]//AbsoluteTiming 

(* 
  {0.67473,{{4,21,17,30,28},{4,21,21,28,26},{4,21,21,30,24},{4,21,21,26,28},<<23635>>,{3,30,28,14,25},{3,30,28,24,15},{3,30,28,21,18}}}

   9266760

   {0.047092,{{4,21,17,30,28}}}
*)

What we see is that the search for all matches is about 1.5x - 2x times slower than the in-memory version, which is similar to what we have seen before. In this case, it is probably mostly due to the tuples generation overhead (even though Carl's function is extremely fast, it is still top level compared to kernel's specialized implementation in C). But of course this is still an order of magnitude more memory efficient.

But for a single match search, we're in for a pleasant finding. Not only is this an order of magnitude more memory efficient than direct built-in, but it is also an order of magnitude faster.

Of course, as I mentioned previously, this isn't really surprising, since lazy computations don't need to generate all tuples if the result can be found within e.g. the first few thousands. And of course one can write imperative code with loops, which would lead to similar or somewhat better performance. But Streaming has all that built in, and allows one to stay within the declarative functional paradigm, while providing decent performance in many such cases.

Summary

I have outlined how problems such as this can be treated within Streaming framework.

The main advantages of Streaming are built-in memory efficiency and laziness, which can automatically lead to efficient computations when only a subset of data is needed at the end, while allowing one to stay within the declarative / functional programming paradigm.

In some cases, Streaming can be not only vastly more memory-efficient than analogous in-memory computations with built-ins, performed in functional style, but also much faster. I have illustrated that in the section on further possible speedups, with entirely compiled Select.

In many cases, when using Streaming, one can control the memory / speed balance by adjusting the chunk size of the computation.

Here's an answer that seems pretty fast and small. It essentially just implements a simple depth-first search, with the current row of s being investigated being at index i, and the current index investigated in each row j being a[j]. It also handles your kmin case, with the default (omitted argument) being k.

It stores some found values, but not many—only one integer per row of s. In particular, it stores a temporary integer p[i] for each i in the range 1, Length[s] representing the partial sum down the current path to that row. One could rewrite this by re-computing the sum down the path each time, though.

It's a bit slower and bigger than your Outer function in the small case, but is faster and uses less memory in even slightly larger cases; consider s = Table[RandomInteger[10, RandomInteger[{3, 10}]], {20}], n = Total[Last /@ s], and compute MaxMemoryUsed[FindSum[s,n]] // AbsoluteTiming. For mine this gives {0.008944, 18672}, and for the Outer one it gives {0.426522, 2449544}. (It also seems that this basic approach apparently outperforms linear programming, which surprised me.)

You might be able to compile it, but you might need to prepare your input by padding each list so it's a uniform array. In that case, you would need to modify the code to get the right values for amax, and thus the right cycling for a. See update below.

I'll update this with a more thorough explanation. Note that using lists for a and p apparently is worse than using the definitions used below, performance-wise, according to tests.

Note also that due to what seems to be a bug causing a recursion limit to be reached with long lists s (I'll be investigating that separately), one needs to define defaults the old-fashioned way—not with :. Or, just play it fast and loose and don't pattern-match the arguments s and n beyond _.

FindSum[s : {{___Integer} ...}, n_Integer] := FindSum[s, n, 0]

FindSum[s : {{___Integer} ...}, n_Integer, kmin0_Integer] :=

 (* Note: consider inserting the MaybeSumQ check mentioned below. *)

 Module[{a, i = 1, k = Length[s], kmin = kmin0, amax = Length /@ s, 
   sum, backincrement, p},
  If[kmin0 == 0, kmin = k];
  Do[a[i] = 1; p[i] = 0, {i, k}];
  p[1] = First@First@s;
  backincrement[] := If[(a[i] = Mod[1 + a[i], amax[[i]], 1]) == 1,
    If[--i == 0, Throw[False], backincrement[]]];
  Catch[While[True,
     sum = p[i];
     Which[
      sum == n && i >= kmin,
      Throw[Table[s[[j, a[j]]], {j, i}]],
      sum >= n || i == k,
      backincrement[];
      p[i] = p[i - 1] + s[[i, a[i]]],
      True,
      i++; p[i] = sum + s[[i, a[i]]];
      ]];
   ]
  ]

There's also a very fast (but admittedly somewhat memory-using—though not as much as the representation of s itself, I don't think) function that can perform a simple check to see that the list at least doesn't fail outright. It does this by ensuring the sum of the minima of the lists doesn't preclude a sum of n by being too high, and that the sum of the maxima of the lists doesn't preclude a sum of n by being too low.

MaybeSumQ[s : {{___Integer} ...}, n_Integer] := MaybeSumQ[s, n, 0]

MaybeSumQ[s : {{___Integer} ...}, n_Integer, kmin_Integer] :=
 If[kmin == 0 || kmin == Length[s],
  (#[[1]] <= n && #[[2]] >= n) &@Total[MinMax /@ s, {1}],
  Total[Min /@ Take[s, kmin]] <= n && Total[Max /@ s] >= n]

Changing the definition of FindSum to If[MaybeSumQ[s, n, kmin0], <body>, False] can save some time on average if a solution is not known to exist and has a high enough chance of not being generatable from s. For a list s containing 10000 lists of length about 7, it took only about 0.005 seconds.

Update: compilation

I tried compiling it, using PadRight[s, Automatic, -1] to pad out the ragged array s with -1s, and using something like

amax = Table[
  With[{si = spadded[[i0]]}, 
   Replace[FirstPosition[si, -1, Length[si]], 
    Dispatch[{{x_} :> x - 1}]]], {i0, k}]

as well as trying handwritten While loops to check each list element one-by-one, and variations of the above with and without Withs and Replaces, and sometimes with /@s. This didn't seem to make much difference.

It compiles, but surprisingly—even excluding the time and memory it takes to pad s—the compiled version seems to run at the same speed, and in fact ever so slightly slower than the uncompiled one, no matter how I wrote amax. It was important to use RepeatedTiming here, since the differences were small and variable. However, it did use less memory, by about a third. To test easily I made functions to generate random lists like this:

(* Helpful functions for generating test cases: *)
Unprotect[rowlength, vrange];
ClearAll[rowlength, vrange]; 
Options[RandomList] = {rowlength -> {3, 10}, vrange -> {1, 10}}; Protect[rowlength, vrange];

RandomList[k_Integer, OptionsPattern[]] := 
 Table[RandomInteger[OptionValue[vrange], 
   RandomInteger[OptionValue[rowlength]]], {k}];

RandomSum[s : {{___Integer} ...}] := Total[RandomChoice /@ s];

s = RandomList[750]; n = RandomSum[s]; s0 = PadRight[s, Automatic, -1];

(One can also use n = Total[Last /@ s] to try to force a sum on a "hard" path.)

For the compiled function, I got that the MaxMemoryUsage[CFindSum0[s0, n, 0]] // RepeatedTiming was {0.3189, 288464}, and for the uncompiled function, I got {0.292063, 463032}

For completeness, here's the compiled version:

CFindSum0 = 
 Compile[{{spadded, _Integer, 2}, {n, _Integer}, {kmin0, _Integer}}, 
  Module[{a, i = 1, k = Length[spadded], kmin = kmin0, amax, j, l, 
    sum, backincrement, p},
   If[kmin0 == 0, kmin = k];
   amax = 
    Table[With[{si = spadded[[i0]]}, 
      Replace[FirstPosition[si, -1, Length[si]], 
       Dispatch[{{x_} :> x - 1}]]], {i0, k}];
   Do[a[i] = 1; p[i] = 0, {i, k}];
   p[1] = First@First@spadded;
   Catch[While[True,
      sum = p[i];
      Which[
       sum == n && i >= kmin,
       Throw[Table[spadded[[j, a[j]]], {j, i}]],
       sum >= n || i == k,
       Label["backincrement"]; 
       If[(a[i] = Mod[1 + a[i], amax[[i]], 1]) == 1,
        If[--i == 0, Throw[False], Goto["backincrement"]]];
       p[i] = p[i - 1] + spadded[[i, a[i]]],
       True,
       i++; p[i] = sum + spadded[[i, a[i]]];
       ]];
    ]
   ],
  CompilationTarget -> "C"]

CFindSum[s : {{___Integer}...}, n_Integer, kmin_Integer] :=
    CFindSum0[PadRight[s, Automatic, -1], n, kmin]

CFindSum[s : {{___Integer}...}, n_Integer] := CFindSum[s, n, 0]

CFindSum can be used with the same argument type and syntax as FindSum, since it packages in the padding and optional arguments.

Here is a recursive implementation. The recursion is stopped when the first result is found using Throw. For more results the recursion can be stopped later.

k is the number of sub lists and kmin is the minimal number of terms in the sum of n`:

s = {{1, 2, 3}, {4, 5, 6, 7}, {1, 3}, {2, 4, 6}};
n = 19; (*searched sum*)
k = Length[s];
kmin = k; (* minimal length of result*)

Clear[step];
step[ind_, sum_] := Module[{len = Length[ind]},
  Which[
   sum > n, Nothing,
   len >= kmin && sum == n, Throw[ Table[s[[i, ind[[i]]]], {i, len}]],
   len < k, 
   Table[step[Append[ind, i], sum + s[[len + 1, i]]], {i, 
     Length[s[[len + 1]]]}]
   , True, Nothing
    ]
  ]
Catch[step[{}, 0]]

(*{3, 7, 3, 6}*)