使用NumPy生成拉丁方阵的高效方法（或在两个轴上独特地随机排列矩阵中的数字）

Question

例如，如果有5个数字1、2、3、4、5，

我希望得到一个随机结果，如下：

[[ 2, 3, 1, 4, 5]
 [ 5, 1, 2, 3, 4]
 [ 3, 2, 4, 5, 1]
 [ 1, 4, 5, 2, 3]
 [ 4, 5, 3, 1, 2]]

确保每个数字在其行和列中都是唯一的。

有没有一种高效的方法来实现这个目标？

我尝试使用 while 循环来生成每次迭代的一行，但是似乎不是很高效。

import numpy as np

numbers = list(range(1,6))
result = np.zeros((5,5), dtype='int32')
row_index = 0
while row_index < 5:
    np.random.shuffle(numbers)
    for column_index, number in enumerate(numbers):
       if number in result[:, column_index]:
           break
       else:
           result[row_index, :] = numbers
           row_index += 1

- KHC

---

我们可以随机选择一行数据。然后将此行[1--(len-1)]列移动以形成下一行。我在下面的回复中放了一些代码。这不是完全随机，但它是以固定格式随机的。~ - Jilin

你所需要的是一种生成拉丁方阵的方法。https://en.wikipedia.org/wiki/Latin_square - norok2

Answer 1

仅供参考，您正在寻找的是生成拉丁方阵的方法。 至于解决方案，它取决于您认为多少随机是“随机”的。

我会设计至少四种主要技术，其中两种已经被提出。 因此，我将简要描述另外两个：

1. 循环遍历所有可能的项目排列，并接受满足行唯一性约束条件的第一个
2. 仅使用循环置换来构建后续行：这些行通过构造在行上满足唯一性约束条件的循环变换进行（循环变换可以向前或向后进行）；为了改善“随机性”，可以对行进行洗牌

假设我们使用标准的Python数据类型，因为我认为使用NumPy没有真正的优点（但如果需要，结果可以轻松转换为np.ndarray），这将是代码（第一个函数只是用来检查解决方案是否正确）：

import random
import math
import itertools

# this only works for Iterable[Iterable]
def is_latin_rectangle(rows):
    valid = True
    for row in rows:
        if len(set(row)) < len(row):
            valid = False
    if valid and rows:
        for i, val in enumerate(rows[0]):
            col = [row[i] for row in rows]
            if len(set(col)) < len(col):
                valid = False
                break
    return valid

def is_latin_square(rows):
    return is_latin_rectangle(rows) and len(rows) == len(rows[0])

# : prepare the input
n = 9
items = list(range(1, n + 1))
# shuffle items
random.shuffle(items)
# number of permutations
print(math.factorial(n))


def latin_square1(items, shuffle=True):
    result = []
    for elems in itertools.permutations(items):
        valid = True
        for i, elem in enumerate(elems):
            orthogonals = [x[i] for x in result] + [elem]
            if len(set(orthogonals)) < len(orthogonals):
                valid = False
                break
        if valid:
            result.append(elems)
    if shuffle:
        random.shuffle(result)
    return result

rows1 = latin_square1(items)
for row in rows1:
    print(row)
print(is_latin_square(rows1))


def latin_square2(items, shuffle=True, forward=False):
    sign = -1 if forward else 1
    result = [items[sign * i:] + items[:sign * i] for i in range(len(items))]
    if shuffle:
        random.shuffle(result)
    return result

rows2 = latin_square2(items)
for row in rows2:
    print(row)
print(is_latin_square(rows2))

rows2b = latin_square2(items, False)
for row in rows2b:
    print(row)
print(is_latin_square(rows2b))

为了对比，我们还提供了一种尝试随机排列并接受有效排列的实现方法（基本上就是@hpaulj提出的方法）。

def latin_square3(items):
    result = [list(items)]
    while len(result) < len(items):
        new_row = list(items)
        random.shuffle(new_row)
        result.append(new_row)
        if not is_latin_rectangle(result):
            result = result[:-1]
    return result

rows3 = latin_square3(items)
for row in rows3:
    print(row)
print(is_latin_square(rows3))

我还没有时间实现另一种方法（使用来自@ConfusedByCode的回溯数独解决方案）。

对于n = 5的时间记录：

%timeit latin_square1(items)
321 µs ± 24.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit latin_square2(items)
7.5 µs ± 222 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit latin_square2(items, False)
2.21 µs ± 69.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit latin_square3(items)
2.15 ms ± 102 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

...并且对于n = 9：

%timeit latin_square1(items)
895 ms ± 18.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit latin_square2(items)
12.5 µs ± 200 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit latin_square2(items, False)
3.55 µs ± 55.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

%timeit latin_square3(items)
The slowest run took 36.54 times longer than the fastest. This could mean that an intermediate result is being cached.
9.76 s ± 9.23 s per loop (mean ± std. dev. of 7 runs, 1 loop each)

所以，解决方案1提供了相当数量的随机性，但速度不是特别快（并且与O(n!)缩放），解决方案2（和2b）要快得多（缩放为O(n)），但不像解决方案1那样随机。解决方案3非常慢，性能可能会有很大差异（可以通过计算而不是猜测来加速最后一次迭代）。

更多技术细节，请参考以下高效算法：

Jacobson, M. T. and Matthews, P. (1996), Generating uniformly distributed random latin squares. J. Combin. Designs, 4: 405-437. doi:10.1002/(SICI)1520-6610(1996)4:6<405::AID-JCD3>3.0.CO;2-J

- norok2

Answer 2

这可能看起来很奇怪，但您基本上描述了生成随机n维数独谜题的过程。引用自Daniel Beer的博客文章：

解决数独谜题的基本方法是通过对每个单元格的候选值进行回溯搜索。一般的步骤如下：

1. 为每个单元格生成一个候选值列表，方法是从所有可能的值集合开始，并消除在检查的单元格所在行、列和盒子中出现的值。

2. 选择一个空单元格。如果没有可用的，则谜题已解决。

3. 如果单元格没有候选值，则谜题无法解决。

4. 对于该单元格中的每个候选值，将该值放入单元格并尝试递归解决谜题。

有两个优化可以极大地提高此算法的性能：

1. 在选择单元格时，始终选择具有最少候选值的单元格。这减少了分支因子。随着网格中添加值，其他单元格的候选数也会减少。

2. 在分析空单元格的候选值时，从上一步开始分析并通过删除最后修改单元格的行、列和盒子中的值来修改它要快得多。这是谜题大小的O(N) ，而从头开始分析是O(N3)。

在您的情况下，“无法解决的谜题”是一个无效的矩阵。在可解的谜题中，矩阵中每个元素在两个轴上都是唯一的。

- user9478968

---

我认为这里谈论的是解决难题，这可能比生成难题本身更复杂。 - norok2

提议的步骤依然可行，只是您需要使用空白的网格开始。随机化第二步，结果将会是一份随机数独谜题。 - user9478968

解决数独谜题就是填写空白方格。在这种情况下，所有的方格都是空白的。这两个任务的复杂度是相同的。 - user9478968

我并不是在争论这个方法一定会让你得到解决方案，但解决方案的复杂度和随机性很大程度上取决于实现。因此，我认为这种宽泛的描述可能暗示着某些解决方案，但实际上并没有太多帮助。 - norok2

Answer 3

编辑：以下是norok2答案中第二种解决方案的实现。

编辑：我们可以再次洗牌生成的正方形，使其真正随机。 因此，解决函数可以进行修改：

def solve(numbers):
    shuffle(numbers)
    shift = randint(1, len(numbers)-1)
    res = []

    for r in xrange(len(numbers)):
        res.append(list(numbers))
        numbers = list(numbers[shift:] + numbers[0:shift])

    rows = range(len(numbers))
    shuffle(rows)

    shuffled_res = []
    for i in xrange(len(rows)):
        shuffled_res.append(res[rows[i]])

    return shuffled_res

编辑：我之前误解了问题。 因此，这里有一种“快速”方法可以生成“在某种程度上”的随机解决方案。 基本思路是，

    a, b, c
    b, c, a
    c, a, b

我们可以通过固定步长移动一行数据来形成下一行，这将符合我们的限制要求。

所以，这是代码：

from random import shuffle, randint


def solve(numbers):
    shuffle(numbers)
    shift = randint(1, len(numbers)-1)
    res = []

    for r in xrange(len(numbers)):
        res.append(list(numbers))
        numbers = list(numbers[shift:] + numbers[0:shift])

    return res


def check(arr):
    for c in xrange(len(arr)):
        col = [arr[r][c] for r in xrange(len(arr))]
        if len(set(col)) != len(col):
            return False
    return True


if __name__ == '__main__':
    from pprint import pprint
    res = solve(range(5))
    pprint(res)
    print check(res)

---

这是使用itertools的可能解决方案，如果您不坚持使用我不熟悉的numpy：

import itertools
from random import randint
list(itertools.permutations(range(1, 6)))[randint(0, len(range(1, 6))]

# itertools returns a iterator of all possible permutations of the given list.

- Jilin

---

我不确定在这里使用“所有可能的排列”是否可行，因为OP正在寻找每行/列中每个数字仅出现一次的结果。 - Jonathan Dayton

是的，这个方法不能确保每个数字在其行或列中都是唯一的，感谢指出。 - Jilin

这是我提出的第二种方法的实现是否更为特定且速度较慢，还是还有其他因素？ - norok2

没错，你之前已经提出了这个解决方案。 - Jilin

Answer 4

我尝试了一种暴力随机选择的方法。生成一行，如果有效，则添加到累积的行中：

def foo(n=5,maxi=200):
    arr = np.random.choice(numbers,n, replace=False)[None,:]
    for i in range(maxi):
        row = np.random.choice(numbers,n, replace=False)[None,:]
        if (arr==row).any(): continue
        arr = np.concatenate((arr, row),axis=0)
        if arr.shape[0]==n: break
    print(i)
    return arr

一些示例运行：

In [66]: print(foo())
199
[[1 5 4 2 3]
 [4 1 5 3 2]
 [5 3 2 1 4]
 [2 4 3 5 1]]
In [67]: print(foo())
100
[[4 2 3 1 5]
 [1 4 5 3 2]
 [5 1 2 4 3]
 [3 5 1 2 4]
 [2 3 4 5 1]]
In [68]: print(foo())
57
[[1 4 5 3 2]
 [2 1 3 4 5]
 [3 5 4 2 1]
 [5 3 2 1 4]
 [4 2 1 5 3]]
In [69]: print(foo())
174
[[2 1 5 4 3]
 [3 4 1 2 5]
 [1 3 2 5 4]
 [4 5 3 1 2]
 [5 2 4 3 1]]
In [76]: print(foo())
41
[[3 4 5 1 2]
 [1 5 2 3 4]
 [5 2 3 4 1]
 [2 1 4 5 3]
 [4 3 1 2 5]]

需要尝试的次数因情况而异，有些超过了我的迭代限制。

不涉及任何理论，快速生成二维排列和生成在某种意义上最大随机性的排列之间存在差异。我怀疑我的方法比更系统和高效的方法更接近这个随机目标（但我无法证明）。

---

def opFoo():
    numbers = list(range(1,6))
    result = np.zeros((5,5), dtype='int32')
    row_index = 0; i = 0
    while row_index < 5:
        np.random.shuffle(numbers)
        for column_index, number in enumerate(numbers):
            if number in result[:, column_index]:
                break
            else:
                result[row_index, :] = numbers
                row_index += 1
        i += 1
    return i, result

In [125]: opFoo()
Out[125]: 
(11, array([[2, 3, 1, 5, 4],
        [4, 5, 1, 2, 3],
        [3, 1, 2, 4, 5],
        [1, 3, 5, 4, 2],
        [5, 3, 4, 2, 1]]))

我的代码比原作者的慢一些，但是我的结果是正确的。

---

这是对我的改进（速度提高了2倍）：

def foo1(n=5,maxi=300):
    numbers = np.arange(1,n+1)
    np.random.shuffle(numbers)
    arr = numbers.copy()[None,:]
    for i in range(maxi):
        np.random.shuffle(numbers)
        if (arr==numbers).any(): continue
        arr = np.concatenate((arr, numbers[None,:]),axis=0)
        if arr.shape[0]==n: break
    return arr, i

---

为什么翻译后的数独求解器比原版慢？

我发现使用Python列表比使用numpy数组更快，这是在翻译Java数独求解器时发现的。

明天我可能会尝试将那个脚本适应到这个问题上。

- hpaulj

---

这与OP算法的效率相比如何？我怀疑这实际上甚至更不高效。 - norok2

我的速度明显较慢。OP的程序在尝试10行左右时就退出了。但是我的程序是正确的。 - hpaulj

在我的回答中所描述的回溯算法会更快且仍然正确。回溯只在必要时重新计算一个单元格，而不是随机猜测一行。 - user9478968

我的链接数独求解器是为9x9的正方形（带有3x3子正方形）设计的，但我已经适应它可以在6x6上工作（'空白'除了一个随机的第一行）。 它比这个foo1慢（3倍），但是foo1不具有良好的可扩展性。 对于n=9，它需要更多的迭代次数。 @ConfusedByCode - hpaulj

Answer 5

手机不能输入代码，以下是伪代码:

1. 创建一个比目标矩阵（3 d）多一维的矩阵

2. 用数字1到5初始化25个元素

3. 遍历这25个元素。

4. 从元素列表中选择第一个元素的随机值（其中包含数字1到5）

5. 从行和列中删除所有元素中随机选择的值。

6. 对所有元素重复步骤4和5。

- saga