python - 如何将特定元素的数量放入矩阵的特定行和列约束中？

Question

我基本上必须创建一个 NxN 的对称矩阵。其中随机填充了 1 和 0。然而，唯一的限制是我只需要在任何行和任何列中有一个“1”。

我写了一个代码来生成矩阵，但它在任何行或列中都有多个“1”。我需要遵循上面提到的约束，我该如何修改我的代码？

import numpy as np
N = int(input("Enter the number of row and col:"))
my_matrix = np.random.randint(2,size=(N,N))
print(my_matrix)

Answer 1

长话短说

每个结果均等概率生成并以O(n)时间复杂度运行:

import random

_prob_cache = [1, 1]


def prob(n):
    try:
        return _prob_cache[n]
    except IndexError:
        pass

    for i in range(len(_prob_cache) - 1, n):
        _prob_cache.append(1 / (i * _prob_cache[-1] + 1))

    return _prob_cache[-1]


def symmetric_permutation(n):
    res = np.zeros((n, n), int)
    remain = list(range(n))

    while remain:
        m = len(remain)
        diag_prob = prob(m)
        row = remain.pop()
        rnd = random.random()
        if rnd < diag_prob:
            col = row
        else:
            nondiag_prob = (1 - diag_prob) / (m - 1)
            idx = int((rnd - diag_prob) / nondiag_prob)
            remain[idx], remain[-1] = remain[-1], remain[idx]
            col = remain.pop()
        res[row, col] = res[col, row] = 1

    return res

---

长答案

从一些推导开始:

设f(n)为n * n矩阵所有设置方案的个数。显然，我们有:

f(1) = 1

然后约定一下:

f(0) = 1

对于 n > 1，我可以从任意行中提取一个位置并将其设置为 1。有两种情况:

1. 如果1在对角线上，我们可以去掉这个1的行和列，在剩下的(n - 1) * (n - 1)矩阵上继续设置，所以数剩余的设置方案是 f(n - 1)。
2. 如果1不在对角线上，那么对称部分也需要置1。那么我们就可以去掉这两个1所在的行和列。我们需要继续设置剩余的(n - 2) * (n - 2)矩阵。因此，剩余设置方案的数量为f(n - 2)。

因此，我们可以推断:

f(n) = f(n - 1) + (n - 1) * f(n - 2)

根据上述策略，如果我们想让每一个设置方案都以相同的概率出现，那么我们在选择索引时应该对对角线索引和其他索引赋予不同的权重。对角线索引的权重应该是:

p(n) = f(n - 1) / f(n)

因此:

   f(n) = f(n - 1) + (n - 1) * f(n - 2)
     f(n)         (n - 1) * f(n - 2)
=> -------- = 1 + ------------------
   f(n - 1)           f(n - 1)
     1
=> ---- = 1 + (n - 1) * p(n - 1)
   p(n)
                  1
=> p(n) = ------------------
          (n - 1) * p(n - 1)

概率函数代码如下:

_prob_cache = [1, 1]


def prob(n):
    """
    Iterative version to prevent stack overflow caused by recursion.
    Old version:
    @lru_cache
    def prob(n):
        if n == 1:
            return 1
        else:
            return 1 / ((n - 1) * prob(n - 1) + 1)
    """
    try:
        return _prob_cache[n]
    except IndexError:
        pass

    for i in range(len(_cache) - 1, n):
        _prob_cache.append(1 / (i * _prob_cache[-1] + 1))

    return _prob_cache[-1]

非对角线索引的权重为:

f(n - 2)   f(n - 2)   f(n - 1)
-------- = -------- * -------- = p(n - 1) * p(n)
  f(n)     f(n - 1)     f(n)

or

f(n - 2)   1 - p(n)
-------- = --------
  f(n)      n - 1

_{这里我选择使用后者来少调用一次函数。}

具体实现:

我们使用一个列表来存储仍然可以使用的索引。在每次循环中，我们将列表的最后一个元素作为行索引(不像前面说的选择第一个元素，这样可以加快从列表中移除元素的速度)，计算两种情况的权重并得到列随机索引，设置相应位置的值，并从列表中删除使用的索引，直到列表为空:

import random
import numpy as np


def symmetric_permutation(n):
    res = np.zeros((n, n), int)
    remain = list(range(n))

    while remain:
        m = len(remain)
        diag_prob = prob(m)
        row = remain.pop()
        rnd = random.random()
        if rnd < diag_prob:
            col = row
        else:
            nondiag_prob = (1 - diag_prob) / (m - 1)
            col = remain.pop(int((rnd - diag_prob) / nondiag_prob))
        res[row, col] = res[col, row] = 1

    return res

优化到 O(n) 时间复杂度:

如果我们不考虑零矩阵的创建，上述策略的时间复杂度是O(n^2)，因为每次我们都有很高的概率从中移除一个索引列表。

但是，暴力移除是不必要的。我们对其余索引的顺序没有要求，因为行索引的选择不会影响列索引的随机性。因此，一种更便宜的解决方案是用最后一个元素覆盖选定的列索引，然后删除最后一个元素。这使得去除中间元素的O(n)操作变成了O(1)操作，所以时间复杂度变成了O(n):

def symmetric_permutation(n):
    res = np.zeros((n, n), int)
    remain = list(range(n))

    while remain:
        m = len(remain)
        diag_prob = prob(m)
        row = remain.pop()
        rnd = random.random()
        if rnd < diag_prob:
            col = row
        else:
            nondiag_prob = (1 - diag_prob) / (m - 1)
            idx = int((rnd - diag_prob) / nondiag_prob)
            remain[idx], remain[-1] = remain[-1], remain[idx]
            col = remain.pop()
        res[row, col] = res[col, row] = 1

    return res

概率测试:

这里我们准备了另一个函数来计算 f(n) 用于以下测试:

def f(n):
    before_prev, prev = 1, 1
    for i in range(1, n):
        before_prev, prev = prev, prev + before_prev * i

    return prev

接下来是概率检验，验证结果是否足够均匀。这里我取n=8构建矩阵500_000次，每行1的列索引作为每条结果的标识，绘制了折线图和直方图每个结果出现的次数:

from collections import Counter

import matplotlib.pyplot as plt


random.seed(0)

n = 8
times = 500_000
n_bin = 30

cntr = Counter()
cntr.update(tuple(symmetric_permutation(n).nonzero()[1]) for _ in range(times))

assert len(cntr) == f(n)

plt.subplot(2, 1, 1).plot(cntr.values())
plt.subplot(2, 1, 2).hist(cntr.values(), n_bin)
plt.show()

从子图1可以看出每个结果的出现次数大致在650±70的范围内，从子图2可以看出每个结果出现次数的分布为接近高斯分布:

对于@AndrzejO的回答，这里使用同样的代码测试，他的解决方案速度更快(经过优化，现在两者的速度差不多了)，但是每个结果的概率似乎并不相等(请注意，各种结果也出现在这里):