计算所有子矩阵的满秩矩阵数

Question

我想计算有多少个元素为1或-1的m×n矩阵具有其所有floor(m/2)+1 by n子矩阵具有满秩的性质。我当前的方法既幼稚又缓慢，使用的是下面的python/numpy代码。它只需迭代所有矩阵并测试所有子矩阵。在

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)

有更好/更快的方法吗？我想对n和m做这个计算，直到20，但任何重大的改进将是伟大的。在

上下文。我有兴趣得到https://math.stackexchange.com/questions/640780/probability-that-every-vector-is-not-orthogonal-to-half-of-the-others的精确解。在

---

作为比较实现的数据点。n,m = 4,4应该输出26880。n,m=5,5太慢，我无法运行。对于n =2和m = 2,3,4,5,6，输出应该是8, 0, 96, 0, 1280。在

---

现状2014年2月2日：

• 李旺忠的回答很快，但对m>；n并不正确。李旺忠正在考虑如何修复。在
• Hooked的答案不适用于m>；n。在

Answer 1

（现在是n=m//2+1的部分解决方案，以及请求的代码。）

设k：=m//2+1

这在某种程度上相当于这样一个问题：“有多少个{-1,1}的mn维向量集合没有大小为min（k，n）的线性相关集？”在

对于这些矩阵，我们知道或可以假设：

• 每个向量的第一个条目是1（如果不是，则将整数乘以-1）。这使计数减少了2**m
• 列表中的所有向量都是不同的（如果不是，则具有两个相同向量的任何子矩阵都具有非满秩）。这消除了很多。存在不同向量的选择矩阵（2**m，n）。在
• 向量列表按字典顺序排序（秩不受排列的影响）。所以我们考虑的是向量集而不是列表。这使计数减少了一倍！（因为我们要求清晰）。在

这样，我们得到了n=4，m=8的解。只有八个不同的向量具有第一个条目为正的属性。只有一个组合（排序列表）的8个不同的向量从8个不同的向量。在

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个4码组合的排名为3。所以有0个矩阵具有这个性质。在

---

对于更一般的解决方案：

注意，有第一个坐标为-1的2**(n-1)向量和choose(2**(n-1),m)矩阵要检查。对于n=8和m=8，有128个向量和1.4297027e+12个矩阵。回答“对于i=1，…，k，有多少个组合具有秩i？”在

或者，“什么样的矩阵（在上述假设下）的秩小于全秩？” 我认为答案是，一个充分条件是，“两列是彼此的倍数”。 我觉得这是真的，我对所有4x4、5x5和6x6矩阵进行了测试。（一定把测试搞砸了）因为第一列被选为齐次的，而且所有齐次向量都是彼此的倍数，任何大小为k的子矩阵的齐次列不是第一列的，其秩都将小于k。在

不过，这不是必要条件。下面的矩阵是奇异的（第一个加上第四个等于第三个加上第二个）。在

^{pr2}$

因为只有两个可能的值（-1和1），所有的mxn矩阵，其中m>2, k := m//2+1, n = k和第一列-1在每列中都有一个多数成员（即至少k个成员是相同的）。所以对于n=k，答案是0。在

---

对于n<；=8，下面是生成向量的代码。在

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
    )

矩阵发生器：

#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))

下面是一个函数，用于检查长度为maxrank的所有子矩阵是否具有满秩。如果任何一个值小于maxrank，它将停止，而不是检查所有组合。在

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))

生成所有具有全秩半矩阵的矩阵 定义尼斯矩阵（m，n）： 最大秩=min（m//2+1，n） return（matr for matr in matr in matr\g（n，m）if halfrank（matr，maxrank））

综合起来：

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)大约需要15秒，其中绝大多数是由{}函数占用的。在

Answer 2

你需要从数学的角度重新思考这个问题。也就是说，即使使用暴力，也有一些编程技巧可以用来加快进程（编程站点也是如此）。像不重新计算int(min(n,rowstochoose))和itertools.combinations(range(m), int(rowstochoose))这样的小技巧可以节省几个百分点——但真正的好处来自记忆。其他人也提到过，但我认为有一个完整、有效的代码示例可能会有用：

import numpy as np
from scipy.misc import comb
import itertools, hashlib

m,n = 4,4

rowstochoose = int(np.floor(m/2)+1)
maxnumber    = comb(m, rowstochoose, exact = True)
combo_itr    = (x for x in itertools.product([-1,1], repeat = m*n))
matrix_itr   = (np.array(x,dtype=np.int8).reshape((n,m)) for x in combo_itr)

sub_shapes = map(list,(itertools.combinations(range(m), int(rowstochoose))))
required_rank = int(min(n,rowstochoose))

memo = {}

no_found = 0
for A in matrix_itr:
    check = True
    for s in sub_shapes:
        view  = A[s].view(np.int8)
        h    = hashlib.sha1(view).hexdigest()

        if h not in memo: 
            memo[h] =  np.linalg.matrix_rank(view)

        if memo[h] != required_rank:
            check = False
            break
    if check: no_found+=1   

print no_found, 2**(m*n)

这使4x4的速度提高了近10倍-如果您愿意等待足够长的时间，您将看到更大矩阵的实质性改进。对于较大的矩阵，秩的成本按比例更昂贵，您可以在哈希运算中提前对矩阵排序：

^{pr2}$

对于4x4的情况，这会使情况稍差，但我希望5x5的情况会反过来。在

Answer 3

既然还没有人回答，这里有一个没有代码的答案。我看到的有用的对称性如下。在

1. 将一行乘以-1。在
2. 将一列乘以-1。在
3. 排列行。在
4. 排列列。在

我将通过exhaustively generating the non-isomorphs来解决这个问题，过滤它们，并求出它们的轨道大小。nauty对于第一步和第三步非常有用。假设大多数矩阵几乎没有对称性（对于nLarge来说这无疑是一个很好的假设，但是事先并不明显有多大），我预计8x8是可行的，9x9是临界的，10x10是遥不可及的。在

扩展伪代码：

1. 生成（m-1）乘（n-1）0-1矩阵的每个轨道的一个代表，由行和列置换生成的群所作用的矩阵，以及轨道的大小（=（m-1）！（n-1）！/自同构组的大小）。也许Tim linked论文的作者愿意分享他的代码；否则，请参见下文。

2. 对于每个矩阵，用（-1）^x替换条目x。添加一行和一列1s。将其轨道大小乘以2^（m+n-1）。这处理符号变化的对称性。

3. 过滤矩阵并求出剩余矩阵的轨道大小。您可以通过自己实现Gram-Schmidt来节省一点计算时间，这样当您按照字典顺序尝试所有组合时，就有机会为共享前缀重用部分结果。

无同构枚举：

McKay的模板可用于从m×n0-1矩阵的代表中生成（m+1）乘n0-1矩阵的代表，其方式适合于深度优先搜索。对于每个m×n0-1矩阵，将二分图与m个黑色顶点、n个白色顶点以及每个1个条目的相应边相关联。对每个m×n代表执行以下操作。在

1. 对于每一个长度为n的向量，构造由代表向量和新向量组成的（m+1）乘n矩阵的图，并运行nauty得到正则标号和顶点轨道。

2. 过滤掉与新向量对应的顶点与数量最少的黑色顶点在不同轨道上的可能性。

3. 用重复的规范标签过滤掉可能性。

nauty还计算自同构群的阶数。在