(现在是对于n = m//2+1的部分解决方案和所需代码)
令k := m//2+1
这有点相当于问："有多少个集合包含m个n维{-1,1}向量，且其中没有大小为min(k,n)的线性相关组？"
对于那些矩阵，我们知道或可以假设:
- 每个向量的第一个条目都为1（如果不是，则将整个向量乘以-1）。这将使计数减少2 ** m倍。
- 列表中的所有向量都是唯一的（如果不是，则具有两个相同向量的任何子矩阵都具有非满秩）。这消除了很多问题。有choose（2 ** m，n）个不同向量的矩阵。
- 向量列表按字典顺序排序（排名不受排列影响）。因此，我们实际上是在思考向量而不是列表的集合。这将使计数减少m！倍（因为我们要求不同）。
有了这个，我们就有了n=4，m = 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。
对于更一般的解决方案：
请注意，有2 ** (n-1)个第一个坐标为-1的向量，以及choose(2 ** (n-1),m)个要检查的矩阵。对于n=8和m=8，有128个向量和1.4297027e + 12个矩阵。回答“对于i=1，...，k，有多少组合的排名为i？”可能会有所帮助。
或者，“什么样的矩阵（在上述假设下）具有不到满秩？”我认为答案是准确的，一个充分的条件是，“两列相互倍增”。 我感觉这是正确的，并且我已经测试了所有4x4、5x5和6x6矩阵。（可能测试出现问题了）由于选择了第一列为齐次列，并且所有的齐次向量都是相互倍增的，任何具有除第一列外其他齐次列的k大小的子矩阵的排名将小于k。
尽管如此，这并不是必要条件。下面的矩阵是奇异的（第一列加上第四列等于第三列加上第二列）。
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)
由于只有两个可能的值(-1和1)，所有当m>2, k := m//2+1, n = k且第一列为-1的mxn矩阵都有每列的多数成员(即至少有k个成员相同)。所以对于n=k，答案为0。
对于n<=8，以下是生成向量的代码。
from numpy import unpackbits, arange, uint8, int8
def nvectors(n):
if n > 8:
raise ValueError
return -1 + 2 * (
unpackbits(arange(2**(n-1),dtype=uint8))
.reshape((-1,8))
[:,-n:]
.view(int8)
)
矩阵生成器：
#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))
生成全部半矩阵均为满秩的矩阵
def nicematrices(m,n):
最大秩数 = min(m//2+1,n)
return (matr for matr in matrix_g(n,m) if halfrank(matr,最大秩数))
将所有内容组合起来:
import numpy as np
from numpy import unpackbits, arange, uint8, int8, array
from itertools import combinations
def nvectors(n):
if n > 8:
raise ValueError
if n==0:
return array([])
return -1 + 2 * (
unpackbits(arange(2**(n-1),dtype=uint8))
.reshape((-1,8))
[:,-n:]
.view(int8)
)
def matrix_g(n,m):
return (array(mat) for mat in combinations(nvectors(n),m))
rankof = np.linalg.matrix_rank
def halfrank(matrix,maxrank):
return all(rankof(submatr) == maxrank for submatr in combinations(matrix,maxrank))
def nicematrices(m,n):
maxrank = min(m//2+1,n)
return (matr for matr in matrix_g(n,m) if halfrank(matr,maxrank))
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秒钟，其中绝大部分时间都被 matrix_rank 函数占用。
min(n, rowstochoose)进行比较，一个JxK矩阵的最大可能排名是min(J, K)，我认为这就是你所说的“满秩”。 - Tim Petersmatrix_g生成了2^32个元素？ - leewz