2

Suppose I have a numpy array

X = np.array([[1,2,3],
              [4,5,6],
              [7,8,9]])

I want to extend this matrix by adding (on the left) the columns resulting by multiplying together all possible pairs of columns. In this example it would become 

X = np.array([[1, 2, 3, 2, 6],
              [4,5,6,20,24,30],
              [7,8,9,56,63,72]])

Where the fourth column is the product of X[:,0] and X[:,1], the fifth column is the product of X[:,0] and X[:,2] and the sixth column is the product of X[:,1] and X[:,2].

My Try

I wanted to use np.hstack for this. However I also know that using loops slows everything down, but I wouldn't know how to do it properly without loops.

for i in range(matrix.shape[1]-1):
    for j in range(matrix.shape[1])[i:]:
        matrix2 = np.hstack((matrix, (matrix[:,i]*matrix[:,j]).reshape(-1,1))).copy()

The problem with this is that it is slow and also I have to use a different matrix, otherwise it will keep on adding columns.. any better idea?

Euler_Salter
  • 3,271
  • 8
  • 33
  • 74

1 Answers1

3

Approach #1

Get the pair-wise column indices with np.triu_indices. Use those to select two sets of blocks that are obtained from column-indexing into the input array. Use those blocks to perform element-wise multiplication and finally stack those as new columns alongside input array with np.concatenate.

Hence, the implementation -

n = X.shape[1]
r,c = np.triu_indices(n,1)
out0 = X[:,r] * X[:,c]
out = np.concatenate(( X, out0), axis=1)

Approach #2

For memory efficiency and hence performance especially with large arrays, one more inspired by this post by looping through groups of pairings -

m,n = X.shape
N = n*(n-1)//2
idx = np.concatenate(( [0], np.arange(n-1,0,-1).cumsum() ))+n
start, stop = idx[:-1], idx[1:]
out = np.empty((m,n+N),dtype=X.dtype)
out[:,:n] = X
for j,i in enumerate(range(n-1)):
    out[:, start[j]:stop[j]] = X[:,i,None]*X[:,i+1:]

Runtime test

In [403]: X = np.random.randint(0,9,(10,100))

In [404]: %timeit app1(X)
     ...: %timeit app2(X)
     ...: 
1000 loops, best of 3: 277 µs per loop
1000 loops, best of 3: 350 µs per loop

In [405]: X = np.random.randint(0,9,(10,1000))

In [406]: %timeit app1(X)
     ...: %timeit app2(X)
     ...: 
10 loops, best of 3: 68.6 ms per loop
100 loops, best of 3: 12.5 ms per loop

In [407]: X = np.random.randint(0,9,(10,2000))

In [408]: %timeit app1(X)
     ...: %timeit app2(X)
     ...: 
1 loop, best of 3: 311 ms per loop
10 loops, best of 3: 44.8 ms per loop
Divakar
  • 218,885
  • 19
  • 262
  • 358