解数独与解方程？

解数独与解方程？

大家最常见的那种数独游戏是给你一个9行9列的网格，里面有一些格子填上了数字，有一些没有，要你填入1到9的数字使每行每列甚至每对角线和每子区块没有重复的数字，对这种类型的简单的数独问题，我们尝试用数学方法来对其进行形式化，用数学的方式来讨论和研究数独问题.

上面这张图片呢给出了一个数独问题，你按照一定的规则在空白格子处填入数字填完所有的空白格子，就算是通关了，开局给出的数字越多越容易，开局给出的数字越少越难因为有更多的数字要去．在本文中我们不讨论求解数独问题的诀窍，最后我们会写出一个定义在 {0,1} 上的线性方程组，通过求解这个方程组来求解这个数独问题.

决定要在一个格子上填0到9范围的哪个数字算是做一个决策，为此我们设$x_{ijk}，x_{ijk} \in \{0,1\}$为决策变量，具体地，如果 xijk=1 就表示第 i 行第 j 列的格子填入了数字 k ，否则如果 xijk=0 就表示第 i 行第 j 列的格子没有填入数字 k ，稍后你会看到，我们按照这种方式设置决策变量，按照这种方式把「决策」用数学符号表示，是为了能够简明易懂地描述约束条件：行不重复列不重复块不重复斜线不重复等.

以上图这个3乘3的微型数独为例，第一行填入了3，1，2，第二行填入了2，3，1，第三行填入了1，2，3，我们可以看成是每个小格子里有3个决策变量，或者看成是每个小格子里有3个「灯泡」💡，哪个灯泡💡点亮，那取值就是多少，比如说你看到第1行1列格子第3个决策变量涂了，那第1行1列格子就取3，第1行2列的第1个决策变量涂了，那第1行2列的格子就取1，以此类推到所有任意i行k列的格．然后我们就可以用决策变量的取值来描述这个数独的各个格子的数字了.

以我们一开始展示的那个专家级别难度的数独问题为例，对于有数值的格子，例如上图画出来的那两个区块， x157=1 表示第1行第5列的格子取7， x169=1 表示第1行第6列的格子取9， x917=1 表示第9行第1列的格子取7， x923=1 表示第9行第2列的格子取3， x838=1 表示第8行第3列的格子取8. 另外你可以看到每个小格子里有9个决策变量，而这9个决策变量只有一个取1剩下的全取0，这是因为一个格子只能填入一个1到9的数字，而不能同时填入两个或多个数字.

刚才我们提到了——一个小格子只能填入一个1到9的数字而不能填入两个或多个，这实际上就是一条最基本的约束条件，怎么描述它呢？其实很简单：

xij1+xij2+xij3+xij4+xij5+xij6+xij7+xij8+xij9=1∀i,j∈{1,2,3,4,5,6,7,8,9}xijk∈{0,1}

有了这样一条约束，每一个小格子就只能填1个数字了——因为为了使式子成立 xij1,⋯,xij9 之间只有一个能取1，也必须要有一个取1. 再结合决策变量的现实意义，就是每个小格子只能填入一个数字且必须填入一个数字.

上面的约束条件保证了在任意单个小格子内数字不漏不重，而一般的数独问题还要求行不重，列不重，块不重，甚至斜线也不重，该怎么用决策变量的方程组来表示这些约束呢？

首先来表示「行不重」这个约束条件，它具体是说：对任意一行，每一行当中，1到9这9个数字要全部用掉且每个数字只能用1次，刚好对应每一行都有9个格子，每行的9个格子就是这么填的，怎么形式化呢？以第一行为例，第一行只能有一个1，于是

x111+x121+x131+x141+x151+x161+x171+x181+x191=1

上边的这个方程中的每一项实际上是第一行的每一个小格子的第一个决策变量，这些决策变量中只有一个取1，也就是说第一行只能有1个小格子取1，例如说 x141=1 就表示第1行第4列的小格子取了1，那其他的 x1j1,j=1,2,3,5,6,7,8,9 就不能取1了，也就说第一行的其他格子就不能再填入1了，这种约束条件是这样实现．写简洁一点，这个「第一行的每一个格子只能有一个填1」的约束条件可以写为

∑j=19x1j1=1

类似的，像是「第一行的每一个格子只能有一个填2」这样的约束条件可以表示为

∑j=19x1j2=1

推广到剩余的其他数字，是这样的

∑j=19x1jk=1,k=1,2,⋯,9

再推广到每一行，是这样的：

∑j=19xijk=1,k=1,2,⋯,9,i=1,2,⋯,9

上面这个式子的意思是，对 i=1,2,⋯ 行，该行不能重复$k=1,2,\cdots,9．「行不得出现重复」的约束就这样表示出来了.

类似的，也可以很轻松很简洁地表述「列不得出现重复」这样的约束：

∑i=19xijk=1,k=1,2,⋯,9,j=1,2,⋯,9

所谓「宫」就是整个9乘9的数独分成的9个3乘3的子数独，有9个宫，如下图所示

现在要表示宫(1,1)中的9个格子只能分别填入1到9，即不得重复，这样的约束条件怎么表示呢？首先，数字1不得在宫中有重复，那么我们要把宫中9个对应数字1的决策变量加起来确保它们等于1：

x111+x121+x131+x211+x221+x231+x311+x321+x331=1

类似地，在宫(1,1)中，1,2,...,9这9个数字也都是如此：

x11k+x12k+x13k+x21k+x22k+x23k+x31k+x32k+x33k=1,k=1,2,⋯,9

类似地，对每个宫，我们先找数字1对应的决策变量，把这9个决策变量加起来，固定加起来得的和为1，再找数字2对应的9个决策变量，加起来，和固定为1，再找数字3对应的9个决策变量，加起来，和固定为1，以此类推到全部9个．所以，对于宫(1,2)，应有这样的约束条件

x14k+x15k+x16k+x24k+x25k+x26k+x34k+x35k+x36k=1,k=1,2,⋯,9

我们把余下7个宫的约束条件一并写出

x17k+x18k+x19k+x27k+x28k+x29k+x37k+x38k+x39k=1,x41k+x42k+x43k+x51k+x52k+x53k+x61k+x62k+x63k=1,x44k+x45k+x46k+x54k+x55k+x56k+x64k+x65k+x66k=1,x47k+x48k+x49k+x57k+x58k+x59k+x67k+x68k+x69k=1,x71k+x72k+x73k+x81k+x82k+x83k+x91k+x92k+x93k=1,x74k+x75k+x76k+x84k+x85k+x86k+x94k+x95k+x96k=1,x77k+x78k+x79k+x87k+x88k+x89k+x97k+x98k+x99k=1,k=1,2,⋯,9.

这样就表示了行约束列约束和宫约束，斜线约束游戏中并没有提到，所以我们就不在列出了，如果要列出的话也不难，找到每一条45度或135度的斜线，对每个数字找出对应的决策变量，相加，固定和为1，这样即可表示斜线的约束.

我们把上面零散的方程片段汇总到这里：

∑k=19xijk=1,i=1,2,⋯,9,j=1,2,⋯,9

∑j=19xijk=1,k=1,2,⋯,9,i=1,2,⋯,9

∑i=19xijk=1,k=1,2,⋯,9,j=1,2,⋯,9

∑i=13∑j=13xijk=1∑i=13∑j=46xijk=1∑i=13∑j=79xijk=1∑i=46∑j=13xijk=1∑i=46∑j=46xijk=1∑i=46∑j=79xijk=1∑i=79∑j=13xijk=1∑i=79∑j=46xijk=1∑i=79∑j=79xijk=1k=1,2,⋯,9.

xijk∈{0,1},∀i,j,k∈{1,2,3,4,5,6,7,8,9}

上面的方程的所有的解对应所有的有行约束列约束和宫约束的数独的．

通过求解数学方程组来求解数独

事实上已经有专门的数独求解算法而非通过求解方程来求解数独，但是为了验证我们对数独问题建立的数学模型的正确性，我们还是打算通过求解这个方程组来求解一个数独问题.

首先呢，我们选择用Wolfram Mathematica来处理这个方程，Mathematica呢它有强大的符号处理能力，解方程之类的，不在话下，但是我们不打算直接用Mathematica来解方程，我们是打算呢，用Mathematica把这个线性方程组啊转化成这么一个矩阵的形式，什么意思呢？例如说有下面这样一个线性方程组

{x+2y+z=12x+y+3z=3

回想起我们大一学过的线性代数（也叫高等代数）的知识，根据矩阵的乘法运算法则，我们呢，可以把上面这个线性方程组写成这样一个等价的形式：

[121213][xyz]=[13]

类似的道理，我们上面列出的那一堆式子组成的方程组呢，其实也是左边是 xijk 的线性组合，而右边全是1，我们想把上面那堆方程组写成

Ax=1324×1

其中的 x 其实是一个列向量，它包含了全部 9×9×9=729 个变量，而 A 是线性方程组的系数矩阵，因为前面那堆方程组总共有324个式子，所有这个系数矩阵 A 它是324行（方程组的个数）729列（变量的个数）的，而等式右边是一个324行的全1矩．但是注意到，这个数独已经给出了一些提示，也就是一些已知信息，例如第一行第二列的格子等于5，也就是 x125=1 ，还有更多的这类信息，那么输入的信息实际上可用表示为

Bx=123×1

其中 B 照样是729列的，也就是 B 这个矩阵的每一行它有729个数，要么取1要么取0，取1的就对应 xijk=1 的系数1， B 有23行，是因为总共有23个已知数．写起来，我们将要求解这么样的一个方程组

[AB]x=1327×1

而我们使用Mathematica，于求解原来的数独问题的意义，就是将前面列出来的

∑xijk=1

这种形式的方程组转换为

Ax=1

这种矩阵形式表示的方程组.

首先要表示的是

∑k=19xijk=1,i=1,2,⋯,9,j=1,2,⋯,9

它对应Mathematica代码

mat1 = Coefficient[#, variables] & /@ Flatten[Table[
     Sum[x[i][j][k], {k, 1, 9}], {i, 1, 9}, {j, 1, 9}
]];

∑j=19xijk=1,k=1,2,⋯,9,i=1,2,⋯,9

对应的Mathematica代码是

mat2 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[x[i][j][k], {j, 1, 9}], 
    {k, 1, 9}, {i, 1, 9}
]];

∑i=19xijk=1,k=1,2,⋯,9,j=1,2,⋯,9

同样转化为Mathematica代码

mat3 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[x[i][j][k], {i, 1, 9}], {k, 1, 9}, {j, 1, 9}
]];

下面是每个宫的9个约束：

∑i=13∑j=13xijk=1∑i=13∑j=46xijk=1∑i=13∑j=79xijk=1∑i=46∑j=13xijk=1∑i=46∑j=46xijk=1∑i=46∑j=79xijk=1∑i=79∑j=13xijk=1∑i=79∑j=46xijk=1∑i=79∑j=79xijk=1k=1,2,⋯,9.

转化为Mathematica代码就是：

mat4 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 1, 3}], {i, 1, 3}], 
    {k, 1, 9}
]];

mat5 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 4, 6}], {i, 1, 3}], 
    {k, 1, 9}
]];

mat6 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 7, 9}], {i, 1, 3}], 
    {k, 1, 9}
]];

mat7 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 1, 3}], {i, 4, 6}], 
    {k, 1, 9}
]];

mat8 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 4, 6}], {i, 4, 6}], 
    {k, 1, 9}
]];

mat9 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 7, 9}], {i, 4, 6}], 
    {k, 1, 9}
]];

mat10 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 1, 3}], {i, 7, 9}], 
    {k, 1, 9}
]];

mat11 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 4, 6}], {i, 7, 9}], 
    {k, 1, 9}
]];

mat12 = Coefficient[#, variables] & /@ Flatten[Table[
    Sum[Sum[x[i][j][k], {j, 7, 9}], {i, 7, 9}], 
    {k, 1, 9}
]];

以及还有已知信息：

inputSudoku = {
   {0, 5, 0, 0, 7, 9, 0, 0, 0},
   {0, 0, 0, 0, 0, 0, 5, 0, 0},
   {0, 9, 2, 0, 0, 0, 0, 6, 0},
   {0, 8, 0, 0, 0, 0, 4, 0, 7},
   {0, 2, 0, 6, 0, 0, 0, 1, 0},
   {0, 7, 0, 2, 5, 0, 0, 8, 0},
   {0, 0, 0, 0, 4, 0, 0, 0, 0},
   {0, 0, 8, 0, 2, 0, 0, 0, 0},
   {7, 3, 0, 0, 0, 1, 0, 0, 0}
};

inputValues = Flatten[Table[
    If[inputSudoku[[i]][[j]] == k, x[i][j][k], Nothing],
    {i, 1, 9}, {j, 1, 9}, {k, 1, 9}
]];

mati = Coefficient[#, variables] & /@ inputValues;

我们把上面所有小矩阵合并成一个大矩阵

matn = Join[
   mat1,
   mat2,
   mat3,
   mat4,
   mat5,
   mat6,
   mat7,
   mat8,
   mat9,
   mat10,
   mat11,
   mat12,
   mati
];

这时如果我们查看总系数矩阵matn的维数

Dimensions[matn]

{347, 729}

接下来我们把这个总系数矩阵导出为.csv后缀名的CSV文件，供MATLAB求解，因为我不知道怎么用Mathematica快速高效地求解二元整数变量组成的方程组：

Export["~/Desktop/matn.csv", matn]

现在matn.csv这个文件就包含了足够清楚地描述这个9乘9数独问题的全部信息，matn.csv就是那个系数矩阵，我们把它导入到MATLAB中并求解，下面是MATLAB求解数独问题的代码

clear
matData = readtable("~/Desktop/matn.csv");
aeq = matData{:, :};
beq = ones([347, 1]);
intcon = 1:729;
c = zeros([1, 729]);
lb = zeros([729, 1]);
ub = ones([729, 1]);
sol = intlinprog(c, intcon, [], [], aeq, beq, lb, ub);
deciders = reshape(sol, 9, 9, 9);
for i = 1:9
    for j = 1:9
        temp = deciders(i, j, :);
        sudoku(i, j) = [1:9]*temp(:);
    end
end
sudoku
writetable( ...
    table(sudoku), ...
    "~/Desktop/sudoku.csv", ...
    "WriteVariableNames", false ...
);

看一下输出结果也无妨：

可以很容易地验证，这个解的每一行当中没有重复的数字，每一列当中没有重复的数字，每一个宫中也没有重复的数字.

额外的讨论

我们算出的系数矩阵是347行729列的，也就是说至少有 729−347=382 个自由变量，因为该系数矩阵的秩最多会是347，解空间实际上至少是个382维的空间，而要在这个382维的空间寻找一个729维的坐标，或者729个元素的向量，这个坐标（向量）的每个分量都只能是0或1，就是依方程组方法求解数独问题的几何解．而在这个347维的空间包含了多少个这样的729元的向量，那么这个数独问题就有多少个解.

我们先用线性方程组，对一类简单的数独问题进行了建模，准确的说是用 {0,1} 二元变量来表示最基本的决策单元，然后我们又用Mathematica将列出的方程组转化为系数矩阵，然后我们用MATLAB在 {0,1} 的定义域内求解这个系数矩阵所表示的方程组，这样就完成了这类简单数独问题的求解.