numpy to find products of all combinations of pairs of numbers in a matrix row

user_14405552 0
Viewed 287

I'm given an n x m matrix and my goal is to find the "cross-product" of all the features, specifically each row in the product matrix is of the form xijjxij', j < j', j = 1, . . . , m, j' = (j + 1), . . . , m so that the resulting row is a the product of all combinations of pairs in that row. Is there an elegant way to do this using numpy functions rather than python loops?

Edit: example

[1, 2, 3, 4]

should become

[1*2, 1*3, 1*4, 2*3, 2*4, 3*4]

which give:

[2, 3, 4, 6, 8, 12]

matrix numpy python
1 Answers

Basically you want to generate all the possible subsets containing 2 elements of your original set.

Short answer:

# With m = 4                            
c = np.multiply(*np.add(np.triu_indices(4,1),1))

General solution for any input array:

If using itertools is an option then you can use:

import numpy as np
import itertools

x = list(itertools.combinations([1,2,3,4], 2))
c = np.prod(x,-1)

c output:

array([ 2,  3,  4,  6,  8, 12])

From the doc:

itertools.combinations(iterables,r) : return r-length tuples in sorted order with no repeated elements.

And the number of elements in c correspond to the binomial coefficient C(n,k): n choose k, where n = len([1,2,3,4]) and k = 2.

Noticed that itertools.combinations() only hide the for loops, but since there is no closed-form formula for this problem a for loop is inevitable.

Numpy only solution:

In your specific case, where your iterable are the suit of n positive integers [1,2,3,4,...,n] then you can noticed that the positive indice of an upper triangle 2D matrice of length n-1 will produce the same result as combinations so:

# Number of elements in your array
n = 4
# Upper triangular matrice 
x = np.triu(np.ones([n-1,n-1]))
# Get the result
c = np.prod(np.argwhere(x)+np.arange(1,3),-1)

And again c output:

array([ 2,  3,  4,  6,  8, 12])

Or (with the help of @Nachikel, I wasn't aware of the existence of np.triu_indices()) the one liner:

c = np.multiply(*np.add(np.triu_indices(4,1),1))

Benchmarking:

enter image description here

and also with itertools:

enter image description here

The code:

import numpy as np
import itertools
import timeit
import matplotlib.pyplot as plt

def itertools1(m):
    x = list(itertools.combinations(np.arange(1,m+1), 2))
    np.prod(x,-1)
    
    
def numpy1(m):
    n = m-1
    x = np.triu(np.ones([n,n]))
    np.prod(np.argwhere(x)+np.arange(1,3),-1)
    
def numpy2(m):
    np.multiply(*np.add(np.triu_indices(m,1),1))
    



   
def benchmark_time(m): 
    SETUP_CODE = ''' 
from __main__ import numpy1 
from __main__ import numpy2
from __main__ import itertools1
'''
    x = np.zeros([3,len(m)])
    for ind, m in enumerate(m):
        print('For m = {}'.format(m))
        TEST_CODE = ''' 
itertools1({})
        '''.format(m)
        # timeit.repeat statement 
        times = timeit.repeat(setup = SETUP_CODE, 
                              stmt = TEST_CODE, 
                              repeat = 10, 
                              number = 50) 
        x[0,ind] = np.average(times)
        print('Itertools1 give:\t{} s'.format(np.round(np.average(times),3)))
        
        TEST_CODE = ''' 
numpy1({})
        '''.format(m)
        
        times = timeit.repeat(setup = SETUP_CODE, 
                              stmt = TEST_CODE, 
                              repeat = 10, 
                              number = 50) 
        x[1,ind] = np.average(times)
        print('Numpy1 give:\t\t{} s'.format(np.round(np.average(times),3)))
        
        TEST_CODE = ''' 
numpy2({})
        '''.format(m)
        
        times = timeit.repeat(setup = SETUP_CODE, 
                              stmt = TEST_CODE, 
                              repeat = 10, 
                              number = 50) 
        x[2,ind] = np.average(times)
        print('Numpy2 give:\t\t{} s\n'.format(np.round(np.average(times),3)))
    return x

m = np.arange(10,150,10)
x = benchmark_time(m)

plt.plot(m,x.T)
plt.legend(('itertools', 'numpy triu', 'numpy triu_indices'))
plt.xlabel('m')
plt.ylabel('sec')
plt.show()
user_43638640
Related