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

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

TL；博士

每个结果均以相等的概率生成，并以

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(_prob_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的回答，这里使用了同样的代码测试，他的解决方案更快（优化后，现在两者的速度几乎相同），但是每个结果的概率似乎并不相等（注意各种结果也出现在这里）：

Answer 2

创建一个包含零的矩阵。然后，您需要随机获取 N 个不重复的行号，以及随机的 N 个不重复的列号。您可以使用

random.sample

来实现此目的。然后在行/列位置上放1。

import numpy as np
from random import sample

N = int(input("Enter the number of row and col:"))
my_matrix = np.zeros((N,N), dtype='int8')
rows = sample(range(N), N)
cols = sample(range(N), N)
points = zip(rows, cols)
for x, y in points:
    my_matrix[x, y] = 1

print(my_matrix)

如果你想要一个对称矩阵：如果N是偶数，我会从N中取出N个随机数，一半作为x，一半作为y；在 (x,y) 和 (y,x) 两个位置上都加 1。如果 N 不均匀，则需要在对角线上的随机位置额外加 1。

import numpy as np
from random import sample, randint

N = int(input("Enter the number of row and col:"))
even = N%2 == 0

my_matrix = np.zeros((N,N), dtype='int8')
N_range = list(range(N))

if not even:
    diagonal = randint(0, N-1)
    N_range.remove(diagonal)
    my_matrix[diagonal, diagonal] = 1
    N = N - 1

rowcol = sample(N_range, N)

rows = rowcol[:N//2]
cols = rowcol[N//2:]

for x, y in zip(rows,cols):
    my_matrix[x, y] = 1
    my_matrix[y, x] = 1

这是一个更好的版本。取第一个空闲行，随机获取一个空闲列，在 (row,col) 和 (col,row) 上放置 1。删除使用的列/行。重复直到使用所有数字 0-(N-1)。

import numpy as np
import random

N = int(input("Enter the number of row and col:"))

my_matrix=np.zeros((N,N))
not_used_number = list(range(N))

while len(not_used_number) != 0:
    current_row = not_used_number[0]
    random_col = random.choice(not_used_number)
    my_matrix[current_row, random_col] = 1
    my_matrix[random_col, current_row] = 1
    
    not_used_number.remove(current_row)
    if current_row != random_col:
        not_used_number.remove(random_col)