优化杂记（二）：Munkres算法与分配问题
俗话说得好，追根结底，所有问题都是分配问题，所以今天来看看运筹学中的分配问题（Assignment Problem），也称为指派问题。情景为：假设有N个大学生，有N个岗位， 第 i 个大学生做第j份工作要价为C_{ij}, 这些C_{ij}组成了N\times N的花费矩阵（Cost Matrix）。要把这N个大学生安排到这N个岗位中，使得总花费最小，其数学优化模型A为：
\begin{align*} \min&~\sum_{1\leq i,j\leq N} C_{ij}x_{ij}\\ s.t.&~ x_{ij}\in\{0,1\}\\ &\sum_{i=1}^N x_{ij}=1,\forall 1\leq j\leq N\\ &\sum_{j=1}^N x_{ij}=1,\forall 1\leq j\leq N \end{align*} \\
[*]x_{ij}是二值变量，值为1表示将工作j分给了大学生i；
[*]\sum_{i=1}^N x_{ij}=1表示工作j安排且只安排了一个大学生；
[*]\sum_{j=1}^N x_{ij}=1表示大学生i安排且只安排了一项工作。
上面是个混合整数规划问题（Mixed Integer Programming），典型组合优化，爆破需要O(N!)，然而将其松弛（relax）一下，变成模型B：
\begin{align*} \min&~\sum_{1\leq i,j\leq N} C_{ij}x_{ij}\\ s.t.&~ 0\leq x_{ij}\leq 1\\ &\sum_{i=1}^N x_{ij}=1,\forall 1\leq j\leq N\\ &\sum_{j=1}^N x_{ij}=1,\forall 1\leq j\leq N \end{align*} \\
注意这个问题是线性规划问题，其可行集是个凸多面体，而每个顶点都是原问题A的解，由于目标函数也是线性的，求解这个线性规划模型B，某个顶点一定在其解集内，那么可以认为得到的解就是模型A的解。
把一个混合整数规划这种典型组合优化问题转化为了线性规划这种凸问题，是非常神奇的，假设花费矩阵为C：
\begin{bmatrix} 5&9&1\\ 10& 3& 2\\ 8& 7& 4\\ \end{bmatrix} \\
那么可以用google的ortools或者别的cvxpy来求解一下：
from ortools.linear_solver import pywraplp
cost = [
,
,
]
# Create the mip solver with the SCIP backend.
solver = pywraplp.Solver.CreateSolver('SCIP')
N = 3
x = {}
for j in range(N):
for k in range(N):
x[(j, k)] = solver.NumVar(0, 1, f'x[{j},{k}]')
print('Number of variables =', solver.NumVariables())
for j in range(N):
solver.Add(sum( for k in range(N)])==1)
for k in range(N):
solver.Add(sum( for j in range(N)])==1)
objective = solver.Objective()
for j in range(N):
for k in range(N):
objective.SetCoefficient(x[(j, k)], cost)
objective.SetMinimization()
status = solver.Solve()
结果为
\begin{bmatrix} 0&0&1\\ 0& 1& 0\\ 1& 0& 0\\ \end{bmatrix} \\
最优值为12。虽然能求解，但注意其变量数为N^2，那么使用线性规划总复杂度差不多是O(N^5)或者稍微再少点. 规模大一点估计就不太行了。于是又有了针对该问题的巧妙解法，比如Hungarian 方法，或者叫Munkres算法，其复杂度为O(N^3).
主要有六大步骤，写得繁琐，其实很简单，算法出口在第3步：
[*]对于每一行，找到最小值，从每一行中减去该最小值，这样每行至少有一个零，这些零的位置代表了让每个大学生都有工作的安排的条件下，最小的消耗。但可能有的工作没人做，后面几个步骤就是不停调整解决这个问题。执行第2步；
[*]逐个遍历矩阵C的零元素，如果该零元素的行与列中没有被“标记”的零元素，就将其“标记”。执行第3步；
[*]对于C的每一列，如果有一个“标记”的零元素，则“覆盖”该列。如果所有列都被“覆盖”，就意味着每一列都有一个被“标记”的零元素，也就代表了每一行刚好也有一个被“标记”的零元素，那么结束了，得到结果了；否则，执行第4步；
[*]尝试在没被“覆盖”的位置上找到一个零元素，如果没有，保存最小的未被覆盖的零元素，执行第6步；否则，“刻记”该元素，如果该元素所在行没有“标记”元素，执行第5步；如果该元素所在行有“标记”元素，将“标记”元素的列取消“覆盖”，并将此行“覆盖”，重复执行第4步；
[*]这步比较神奇，注意这步是从第4步过来，从该“刻记”元素出发，尝试找到同列的“标记”元素，再找“同行”的“刻记”元素。。。直到某个“刻记”元素找不到下个一个“标记”元素。清除矩阵中的所有“标记”元素，然后将序列中的所有“刻记”元素转化为“标记”元素，清除所有覆盖，执行第3步，（“标记”元素马上又会被覆盖）。
[*]未覆盖区域减去未覆盖区域中最小的元素，执行第4步，注意此时产生了一个未覆盖的零元素；
python的scipy已经给出了实现，直接执行即可，没啥观赏性，就不重造轮子了。好奇的话可以查看一个叫munkres的python包看到具体代码。有了指派问题求解模型，就可以回答下面一个问题：
平面上有一组N个点
A=\{P_1,P_2,\cdots,P_N\} \\
要将其挪到指定的N个位置：
B=\{Q_1,Q_2,\cdots,Q_N \} \\
具体的情景可以是N个货物发往N个地点，代价如何最小，这就可以建模为典型的指派问题。
好，这里给定点集A为一组随机点（理解为货物）：
位置B为某些给定的位置（理解为收货地）：
(B的生成作图参考另一篇专栏文章)
代价为对应点与位置之间距离的平方和，其指派方案做成动图为：
代码为：
from scipy.optimize import linear_sum_assignment
points2 = np.load("covid.npz")["xy"]
N = len(points2)
points1 = np.random.rand(N, 2)
plt.scatter(points1[:, 0], points1[:, 1])
plt.scatter(points2[:, 0], points2[:, 1], s=1)
cost_matrix = points1.reshape((N, 1, 2)) - points2.reshape((1, N, 2))
cost_matrix = (cost_matrix ** 2).sum(-1)
indexes = linear_sum_assignment(cost_matrix)
indexes = list(zip(*indexes))
total = 0
for row, column in indexes:
value = cost_matrix
total += value
今天了解了一下指派问题，用在了计算生成艺术上，就这样吧~
参考资料：
[*]https://software.clapper.org/munkres/
[*]https://brc2.com/the-algorithm-workshop/
[*]https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.optimize.linear_sum_assignment.html
[*]https://developers.google.com/optimization/assignment/assignment_example
页:
[1]