I can't think of an elegant solution but here's a not-so-elegant pure numpy one:
import numpy as np
def combination_matrices(K):
# get combination indices
i, j = np.indices((K, K))
comb_indices = np.transpose((i < j).nonzero()) # (num_combs, 2) array where ones are
num_combs = comb_indices.shape[0] # K*(K-1)/2
# create a matrix of the desired shape, first axis enumerates combinations
matrices = np.zeros((num_combs, K, 2), dtype=int)
# broadcasting assignment of ones
comb_range, col_index = np.ogrid[:num_combs, :2]
matrices[comb_range, comb_indices, col_index] = 1
return matrices
This first uses the indices of a (K, K)
-shaped array to find the index pairs for every combination (these are indices that encode the upper triangle of the array, excluding the diagonal). Then we use a bit tricky broadcasting assignment (heavy fancy indexing) to set each corresponding element of the pre-allocated output array to 1.
Note that I put the K*(K-1)/2
-sized axis first, because this makes the most sense in numpy with C-contiguous memory layout. This way when you take the matrix for combination index 3
, arr[3, ...]
will be a contiguous chunk of memory of shape (K, 2)
that's fast to work with in vectorised operations.
The output for K = 4
:
[[[1 0]
[0 1]
[0 0]
[0 0]]
[[1 0]
[0 0]
[0 1]
[0 0]]
[[1 0]
[0 0]
[0 0]
[0 1]]
[[0 0]
[1 0]
[0 1]
[0 0]]
[[0 0]
[1 0]
[0 0]
[0 1]]
[[0 0]
[0 0]
[1 0]
[0 1]]]