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.

Here's the code I have made to do the task:

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.
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.