# Code to find the sum of the product of all size-n combinations of an array

The title is a bit of a mouthful, so here's a quick example:

You've got an array, maybe just the digits 1-3, and you want to get the sum of product of all size-2 combinations of the array. Those combinations are (1,2), (1,3), and (2,3), so the answer I'm looking for is this:

$$1\cdot2+1\cdot3+2\cdot3=11$$

Except now, the array is about 3-4 orders of magnitude in size, and the problem is kinda intractable beyond size-2 combinations.

import numpy as np
from itertools import combinations as comb

#A is the target array, in this case, a random array of 1000 elements.
#B is a iterator of length-2 tuples which should act as indices for reading A.

n = 2
A = np.random.rand(1000)
B = comb(range(1000),n)

print np.sum(np.array([np.prod(A[i,]) for i in B]))


The code is fairly condensed, but is kinda simple. It first makes an iterator object containing all the unique combinations of the indices of A, and then uses those tuples in an (ab)use of numpy's array slicing notation to get reference those indices, find the product there, and then sum it all together.

The greatest weakness of this code, I would guess, is the list comprehension, since it doesn't make use of numpy functions. As it stands, the code runs fast enough for size-2 combinations, but falls flat for anything higher. I would like it to at least be tractable for size-3 combinations, but as it stands, this code doesn't make the cut.

Any suggestions or things to look into would be greatly appreciated.

• This is a notoriously non-trivial problem. Start here – vnp Sep 4 '18 at 20:17

It first makes an iterator object containing all the unique combinations of the indices of A

Brute force is the first step in "Make it work; then make it fast". Since it's the bottleneck, you need to take the second step: exploit the structure of the problem. Here there's plenty of algebraic structure to exploit (as long as numerical stability isn't a major issue).

In particular, think about divide and conquer. Rather than compute all $\binom{S}{2} = \Theta(S^2)$ pairs separately, take a running total and multiply each element by each of the previous elements simultaneously:

result, running_total = 0, 0
for a in A:
result += running_total * a
running_total += a


and you have $O(S)$ time. This can be extended for larger values of $n$ to work in $O(Sn)$ time and $O(n)$ space.

• In the case of a size-3 combinations, wouldn't I still have to deal with $O(S^2)$ code? I'm not exactly sure how I would divide and conquer here - are you suggesting that I group the multiplications necessary for more efficiency? – rp0 Sep 4 '18 at 16:24
• The logic for pairs is that we can group the pairs by their last element, so we multiple the last element by the sum of one-element subsets drawn from the earlier ones. Exactly the same logic applies for triples: multiply the last element by the sum of two-element subsets drawn from the earlier ones. We're already calculating that sum of two-element subsets, although it might need renaming... – Peter Taylor Sep 4 '18 at 16:30
• That is much better! – Here is a generalization to combinations of n elements: gist.github.com/martinr448/ebc3e81eb394eace14df0d0846b9e5d5 (written mainly in order to understand the algorithm). – Martin R Sep 4 '18 at 18:06
• @MartinR Would you mind if I take and adapt this code? I'll need to understand and numpy-ify it, but this is exactly what I need. – rp0 Sep 4 '18 at 19:26
• @rp0: Not at all. It is just a generalization of Peter's excellent approach (and can probably be further improved). Feel free to use it as you like. – Martin R Sep 4 '18 at 19:29

Some suggestions to simplify the code:

• Instead of creating all possible combinations of indices and using these as subscript into the list, comb(A, n) can be used directly to get all n-element combinations.
• There is no need to create an intermediate list in order to pass it to np.sum().
• Variable names should be lowercase.
• Put the code into a function so that it can be easily tested and reused.

Generally you might consider to switch to Python 3, or at least make it Python 3 compatible (compare Porting Python 2 Code to Python 3).

Then it would look like this:

from __future__ import print_function
import numpy as np
from itertools import combinations as comb

def sum_of_products(a, n):
return np.sum(np.prod(x) for x in comb(a, n))

if __name__ == "__main__":
n = 2
a = np.random.rand(1000)
print(sum_of_products(a, n))


This is not necessarily faster. There are 166167000 ways to pick three elements out of 1000 – I am curious if someone comes up with a better approach to compute that in reasonable time!

For posterity, here's how this code can be done without a list comprehension:

import numpy as np
from itertools import combinations as comb

subsetSize = 4
elements = 100

targetArray = np.random.rand(elements)

arrayCombs = np.fromiter(comb(targetArray,subsetSize),'f,f,f')
arrayCombs = np.array((arrayCombs['f0'],arrayCombs['f1'],arrayCombs['f2']))

print np.sum(np.prod(arrayCombs,0))


I made use of the np.fromiter function, which creates an array based on an iterator. It only works to make a 1D array, though, so I had to make a data type that could hold n elements. Unfortunately, the result is still a 1D array, so further effort needs to be added to change the array to be something that can be easily broadcast to.

This isn't generalized to any n, and would have to take some fiddling to always work. It is mountains faster than the original code, but still suffers from a $O(n^3)$ growth rate. Modifying the algorithm itself will be necessary for any further improvements.