Count how many matrices have full rank for all submatrices
Question: I would like to count how many m by n matrices whose elements are 1 or -1 have
the property that all its `floor(m/2)+1 by n` submatrices have full rank. My
current method is naive and slow and is in the following python/numpy code. It
simply iterates over all matrices and tests all the submatrices.

    
    
    import numpy as np
    import itertools
    from scipy.misc import comb
    
    m = 8
    n = 4
    
    rowstochoose = int(np.floor(m/2)+1)
    
    maxnumber = comb(m, rowstochoose, exact = True)
    
    matrix_g=(np.array(x).reshape(m,n) for x in itertools.product([-1,1], repeat = m*n))
    
    nofound = 0
    for A in matrix_g:
        count = 0
        for rows in itertools.combinations(range(m), int(rowstochoose)):
           if (np.linalg.matrix_rank(A[list(rows)]) == int(min(n,rowstochoose))):
               count+=1
           else:
               break
        if (count == maxnumber):
             nofound+=1   
    print nofound, 2**(m*n)
    

Is there a better/faster way to do this? I would like to do this calculation
for n and m up to 20 but any significant improvements would be great.

**Context.** I am interested in getting some exact solutions for
<http://math.stackexchange.com/questions/640780/probability-that-every-vector-
is-not-orthogonal-to-half-of-the-others> .

* * *

As a data point to compare implementations. `n,m = 4,4` should output 26880 .
`n,m=5,5` is too slow for me to run. For `n =` 2 and `m = 2,3,4,5,6` the
outputs should be `8, 0, 96, 0, 1280`.

* * *

**Current status Feb 2, 2014:**

  * The answer of leewangzhong is fast but is not correct for m > n . leewangzhong is considering how to fix it.
  * The answer of Hooked does not run for m > n .


Answer: (Now a partial solution for n = m//2+1, and the requested code.)

Let k := m//2+1

This is somewhat equivalent to asking, "How many collections of m
n-dimensional vectors of {-1,1} have no linearly dependent sets of size
min(k,n)?"

For those matrices, we know or can assume:

  * The first entry of every vector is 1 (if not, multiply the whole by -1). This reduces the count by a factor of 2**m.
  * All vectors in the list are distinct (if not, any submatrix with two identical vectors has non-full rank). This eliminates a lot. There are choose(2**m,n) matrices of distinct vectors.
  * The list of vectors are sorted lexicographically (rank isn't affected by permutations). So we're really thinking about sets of vectors instead of lists. This reduces the count by a factor of m! (because we require distinctness).

With this, we have a solution for n=4, m=8. There are only eight different
vectors with the property that the first entry is positive. There is only one
combination (sorted list) of 8 distinct vectors from 8 distinct vectors.

    
    
    array([[ 1,  1,  1,  1],
           [ 1,  1,  1, -1],
           [ 1,  1, -1,  1],
           [ 1,  1, -1, -1],
           [ 1, -1,  1,  1],
           [ 1, -1,  1, -1],
           [ 1, -1, -1,  1],
           [ 1, -1, -1, -1]], dtype=int8)
    

100 size-4 combinations from this list have rank 3. So there are 0 matrices
with the property.

* * *

For a more general solution:

Note that there are `2**(n-1)` vectors with first coordinate -1, and
`choose(2**(n-1),m)` matrices to inspect. For n=8 and m=8, there are 128
vectors, and 1.4297027e+12 matrices. It might help to answer, "For i=1,...,k,
how many combinations have rank i?"

Alternatively, "What kind of matrices (with the above assumptions) have less
than full rank?" ~~And I think the answer is exactly,~~ A sufficient condition
is, "Two columns are multiples of each other". ~~I have a feeling that this is
true, and I tested this for all 4x4, 5x5, and 6x6 matrices.~~(Must've screwed
up the tests) Since the first column was chosen to be homogeneous, and since
all homogeneous vectors are multiples of each other, any submatrix of size k
with a homogeneous column other than the first column will have rank less than
k.

This is not a necessary condition, though. The following matrix is singular
(first plus fourth is equal to third plus second).

    
    
    array([[ 1,  1,  1,  1,  1],
           [ 1,  1,  1,  1, -1],
           [ 1,  1, -1, -1,  1],
           [ 1,  1, -1, -1, -1],
           [ 1, -1,  1, -1,  1]], dtype=int8)
    

Since there are only two possible values (-1 and 1), all mxn matrices where
`m>2, k := m//2+1, n = k` and with first column -1 have a majority member in
each column (i.e. at least k members are the same). So for n=k, the answer is
0.

* * *

For n<=8, here's code to generate the vectors.

    
    
    from numpy import unpackbits, arange, uint8, int8
    
    #all distinct n-length vectors from -1,1 with first entry -1
    def nvectors(n):
        if n > 8:
            raise ValueError #is that the right error?
        return -1 + 2 * (
            #explode binary numbers to arrays of 8 zeroes and ones
            unpackbits(arange(2**(n-1),dtype=uint8)) #unpackbits only takes uint
                .reshape((-1,8)) #unpackbits flattens, so we need to shape it to 8 bits
                [:,-n:] #only take the last n bytes
                .view(int8) #need signed
        )
    

Matrix generator:

    
    
    #generate all length-m matrices that are combinations of distinct n-vectors
    def matrix_g(n,m):
        return (array(mat) for mat in combinations(nvectors(n),m))
    

The following is a function to check that all submatrices of length maxrank
have full rank. It stops if any have less than maxrank, instead of checking
all combinations.

    
    
    rankof = np.linalg.matrix_rank
    #all submatrices of at least half size have maxrank
    #(we only need to check the maxrank-sized matrices)
    def halfrank(matrix,maxrank):
    return all(rankof(submatr) == maxrank for submatr in combinations(matrix,maxrank))
    

Generate all matrices that have all half-matrices with full rank def
nicematrices(m,n): maxrank = min(m//2+1,n) return (matr for matr in
matrix_g(n,m) if halfrank(matr,maxrank))

Putting it all together:

    
    
    import numpy as np
    from numpy import unpackbits, arange, uint8, int8, array
    from itertools import combinations
    
    #all distinct n-length vectors from -1,1 with first entry -1
    def nvectors(n):
        if n > 8:
            raise ValueError #is that the right error?
        if n==0:
            return array([])
        return -1 + 2 * (
            #explode binary numbers to arrays of 8 zeroes and ones
            unpackbits(arange(2**(n-1),dtype=uint8)) #unpackbits only takes uint
                .reshape((-1,8)) #unpackbits flattens, so we need to shape it to 8 bits
                [:,-n:] #only take the last n bytes
                .view(int8) #need signed
        )
    
    #generate all length-m matrices that are combinations of distinct n-vectors
    def matrix_g(n,m):
        return (array(mat) for mat in combinations(nvectors(n),m))
    
    rankof = np.linalg.matrix_rank
    #all submatrices of at least half size have maxrank
    #(we only need to check the maxrank-sized matrices)
    def halfrank(matrix,maxrank):
        return all(rankof(submatr) == maxrank for submatr in combinations(matrix,maxrank))
    
    #generate all matrices that have all half-matrices with full rank
    def nicematrices(m,n):
        maxrank = min(m//2+1,n)
        return (matr for matr in matrix_g(n,m) if halfrank(matr,maxrank))
    
    #returns (number of nice matrices, number of all matrices)
    def count_nicematrices(m,n):
        from math import factorial
        return (len(list(nicematrices(m,n)))*factorial(m)*2**m, 2**(m*n))
    
    for i in range(0,6):
        print (i, count_nicematrices(i,i))
    

`count_nicematrices(5,5)` takes about 15 seconds for me, the vast majority of
which is taken by the `matrix_rank` function.

