寻找更好的矩阵计数方法

8

我想统计只有1和0条目的二维数组数量,其中有一对不相交的行具有相等的向量和。 对于4x4矩阵,以下代码通过迭代所有矩阵并依次测试每个矩阵来实现此目的。

import numpy as np
from itertools import combinations
n = 4
nxn = np.arange(n*n).reshape(n, -1)
count = 0
for i in xrange(2**(n*n)):
   A = (i >> nxn) %2
   p = 1
   for firstpair in combinations(range(n), 2):
       for secondpair in combinations(range(n), 2):
           if firstpair < secondpair and not set(firstpair) & set(secondpair):
              if (np.array_equal(A[firstpair[0]] + A[firstpair[1]], A[secondpair[0]] + A[secondpair[1]] )):
                  if (p):
                      count +=1
                      p = 0
print count

输出结果为3136。

问题在于它使用了2^(4^2)次迭代,而我想要运行它直到n等于8。有没有更聪明的方法来计算这些,而不是遍历所有矩阵?例如,反复创建相同矩阵的排列似乎毫无意义。


你想要“两个不相交的行对,它们的向量和相等”,还是这是一个打字错误,你的意思是“两个不相交的行对,它们的向量和相等”? - Gastón Bengolea
矩阵是任何2D矩阵或类似于您在那里创建的那个吗? - Gastón Bengolea
我正在考虑所有元素为1或0的n×n矩阵。我的代码仅针对n = 4进行迭代。我希望它适用于n = 5、6、7、8,但目前运行时间太长。 - marshall
在今天的硬件上,为每个8x8矩阵执行单个CPU周期大约需要10000年;正如我们之前发现的那样,您需要数百万个周期来检查一个矩阵。因此,您确实需要想出比尝试所有这些更好的方法。即使进行行/列置换取模,我的直觉告诉我您仍然会遇到很大的麻烦。您可以通过随机抽样所有矩阵空间来获得这些矩阵的总数的良好估计吗? - Eelco Hoogendoorn
例如,100,000个8x8矩阵中有2276个测试结果为阳性,可以在十几秒内计算出来。这种统计数据至少足以描述该属性随着矩阵大小的变化而表现出的行为。(编辑:这是针对三重求和的;对于双行求和,它是7255/100,000) - Eelco Hoogendoorn
显示剩余9条评论
3个回答

8

使用CPython 3.3在我的电脑上计算大约花费一分钟:

4 3136
5 3053312
6 7247819776
7 53875134036992
8 1372451668676509696

基于记忆化的容斥原理代码:

#!/usr/bin/env python3
import collections
import itertools

def pairs_of_pairs(n):
    for (i, j, k, m) in itertools.combinations(range(n), 4):
        (yield ((i, j), (k, m)))
        (yield ((i, k), (j, m)))
        (yield ((i, m), (j, k)))

def columns(n):
    return itertools.product(range(2), repeat=n)

def satisfied(pair_of_pairs, column):
    ((i, j), (k, m)) = pair_of_pairs
    return ((column[i] + column[j]) == (column[k] + column[m]))

def pop_count(valid_columns):
    return bin(valid_columns).count('1')

def main(n):
    pairs_of_pairs_n = list(pairs_of_pairs(n))
    columns_n = list(columns(n))
    universe = ((1 << len(columns_n)) - 1)
    counter = collections.defaultdict(int)
    counter[universe] = (- 1)
    for pair_of_pairs in pairs_of_pairs_n:
        mask = 0
        for (i, column) in enumerate(columns_n):
            mask |= (int(satisfied(pair_of_pairs, column)) << i)
        for (valid_columns, count) in list(counter.items()):
            counter[(valid_columns & mask)] -= count
    counter[universe] += 1
    return sum(((count * (pop_count(valid_columns) ** n)) for (valid_columns, count) in counter.items()))
if (__name__ == '__main__'):
    for n in range(4, 9):
        print(n, main(n))

如果您对n>8的值感兴趣,仍然有相当数量的对称性可以通过合并在行置换下同构的“valid_columns”项来利用。 - David Eisenstat
聪明;我曾经考虑过这样的解决方案,但卡住了。 - Eelco Hoogendoorn
非常好。通过修改您的代码以处理N行x M列的情况,我认为我已经得到了(经验性的)生成函数,用于N = 4,5,6,7的数据。如果我能破解模式,我至少可以推测出更高的N,M的值。 - DSM
太棒了!这很流畅 :-) - Tim Peters
谢谢!虽然这已经超出了 SO 的要求,但你能否用数学术语解释一下你的代码正在计算什么?例如,它可以写成一个简单的总和吗?另外,你知道它的时间复杂度是多少吗? - marshall
@marshall 它是通过容斥原理计算并集的基数,其中并集中的每个集合都是满足特定行方程(以及可能的其他方程)的二进制矩阵集合。满足一组特定行方程的n乘m二进制矩阵的数量是满足该集合的列数的m次幂。我不知道这个算法的复杂度;它运行得比2 **(2 ** n)快,只要有比那更少的不同交集。 - David Eisenstat

3
您可以将此归类为“好过没有”;-) 这是一个纯Python3代码,它重新思考了这个问题。也许numpy技巧可以极大地加速它,但难以看出如何做到这一点。
  1. 这里的“一行”是range(2**n)中的整数。因此,数组只是由整数组成的元组。
  2. 因此,通过combinations_with_replacement()轻松生成所有在行置换下唯一的数组。这将外部循环的旅行次数从2**(n**2)减少到(2**n+n-1)-choose-n)。这是一个巨大的减少,但仍然...
  3. 预计算的字典将行对(这意味着这里是整数对!)映射到它们的向量和作为元组。因此,在测试时不需要进行任何数组操作,除了测试元组是否相等。通过一些更巧妙的方法,元组可以编码为(比如)基于3的整数,将内层循环测试减少到从一对字典查找检索的两个整数之间进行比较。
  4. 用于预计算的时间和空间是相对微不足道的,因此没有尝试加速该部分。
  5. 内循环每次选择4个行索引,而不是您的两个循环每次选择两个索引对。一次性做完所有4个操作更快,这在很大程度上是因为没有必要除去带有重复索引的对。
以下是代码:
def calc_row_pairs(n):
    fmt = "0%db" % n
    rowpair2sum = dict()
    for i in range(2**n):
        row1 = list(map(int, format(i, fmt)))
        for j in range(2**n):
            row2 = map(int, format(j, fmt))
            total = tuple(a+b for a, b in zip(row1, row2))
            rowpair2sum[i, j] = total
    return rowpair2sum

def multinomial(n, ks):
    from math import factorial as f
    assert n == sum(ks)
    result = f(n)
    for k in ks:
        result //= f(k)
    return result

def count(n):
    from itertools import combinations_with_replacement as cwr
    from itertools import combinations
    from collections import Counter
    rowpair2sum = calc_row_pairs(n)
    total = 0
    class NextPlease(Exception):
        pass
    for a in cwr(range(2**n), n):
        try:
            for ix in combinations(range(n), 4):
                for ix1, ix2, ix3, ix4 in (
                       ix,
                       (ix[0], ix[2], ix[1], ix[3]),
                       (ix[0], ix[3], ix[1], ix[2])):
                    if rowpair2sum[a[ix1], a[ix2]] == \
                       rowpair2sum[a[ix3], a[ix4]]:
                        total += multinomial(n, Counter(a).values())
                        raise NextPlease
        except NextPlease:
            pass
    return total

这就足以通过n = 6找到结果了,虽然最后一个结果需要很长时间才能完成(具体多久不知道-没有计时 - 大约需要一小时左右,“很长时间”是相对的;-)):

>>> count(4)
3136
>>> count(5)
3053312
>>> count(6)
7247819776

编辑 - 删除了一些不必要的索引

通过将主函数更改为以下内容,可以获得较快的速度:

def count(n):
    from itertools import combinations_with_replacement as cwr
    from itertools import combinations
    from collections import Counter
    rowpair2sum = calc_row_pairs(n)
    total = 0
    for a in cwr(range(2**n), n):
        for r0, r1, r2, r3 in combinations(a, 4):
            if rowpair2sum[r0, r1] == rowpair2sum[r2, r3] or \
               rowpair2sum[r0, r2] == rowpair2sum[r1, r3] or \
               rowpair2sum[r0, r3] == rowpair2sum[r1, r2]:
                total += multinomial(n, Counter(a).values())
                break
    return total

编辑 - 加快求和测试速度

虽然这只是个小问题,但由于目前似乎这是最精确的方法之一,因此可以尝试提高一些性能。 正如之前所指出的,由于每个和都在range(3)中,因此每个和的元组可以用一个整数代替(将元组视为给出基于3的整数的数字)。 将calc_row_pairs() 替换为以下内容:

def calc_row_pairs(n):
    fmt = "0%db" % n
    rowpair2sum = dict()
    for i in range(2**n):
        row1 = list(map(int, format(i, fmt)))
        for j in range(2**n):
            row2 = map(int, format(j, fmt))
            total = 0
            for a, b in zip(row1, row2):
                t = a+b
                assert 0 <= t <= 2
                total = total * 3 + t
            rowpair2sum[i, j] = total
    return rowpair2sum

我相信numpy一定有更快的方法来完成这个任务,但是calc_row_pairs()花费的时间微不足道,所以为什么要费心呢?顺便说一句,这样做的好处是内部循环中的==测试由比较元组变成了比较小整数。这对于普通的Python来说是有益的,但我敢打赌pypy能够获得更多的好处。


出于兴趣,pypy能加速这个吗? - marshall
不了解pypy - 目前还没有任何版本的pypy与Python3兼容。而且,虽然可能看起来很奇怪,但是当我处理“困难”问题时,我坚持使用纯标准的Python(没有pypy,没有numpy,没有外部包)。为什么?因为它非常慢,这迫使你考虑更聪明的方法 :-) - Tim Peters
实际上,那对我来说很有道理 :) - marshall
1
只是提醒一下,你的代码似乎也是在Python 2.x中编写的,并且可以在pypy上正常运行,n = 6时需要17分钟。 - marshall
1
@TimPeters 哇,速度提升得不错。Marshall 验证了 n=5 或 n=6 的答案吗?顺便说一句,我可能会看看 numpy 或 numba 的解决方案。无法想象使用“高级索引”不会加速。这只是说说而已,除非我发布一个,所以干得好! - Phil Cooper
@PhilCooper,在最后一次编辑之后,rowpair2sum 字典现在以整数对为索引并返回一个整数 - 它非常适合用于 2D numpy 数组。但是,巨大的速度提升需要不同的方法。我试图将其视为一个数学包含-排除问题,但我的头总是被扭成了一个结;-) - Tim Peters

2
并非直接回答你的问题,但正如我所指出的,我认为您可以安全地忘记对任何重要的n测试所有矩阵。但这个问题很适合使用随机特征描述。有趣的是,在某些条件下,三倍和的情况比双倍和更常见!虽然击中的可能性似乎是n和m的一个相当简单(单调)函数,但这并不令人惊讶。“double “triple

你能稍微解释一下你的图片吗?另外,当你说“在某些条件下,三重求和比双重求和更常见”时,这似乎令人惊讶。是什么样的条件呢? - marshall
某种程度上图像标题未显示。顶部图片为双重求和,底部图片为三重求和,每个单元中有10,000个样本。垂直轴为n,水平轴为m,范围从4到20(含)。红色表示100%点击率,蓝色表示0%点击率。在最后一行(n=20)中,我们可能会看到三重求和比双重求和更多地出现;至少对于适度的m而言是如此。随着m的增大,三重求和似乎更快地下降,因此这可能并不适用于所有的m。 - Eelco Hoogendoorn
在n=300,m=20周围的区域内,您几乎肯定会有三重碰撞,并且肯定不会有双重碰撞;而这种对比可以被推到任意高。相反的情况似乎在参数空间中没有出现过。当n<=m时,双重碰撞更常见,但它们从未以同样的方式占主导地位。 - Eelco Hoogendoorn

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接