All Hail Scarlet!

0%

N Queens: What's the Limit?

Posted on Edited on
N-皇后问题求解的一系列尝试和优化方法，并介绍了一种构造性解法。

问题背景

八皇后问题是搜索算法中一个著名的问题。

（本段资料来自百度百科）该问题是国际西洋棋棋手马克斯·贝瑟尔于1848年提出：在 8\times8 格的国际象棋棋盘上摆放八个皇后，使其不能互相攻击，即任意两个皇后都不能处于同一行、同一列或同一斜线上，问有多少种摆法。

用计算机不难验证答案是 92 种。那么用计算机在合理时间内（数秒），最多可以求解 N 为多少的 N 皇后问题呢？

我们知道，随着 N 的增加，枚举所有可能摆法所需要的时间是指数级增加的。如果用程序去枚举随 N 增加的解的数量变化规律，不难发现解的数量也是指数级增加的。例如，在 N = 25 的时候，解的数量已经高达 10^{15} 数量级。

也就是说，就算有一个算法可以快速找到所有的解，它所消耗的时间仍然随着 N 的增加而指数级增加。

因而我们转而考虑这样一个 N 皇后问题：找到一种方法，在 N \times N 格的棋盘上摆放 N 个皇后，使其不能互相攻击。

盲目搜索

第一类解决方法，就是用以解决八皇后问题的盲目搜索算法。这里给出博客中代码共用的Code Base。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import queue


def queen_data(x, y):
    return (x, y, x + y, x - y)


def collides(qdata_p, qdata_q):
    return any(p == q for p, q in zip(qdata_p, qdata_q))


def evaluate(queens):
    judge = [collides(q1, q2) for q1 in queens for q2 in queens if q1 != q2]
    return judge.count(True)


def build_matrix(queens, N):
    output = [[0] * N for _ in range(N)]
    for queen in queens:
        qx, qy, _, _ = queen
        output[qx][qy] = 1
    return output


def print_state(queens, N):
    output = build_matrix(queens, N)
    for i in range(N):
        print(output[i])
    print()

首先，我们来看看最最最朴素，且不使用任何剪枝等优化技术的DFS搜索可以做得怎么样。（本机环境：Core i7 4700MQ）

1
2
3
4
5
6
7
8
9
10
11
12
def uninformed_search(N):
    dfq = queue.deque()
    dfq.append([])
    while len(dfq) > 0:
        cur = dfq.pop()
        if len(cur) == N:
            if evaluate(cur) == 0:
                return cur
        else:
            nx = len(cur)
            for ny in range(N):
                dfq.append(cur + [queen_data(nx, ny)])
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
$ time print_state(uninformed_search(4), 4)
[0, 0, 1, 0]
[1, 0, 0, 0]
[0, 0, 0, 1]
[0, 1, 0, 0]

Wall time: 2 ms

$ time print_state(uninformed_search(7), 7)
[0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 1, 0, 0]
[0, 0, 1, 0, 0, 0, 0]
[1, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 1, 0]
[0, 0, 0, 1, 0, 0, 0]
[0, 1, 0, 0, 0, 0, 0]

Wall time: 2.26 s

嗯，性能意料之内的糟糕。8以上基本就跑不动了。（注：你或许会说创建每个新状态用了 O(N) 的时间，以及python语言的执行效率等等——We don’t care about this here）
加上简单的剪枝：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
def uninformed_search_cut(N):
    dfq = queue.deque()
    dfq.append([])
    while len(dfq) > 0:
        cur = dfq.pop()
        if len(cur) == N:
            if evaluate(cur) == 0:
                return cur
        else:
            nx = len(cur)
            for ny in range(N):
                new_state = cur + [queen_data(nx, ny)]
                if evaluate(new_state) == 0:
                    dfq.append(new_state)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
$ time print_state(uninformed_search_cut(7), 7)
[0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 1, 0, 0]
[0, 0, 1, 0, 0, 0, 0]
[1, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 1, 0]
[0, 0, 0, 1, 0, 0, 0]
[0, 1, 0, 0, 0, 0, 0]

Wall time: 2 ms

$ time print_state(uninformed_search_cut(15), 15)
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0]
[0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0]
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0]

Wall time: 2.55 s

N = 16 已经需要超过20s的运行时间。简单的剪枝让可以计算的 N 值翻了一倍。
下面我们稍微换一下思路：仍然是盲目搜索，但是我们通过“把其他行可能放置皇后的位置不断剔除”的方式来判断是否还有可能解进行剪枝。大体效率应当没有提升。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def determine_next(var_domains, vx):
    if vx == len(var_domains):
        return [v for v, in var_domains]
    vdom = var_domains[vx]
    for v in vdom:
        ndom = []
        for vsx, vsdom in enumerate(var_domains):
            if vsx == vx:
                ndom.append([v])
                continue
            vndom = [vs for vs in vsdom
                     if not collides(queen_data(vx, v),
                                     queen_data(vsx, vs))]
            if len(vndom) > 0:
                ndom.append(vndom)
            else:
                ndom = None
                break
        if ndom is not None:
            result = determine_next(ndom, vx + 1)
            if result is not None:
                return result


def domain_reduction(N):
    init_var_domains = [list(range(N)) for _ in range(N)]
    result = determine_next(init_var_domains, 0)
    queens = [queen_data(x, y) for x, y in enumerate(result)]
    assert evaluate(queens) == 0
    return queens
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
$ time domain_reduction(15)
Wall time: 44 ms
$ time domain_reduction(16)
Wall time: 326 ms
$ time domain_reduction(17)
Wall time: 197 ms
$ time domain_reduction(18)
Wall time: 1.44 s
$ time domain_reduction(19)
Wall time: 104 ms
$ time domain_reduction(20)
Wall time: 7.78 s
$ time domain_reduction(21)
Wall time: 346 ms
$ time domain_reduction(23)
Wall time: 1.1 s
$ time domain_reduction(25)
Wall time: 2.28 s

由于解的输出太过庞大，就不打出来了。可以看到一件有趣的事情：对奇数的情况远比偶数的情况要快。因而，找到解的效率可能和具体结构有很大关联。为了屏蔽这种关联性，我们将筛选值域的顺序随机化。

1
2
3
4
5
6
7
8
9
def domain_reduction_random(N):
    init_var_domains = [list(range(N)) for _ in range(N)]
    import random
    for dom in init_var_domains:
        random.shuffle(dom)
    result = determine_next(init_var_domains, 0)
    queens = [queen_data(x, y) for x, y in enumerate(result)]
    assert evaluate(queens) == 0
    return queens
1
2
3
4
5
6
7
8
9
10
$ time domain_reduction_random(25)
Wall time: 9 ms
$ time domain_reduction_random(40)
Wall time: 193 ms
$ time domain_reduction_random(80)
Wall time: 247 ms
$ time domain_reduction_random(160)
Wall time: 1.66 s
$ time domain_reduction_random(320)
Wall time: 13.3 s

随机化的效率提升是不是有些出乎意料？可见八皇后问题的解是有其特定结构的，而自然顺序是一个较为糟糕的顺序。

终极解法：构造法

构造法利用解的特定结构，直接构造一组解来解决问题。

N 皇后问题解不唯一，构造的方式也有很多。下面介绍一种由E. J. Hoffman, J. C. Loessi和R. C. Moore提出的构造方式。

首先用两种方式构造边长为偶数格的情形。
A. 若 N = 2m,\ m \neq 3k + 1 ，则从第一行第二列开始，每行向右2格。到一行某尾后下一行从第一列开始，随后继续每行向右两格。例如对 12\times12 ，构造的结果为：
12 by 12 construction
B. 若 N = 2m,\ m \neq 3k ，则对前 m 行，在 (i, 1 + ((2i + m - 3)\ mod\ N)) 上放置皇后。随后沿棋盘中心点作中心对称图形。例如对 14\times14 ，构造的结果为：
14 by 14 construction
随后解决边长为奇数格的情况。
C. 若 N 为奇数，在棋盘右下角放一个皇后，随后对左上角的棋盘使用构造方法A或B。

这样我们就得到了一组 N 皇后问题的解的构造法。

一点花絮

局部搜索和基于Genetic-Programming的搜索效率蜜汁低下，且容易陷入局部最小的困境。本来以为可以比盲目搜索更快的来着……可能“冲突数”在这个问题中不是一个很好的评估吧。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
def gp(N):
    import random as R
    import math
    import functools

    def mutation(indiv, n_mutations=N//10):
        newgen = indiv.copy()
        for i in range(n_mutations):
            x = R.randint(0, N - 1)
            newgen[x] = queen_data(x, R.randint(0, N - 1))
        return newgen

    def crossover(indiv_p, indiv_q):
        split = R.randint(1, N - 1)
        return indiv_p[:split] + indiv_q[split:]

    def fitness(indiv):
        return -evaluate(indiv)

    def search_up(individuals, afit, rs):
        prefix = 0.0
        result = []
        require = len(rs)
        rs = iter(rs)
        r = next(rs)
        for i, fit in enumerate(afit):
            if prefix < r and prefix + fit >= r:
                result.append(individuals[i])
                if len(result) == require:
                    return result
                r = next(rs)
            prefix += fit
        while len(result) < require:
            result.append(individuals[-1])
        return result

    def choose_parent(individuals, fitnesses, count=1, kdiv=10.0):
        mina, maxa = min(fitnesses), max(fitnesses)
        afit = list(map(lambda x: math.exp(kdiv * (x - mina) / (maxa - mina)),
                        fitnesses))
        suma = sum(afit)
        rs = sorted(R.random() * suma for _ in range(count))
        return search_up(individuals, afit, rs)

    def random_indiv():
        return [queen_data(x, R.randint(0, N - 1)) for x in range(N)]

    def next_generation(individuals,
                        kcopy=0.3,
                        kmutation=0.03):
        fitnesses = [fitness(indiv) for indiv in individuals]
        for i, fit in enumerate(fitnesses):
            if fit == 0:
                return [individuals[i]], 0
        parents = functools.partial(choose_parent, individuals, fitnesses)
        prec = parents(count=int(len(individuals) * kcopy))
        while len(prec) < len(individuals):
            prec.append(crossover(*parents(count=2)))
        for i in range(len(prec)):
            if R.random() < kmutation:
                prec[i] = mutation(prec[i])
        return prec, max(fitnesses)

    population = [random_indiv() for _ in range(500)]
    for gen in range(100):
        population, bestfit = next_generation(population)
        if bestfit == 0:
            indiv, = population
            assert evaluate(indiv) == 0
            return indiv
        print("Generation", gen, "best fitness:", bestfit)