Fastest way for working with itertools.combinations I need to speed up the function below:

import numpy as np
import itertools

def combcol(myarr):
    ndims = myarr.shape[0]
    solutions = []
    for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(np.arange(ndims), 6):
        c1, c2, c3, c4, c5, c6 = myarr[idx1,1], myarr[idx2,2], myarr[idx3,1], myarr[idx4,2], myarr[idx5,1], myarr[idx6,2]
        if c1-c2&gt;0 and c2-c3&lt;0 and c3-c4&gt;0 and c4-c5&lt;0 and c5-c6&gt;0 :
            solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
    return solutions

X = np.random.random((20, 10))  
Y = np.random.random((40, 10))  

if __name__=='__main__':
    from timeit import Timer
    t = Timer(lambda : combcol(X))
    t1 = Timer(lambda : combcol(Y))
    print('t : ',t.timeit(number=1),'t1 : ',t1.timeit(number=1))


Results:

t :  0.6165180211451455 t1 :  64.49216925614847 


The algorithm is too slow for my standard use (myarr.shape[0] = 500).  Is there a NumPy way to decrease execution time of this function (without wasting too much memory)?  Is it possible to implement the problem in Cython?
Your algorithm is of time complexity O(n^6) in ndims. If your current code needs 65 s for ndims = 40, then one could estimate the time for ndims = 500 to be roughly 8 years. Even an improvement of the most inner block's execution time by a factor of 1000 would yield an execution time of about 3 days. I guess this is not good enough.

So you should better try to improve the algorithm itself. However I do not really understand the problem the function is trying to solve and I couldn't come up with a better solution on the spot.

Also I haven't looked at numpy approaches, because you seem to have no vectorizable calculations in there.

Nonetheless I tried to make improvements to the code as is. I came up with the following:


You use only the a small fraction of myarr, so extract these parts first

myarr1 = myarr[:,1]
myarr2 = myarr[:,2]

Since you do many index lookups on these arrays turn them into lists. Numpy is not good at single index lookups.

myarr1 = myarr[:,1].tolist()
myarr2 = myarr[:,2].tolist()

Instead of substracting and comparing to 0, compare the variables directly saving one operation:

if c1&gt;c2 and c2&lt;c3 and c3&gt;c4 and c4&lt;c5 and c5&gt;c6:

Instead of using and, you can also connect the comparisons all together, potentially saving the double variable lookups:

 if c1 &gt; c2 &lt; c3 &gt; c4 &lt; c5 &gt; c6:

Instead of assigning the array elements to c1, c2, etc., use them directly. Since later comparisons are not executed if earlier one already return False this might save array lookups:

if myarr1[idx1] &gt; myarr2[idx2] &lt; myarr1[idx3] &gt; myarr2[idx4] &lt; myarr1[idx5] &gt; myarr2[idx6]:



All in all I get the following code:

import numpy as np
import itertools

def combcol(myarr):
    ndims = myarr.shape[0]
    solutions = []
    for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(np.arange(ndims), 6):
        c1, c2, c3, c4, c5, c6 = myarr[idx1,1], myarr[idx2,2], myarr[idx3,1], myarr[idx4,2], myarr[idx5,1], myarr[idx6,2]
        if c1-c2&gt;0 and c2-c3&lt;0 and c3-c4&gt;0 and c4-c5&lt;0 and c5-c6&gt;0 :
            solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
    return solutions

def combcol2(myarr):
    ndims = myarr.shape[0]
    myarr1 = myarr[:,1].tolist()
    myarr2 = myarr[:,2].tolist()
    solutions = []
    for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(range(ndims), 6):
        if myarr1[idx1] &gt; myarr2[idx2] &lt; myarr1[idx3] &gt; myarr2[idx4] &lt; myarr1[idx5] &gt; myarr2[idx6]:
            solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(myarr1[idx1], myarr2[idx2], myarr1[idx3], myarr2[idx4], myarr1[idx5], myarr2[idx6])))
    return solutions

X = np.random.random((40, 10))
Y = [np.random.random((10, 10)) for i in range(100)]

if __name__=='__main__':
    from timeit import Timer
    for test_case in Y: assert combcol(test_case) == combcol2(test_case)
    t = Timer(lambda : combcol(X))
    t1 = Timer(lambda : combcol2(X))
    print('t : ',t.timeit(number=1),'t1 : ',t1.timeit(number=1))


With results like these:

t :  12.221544027328491 t1 :  0.8104081153869629
t :  10.69440507888794 t1 :  0.6592199802398682
t :  12.91178011894226 t1 :  0.8727250099182129
t :  11.21804690361023 t1 :  0.6851639747619629


The improvement is not even close to what you need. Also I tried to run the algorithm with ndims = 80 and got memory problems. The reason is the number of solutions: There were 18679356 solutions taking up 5 GB of RAM. With other values for ndims I found that around 5% - 25% of all possible solutions are valid! So it seems you need to change the way you output your solutions if you want to get beyond ndims = 100 on any realistic machine.

Edit: I just notice that 5 GB is a bit too much for the number of solutions. However this was the memory allocated by the python process. A run with ndims = 100 was unsuccessful due to lack of available memory.

Edit2: There is actually something far more efficient you can do, by interweaving your solution checks with the combination generation (abandoning itertools.combinations). This doesn't look so nice however:

import numpy as np
import itertools


def combcol(myarr):
    ndims = myarr.shape[0]
    solutions = []
    for idx1, idx2, idx3, idx4, idx5, idx6 in itertools.combinations(np.arange(ndims), 6):
        c1, c2, c3, c4, c5, c6 = myarr[idx1,1], myarr[idx2,2], myarr[idx3,1], myarr[idx4,2], myarr[idx5,1], myarr[idx6,2]
        if c1-c2&gt;0 and c2-c3&lt;0 and c3-c4&gt;0 and c4-c5&lt;0 and c5-c6&gt;0 :
            solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
    return solutions


def combcol2(myarr):
    ndims = myarr.shape[0]
    myarr1 = myarr[:,1].tolist()
    myarr2 = myarr[:,2].tolist()
    solutions = []

    for idx1 in range(ndims):
        c1 = myarr1[idx1]
        for idx2 in range(idx1+1, ndims):
            c2 = myarr2[idx2]
            if c1 &gt; c2:
                for idx3 in range(idx2+1, ndims):
                    c3 = myarr1[idx3]
                    if c2 &lt; c3:
                        for idx4 in range(idx3+1, ndims):
                            c4 = myarr2[idx4]
                            if c3 &gt; c4:
                                for idx5 in range(idx4+1, ndims):
                                    c5 = myarr1[idx5]
                                    if c4 &lt; c5:
                                        for idx6 in range(idx5+1, ndims):
                                            c6 = myarr2[idx6]
                                            if c5 &gt; c6:
                                                solutions.append(((idx1, idx2, idx3, idx4, idx5, idx6),(c1, c2, c3, c4, c5, c6)))
    return solutions


X = np.random.random((40, 10))
Y = [np.random.random((10, 10)) for i in range(100)]

if __name__=='__main__':
    from timeit import Timer
    for test_case in Y: assert combcol(test_case) == combcol2(test_case)
    t = Timer(lambda : combcol(X))
    t1 = Timer(lambda : combcol2(X))
    print('t : ',t.timeit(number=1),'t1 : ',t1.timeit(number=1))


The result:

t :  11.706523180007935 t1 :  0.15908312797546387
t :  12.531341075897217 t1 :  0.16385602951049805
t :  11.581043004989624 t1 :  0.14388418197631836


So with that we already have nearly a factor of 100x speed-up.