🟠 七桥问题与图论\#
哥尼斯堡七桥问题算是图论问题的开端了，大家在入门图论的时候应该都是从这个问题开始的吧。在本文中，我将淡化这个问题中的数学部分，更多地在计算机领域呈现这个问题，毕竟图论也属于计算机中的一大研究领域
问题抽象\#
七桥问题讲的是哥尼斯堡有七座桥按如下方式连接，现在要一次不重复地走过全部这七座桥
如果我们想将这个问题当作数学问题来解决，我们需要对问题中出现的元素进行一些抽象和简化，以更好地形成数学概念并且应用数学方法解决
欧拉在思考这个问题时，认为七桥问题与岛的形状、大小没有关系，与河岸的形状、长短没有关系，与桥的形状、长短也没有关系，重要的、有关系的是桥与桥、桥与河岸、桥与岛、岛与河岸的位置关系
于是欧拉首先把岛和岸都抽象成点，把桥抽象成线，原来的问题就转化为“不重复地走遍这七条线”
问题可以抽象为点之间这样的连接关系
然后欧拉把哥尼斯堡七桥问题抽象成“一笔画问题”：笔尖不离开纸面，一笔画出给定的点线图，不允许重复过任何一条线（可以重复过点），这样的点线图简称为“一笔画”。现在要解决的问题是：找到“一个点线图是一笔画”的充分必要条件，并对是一笔画的点线图给出一笔画的方法
一笔画问题\#
欧拉发现，能够一笔画的点线图都必须是“连通”的，然后注意到每个点都是若干条线的端点，他把点线图上的点分成两类：
- 如果以某点为端点的线有偶数条，就称此点为偶节点
- 如果以某点为端点的线有奇数条，就称此点为奇节点
因此如果我们要不重复地走完七座桥，那么展开的拓扑关系应该是这样的：
我们可以把其中的点分为两类：端点和中间点
对于端点，我们可以发现，连接数有两种情况：作为起点时流出，作为终点时流入；同时作为起点和终点时，流入和流出。因此端点可以存在有奇数个点连接的情况，但是端点只能存在两个
对于中间点，我们发现必然有流入和流出，因此中间点连接的点必然为偶数个
要想不重复地一笔画出某点线图，除去起始点和终止点两个点以外，其余每个点，如果画进去一条线，就一定要画出来一条线，从而都必须时偶节点，于是一笔画的必要条件是：点线图中的奇结点个数为0或2（当起始点和终止点重合时为0）；反之，如果点线图中的奇结点个数为0或2时，就一定能完成一笔画。当点线图中有两个奇结点时，以其中一个为起始点，另一个为终止点，就能完成一笔画。当点线图中没有奇结点时，从任何一个点起始都可以完成一笔画。
这样，欧拉就得到了点线图是一笔画的充要条件：点线图中的奇结点个数为0或2
现在我们再看哥尼斯堡七桥问题，我们发现图中的4个结点全都是奇结点，显然不符合一笔画的条件
欧拉公式\#
对于平面，有： $$ V-E+F=1 $$ 其中，\(V\)为网络的结点个数，\(E\)为网络的边数，\(F\)为面数
对于一个网络，我们可以考虑从边缘入手，假如有这样一个外边\(AB\)，我们将它去掉，显然这时\(E\)减少了\(1\)，\(F\)也减少了\(1\)，而\(V\)则保持不变，因此我们可以反复进行这样的操作，因为外边总是有的，我们可以发现\(V-E+F\)的值保持不变，最后直到得到一张只有一个结点的网络，这时\(V =1, E=F=0\)，我们就得到了这样一个等式
同理我们还可以得到多面体中\(V-E+F =2\)
继续挖掘下去我们还能发现更多有关图的性质
传球问题\#
在正式介绍哈密顿问题之前，我们先来看一看传球问题：
A, B, C, D四个人传球6次，从A开始，最终回到A手里，有多少种传法？
先用动态规划的思路来看看
动态规划\#
定义一个递推数组\(dp[i][j]\)，这个数组用来表示在第\(i\)次传球后，球传到第\(j\)个人手里的不同传球方式的数量，\(j\)的取值可以为\(A,B,C,D\)，这里我们为了方便，用数字\(1,2,3,4\)来表示
为了方便理解，我们用更数学的语言来描述一下我们现在有的条件：
- 如果球在某个人手里，下一步可以选择传给其他任意三个人中的一个
- 初始状态：从\(A\)开始传
- 最终结果：第6次传球后球回到\(A\)手中，求传法的数量
初始化：
-
\(dp[0][1] = 1\)
-
\(dp[0][1],dp[0][2],dp[0][3]=0\)
状态转移：
对于每一次传球\(i\)，每一个人\(j\)，如果球在\(j\)手中，他就可以传给\(k\)，且\(k \ne j\)： $$ dp[i][j]+=dp[i-1][j] $$ 有了递推公式，我们可以初始化一个递推表：
# 初始化dp数组，因为传球6次+初始状态，总共有7个阶段，4个状态(A,B,C,D)
dp = [[0]* 4 for _ in range(7)]
dp[0][0] = 1
# 填充dp表
for i in range(1,7):
for j in range(4):
for k in range(4):
if j != k:
dp[i][k] += dp[i-1][j]
print("\n----------------------\n",dp[6][0])
邻接矩阵\#
上面的方法奏效，但是不是很直观，并且让人来算递推有些折磨，采用邻接矩阵的方式可以更加直观地计算
对于4人传球问题，我们可以视为一个无向图：
对于这个无向图，我们可以采用矩阵的方式来描述： $$ A = \begin{pmatrix} 0&1&1&1 \ 1&0&1&1 \ 1&1&0&1 \ 1&1&1&0 \end{pmatrix} $$ 像这样，其中所有非对角线元素为1，对角线元素为0的矩阵，被称为完全图的邻接矩阵，我们希望将这个矩阵连续自乘来进行状态转移，或者也可以叫迭代
特征值分解\#
对于这类特定的矩阵，其特性可以用于简化矩阵次方的计算，根据特征方程： $$ \det (A -\lambda I) = 0 $$ 求解得到特征值：\(\lambda_1 = 3,\lambda_2 =\lambda_3=\lambda_4 -1\)，至此我们得到特征值矩阵\(D\) $$ D = \begin{pmatrix} -3&0&0&0 \ 0&1&0&0 \ 0&0&1&0 \ 0&0&0&1 \end{pmatrix} $$ 根据特征向量的定义 $$ A\overrightarrow{v} = \lambda \overrightarrow{v} \\ (A - \lambda I)\boldsymbol{v} = 0 $$
我们先计算\(\lambda = 3\)时的特征向量：
根据对称性，很容易得到特征向量\(x=y=z=w\)，因此特征向量可以取任意倍数的\((1,1,1,1)^T\)，在这里我们需要将它规范化成单位向量\(\frac{1}{2}(1,1,1,1)^T\)
对于特征值\(\lambda = -1\)，带入我们发现矩阵所有元素全都是\(1\)，显然这是一个有3个自由度的矩阵，特征向量可以选择任何正交于\((1,1,1,1)^T\)的向量：
- \((1,-1,0,0)\)
- \((1,0,-1,0)\)
- \((1,0,0,-1)\)
现在我们就得到了一个特征矩阵： $$ P = \begin{pmatrix} \frac{1}{2}&1&1&1 \ \frac{1}{2}&-1&0&0 \ \frac{1}{2}&0&-1&0 \ \frac{1}{2}&0&0&-1 \end{pmatrix} $$
我们使用特征值分解：
至此，我们就可以快速算出\(A^k\)了
P = np.array([[1/2,1,1,1],
[1/2,-1,0,0],
[1/2,0,-1,0],
[1/2,0,0,-1]])
P_inv = np.linalg.inv(P)
D = np.diag([3,-1,-1,-1])
D_6 = np.diag([3**6,(-1)**6,(-1)**6,(-1)**6])
A_6_P = P @ D_6 @ P_inv
print(A_6_P)
虽然我们现在能够得到正确的结果，但是这个算法貌似还有可以优化的余地，并且我们还可以更深入地挖掘算法中的细节
正交化和归一化\#
我们先来看一下正交矩阵(Orthogonal Matrix)，首先正交矩阵是一个方阵 $$ Q^TQ = QQ^T = I_{n\times n} $$ Where \(Q^T\) is the transpose of \(Q\) and \(I\) is the identity matrix, so a matrix \(Q\) is orthogonal if its transpose is equal to its inverse: $$ Q^T = Q^{-1} $$ 因此如果我们一开始就将向量正交化，在大矩阵中就不必计算特征矩阵的逆了(正交矩阵还在SVD奇异值分解中有重要作用) $$ A^k = PDkPT $$ 那么什么是向量正交呢？我在学习傅立叶级数的时候就碰到过三角函数的正交性，其中也用到了向量的正交，在几何上，正交的向量可以认为是相互垂直的，即点乘为\(0\)，这也表示正交向量是线性无关的，各个向量在各自的方向上不互相依赖，在更高维的空间中，正交向量提供了一种清晰地描述空间和分析空间结构的方法。它们定义了空间中的“轴”，使得不同的轴是互相独立的，从而简化了很多的数学上的分析
施密特正交化 $$ 不正交的{\boldsymbol{v_1},\boldsymbol{v_2},\cdots,\boldsymbol{v_n}} \Longrightarrow 正交归一的{\boldsymbol{e_1},\boldsymbol{e_2},\cdots, \boldsymbol{e_n}} $$ 而施密特正交化的这个变换，仅仅相当于对矩阵进行了线性变换，前后向量表达的含义是等价的
我们首先进行正交化：
对于第一个向量\(\boldsymbol{u_1} = \boldsymbol{v_1}\)，然后我们对后面的向量\(\boldsymbol{v_k}\)进行正交化，这里我们得注意\(2<k<n\)
首先得计算\(\boldsymbol{v_k}\)在每个\(\boldsymbol{u_i}\)上的投影: $$ \mathrm{projection}{\boldsymbol{u_i}}(\boldsymbol{v_k}) = \frac{\boldsymbol{u_i}·\boldsymbol{v_k}}{\boldsymbol{u_i}·\boldsymbol{u_i}}\boldsymbol{u_i} $$ 更新\(\boldsymbol{v_k}\)：我们从\(\boldsymbol{v_k}\)中减去这些投影 $$ \boldsymbol{u_k} = \boldsymbol{v_k} - \sum $$ 这里的重点是理解正交化是怎么实现的，我们发现更新的向量，它们会带上在它们之前计算的向量的投影，我们可以看看： $$ \begin{align} \boldsymbol{u_2} &= \boldsymbol{v_2} -\mathrm{proj}\boldsymbol{v_1}\tag{1}\ \boldsymbol{u_3} &= \boldsymbol{v_3} -\mathrm{proj}\boldsymbol{v_1}-\mathrm{proj}\boldsymbol{v_2} \tag{2} \end{align} $$ 对于(1)，等号右边是向量减法，如果我们将这两个向量看成是二维向量，我们很容易得到做了向量减法后的向量必然垂直于原向量}^{k-1}\mathrm{projection}_{\boldsymbol{u_i}}\boldsymbol{v_k\(\boldsymbol{u_1}\)，而二维平面中的垂直关系，在三维空间也显然成立，推广到无穷维也都是成立的
对于(2)，我们可以认为这三个向量是三维向量，那么\(\boldsymbol{v_3} -\mathrm{proj}\boldsymbol{v_1}\)只能保证\(\boldsymbol{u_3}\)与\(\boldsymbol{u_1}\)的垂直关系，不足以保证\(\boldsymbol{u_3}\)与\(\boldsymbol{u_2}\)的垂直关系，而只需要对投影做向量减法就能保证向量的垂直…因此我们对后来的向量，不断减掉前面向量投影的累加即可一直保证新的向量能够垂直于前面的所有向量，至此，我们就完成了所有向量的正交化
归一化
对于每个\(\boldsymbol{u}_k\)，让它们的模长都变成\(1\) $$ \boldsymbol{e}_k = \frac{\boldsymbol{u}_k}{||\boldsymbol{u}_k||} $$ 这样做有一个好处，当我们在进行计算的时候，正交化让所有向量的长度都相等，因此出现的误差带来的影响也是一样的，以免误差分布不均
在python中进行施密特正交化：
import numpy as np
# 重新定义原始矩阵 P
P = np.array([
[0.5, 1, 1, 1],
[0.5, -1, 0, 0],
[0.5, 0, -1, 0],
[0.5, 0, 0, -1]
])
# 定义一个矩阵来存储正交化后的向量
U = np.zeros_like(P, dtype=float)
# 对第一个向量不做改变，直接归一化
U[:, 0] = P[:, 0] / np.linalg.norm(P[:, 0])
# 施密特正交化过程
for i in range(1, P.shape[1]):
orthogonal = P[:, i]
for j in range(i):
orthogonal -= np.dot(U[:, j], P[:, i]) / np.dot(U[:, j], U[:, j]) * U[:, j]
U[:, i] = orthogonal / np.linalg.norm(orthogonal)
print(U)
运行我们可以得到：
[[ 0.5 0.70710678 0.40824829 0.28867513]
[ 0.5 -0.70710678 0.40824829 0.28867513]
[ 0.5 0. -0.81649658 0.28867513]
[ 0.5 0. 0. -0.8660254 ]]
至此，对于邻接矩阵的特征值分解计算，我们基本上已经找到了比较理想的算法了，我们将循环重复不断计算大矩阵的运算变成了只需计算一次大矩阵，对角矩阵和向量计算相较于一般矩阵，复杂度大大降低
哈密顿路径问题\#
哈密顿问题属于NP完全问题，，可以在多项式时间内验证所提出的解决方案，是旅行商问题的一种特例
我们来看一个具体的哈密顿路径问题，现在在一个\(6\times 3\)的网络中，我们要不重复地经过所有的结点，那么有多少种画法？这个问题跟七桥问题不同的地方在于，七桥问题是可以重复经过结点，而哈密顿问题则不允许重复经过结点，我们对比刚才的图，可以猜到，哈密顿问题的条件更为严格，，因为对于两个结点，你不能用两条及以上的边同时连接这两个结点
对于哈密顿问题，我们也可以采用类似的方法来构建邻接矩阵，由于上面的性质，在哈密顿问题中，邻接矩阵的元素也只有0和1
import time
def construct_adjacency_matrix(rows, cols):
num_vertices = rows * cols
adj_matrix = [[0] * num_vertices for _ in range(num_vertices)]
for i in range(rows):
for j in range(cols):
index = i * cols + j
# 上
if i > 0:
adj_matrix[index][index - cols] = 1
# 下
if i < rows - 1:
adj_matrix[index][index + cols] = 1
# 左
if j > 0:
adj_matrix[index][index - 1] = 1
# 右
if j < cols - 1:
adj_matrix[index][index + 1] = 1
return adj_matrix
def find_paths(matrix, path, results):
if len(path) == len(matrix):
results.append(path.copy())
return
if not path:
start_points = range(len(matrix))
else:
start_points = [i for i, val in enumerate(matrix[path[-1]]) if val == 1 and i not in path]
for next_point in start_points:
path.append(next_point)
find_paths(matrix, path, results)
path.pop()
def main():
time1 = time.time()
rows, cols = 8, 4
adj_matrix = construct_adjacency_matrix(rows, cols)
results = []
for start in range(rows * cols):
find_paths(adj_matrix, [start], results)
time2 = time.time()
print(f"Time taken: {time2 - time1} seconds")
print(f"Total Hamiltonian paths found: {len(results)}")
return results
results = main()
在哈密顿问题中，复杂度是一个很困难的问题，用上面的python代码寻找一个\(8\times 4\)的网络的总路线：
可以发现，8w条不到的路线，竟然花了87s，网络复杂度每增加1，花费的时间都以指数级增长