How can one generate a sequence of unique k-sparse matrices without rejection sampling for an arbitrary k efficiently?

Question

I would like to efficiently generate $k$-sparse matrices, i.e., matrices that have only $k$ nonzero entries. The catch is that all such matrices must be different. Matrices are filled with 1 at the $k$ nonzero entries and zero everywhere else. Unique means exactly identical.

• I can't randomly fill each new matrix at $k$ random indices, because if I sample a lot, there will be some indices that are not unique.
• And I don't want to have to check each time if I have already sampled a particular matrix, because this will get tedious towards the end, where almost all possible k-sparse matrices have been sampled already (it will take forever to get the last cases by random trial).
• I also don't want to precompute in an orderly fashion all the cases and then draw from this, because it will be way too large in memory for larger matrices.

Also, sometimes I only need a relatively short sequence. For the case of 1-sparse matrices it is simple to create a generator function that efficiently only returns unique matrices. But I need to go up to $k=10$ at least. Dimensions of the matrices can vary from $10\times 10$ to $100,000 \times 100,000$.

Answer 1

Assume you want a uniform random sample from the set of $m\times n$ matrices with $k$ entries with $1$, and $mn-k$ entries with $0$, without replacement. With $N=mn$, the problem is equivalent to drawing uniformly at random without replacement from the set of $N$-dimensional vectors with $k$ coordinates which are $1$, and $N-k$ coordinates which are $0$.

Conversion from this linear indexing into row and column indices of a matrix can be done with Euclidean division by the number of rows or by the number of columns. The quotient and remainder give the column and row indices, and vice versa, respectively. Matlab even allows you to index into a matrix with linear indexing, but also provides the routine ind2sub, and the inverse sub2ind. In NumPy for python, all you need is the numpy.reshape routine applied to the vector. But you can convert linear indexing into matrix coordinates with numpy.unravel_index. Its inverse is numpy.ravel_multi_index.

And that problem is equivalent to drawing uniformly at random without replacement $k$-element subsets of an $N$-element set.

1. Generate a random sample of the required size from $\left\{1,2,\dots,{{N}\choose{k}}\right\}$ with $N=mn$ without replacement.

   • For this, call the appropriate sampling function, or possibly the random permutation of $\left\{1,2,\dots,{{N}\choose{k}}\right\}$ and select the required number of elements from the beginning. However, the random permutation can run out of memory if ${{N}\choose{k}}$ is too large.

   • Alternatively, draw random samples from $\left\{1,2,\dots,{{N}\choose{k}}\right\}$ and reject proposals if they have already been drawn. Using the right data structure (binary search tree?), the comparison of proposals with already drawn samples is fast.

2. To each of these numbers, associate a $k$-element subset of the $N$-element set using the combinatorial number system. The combinatorial number system is a one-to-one mapping between $\left\{1,2,\dots,{{N}\choose{k}}\right\}$ and $k$-element subsets of the $N$-element set (their binary encoding), based on lexicographic ordering. It is a coordinate system that allows you to draw (and reject) simple nonnegative integers that are used for the indexing.

Here is a python implementation.

import numpy as np
import random, itertools
import scipy.special    # https://docs.scipy.org/doc/scipy/reference/special.html
import sortedcontainers        # http://www.grantjenks.com/docs/sortedcontainers/

nrows=3
ncols=nrows
k=2     # number of 1s
sample_size=5

def generate_matrix(nrows, ncols, k, idx):
    # Use the combinatorial number system to map from idx \in {1,2,..., N_choose_k} to corresponding N-dimensional binary vector with k 1s.
    # https://en.wikipedia.org/wiki/Combinatorial_number_system#Finding_the_k-combination_for_a_given_number
    mat=np.zeros((nrows*ncols,), dtype=np.uint8)
    for i in range(k,0,-1): # i means how many entries still to set to 1
        j=i-1
        found=False
        while found==False:
            testvalue=scipy.special.comb(j, i, exact=True) # Its first value is 0==((i-1) choose i)
            if testvalue==idx:
                found=True
                mat[j]=1
                idx=idx-testvalue
            elif testvalue>idx:
                found=True
                mat[j-1]=1
                idx=idx-testvalue_prev
            else:
                j=j+1
                testvalue_prev=testvalue
    return mat.reshape((nrows,ncols))

assert nrows*ncols >= k

Nchoosek = scipy.special.comb(nrows*ncols, k, exact=True)
assert Nchoosek >= sample_size

random.seed(1000)

# Option A
samples=random.sample(range(Nchoosek), sample_size) # This gives "OverflowError: Python int too large to convert to C ssize_t" already for 100x100 with k=50: Nchoosek == 2.9e+135.

'''
# Option B
samples=sortedcontainers.SortedList()
while len(samples)<sample_size:
    proposed=random.randrange(Nchoosek)
    if not proposed in samples:
        samples.add(proposed)
'''
for i in samples:
    mat = generate_matrix(nrows, ncols, k, i)
    print('\n',mat)

The weak point of my code is finding j in generate_matrix(). This code increments j one by one until ${{j}\choose{i}}$ first reaches idx. There should be a better way to guess and home on what the largest $j$ is such that ${{j}\choose{i}}\le I$ with given $i,I$.

For instance, this can be achieved by observing $${{j}\choose{i}} = \frac{j!}{i!(j-i)!} \le \frac{j^i}{i!} = \frac{j^i}{\Gamma(i+1)}.$$ Define $I=$idx. If you now rearrange $\frac{j^i}{\Gamma(i+1)} = I$, then you get $$j=\exp\left(\frac{\log I + \log \Gamma(i+1)}{i}\right),$$ which is a much better starting point for incrementing $j$.

Therefore I recommend to include

import math

and replace the line j=i-1 with

j=int(math.floor((1-1e-10)*math.exp((math.log(idx)+scipy.special.gammaln(i+1))/i)))

I included (1-1e-10)* to compensate for possible rounding errors but I haven't tested this much so it might be too tight and you should use e.g. (1-1e-9)*.

With this code I manage to generate a $100,\!000\times 100,\!000$ matrix with $k=500$ entries of $1$ in 18 sec.

You could also do a binary search for j. The other side of the estimate is $${{j}\choose{i}} = \frac{j!}{i!(j-i)!} \ge \frac{(j-i)^i}{i!} = \frac{(j-i)^i}{\Gamma(i+1)},$$ giving the upper end point for j: $$j=i+\exp\left(\frac{\log I + \log \Gamma(i+1)}{i}\right).$$

Answer 2

First, I would turn the $n\times n$ matrices into vectors of length $n^2$ by linear indexing because the tabular arrangement is irrelevant here. (Once you have got a solution for vectors, you can go back to matrix indexing by using Euclidean division of the linear index by the number of rows in the matrix: the quotient is the column index, the remainder is the row index.)

A trivial deterministic solution to the question is to use lexicographic ordering and (e.g. for $k=5$) to go through
$000000011111$
$000000101111$
$000000110111$
etc.
where each row is one matrix with matrix entries in linear indexing.

From now on, I assume you actually want a random sample, and specifically, one from the uniform distribution. I haven't got an exact solution but I suggest an approximation to the uniform which is random and should be good enough until a better is found.

It's based on the totally asymmetric simple exclusion process (TASEP) with $k$ balls. Start for $k=4$ with:
$111100000000\dots$
and in each step, one ball ($1$) jumps right if the right neighbour is empty ($0$)‎; then that $0$ becomes $1$, and the $1$ becomes $0$. (I'm afraid there is rejection sampling in each step.) Here is a realisation:
$111100000000\dots$
$111010000000\dots$
$110110000000\dots$
$110101000000\dots$

‎If you want to make this 'more random', then start the whole thing with uniformly randomly permuting the linear indexing of the matrix and fixing the indexing that way. So one random permutation and then steps of TASEP.
If the random permutation were too costly, then it may be substituted by a sampling without replacement of many elements. (Depending on the relationship between $n$, $k$ and the number of samples needed, the TASEP may not need to go very far from its starting position, and you may not need to reach far into the permuted, vectorised matrix.)

Other improvements: use burn-in for the TASEP, i.e. discard the first many steps to achieve that the balls are already spread out when you start sampling. And possibly discard many steps inbetween sampling two matrices. This helps to reduce dependence between consecutive matrices.