Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Count how many matrices have full rank for all submatrices

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 https://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 .
like image 651
marshall Avatar asked Jan 22 '14 10:01

marshall


1 Answers

(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.

like image 101
leewz Avatar answered Sep 30 '22 10:09

leewz