[图] - 所有结点对的最短路径问题(Floyd Warshall)

单源最短路径

我们可以用已经知道的算法，来对所有点计算一次最短路径

无权重的图

使用$V$次BFS算法，时间复杂度是$O(V^2+VE)$

非负权重的图

使用$V$次$Dijkstra$算法

• 线性数组：$O(V^3+VE)=O(V^3)$
• 二叉堆：$O(VElgV)$
• 斐波那契堆：$O(V^2lgV+VE)$

含有负权重的图

使用$V$次$Bellman-Ford$算法，时间复杂度是$O(V^4)$

看时间复杂度，以上的方法都不是很好，所以我们需要更好的方法

Floyd-Warshall

Floyd-Warshall的本质是使用动态规划，其步骤为

1. 分析最优解结构
2. 递归定义最优解的值
3. 自底向上计算最优解的值
4. 根据计算出的最优解的值构造最优解

此算法可以计算有权重为负的边，但是不能计算存在负权重的环路。

流程

我们假定一个图$G=(E, V)$，其中$V=\lbrace1,2,3,\cdots,n\rbrace$, 我们选取任意一个子集$V'$，使得

$V' = \lbrace1,2,3,...,k\rbrace$

$V' \in V$

如图，蓝色的集合就是$V'$

对于任意结点$i,j \in V$，我们考虑所有从结点$i$到结点$j$的所有中间结点都$\in V'$的路径
这里的中间结点指的是：简单路径$p=\lbrace n_1,n_2,n_3,\cdots,n_k\rbrace$上除结点$v_1, v_k$结点以外点任意结点。

这里，我们设$p$为这些路径中权重最小的一条边，我们需要通过$p$和集合$\lbrace1,2,3,\cdots,k-1\rbrace$的关系，来寻找对于所有结点点最短路径。

总共分2种情况，我们先定义一个变量$d^{(k)}_{ij}$：表示从结点$i$到$j$的所有中间结点均$\in \lbrace1,2,3,\cdots,k\rbrace$的一条最短路径的权重。

结点$k$不是路径$p$的中间结点

现在，路径$p$的中间结点都属于集合$\lbrace1,2,3,\cdots,k-1\rbrace$，就可以推出

$d^{(k)}_{ij} = d^{(k-1)}_{ij}$

结点$k$是路径$p$的中间结点

现在，路径$p$的被中间结点$k$划分为$p_1, p_2$ 2条路径，其中每条路径的结点都$\in \lbrace 1, 2, 3, \cdots, k - 1 \rbrace$，就可以推出

$d^{(k)}_{ij} = d^{(k-1)}_{ik} + d^{(k-1)}_{kj}$

递归定义最优解

根据上方，我们可以定义出$d^{(k)}_{ij}$与$d^{(k-1)}_{ij}$的关系，从而得到一个递归解

$d^{(k)}_{ij} = \begin{cases} \varpi_{i,j} & k = 0 \\ min(d^{(k-1)}_{ij}, d^{(k-1)}_{ik} + d^{(k-1)}_{kj}) & k >= 1 \\ \end{cases}$

我们还规定

$\varpi_{i,j} = \begin{cases} 0 & i = j \\ \varpi_{i,j} & i \neq j \\ \infty & i \text{与} j \text{不连通} \\ \end{cases}$

这里$k >= 1$的情况，通俗点说明就是当路径$p$的被中间结点$k$划分为$p_1, p_2$的2条路径时，且这2条路径加起来比$d^{(k-1)}_{ij}$的权重小的话，那么就选$d^{(k-1)}_{ik} + d^{(k-1)}_{kj}$，大的话就选$d^{(k-1)}_{ij}$。要注意，这里的$k$，我们是不知道在哪里的，需要对集合$k \in V$进行遍历。

C语言代码

int **floyd_warshall(int graph[][MAX_NODE], int length)
{
    int n = length;
    int **L = (int **)malloc(sizeof(int *) * n);
    for (int i = 0; i < n; i++)
        L[i] = (int *)malloc(sizeof(int) * n);

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            L[i][j] = graph[i][j];

    for (int k = 0; k < n; k++)
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                L[i][j] = min(L[i][j], L[i][k] + L[k][j]);
    return L;
}

代码的3-6行，我们创建一个二维数组L[i][j]，来记录当前$i$与$j$的路径的长度，而代码的12-15行，我们对每一个k结点进行松弛（relax）操作，如果通过k结点到达i更短，那么就更新当前路径权重。

关于k的循环的顺序问题

Floyd-Warshall里面有个定理：

定理： 假设$i$和$j$之间的最短路径上的结点集里(不包含$i$,$j$)，编号最大的一个是$x$.那么在外循环$k=x$时,L[i][j]肯定得到了最小值。

我们可以通过这张图来看，

1. 首先是k=1，通过点1并不会有什么改变，所以什么都不会变
2. 然后是k=2，通过点2，我们会发现L[1][3]取得了最小值，此时L[1][2]和L[2][3]已经是最小值了。
3. 然后是k=3，通过点3，我们会发现L[1][4]取得了最小值=此时L[1][2]和L[2][3]和L[3][4]已经是最小值了。
4. 然后是k=4，通过点4，我们会发现L[1][5]取得了最小值，此时L[1][2]和L[2][3]和L[3][4]和L[4][5]已经是最小值了。

即使把5和2或者其他任意的顺序交换也会得到同样的结果。所以对k的循环必须在最外层。这也就是为什么本来是3层循环，应该是三维数组来储藏，但是最后却只用了二维数组就做到了的原因，因为上一步计算出点最小值，被用于了下一步，而不需要重新再计算。这也就是动态规划的基本。

求最短路径

上面我们只是把最短路径点权重计算出来了，但是我们还不知道最短路径具体是什么，通过关于k的循环的顺序问题，我们知道了，当我们在本次通过x取得L[1][k]的最小值的时候，L[1][x]已经取得了最小值。所以我们可以记录前驱结点，也就是k的前驱结点就是x，同理，x也有自己对应的前驱结点，这样递归点找下去，即可找到最短路径。

于是

$k=0$

$\pi^{(0)}_{i,j} = \begin{cases} NIL & i = j\ or\ \varpi_{i,j} = \infty \\ i & i \neq j\ and\ \varpi_{i,j} \neq \infty \\ \end{cases}$

$k>0$

$\pi^{(k)}_{i,j} = \begin{cases} \pi^{(k-1)}_{i,j} & d^{(k)}_{ij} = d^{(k-1)}_{ij} \\ \pi^{(k-1)}_{k,j} & d^{(k)}_{ij} = d^{(k-1)}_{ik} + d^{(k-1)}_{kj} \\ \end{cases}$

求最短路径的C语言代码

int **floyd_warshall(int graph[][MAX_NODE], int length)
{
    int n = length;
    int **L = (int **)malloc(sizeof(int *) * n);
    int **P = (int **)malloc(sizeof(int *) * n);
    for (int i = 0; i < n; i++)
    {
        L[i] = (int *)malloc(sizeof(int) * n);
        P[i] = (int *)malloc(sizeof(int) * n);
    }

    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            L[i][j] = graph[i][j];
            P[i][j] = j;
        }
    }
    for (int k = 0; k < n; k++)
    {
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                if (L[i][j] > L[i][k] + L[k][j])
                {
                    L[i][j] = L[i][k] + L[k][j];
                    P[i][j] = k;
                }
            }
        }
    }
    return P;
}

矩阵乘法

我们来看看矩阵乘法的代码

void multiply(int mat1[][N], int mat2[][N], int res[][N])
{
    int i, j, k;
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            res[i][j] = 0;
            for (k = 0; k < N; k++)
                res[i][j] += mat1[i][k] * mat2[k][j];
        }
    }
}

可以看出，上方的求最短路径点代码和矩阵点乘法非常相似。那么矩阵乘法也可以得到应用。假设我们有了恰好经过$k-1$条边的方案数矩阵$B \in \mathbb{R}^{n\times n}$，还有图的临接矩阵$A \in \mathbb{R}^{n\times n}$，那么$B \times A$就是恰好经过 K 条边互相可达的方案数矩阵

定理： 设$G = \langle V, E,\psi \rangle$是一简单有向图，并且$A$是$G$的临接矩阵，对于$m=1,2,3,\cdots$来说，矩阵$A^m$中的元素$a_{ij}$的值，等于从$v_i$到$v_j$长度为m的路径数目。

证：对于$m$进行归纳证明。当$m=1$时，很明显$A=A$
设矩阵$A^l$中的元素$(i, j)$值是$a^l_{ij}$，且对于$m=1$来说结论为真。
因为$A^{l+1} = A^lA$，所以应有

$a^{l+1}_{ij}=\sum^n_{k=1}a^l_{ik}a_{kj}$

这里$v_k$是倒数第二个结点，因此，$a^{l+1}_{ij}$应是从结点$v_i$出发，经过任意的倒数第二个结点到$v_j$的长度为$l+1$的路径总数。

我们把矩阵乘法用在求恰好经过$K$条边最短路径上

void multiply(int A[][N], int Al[][N], int res[][N])
{
    int i, j, k;
    for (i = 0; i < N; i++)
    {
        for (j = 0; j < N; j++)
        {
            res[i][j] = 0;
            for (k = 0; k < N; k++)
                res[i][j] = Al[i][k] + A[k][j];
        }
    }
}

其中Al[i][k]就相当与$a^{l}_{ik}$，那么$a^{l}_{ik} + a_{kj} = a^{l+1}_{ij}$，但是请注意，这里是从假设我们有了恰好经过$k-1$条边的方案数矩阵$B \in \mathbb{R}^{n\times n}$开始的，而当前这一步的时间复杂度就达到了$O(n^3)$，由于要求$n$个点，所以总的时间复杂度达到$O(n^4)$。

重复平方

我们对$A^m, m \in R$的结果并不感兴趣，因为对于$n$个结点点图，两点之间的最短距离最大只有$n-1$条边，所以我们真正关心的是$A^{n-1}$，我们可以通过重复平方的方法来改良当前的算法

$\begin{aligned} A &= A \\ A^2 &= A * A \\ A^4 &= A^2 * A^2 \\ A^8 &= A^4 * A^4 \\ &\vdots \\ A^{2lg(n-1)} &= A^{2lg(n-1)} * A^{2lg(n-1) - 1} \end{aligned}$

由于$2lg(n-1) \geq n-1$，所以$A^{2lg(n-1)} = A^{n-1}$，这样优化后时间复杂度从$O(n^4)$降到$O(n^3lgn)$，和Floyd Warshall（$O(N^3)$）还是有很大的差距的。

总结

Floyd Warshall利用动态规划，可以在$O(N^3)$就能求得所有结点对的最短路径。
而通过临接矩阵相乘来计算，时间复杂度比Floyd Warshall高得多，既是优化后也提升不大，但是其中的思想还是可以学习的。