A*（念做：A Star）算法是一种很常用的路径查找和图形遍历算法。它有较好的性能和准确度。本文在讲解算法的同时也会提供Python语言的代码实现，并会借助matplotlib库动态的展示算法的运算过程。

• 算法介绍
• 广度优先搜索
• Dijkstra算法
• 最佳优先搜索
• A*算法
  • 启发函数
  • 关于距离
    • 曼哈顿距离
    • 对角距离
    • 欧几里得距离
• 算法实现
  • 坐标点与地图
  • 算法主体
  • 测试入口
• 算法变种
  • ARA*
  • D*
  • Field D*
  • Block A*
• 参考资料与推荐读物

算法介绍

A*算法最初发表于1968年，由Stanford研究院的Peter Hart, Nils Nilsson以及Bertram Raphael发表。它可以被认为是Dijkstra算法的扩展。

由于借助启发函数的引导，A*算法通常拥有更好的性能。

广度优先搜索

为了更好的理解A*算法，我们首先从广度优先（Breadth First）算法讲起。

正如其名称所示，广度优先搜索以广度做为优先级进行搜索。

从起点开始，首先遍历起点周围邻近的点，然后再遍历已经遍历过的点邻近的点，逐步的向外扩散，直到找到终点。

这种算法就像洪水（Flood fill）一样向外扩张，算法的过程如下图所示：

在上面这幅动态图中，算法遍历了图中所有的点，这通常没有必要。对于有明确终点的问题来说，一旦到达终点便可以提前终止算法，下面这幅图对比了这种情况：

在执行算法的过程中，每个点需要记录达到该点的前一个点的位置 – 可以称之为父节点。这样做之后，一旦到达终点，便可以从终点开始，反过来顺着父节点的顺序找到起点，由此就构成了一条路径。

Dijkstra算法

Dijkstra算法是由计算机科学家Edsger W. Dijkstra在1956年提出的。

Dijkstra算法用来寻找图形中节点之间的最短路径。

考虑这样一种场景，在一些情况下，图形中相邻节点之间的移动代价并不相等。例如，游戏中的一幅图，既有平地也有山脉，那么游戏中的角色在平地和山脉中移动的速度通常是不相等的。

在Dijkstra算法中，需要计算每一个节点距离起点的总移动代价。同时，还需要一个优先队列结构。对于所有待遍历的节点，放入优先队列中会按照代价进行排序。

在算法运行的过程中，每次都从优先队列中选出代价最小的作为下一个遍历的节点。直到到达终点为止。

下面对比了不考虑节点移动代价差异的广度优先搜索与考虑移动代价的Dijkstra算法的运算结果：

当图形为网格图，并且每个节点之间的移动代价是相等的，那么Dijkstra算法将和广度优先算法变得一样。

最佳优先搜索

在一些情况下，如果我们可以预先计算出每个节点到终点的距离，则我们可以利用这个信息更快的到达终点。

其原理也很简单。与Dijkstra算法类似，我们也使用一个优先队列，但此时以每个节点到达终点的距离作为优先级，每次始终选取到终点移动代价最小（离终点最近）的节点作为下一个遍历的节点。这种算法称之为最佳优先（Best First）算法。

这样做可以大大加快路径的搜索速度，如下图所示：

但这种算法会不会有什么缺点呢？答案是肯定的。

因为，如果起点和终点之间存在障碍物，则最佳优先算法找到的很可能不是最短路径，下图描述了这种情况。

A*算法

对比了上面几种算法，最后终于可以讲解本文的重点：A*算法了。

下面的描述我们将看到，A*算法实际上是综合上面这些算法的特点于一身的。

A*算法通过下面这个函数来计算每个节点的优先级。

\[f(n) = g(n) + h(n)\]

其中：

• $f(n)$ 是节点n的综合优先级。当我们选择下一个要遍历的节点时，我们总会选取综合优先级最高（值最小）的节点。
• $g(n)$ 是节点n距离起点的代价。
• $h(n)$ 是节点n距离终点的预计代价，这也就是A*算法的启发函数。关于启发函数我们在下面详细讲解。

A*算法在运算过程中，每次从优先队列中选取$f(n)$值最小（优先级最高）的节点作为下一个待遍历的节点。

另外，A*算法使用两个集合来表示待遍历的节点，与已经遍历过的节点，这通常称之为open_set和close_set。

完整的A*算法描述如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

* 初始化open_set和close_set；
* 将起点加入open_set中，并设置优先级为0（优先级最高）；
* 如果open_set不为空，则从open_set中选取优先级最高的节点n：
    * 如果节点n为终点，则：
        * 从终点开始逐步追踪parent节点，一直达到起点；
        * 返回找到的结果路径，算法结束；
    * 如果节点n不是终点，则：
        * 将节点n从open_set中删除，并加入close_set中；
        * 遍历节点n所有的邻近节点：
            * 如果邻近节点m在close_set中，则：
                * 跳过，选取下一个邻近节点
            * 如果邻近节点m也不在open_set中，则：
                * 设置节点m的parent为节点n
                * 计算节点m的优先级
                * 将节点m加入open_set中

启发函数

上面已经提到，启发函数会影响A*算法的行为。

• 在极端情况下，当启发函数$h(n)$始终为0，则将由$g(n)$决定节点的优先级，此时算法就退化成了Dijkstra算法。
• 如果$h(n)$始终小于等于节点n到终点的代价，则A*算法保证一定能够找到最短路径。但是当$h(n)$的值越小，算法将遍历越多的节点，也就导致算法越慢。
• 如果$h(n)$完全等于节点n到终点的代价，则A*算法将找到最佳路径，并且速度很快。可惜的是，并非所有场景下都能做到这一点。因为在没有达到终点之前，我们很难确切算出距离终点还有多远。
• 如果$h(n)$的值比节点n到终点的代价要大，则A*算法不能保证找到最短路径，不过此时会很快。
• 在另外一个极端情况下，如果$h(n)$相较于$g(n)$大很多，则此时只有$h(n)$产生效果，这也就变成了最佳优先搜索。

由上面这些信息我们可以知道，通过调节启发函数我们可以控制算法的速度和精确度。因为在一些情况，我们可能未必需要最短路径，而是希望能够尽快找到一个路径即可。这也是A*算法比较灵活的地方。

对于网格形式的图，有以下这些启发函数可以使用：

• 如果图形中只允许朝上下左右四个方向移动，则可以使用曼哈顿距离（Manhattan distance）。
• 如果图形中允许朝八个方向移动，则可以使用对角距离。
• 如果图形中允许朝任何方向移动，则可以使用欧几里得距离（Euclidean distance）。

关于距离

曼哈顿距离

如果图形中只允许朝上下左右四个方向移动，则启发函数可以使用曼哈顿距离，它的计算方法如下图所示：

计算曼哈顿距离的函数如下，这里的D是指两个相邻节点之间的移动代价，通常是一个固定的常数。

1
2
3
4

function heuristic(node) =
    dx = abs(node.x - goal.x)
    dy = abs(node.y - goal.y)
    return D * (dx + dy)

对角距离

如果图形中允许斜着朝邻近的节点移动，则启发函数可以使用对角距离。它的计算方法如下：

计算对角距离的函数如下，这里的D2指的是两个斜着相邻节点之间的移动代价。如果所有节点都正方形，则其值就是\\sqrt{2} \* D。

1
2
3
4

function heuristic(node) =
    dx = abs(node.x - goal.x)
    dy = abs(node.y - goal.y)
    return D * (dx + dy) + (D2 - 2 * D) * min(dx, dy)

欧几里得距离

如果图形中允许朝任意方向移动，则可以使用欧几里得距离。

欧几里得距离是指两个节点之间的直线距离，因此其计算方法也是我们比较熟悉的：\\sqrt{(p2.x-p1.x)^2 + (p2.y-p1.y)^2}。其函数表示如下：

1
2
3
4

function heuristic(node) =
    dx = abs(node.x - goal.x)
    dy = abs(node.y - goal.y)
    return D * sqrt(dx * dx + dy * dy)

算法实现

虽然前面介绍了很多内容，但实际上A*算法并不复杂，实现起来也比较简单。

下面我们给出一个Python语言的代码示例。

之所以使用Python语言是因为我们可以借助matplotlib库很方便的将结果展示出来。在理解了算法之后，通过其他语言实现也并非难事。

我们的算法演示的是在一个二维的网格图形上从起点找寻终点的求解过程。

坐标点与地图

首先，我们创建一个非常简单的类来描述图中的点，相关代码如下：

1
2
3
4
5
6
7
8
9

# point.py

import sys

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y
        self.cost = sys.maxsize

接着，我们实现一个描述地图结构的类。为了简化算法的描述:

我们选定左下角坐标[0, 0]的点是算法起点，右上角坐标[size - 1, size - 1]的点为要找的终点。

为了让算法更有趣，我们在地图的中间设置了一个障碍，并且地图中还会包含一些随机的障碍。该类的代码如下：

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

# random_map.py

import numpy as np

import point

class RandomMap:
    def __init__(self, size=50): ①
        self.size = size
        self.obstacle = size//8 ②
        self.GenerateObstacle() ③

    def GenerateObstacle(self):
        self.obstacle_point = []
        self.obstacle_point.append(point.Point(self.size//2, self.size//2))
        self.obstacle_point.append(point.Point(self.size//2, self.size//2-1))

        # Generate an obstacle in the middle
        for i in range(self.size//2-4, self.size//2): ④
            self.obstacle_point.append(point.Point(i, self.size-i))
            self.obstacle_point.append(point.Point(i, self.size-i-1))
            self.obstacle_point.append(point.Point(self.size-i, i))
            self.obstacle_point.append(point.Point(self.size-i, i-1))

        for i in range(self.obstacle-1): ⑤
            x = np.random.randint(0, self.size)
            y = np.random.randint(0, self.size)
            self.obstacle_point.append(point.Point(x, y))

            if (np.random.rand() > 0.5): # Random boolean ⑥
                for l in range(self.size//4):
                    self.obstacle_point.append(point.Point(x, y+l))
                    pass
            else:
                for l in range(self.size//4):
                    self.obstacle_point.append(point.Point(x+l, y))
                    pass

    def IsObstacle(self, i ,j): ⑦
        for p in self.obstacle_point:
            if i==p.x and j==p.y:
                return True
        return False

这段代码说明如下：

1. 构造函数，地图的默认大小是50x50；
2. 设置障碍物的数量为地图大小除以8；
3. 调用GenerateObstacle生成随机障碍物；
4. 在地图的中间生成一个斜着的障碍物；
5. 随机生成其他几个障碍物；
6. 障碍物的方向也是随机的；
7. 定义一个方法来判断某个节点是否是障碍物；

算法主体

有了基本的数据结构之后，我们就可以开始实现算法主体了。

这里我们通过一个类来封装我们的算法。

首先实现一些算法需要的基本函数，它们如下：

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

# a_star.py

import sys
import time

import numpy as np

from matplotlib.patches import Rectangle

import point
import random_map

class AStar:
    def __init__(self, map):
        self.map=map
        self.open_set = []
        self.close_set = []

    def BaseCost(self, p):
        x_dis = p.x
        y_dis = p.y
        # Distance to start point
        return x_dis + y_dis + (np.sqrt(2) - 2) * min(x_dis, y_dis)

    def HeuristicCost(self, p):
        x_dis = self.map.size - 1 - p.x
        y_dis = self.map.size - 1 - p.y
        # Distance to end point
        return x_dis + y_dis + (np.sqrt(2) - 2) * min(x_dis, y_dis)

    def TotalCost(self, p):
        return self.BaseCost(p) + self.HeuristicCost(p)

    def IsValidPoint(self, x, y):
        if x < 0 or y < 0:
            return False
        if x >= self.map.size or y >= self.map.size:
            return False
        return not self.map.IsObstacle(x, y)

    def IsInPointList(self, p, point_list):
        for point in point_list:
            if point.x == p.x and point.y == p.y:
                return True
        return False

    def IsInOpenList(self, p):
        return self.IsInPointList(p, self.open_set)

    def IsInCloseList(self, p):
        return self.IsInPointList(p, self.close_set)

    def IsStartPoint(self, p):
        return p.x == 0 and p.y ==0

    def IsEndPoint(self, p):
        return p.x == self.map.size-1 and p.y == self.map.size-1

这里的函数说明如下：

• __init__：类的构造函数。
• BaseCost：节点到起点的移动代价，对应了上文的$g(n)$。
• HeuristicCost：节点到终点的启发函数，对应上文的$h(n)$。由于我们是基于网格的图形，所以这个函数和上一个函数用的是对角距离。
• TotalCost：代价总和，即对应上面提到的$f(n)$。
• IsValidPoint：判断点是否有效，不在地图内部或者障碍物所在点都是无效的。
• IsInPointList：判断点是否在某个集合中。
• IsInOpenList：判断点是否在open_set中。
• IsInCloseList：判断点是否在close_set中。
• IsStartPoint：判断点是否是起点。
• IsEndPoint：判断点是否是终点。

有了上面这些辅助函数，就可以开始实现算法主逻辑了，相关代码如下：

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

# a_star.py
def RunAndSaveImage(self, ax, plt):
    start_time = time.time()

    start_point = point.Point(0, 0)
    start_point.cost = 0
    self.open_set.append(start_point)

    while True:
        index = self.SelectPointInOpenList()
        if index < 0:
            print('No path found, algorithm failed!!!')
            return
        p = self.open_set[index]
        rec = Rectangle((p.x, p.y), 1, 1, color='c')
        ax.add_patch(rec)
        self.SaveImage(plt)

        if self.IsEndPoint(p):
            return self.BuildPath(p, ax, plt, start_time)

        del self.open_set[index]
        self.close_set.append(p)

        # Process all neighbors
        x = p.x
        y = p.y
        self.ProcessPoint(x-1, y+1, p)
        self.ProcessPoint(x-1, y, p)
        self.ProcessPoint(x-1, y-1, p)
        self.ProcessPoint(x, y-1, p)
        self.ProcessPoint(x+1, y-1, p)
        self.ProcessPoint(x+1, y, p)
        self.ProcessPoint(x+1, y+1, p)
        self.ProcessPoint(x, y+1, p)

这段代码应该不需要太多解释了，它就是根据前面的算法逻辑进行实现。为了将结果展示出来，我们在算法进行的每一步，都会借助于matplotlib库将状态保存成图片。

上面这个函数调用了其他几个函数代码如下：

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

# a_star.py
def SaveImage(self, plt):
    millis = int(round(time.time() * 1000))
    filename = './' + str(millis) + '.png'
    plt.savefig(filename)

def ProcessPoint(self, x, y, parent):
    if not self.IsValidPoint(x, y):
        return # Do nothing for invalid point
    p = point.Point(x, y)
    if self.IsInCloseList(p):
        return # Do nothing for visited point
    print('Process Point [', p.x, ',', p.y, ']', ', cost: ', p.cost)
    if not self.IsInOpenList(p):
        p.parent = parent
        p.cost = self.TotalCost(p)
        self.open_set.append(p)

def SelectPointInOpenList(self):
    index = 0
    selected_index = -1
    min_cost = sys.maxsize
    for p in self.open_set:
        cost = self.TotalCost(p)
        if cost < min_cost:
            min_cost = cost
            selected_index = index
        index += 1
    return selected_index

def BuildPath(self, p, ax, plt, start_time):
    path = []
    while True:
        path.insert(0, p) # Insert first
        if self.IsStartPoint(p):
            break
        else:
            p = p.parent
    for p in path:
        rec = Rectangle((p.x, p.y), 1, 1, color='g')
        ax.add_patch(rec)
        plt.draw()
        self.SaveImage(plt)
    end_time = time.time()
    print('===== Algorithm finish in', int(end_time-start_time), ' seconds')

这三个函数应该是比较容易理解的：

• SaveImage：将当前状态保存到图片中，图片以当前时间命名。
• ProcessPoint：针对每一个节点进行处理：如果是没有处理过的节点，则计算优先级设置父节点，并且添加到open_set中。
• SelectPointInOpenList：从open_set中找到优先级最高的节点，返回其索引。
• BuildPath：从终点往回沿着parent构造结果路径。然后从起点开始绘制结果，结果使用绿色方块，每次绘制一步便保存一个图片。

测试入口

最后是程序的入口逻辑，使用上面写的类来查找路径：

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

# main.py

import numpy as np
import matplotlib.pyplot as plt

from matplotlib.patches import Rectangle

import random_map
import a_star

plt.figure(figsize=(5, 5))

map = random_map.RandomMap() ①

ax = plt.gca()
ax.set_xlim([0, map.size]) ②
ax.set_ylim([0, map.size])

for i in range(map.size): ③
    for j in range(map.size):
        if map.IsObstacle(i,j):
            rec = Rectangle((i, j), width=1, height=1, color='gray')
            ax.add_patch(rec)
        else:
            rec = Rectangle((i, j), width=1, height=1, edgecolor='gray', facecolor='w')
            ax.add_patch(rec)

rec = Rectangle((0, 0), width = 1, height = 1, facecolor='b')
ax.add_patch(rec) ④

rec = Rectangle((map.size-1, map.size-1), width = 1, height = 1, facecolor='r')
ax.add_patch(rec) ⑤

plt.axis('equal') ⑥
plt.axis('off')
plt.tight_layout()
#plt.show()

a_star = a_star.AStar(map)
a_star.RunAndSaveImage(ax, plt) ⑦

这段代码说明如下：

1. 创建一个随机地图；
2. 设置图像的内容与地图大小一致；
3. 绘制地图：对于障碍物绘制一个灰色的方块，其他区域绘制一个白色的的方块；
4. 绘制起点为蓝色方块；
5. 绘制终点为红色方块；
6. 设置图像的坐标轴比例相等并且隐藏坐标轴；
7. 调用算法来查找路径；

由于我们的地图是随机的，所以每次运行的结果可能会不一样，下面是我的电脑上某次运行的结果：

算法变种

A*算法有不少的变种，这里我们介绍最主要的几个。

更多的内容请以访问维基百科：A* Variants。

ARA*

ARA* 全称是Anytime Repairing A*，也称为Anytime A*。

与其他Anytime算法一样，它具有灵活的时间成本，即使在它结束之前被中断，也可以返回路径查找或图形遍历问题的有效解决方案。方法是在逐步优化之前生成快速，非最优的结果。

在现实世界的规划问题中，问题的解决时间往往是有限的。与时间相关的规划者对这种情况都会比较熟悉：他们能够快速找到可行的解决方案，然后不断努力改进，直到时间用完为止。

启发式搜索ARA*算法，它根据可用的搜索时间调整其性能边界。它首先使用松散边界快速找到次优解，然后在时间允许的情况下逐渐收紧边界。如果有足够的时间，它会找到可证明的最佳解决方方案。在改进其约束的同时，ARA*重复使用以前的搜索工作，因此，比其他随时搜索方法更有效。

与A*算法不同，Anytime A*算法最重要的功能是，它们可以被停止，然后可以随时重启。该方法使用控制管理器类来处理时间限制以及停止和重新启动A*算法以找到初始的，可能是次优的解决方案，然后继续搜索改进的解决方案，直到达到可证明的最佳解决方案。

关于ARA*的更多内容可以阅读这篇论文：

• ARA* - Anytime A* with Provable Bounds on Sub-Optimality。

D*

D*是Dynamic A*的简写，其算法和A*类似，不同的是，其代价的计算在算法运行过程中可能会发生变化。

D*包含了下面三种增量搜索算法：

• 原始的D*由Anthony Stentz发表。
• Focussed D*由Anthony Stentz发表，是一个增量启发式搜索算法，结合了A*和原始D*的思想。
• D* Lite是由Sven Koenig和Maxim Likhachev基于LPA*构建的算法。

所有三种搜索算法都解决了相同的基于假设的路径规划问题，包括使用自由空间假设进行规划。在这些环境中，机器人必须导航到未知地形中的给定目标坐标。它假设地形的未知部分（例如：它不包含障碍物），并在这些假设下找到从当前坐标到目标坐标的最短路径。

然后机器人沿着路径行进。当它观察到新的地图信息（例如以前未知的障碍物）时，它会将信息添加到其地图中，并在必要时将新的最短路径从其当前坐标重新添加到给定的目标坐标。它会重复该过程，直到达到目标坐标或确定无法达到目标坐标。在穿越未知地形时，可能经常发现新的障碍，因此重新计划需要很快。增量（启发式）搜索算法通过使用先前问题的经验来加速搜索当前问题，从而加速搜索类似搜索问题的序列。假设目标坐标没有改变，则所有三种搜索算法都比重复的A*搜索更有效。

D*及其变体已广泛用于移动机器人和自动车辆导航。当前系统通常基于D* Lite而不是原始D*或Focussed D*。

关于D*的更多内容可以阅读这两篇文章：

• Project “Fast Replanning （Incremental Heuristic Search）”
• Real-Time Replanning in Dynamic and Unknown Environments

Field D*

Field D*扩展了D*和D* Lite，是一种基于插值（ interpolation-based ）的规划算法，它使用线性插值来有效地生成低成本路径，从而消除不必要的转向。

在给定线性插值假设的情况下，路径是最优的，并且在实践中非常有效。该算法目前被各种现场机器人系统使用。

关于Field D*的详细内容可以看下面这篇论文：

• Field D*: An Interpolation-based Path Planner and Replanner

Block A*

Block A*扩展自A*，但它操作是一块（block）单元而不是单个单元。

其open_set中的每个条目都是已到达但尚未扩展的块，或者需要重新扩展的块。

open_set中块的优先级称为其堆值（heap value）。与A*类似，Block A*中的基本循环是删除具有最低堆值的条目并将其展开。在扩展期间使用LDDB来计算正在扩展的块中的边界单元的g值。

LDDB是一种新型数据库，它包含了本地邻域边界点之间的距离。

关于Block A*的更多内容可以看下面这篇论文：

• Block A*: Database-Driven Search with Applications in Any-angle Path-Planning

参考资料与推荐读物

• Stanford: Introduction to A*
• Wikipedia: A* search algorithm
• PythonRobotics: A* algorithm
• ARA* - Anytime A* with Provable Bounds on Sub-Optimality

留言與分享

A*搜索算法简介及保姆级代码解读

• 1. A*算法简单介绍
• • 1.1 A*算法理论基础
  • • 1.1.1 节点计算
    • 1.1.2 由计算得出的小结论
  • 1.2 算法逻辑结构
• 2. 代码解析
• • 2.1 引入地图
  • 2.2 预处理
  • 2.3 定义父节点`parent`
  • 2.4 主循环
  • • 2.4.1
    • 2.4.2
    • 2.4.3
    • 2.4.4
    • 2.4.5
  • 2.5 画图
• 3. 结果

1. A*算法简单介绍

A*算法是一种路径规划算法，和传统的Dijkstra算法有所不同，该算法有选择地进行节点搜索，因此比Dijkstra算法更快、搜索的点更少。
阅读本文，不需要掌握Dijkstra算法的知识，请放心食用。
注意：本文只介绍二维A*算法及相关示例。

PS：本文代码编写参考B站up主
小黎的Ally的路径规划与轨迹跟踪系列算法学习_第4讲_A*算法，讲解详细，本文代码部分是将其代码进行了些许改动并加以解释，在此对up主的辛苦表达感谢！！

PPS：本文所使用的把地图栅格化的函数来自于博客
Matlab中将一幅图片做成栅格地图，本文进行了些许改动，再次一并感谢博主的辛苦！

1.1 A*算法理论基础

A*算法首先将要搜索的区域划分为若干栅格（grid），并有选择地标识出障碍物（Obstacle）与空白区域。一般地，栅格划分越细密，搜索点数越多，搜索过程越慢，计算量也越大；栅格划分越稀疏，搜索点数越少，相应地搜索精确性就越低。

如上图，引入地图信息后画出栅格，该图片采用 100 × 100 100 \times 100 100×100的栅格划分，图中黑色区域为障碍物区域。图中绿点为起始点，红点为终点。

对每个节点，在计算时同时考虑两项代价指标：当前节点与起始点的距离，以及当前节点与目标点的距离：
f = g + h f = g + h f=g+h其中 f f f为总代价， g g g为当前节点距离起始点的距离， h h h为当前节点距离目标点的距离。
而在计算距离时，又可以采用两种方式：
欧氏距离：
L = ( x 1 − x 2 ) 2 + ( y 1 − y 2 ) 2 L = \sqrt{\left( x_1 - x_2 \right) ^2 + \left( y_1 - y_2 \right) ^2} L=(x1​−x2​)2+(y1​−y2​)2 ​曼哈顿距离：
L = ∣ x 1 − x 2 ∣ + ∣ y 1 − y 2 ∣ L = \lvert x_1 - x_2 \rvert + \lvert y_1 - y_2 \rvert L=∣x1​−x2​∣+∣y1​−y2​∣为了计算方便，本文计算 h h h时采用曼哈顿距离，这也是A*算法中的一贯做法。
为了方便计数，在计算每一个节点时，在栅格左上角写 f f f值，左下角写 g g g值，右下角写 h h h值。

1.1.1 节点计算

这里举一个例子。

如图所示，A点为起始点，M点为终点。对于A点来说，对其周边的8个节点进行寻找，A点本身为父节点，周边8个点为子节点。

假设：

• 格子边长为10，这样水平和竖直位移一格为10，而对角位移一格为14；
• 从父节点到子节点可以水平、竖直、对角线位移计算 g g g值，而从子节点到目标点只能使用水平、竖直位移计算曼哈顿距离 h h h值。

对于点B：从A到B只需右移一格，因此B的 g = 10 g=10 g=10；从B到M需要先右移四格，再上移三格，因此B的 h = 40 + 30 = 70 h = 40+30=70 h=40+30=70。这样B的 f = g + h = 10 + 70 = 80 f=g+h=10+70=80 f=g+h=10+70=80。将三者都记在B格中。
对于点C：从A到C只需向右上方平移一格，因此C的 g = 14 g=14 g=14；从C到M需要先右移四格，再上移两格，因此C的 h = 40 + 20 = 60 h=40+20=60 h=40+20=60，继而C的 f = g + h = 74 f=g+h=74 f=g+h=74。
同理可以计算出D的 g = 14 , h = 80 , f = 94 g=14, h=80,f=94 g=14,h=80,f=94，以及其他几个子节点的值。

将这8个子节点进行对比，可以发现，C点的 f f f值最小，因此选取C点为下一步搜寻的父节点， A C ‾ \overline{AC} AC即为路径迈出的第一步。

如图所示，现在将C作为新的父节点。

对于J点：
从C到J只需右移一格，因此J的 g = 10 + g C = 10 + 14 = 24 g=10+g_C =10+14=24 g=10+gC​=10+14=24。需要注意的是，J点的 g g g值是从C到J的 g g g加上从A到C的 g g g，即当前子节点的 g g g值是从起点到该点的所有 g g g累和。
从J到M需要先右移三格，再上移两格，因此J的 h = 30 + 20 = 50 h = 30+20=50 h=30+20=50。这样J的 f = g + h = 24 + 50 = 74 f=g+h=24+50=74 f=g+h=24+50=74。将三者都记在J格中。

对于I点：
从C到I只需向右上方平移一格，同样地，I点的 g g g值是从C到I的 g g g加上从A到C的 g g g，即当前子节点的 g g g值是从起点到该点的所有 g g g累和，因此I的 g = 14 + g C = 14 + 14 = 28 g=14+g_C =14+14=28 g=14+gC​=14+14=28。
从I到M需要先右移三格，再上移一格，因此I的 h = 30 + 10 = 40 h = 30+10=40 h=30+10=40。这样I的 f = g + h = 28 + 40 = 68 f=g+h=28+40=68 f=g+h=28+40=68。将三者都记在I格中。
同样计算出其他点的值。
可以看出，I点的 f f f值最小，因此选择I点作为路径上的下一个点，此时路径变为 A C I ‾ \overline{ACI} ACI。

如此一步步进行迭代，最后找到最优轨迹。

1.1.2 由计算得出的小结论

• 每个节点中需要存储至少3个值： g , h , f g,h,f g,h,f。
• 当子节点迭代到目标点本身时， h = 0 h=0 h=0，即当前节点到目标点的距离为0.
• 在算法实施时， g g g的计算可以取10或14，即从父节点到子节点可以水平/竖直位移，也可以沿对角线位移；而 h h h的计算只能取10的倍数，即从当前子节点到目标点只能水平/竖直位移。前者可以采用欧氏距离，后者只能采用曼哈顿距离。
• 当某个子节点被选中后，就作为下一次搜索的父节点；并且为了避免节点的重复计算和筛选，在下一次搜索时，需要将上一步选中的节点从“可搜索”列表中删除。如上图中，当C为父节点时，选中I作为下一次的父节点；那么当第二步I为父节点时，由于C已经在已选轨迹上了，所以I的子节点实际上并不包括C，只计算7个子节点即可。
• 鉴于上一点，可以预想，算法的具体实施过程中，需要有两个数组保存“待计算子节点”和“已被选中的节点”。当“待计算子节点”中某个点符合 f f f最小时，就将其加入“已选中的节点”，并在“待计算子节点”中删除该点，防止后续重复计算。
• 当算法结束时，保存在“已选中的节点”中的所有点，按照顺序即为找出的路径。
• 算法结束时，最后一个点的 f f f值，即为从起点到终点所用的距离。
• 算法的停止条件：1)当“待计算子节点”中没有点时，即已经没有点可供寻找了；2)当当前父节点恰为目标节点时，即 h = 0 h=0 h=0。

1.2 算法逻辑结构

A*算法的逻辑结构如下：
1）初始化。导入地图信息，设置障碍物区域，设置起点start、终点target、栅格数量 m × n m \times n m×n等。
2）数据预处理。建立“待计算子节点”的数组openlist，“已选中的节点”的数组closelist，保存路径的数组closelist_path。除此之外，还需建立一个当前子节点集合children，用来保存当前父节点周围8个子节点的坐标，以及父节点本身parent；还有保存代价值 g , h , f g, h,f g,h,f的数组openlist_cost和closelist_cost。
3）对子节点们children中的每个节点child：
若该子节点不在“待计算子节点”节点openlist中，则追加进去；
若在，则计算出该child的 g g g值，该 g g g值是从起点到父节点parent的距离加上父节点到该子节点的距离。若该 g g g值小于之前openlist_cost中的 g g g最小值，那么就将openlist_cost中的最小 g g g值更新；
4）由于该代价最小点已经加入了轨迹，因此将该点加入clost_list和closelist_path，并从openlist中剔除；
5）更新openlist中的最小代价值，并以其为父节点开始新一轮搜索。

2. 代码解析

2.1 引入地图

代码块：

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

%% 画地图

% 栅格地图的行数、列数定义
m = 150;
n = 150;
% 地图m行n列

start = [10, 20];        % 起始节点
target = [130, 80];       % 终止节点

obs = TrunToGridMap(m, n);


% 画格子
for i = 0 : 5 : m
    plot([0, n], [i, i], 'k', 'handlevisibility', 'off');
    hold on;
end

for j = 0 : 5 : n
    plot([j, j], [0, m], 'k', 'handlevisibility', 'off');
end

axis equal;
xlim([0, n]);
ylim([0, m]);

% 绘制障碍物、起止点颜色块
scatter(start(1), start(2), 700, 'pg', 'filled');
scatter(target(1), target(2), 700, 'pr', 'filled');

for i = 1 : size(obs, 1) - 1
    temp = obs(i, :);
    fill([temp(1)-1, temp(1), temp(1), temp(1)-1],...
        [temp(2)-1, temp(2)-1, temp(2), temp(2)], 'k', 'handlevisibility', 'off');
end

temp = obs(size(obs, 1), :);
fill([temp(1)-1, temp(1), temp(1), temp(1)-1],...
    [temp(2)-1, temp(2)-1, temp(2), temp(2)], 'k');

绘制一个 150 × 150 150 \times 150 150×150的地图，起点设置为 ( 10 , 20 ) (10,20) (10,20)，终点设置为 ( 130 , 80 ) (130,80) (130,80)。
障碍物的坐标通过函数TrunToGridMap(m, 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

function obs = TrunToGridMap(a, b)
    I=imread('此处放地图图片的文件名');   %读入图片
    I = rgb2gray(I);     %将图片转为灰度图
    I = imrotate(I, -90);

    l=1;    %网格边长
    B = imresize(I,[a/l b/l]);
    J=floor(B/255); 

    axes('GridLineStyle', '-');
    axis equal;

    hold on
    grid on
    axis([0,a,0,b]);
    set(gca,'xtick',0:10:a,'ytick',0:10:b);

    obs = [];

    %障碍物填充为黑色
    for i=1:a/l-1
        for j=1:b/l-1
            if(J(i,j)==0)
                obs(end+1, :) = [i, j];
            end
        end
    end
end

该函数读取一个图片文件，将其转化为灰度图像，并将灰度图中黑色色块所在的坐标返回为障碍物坐标。

2.2 预处理

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

%% 预处理

% 初始化closelist
closelist = start;
closelist_path = {start, start};      % 路径，从自身到自身
closelist_cost = 0;
children = child_nodes_cal(start, m, n, obs, closelist);

% 初始化openlist
openlist = children;


for i = 1 : size(openlist, 1)   % i为第i个节点
    openlist_path{i, 1} = openlist(i, :);   % openlist_path的第i行第1列为第i个节点child
    openlist_path{i, 2} = [start; openlist(i, :)]; % openlist_path的第i行第2列为一个列向量，分别是起始节点和当前child
end

for i = 1 : size(openlist, 1)
    g = norm(start - openlist(i, 1:2));
    h = abs(target(1) - openlist(i, 1)) + abs(target(2) - openlist(i, 2));
    f = g + h;
    openlist_cost(i, :) = [g, h, f];
end

这一部分主要做了以下几步：

• 把起点start设为轨迹的第一个点：closelist = start；
• 路径初始化，即从起点到起点：closelist_path = {start, start}；
• 代价首先置为0：closelist_cost = 0；
• 利用child_nodes_cal函数计算当前父节点（即起点）周围的子节点们；
• 这些子节点们children即为待计算的子节点，也就是openlist；
• 随后进入一个循环，对每一个子节点i，openlist_path中第i行第1个元素储存第i个节点child，第2个元素为一个列向量，分别是第i个child的起点和它本身；
• 第二个循环则是计算代价的循环，对每一个子节点i，计算 g g g（这里使用范数norm）， h h h（曼哈顿距离）和 f f f；之后，在代价数组openlist_cost中储存这三个代价值，因此openlist_cost第i行的元素为一个行向量；

其中用到的子节点计算函数如下：

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89

function child_nodes = child_nodes_cal(parent_node, m, n, obs, closelist)
    child_nodes = [];
    field = [1, 1;
        n, 1;
        n, m;
        1, m];
    
    % 左上子节点
    child_node = [parent_node(1) - 1, parent_node(2) + 1];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        % [in, on] = inpolygon, 返回 in，以指明 xq 和 yq 所指定的查询点是在 xv 和 
        % yv 定义的多边形区% 域的边缘内部还是在边缘上,in为内部，on为边缘上
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    
    % 上子节点
    child_node = [parent_node(1), parent_node(2) + 1];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    
    % 右上子节点
    child_node = [parent_node(1) + 1, parent_node(2) + 1];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    % 左子节点
    child_node = [parent_node(1) - 1, parent_node(2)];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    %右子节点
    child_node = [parent_node(1) + 1, parent_node(2)];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    
    %左下子节点
    child_node = [parent_node(1) - 1, parent_node(2) - 1];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    
    % 下子节点
    child_node = [parent_node(1), parent_node(2) - 1];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    
    % 右下子节点
    child_node = [parent_node(1) + 1, parent_node(2) - 1];
    if inpolygon(child_node(1), child_node(2), field(:, 1), field(:, 2))
        if ~ismember(child_node, obs, 'rows')
            child_nodes = [child_nodes; child_node];
        end
    end
    
    
    %% 排除已经存在于closelist的节点
    delete_idx = [];
    for i = 1 : size(child_nodes, 1)
        if ismember(child_nodes(i, :), closelist, 'rows')
            delete_idx(end+1, :) = i;
        end
    end
    
    child_nodes(delete_idx, :) = [];
end

2.3 定义父节点parent

1
2
3
4

%% 定义父节点
% 从openlist开始搜索移动代价最小的节点
[~, min_idx] = min(openlist_cost(:, 3));    % 看f值最小，min_idx为f最小的那一行
parent = openlist(min_idx, :);% 以min_idx该行的子节点child_node为新的父节点

这一步搜索openlist_cost代价数组中 f f f最小的元素，记下其索引min_idx；
随后该索引对应的openlist中的元素即为“ f f f值最小的待计算子节点”，作为下一步计算的父节点parent。

2.4 主循环

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

%% 进入循环
flag = 1;
while flag
    % 找出父节点的忽略closelist的子节点
    children = child_nodes_cal(parent ,m, n, obs, closelist);
    
    % 判断这些子节点是否在openlist中，若在，则更新；没在，则追加到openlist中
    for i = 1 : size(children, 1)
        child = children(i, :);
        [in_flag, openlist_idx] = ismember(child, openlist, 'rows');
        % in_flag表示该child_node是否在openlist中
        % openlist_idx表示该child_node在openlist中的行数
        
        g = openlist_cost(min_idx, 1) + norm(parent - child);
        % 原来的g加上新的g
        h = abs(child(1) - target(1)) + abs(child(2) - target(2));
        f = g + h;
        
        if in_flag  % 若在，则比较更新g, f
            if g < openlist_cost(openlist_idx, 1)   
                % openlist_cost(openlist_idx,1)指的是openlist_cost中idx这一行（即child_node所在的一行）的第一个坐标，即g
                openlist_cost(openlist_idx, 1) = g;
                openlist_cost(openlist_idx, 3) = f;
                openlist_path{openlist_idx, 2} = [openlist_path{min_idx, 2}; child];
                % openlist_path的第i行第2列为一个列向量，分别是起始节点和当前child，此处定义新的起始节点为openlist_path(min_idx,2)， 而openlist_path(min_idx,2)指第min_idx行所对应的 child在openlist_path中对应的路径，相当于把新的child附加到了路径上，延长了路径
            end
        else
            openlist(end+1, :) = child;
            openlist_cost(end+1, :) = [g, h, f];
            openlist_path{end+1, 1} = child;
            openlist_path{end, 2} = [openlist_path{min_idx, 2}; child];
        end
    end
    
    
    % 从openlist移除代价最小的节点到closelist
    closelist(end+1, :) = openlist(min_idx, :);
    closelist_cost(end+1, :) = openlist_cost(min_idx, 3);
    closelist_path(end+1, :) = openlist_path(min_idx, :);
    
    % 同样地，openlist中少了该节点
    openlist(min_idx, :) = [];
    openlist_cost(min_idx, :) = [];
    openlist_path(min_idx, :) = [];
    
    
    % 重新搜索：从openlist搜索移动代价最小的节点，作为新的父节点
    [~, min_idx] = min(openlist_cost(:, 3));
    parent = openlist(min_idx, :);
    
    % 判断是否搜索到终点
    if parent == target
        closelist(end+1, :) = openlist(min_idx, :);
        closelist_cost(end+1, 1) = openlist_cost(min_idx, 1);
        closelist_path(end+1, :) = openlist_path(min_idx, :);
        flag = 0;
    end

end

2.4.1

对这一部分循环来说，首先利用child_nodes_cal函数得出当前父节点parent周围的子节点们children；
对每一个子节点child，判断其是否在“待计算子节点”openlist列表中——in_flag=1表示在列表中，同时openlist_idx为该子节点在openlist中的索引号。
判断完成后，首先计算该子节点的 g , h , f g,h,f g,h,f值，注意： g g g值为父节点parent的 g g g加上该child子节点的 g g g值。

2.4.2

计算 g , h , f g,h,f g,h,f之后，再来看子节点在openlist中的情况。由于循环中查看了parent周围所有子节点的情况，所以一定会存在某个子节点的 g g g比其他子节点的 g g g都小。如果找到了这样的子节点，那么就更新该子节点在openlist_cost中对应位置的 g , h , f g,h,f g,h,f值。

同时，还要将该子节点加入到路径openlist_path中，即这一句代码

1

openlist_path{openlist_idx, 2} = [openlist_path{min_idx, 2}; child];

这行代码的含义是：
openlist_idx表示当前子节点child在openlist中的索引，min_idx表示之前所有子节点中代价最小的子节点的索引。之前在预处理一节中提到，openlist_path的第i行第2列为一个列向量，分别是父节点和当前child，相当于第2列储存了路径。

而min_idx表示“代价最小的节点”，也就是父节点parent。因此这里openlist_path{min_idx, 2}表示父节点的路径。而[openlist_path{min_idx, 2}; child]则是把新的child附加到这一个列向量上，相当于把新的child附加到路径尾端，把向量长度延长了一个child，延长了路径。

这样，[openlist_path{min_idx, 2}; child]构成了一个“父节点路径-当前child”的新路径。

把这个新路径赋值给openlist_idx索引所对应的openlist_path上，即为该openlist_idx索引对应的子节点的路径。

2.4.3

另一方面，如果该child不在openlist中，那么就把该child添加到“待计算子节点”列表中，其 g , h , f g,h,f g,h,f值添加到openlist_cost代价列表中，其路径[openlist_path{min_idx, 2}; child]添加到路径openlist_path中。

2.4.4

由于移动代价最小的节点已经是路径上的一点了，所以为了避免重复计算，应当把她从“待计算子节点”列表中删除，加入“已计算节点”中，即

1
2
3
4
5
6
7

closelist(end+1, :) = openlist(min_idx, :);
closelist_cost(end+1, :) = openlist_cost(min_idx, 3);
closelist_path(end+1, :) = openlist_path(min_idx, :);
% openlist中少了该节点
openlist(min_idx, :) = [];
openlist_cost(min_idx, :) = [];
openlist_path(min_idx, :) = [];

仔细观察不难发现，min_idx对应的正是这一步的parent节点，直到这里，我们才把它加入到“已计算节点”列表中，在之前它一直呆在openlist中。

2.4.5

父节点加入到了“已计算节点”中了，那么下一步就没有父节点了，所以需要找出新的父节点：

1
2

 [~, min_idx] = min(openlist_cost(:, 3));
parent = openlist(min_idx, :);

同时还需要判断一下，是否已经进行到了程序结束：

1
2
3
4
5
6
7

% 判断是否搜索到终点
if parent == target
    closelist(end+1, :) = openlist(min_idx, :);
    closelist_cost(end+1, 1) = openlist_cost(min_idx, 1);
    closelist_path(end+1, :) = openlist_path(min_idx, :);
    flag = 0;
end

需要注意，即使当当前父节点已经是目标点了，也不要忘了把这个父节点加入到“已计算节点”中，相当于把目标点添加入路径中，形成路径上最后一个点。

2.5 画图

这一步就是将结果绘制出来了：

1
2
3
4
5
6
7
8
9

path_opt = closelist_path{end, 2};
% closelist_path中第二列存放路径，所以path_opt存放的是路径的(x,y)值
path_opt(:, 1) = path_opt(:, 1) - 0.5;
path_opt(:, 2) = path_opt(:, 2) - 0.5;
plot(path_opt(:, 1), path_opt(:, 2), 'm', 'linewidth', 1.5);

title(['Total length of path: ' num2str(closelist_cost(end, 1))]);
legend('Start node', 'Target node', 'Obstacle', 'Path', 'location', 'northwest');
set(gca, 'fontsize', 35, 'fontname', 'times new roman');

3. 结果

这里导入稻妻地图作为路径规划的示例，该地图多为岛屿，路径复杂，障碍物分散且数量多，十分适合作为示例使用。

将其灰度化得到

把地图划分为 m × n = 150 × 150 m \times n = 150 \times 150 m×n=150×150的栅格，并设置起点为 ( 10 , 20 ) (10,20) (10,20)，终点为 ( 130 , 80 ) (130,80) (130,80)：

图中绿色为起点，红色为终点。

路径规划结果：

可以看出，A* 算法可以找到期望路径。

文件所有代码在笔者的资源二维A*算法路径规划matlab代码中可以下载得到。

留言與分享

不是有个叫迪杰斯特拉(Dijkstra)的算法吗？名字难记不说，加权图的最短路径问题在实际开发中也很少需要自己实现，所以总是很快就忘了。。。

迪杰斯特拉算法是什么

迪杰斯特拉算法是一种用于求解图的最短路径的算法。很多人可能只是听过这个名字。虽然它的原理简单易懂，但初次接触可能有些难以理解，不过只要逐步分析就会发现它其实是一个相当直观的算法。
对于没有权重的迷宫搜索等问题，可以用广度优先搜索(BFS)解决，但如果每条边都有权重，就需要计算所有可能的路径。假设每个顶点只经过一次，且有E条边，那么时间复杂度会是O(E!)，计算量会爆炸式增长。
这样计算起来就很不现实。

而迪杰斯特拉算法正是高效解决这类问题的算法。

需要注意的是，各边的成本必须是非负值（大于等于0）。如果包含负数，则需要使用贝尔曼-福特(Bellman-Ford)等算法。

步骤

迪杰斯特拉算法的步骤非常简单：

1. 将起点的最短距离设为0
2. 从未访问的点中选择已知最短距离且距离最小的顶点移动
3. 设置该顶点连接的其他顶点的最短距离。如果可以更新更短的距离，则更新。
4. 重复以上步骤，直到所有顶点的最短距离都确定

光这么说可能有点抽象，下面我们通过一个具体例子来说明。

示例

虽然方法很原始，但这是最能表达清楚的方式…
我们来看下面这张图的最短路径：

绿色表示最短路径已确定并已访问的顶点，红色是起点顶点。每个顶点上的数字表示当前的最短路径：

可以看到，算法从起点开始，每次移动到当前最短路径的顶点，并计算相邻顶点的最短路径。

实现

接下来我们来实现这个算法。
按照上述步骤直接实现：

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

function main(nodes) {
    const start = nodes[0]
    // 记录已访问的顶点
    const visited = new Set()
    const routesFromStart = new Map()
     // 记录从起点出发的距离

    routesFromStart.set(start, {distance: 0})
    for(const n of nodes) {
        if(n != start) {
            // 除起点外所有顶点初始化为无穷大
            routesFromStart.set(n, {distance: Number.MAX_VALUE})
        }
    }
    let current = start
    let routes = new Map()
    while(current != null) {
        visited.add(current)
        for(const edge of current.edges) {
             // 计算相邻顶点的最短距离，如果更小则更新
            if(edge.cost + routesFromStart.get(current).distance < routesFromStart.get(edge.to).distance) {
                routesFromStart.set(edge.to, {distance: edge.cost + routesFromStart.get(current).distance})
                routes.set(current, edge.to)
            }
        }
        let cheapestNodeDistance = Number.MAX_VALUE
        current = null
        // 从已计算最短距离的未访问顶点中选择最小的顶点

        for(const city of routesFromStart.keys()) {
            if(!visited.has(city) && cheapestNodeDistance > routesFromStart.get(city).distance){
                cheapestNodeDistance = routesFromStart.get(city).distance
                current = city
            }
        }
    }
    return routesFromStart.get(nodes[nodes.length - 1]).distance
}

假设顶点数为V，这段代码在最坏情况下需要遍历所有边，并在内部循环中选择最小顶点，因此时间复杂度是O(V² + E)。空间复杂度需要记录每个顶点，所以是O(V)。

关于使用优先队列的实现

细心的读者可能已经发现，这段代码中选取最小顶点的逻辑可以优化，这就是优先队列(Priority Queue)。
优先队列的插入和取出操作需要O(logN)的时间复杂度，但比线性搜索最小顶点更快。

JavaScript没有内置优先队列，以下是Python实现：

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

def dijkstra(nodes):
    start_node = nodes[0]
    routes_from_start = {n: math.inf for n in nodes}

    # 起点距离设为0
    routes_from_start[start_node] = 0
    
    minHeap = []

    # 添加起点
    heappush(minHeap, (0, start_node))

    # 直到堆为空
    while minHeap:
        (cost, current_node) = heappop(minHeap)

        # 检查priority key是否重复
        if cost > routes_from_start[current_node]:
            continue

        for node in current_node.routes:
            price_info = current_node.routes[node]
            if routes_from_start[node] > price_info + routes_from_start[current_node]:
                routes_from_start[node] = price_info + routes_from_start[current_node]
                # 记录更新最短距离的节点
                heappush(minHeap, (price_info + routes_from_start[current_node], node))

    return routes_from_start[nodes[-1]]

优先队列的说明我们以后再讨论。
优化后的算法时间复杂度为O(V + ElogE)，比最初的实现更高效。

记录路径

现在我们已经能求出最短成本，但这是"最短路径"问题，我们自然还想知道具体路径。
改进上面的代码：

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

def dijkstra(nodes):
    start_node = nodes[0]
    routes_from_start = {n: math.inf for n in nodes}

    # 起点距离设为0
    routes_from_start[start_node] = 0
    
    minHeap = []

    # 添加起点
    heappush(minHeap, (0, start_node))
    path = collections.defaultdict(Node)

    # 直到堆为空
    while minHeap:
        (cost, current_node) = heappop(minHeap)

        # 检查priority key是否重复
        if cost > routes_from_start[current_node]:
            continue

        for node in current_node.routes:
            price_info = current_node.routes[node]
            if routes_from_start[node] > price_info + routes_from_start[current_node]:
                routes_from_start[node] = price_info + routes_from_start[current_node]
                # 记录更新最短距离的节点
                path[node.id] = current_node.id
                heappush(minHeap, (price_info + routes_from_start[current_node], node))

    current_node = nodes[-1].id
    path_array = []

    # 从终点回溯记录的最短路径节点
    while current_node:
        path_array.append(current_node)
        if current_node not in path:
            break
        current_node = path[current_node]

    return routes_from_start[nodes[-1]], path_array[::-1]

迪杰斯特拉算法会记录更新最短距离的节点，最后回溯即可。时间复杂度会增加最短路径节点数的计算量。

为什么这样能求出最短路径

看到这里，很多人可能会想：算法本身很简单，实现也不难，但为什么能保证求出最短距离呢？我们简单验证一下：

假设集合L中的顶点到起点S的最短距离已确定，那么从L连接到的最小顶点也应该是S的最短距离。

在集合T中选取最小顶点i，则d[i] = min(T)。对于任意顶点k，最短距离d[k] ≥ d[i]是确定的，因为d[i]是最小值且各边权重非负。
通过归纳法可以证明这一点。

其实这就是一个递推公式：

d[i] = min(k ⊂ T) + i到L中相邻顶点的最短距离

说到递推公式就想到动态规划(DP)。关于DP可以参考这篇文章：
https://qiita.com/drken/items/a5e6fe22863b7992efdb

用DP来看值的更新过程：

纵轴表示迭代次数，横轴表示顶点。原来迪杰斯特拉算法是DP的一种啊。

总结

通过以上分析，迪杰斯特拉算法一旦理解后其实相当简单。以后在算法问题中遇到类似问题时，希望能快速联想到这个算法。
*顺便说一句，我也想写一篇关于DP的文章

讲解视频可以在这里观看：
https://youtu.be/jz8b0q5R1Ss

参考

http://www.lab2.kuis.kyoto-u.ac.jp/~shuichi/algintro/alg-6s.pdf
https://www.youtube.com/watch?v=X1AsMlJdiok

留言與分享

23. 最小生成树

引：在设计电子电路时，我们常常需要将多个组件的针脚连起。若要连接n个针脚，可以用n-1根连线，我们希望使用的最短的连线来连接针脚。

最小生成树问题：
一个连通无向图$G = (V, E)$中，对于每条边$(u, v) \in E$，为其赋予权重$w(u, v)$。找到一个无环子集$T \subseteq E$，使得所有的结点可以互相到达，并且具有最小的权重，即：$w(T) = \sum_{(u, v) \in T}w(u, v)$的值最小。$T$就是图$G$的最小生成树。

下图中阴影部分组成的就是最小生成树。

生成最小生成树

本章将介绍两种生成最小生成树的算法，他们都用了贪心策略。

思想：创建一个边的集合A，首先讲集合A设置为空集，每次循环加入一条安全的边到集合A中，是的A是最小生成树T的子集，直到集合A覆盖了所有结点。

下面是这个贪心策略的伪代码：

如何找出一个安全的边

这个贪心策略的重点就是如何找出安全的边。

设边集合A覆盖的结点集合为S，没覆盖到的结点集合为V-S。想要生成最小生成树，则必须有且仅有一条边连接S和V-S（连接两个集合的意思为边的一头在S中，另一头在V-S中），所有连接S和V-S的边中权重最小的边就是安全的边。

可以用反证法证明这种选择安全的边的方法是正确的。

Kruskal算法

思想：初始时，将$|V|$个结点视为$|V|$棵树并组成森林。所有边的权重按照从小到大排序，从最小的边开始编译，如果这个边能连接森林中的两棵树，则将这条边加入集合A，并且将这两棵树合并为一棵树。直到森林中只剩一棵树。

伪代码和事例如下：

Prim算法

思想：初始时，给所有的结点赋予一个无穷大的key值，并加入优先队列。每次循环时，从优先队列中取出一个key值最小的结点u加入集合A中，并且更新优先队列中所有与结点u相连的结点v的key值（如果$w(u,v) < v.key$，则$v.key = w(u,v)$）。直到优先队列为空。

伪代码和事例如下：

TODO:时间复杂度分析

留言與分享

22. 基本的图算法

本章介绍图的表示和搜索。
许多的图算法在一开始都会先通过搜索来获得图的结构，其他一些图算法则是对基本的搜索加以优化。
可以说，图的搜索技巧是整个图算法领域的核心。

图的表示

$G = (V, E)$

图G由结点的集合V 和 边的集合E 组成

可以用两种方法表示，一种是邻接链表（下图b），一种是邻接矩阵（下图c）。
邻接链表一般用于表示稀疏图（边数$|E|$远小于$|V|^2$），默认都使用邻接链表表示。在稠密图（$|E|$接近$|V|^2$）的情况下，倾向使用邻接矩阵。

无向图：

有向图：

邻接链表表示

由一个链表数组$Adj$组成，数组大小为$|V|$，链表$Adj[u]$存储了结点$u$出发的边的终点。
比如上面无向图中，从结点1出发有两条边，边的终点分别是2和4，反应到邻接链表中就是$Adj[1]$存储了两个节点2和4。

邻接链表需要的存储空间为$\Theta(|V|+|E|)$

邻接矩阵表示

由一个$|V| \times |V|$ 的矩阵$A = (a_{ij})$表示，并满足以下条件：

$a_{ij} = \left \{ \begin{aligned} & 1 & if (i, j) \in E, &\\ & 0 & other, & \end{aligned} \right.$

当边有权重时，$a_{ij}$可以存储权重值。

邻接矩阵需要的存储空间为$\Theta(|V|^2)$

一些术语

一个结点的出度为，从这个结点触发的边的个数。
一个结点的入度为，到达这个结点的边的个数。

广度优先搜索

思路：利用一个队列，首先将头结点插入队列，循环的取出队列中的一个结点，并将于该结点相连的结点加入队列，直到队列变空，搜索结束。

辅助的给每个结点加入三个属性：

• color：白色表示还未被搜索，灰色表示加入队列，黑色表示从队列里出来
• d：该节点在第几轮被发现
• $\boldsymbol{\pi}$：该结点的前驱（结点第一次被发现时，正在查询的结点）

下面是广度优先搜索的伪代码和事例：

广度优先搜索的正确性

TODO

广度优先搜索的性质

最短路径

从某个结点u做广度优先搜索，可以找出这个结点到其他结点的最短路径。

结点v和结点u的最短路径距离为v.d

从结点v开始 循环列举前驱结点，就是结点v和结点u的最短路径中经过的结点。

广度优先树

广度优先搜算中，我们在$\pi$属性中存储了每个结点的前驱结点，另前驱结点作为该结点的父节点，我们可以得到一个广度优先树。发现新结点的边为树边。

深度优先搜索

思路：尽可能的深入。总是对最近发现的结点v出发的边进行探索，直到结点v出发的边全部被发现时，回溯到v的前驱结点继续探索，直到从源节点出发所有可以到达的结点都被发现。这时如果还有未被搜索的结点，则再随机找一个未被搜索的结点作为源节点进行搜索。

因为可能从一个源节点出发不能搜索到所有结点，所以深度优先搜索会产生多个深度优先树组成深度优先深林。

辅助的给每个结点加入四个属性：

• color：白色表示还未被搜索，灰色表示已被发现但出发的边还没搜索完，黑色表示已被发现并且出发的边已被搜索完
• d：该节点的发现时间
• f：该结点出发搜索结束的时间
• $\boldsymbol{\pi}$：该结点的前驱（结点第一次被发现时，正在查询的结点）

下面是深度优先搜索的伪代码和事例：

结点中的x/y表示该结点的发现时间为x，完成时间为y
图中树边被标记为灰色，向前，横向，向前的边分别被标记为B, C, F

深度优先搜索的性质

括号花定理

如上图b所示，每个结点的发现时间和完成时间组成一个括号花结构。

对于两个结点u和v，如果区间[u.d, u.f] 和 [v.d, v.f] 不相交，则u不是v的后代，v也不是u的后代；
如果[u.d, u.f] 包含在 [v.d, v.f] 内，则u是v后代，反之亦然。

白色路径定理

如果结点v是u的后代，当u.d时间时，存在一条从u到v全部由白色结点组成的路径。

拓扑排序

对于有向无环图$G = (V, E)$，其拓扑排序是G种所有结点的一种线性次序，该次序满足如下条件：如果G包含边(u, v)，则结点u在拓扑排序中处于结点v的前面。

许多实际应用都需要使用有向无环图来说明事件的优先次序。

下图是一个穿衣服的实例，图a中表示了穿戴某一件衣物的依赖关系，图b是图a拓扑排序的结果，按照图b相反的顺序穿衣服就是没问题的了。

拓扑排序的思想为：对图G做DFS，然后按照结点搜索的完成时间逆序排序。伪代码如下：

拓扑排序的正确性：在向无环图中，如果有变(u, v)，那么v肯定比u先完成搜索。

强联通分量

有向图$G = (V, E)$的强连通分量是一个最大结点集合$C \subseteq V$，对于该集合的任意一对结点可以互相到达。

下图a的灰色覆盖区域就是图的强联通分量。

算法伪代码如下：

1. 对G做DFS，并记下每个结点的完成时间u.f
2. 转置G，得到$G^T$，如上图b
3. 按照u.f大小，从大到小对$G^T$做DFS
4. 第三步得到的深度优先深林中，每棵深度优先树的结点都组成为G的一个强联通分量

思想：
  对于分量图$G^{SCC} = (V^{SCC}, E^{SCC})$，$V^{SCC} = {v_1, v_2, ... , v_k}$，图$G$中的强联通分量$C_i$集合代表分量图中的结点$v_i$：如果有$x\in C_i, y\in C_j$，且图$G$中有边$(x, y)$，则边$(v_i, v_j) \in E^{SCC}$，且图$G^{SCC}$是有向无环图。

证明：
  引理1：$C$和$C^\prime$是图$G$的两个强连通分量。当$u, v \in C; u^\prime, v^\prime \in C^\prime$， 若有$u \rightsquigarrow u^\prime$，则不可能有$v^\prime \rightsquigarrow v$。
     可用反证法证明，若有$u \rightsquigarrow u^\prime$且$v^\prime \rightsquigarrow v$，则$C$和$C^\prime$将成为同一个强连通分量。
  引理2：$C$和$C^\prime$是图$G$的两个强连通分量。当$(u, v) \in E$且$u \in C, v \in C^\prime$，则$f(C) > f(C^\prime)$
     当$d(C) < d(C^\prime)$时，若$x$为$C$中最早被发现的结点，则$C$和$C^\prime$中的其它结点都是$x$的子孙，也就是所有结点的完成时间都比$x$早，所以$x.f=f(C)>f(C^\prime)$
     当$d(C) > d(C^\prime)$时，根据引理1，没有从$C$到$C^\prime$的边，也就是$C^\prime$完成之后不会到达任意$C$中的结点，所以$f(C)>f(C^\prime)$
  引理3：$C$和$C^\prime$是图$G$的两个强连通分量。当$(u, v) \in E^T$且$u \in C, v \in C^\prime$，则$f(C) < f(C^\prime)$
     $(u, v) \in E^T \Longrightarrow (v, u) \in E$，由引理2可直接得知$f(C) < f(C^\prime)$
  证明：根据数学归纳法，假设算法第三行所生成的前$k$棵树都是强连通分量，现在需要考虑第$k+1$棵树是否为强连通分量：
     设第$k+1$棵树是$C$，$C$的根节点为$u$，$u$所在的强连通分量为$C^{\prime\prime}$。
     若$G^T$中从$C^{\prime\prime}$到其他强连通分量$C^{\prime\prime\prime}$有边，那么根据引理3可知$f(C^{\prime\prime}) < f(C^{\prime\prime\prime})$。
     而算法对$G^T$做DFS时，是按照$f$大小逆序进行，所以从$C^{\prime\prime}$出来的变都是到前$k$棵树的。
     所以$C=C^{\prime\prime}$，也就是说第$k+1$棵树为强连通分量。

留言與分享

17. 摊还分析

在摊还分析中，我们求数据结构的一个操作序列中所执行的所有操作的平均时间，来评价操作的代价。这样，我们就可以说明一个操作的平均代价是很低的，即使序列中某个单一操作的代价很高。摊还分析不同于平均情况分析，它并不涉及概率，它可以保证最坏情况下每个操作的平均性能。

本章介绍三种摊还分析方法：

• 聚合分析
• 核算法
• 势能法

首先我们先看两个小问题

栈操作

对传统的栈操作MULTIPOP(S, k)扩展，伪代码如下，相当于做k次POP(S)

1
2
3
4

MULTIPOP(S, k)
1   while not STACK-EMPTY(S) and k > 0
2       POP(S)
3       k = k - 1

如果在一个包含s个对象的栈上做MULTIPOP(S, k)操作，POP的次数为min(s, k)。

我们来分析一个由n个PUSH/POP/MULTIPOP组成的操作序列在一个空栈上的执行情况。
MULTIPOP的最坏代价为$O(n)$，那么n个操作的最坏代价是$O(n^2)$吗？

二进制计数器递增

k位二进制计数器递增，初始值为0。用一个位数组A[0…k-1]作为计数器。当计数器保存二级制数x时，x的最低位保存在A[0]中，最高位保存在A[k-1]中，因此$x = \sum\limits_{i=0}^{k-1}A[i]\cdot 2^i$。INCREMENT的伪代码如下：

1
2
3
4
5
6
7

INCREMENT(A)
1   i = 0
2   while i < A.length and A[i] == 1
3       A[i] = 0
4       i = i + 1
5   if i < A.length
6       A[i] = 1

可以想象，最坏情况当k位全部是1时，INCREMENT操作要做k次翻转操作，也就是要花费$\Theta(k)$时间，所以对初值为0的计数器，做n次操作的最坏情况花费$O(nk)$次操作？

从上两个问题，我们都根据一个操作最坏情况的时间消耗推测，整个操作序列的时间消耗，显然不是一个确界。下面我们通过三种摊还分析方法来分析这两个问题中每个操作的平均代价（摊还代价）。

聚合分析

聚合分析——确定n个操作最坏情况的总代价$T(n)$，所以每个操作的平均代价为$T(n)/n$。

栈操作：
虽然单独的MULTIPOP操作可能代价很高，但在一个空栈上实行n个操作，POP和MULTIPOP操作的总消耗最多为O(n)，因为只有push进去数据，POP和MULTIPOP才会有消耗，而且PUSH操作每次只能push一个数据进来。所以n个操作最坏情况的总代价为O(n)，除以n得到每个操作的平均代价为O(1)。

二进制计数器递增：
下图展示了二进制计数器从0递增到16的过程，阴影部分表示递增后需要翻转的位。我们可以发现A[0]每次递增时都发生翻转，A[1]每两次递增发生一次翻转，A[2]每四次递增发生一次翻转，以此类推第i位每$2^i$次递增发生一次翻转。也就是说n次INCREMENT操作后，第i位会翻转$\lfloor n/2^i\rfloor$ 次。
执行n次INCREMENT操作产生的总翻转次数为：

$\sum_{i=0}^{k-1}\lfloor n/2^i\rfloor < n\sum_{i=0}^\infty \frac{1}{2^i} = 2n$

所以，n个INCREMENT操作的最坏情况用时为O(n)，除以n得到每个操作的平均代价为O(1)。

核算法

核算法——对不同的操作赋予不同的费用（可以大于或者小于实际的代价，但不能为负值），称其为摊还代价。当一个操作摊还代价大于实际代价时，我们将差额存起来，称其为信用/存款。后续操作如果有摊还代价小于实际代价时，我们将之前存的信用拿出一部分来抵消这里的差值。如果用$c_i$表示第i个操作的实际代价，用$\hat{c}_i$表示第i个操作的摊还代价，则对任意n个操作的序列，要求：

$\sum_{i=1}^n \hat{c}_i \ge \sum_{i=1}^n c_i$

信用总能刚好抵债或者有剩余，从而使得总的摊还代价为总实际代价的上界。

栈操作：
栈操作的实际代价如下：

1
2
3

PUSH       1
POP        1
MULTIPOP   min(k, s)

我们将栈操作的摊还代价设置为：

1
2
3

PUSH       2
POP        0
MULTIPOP   0

每次PUSH操作支付1的实际代价，剩下的1存起来。存起来的钱和堆栈里的数据量相同，而POP或MULTIPOP操作的实际代价为要pop的数据量，所以存款总能保证大于等于零。
每个操作的摊还代价都是常数，所以总的摊还代价为O(n)。

二进制计数器递增：
将INCREMENT中的操作分为置位操作（第6行，0->1）和复位操作（while循环中的操作，1->0）。令置位操作的摊还代价为2，复位操作的摊还代价为0。每次置位操作可以多出1存起来（存起来的钱和计数器中1的个数相同），复位操作需要预支存款，但复位操作只有当某一位是1的时候才会执行，由于计数器中1的个数永远不会为负，所以存款总能保证大于等于零。
每个操作的摊还代价都是常数，所以总的摊还代价为O(n)。

势能法

势能法将数据结构和势能关联起来——对一个初始化数据结构$D_0$执行n个操作。$c_i$为第i个操作的实际代价，$D_i$为第i次操作后的数据结构。势函数$\Phi$将每个数据结构$D_i$映射到一个实数$\Phi(D_i)$，为数据结构$D_i$的势。第i个操作的摊还代价$\hat{c}_i$用势函数$D_0$定义为：

即，每个操作的摊还代价为实际代价加上操作引起的势变化。所以总摊还代价为：

根据公式1可以得到每个步骤的摊还代价，从而算出总的摊还代价，然后根据公式2可以得出总实际代价的上界/下界（$\Phi(D_n) - \Phi(D_{0})$为正数则得到上界，$\Phi(D_n) - \Phi(D_{0})$为负数则得到下界）。

栈操作：
将势函数定义为当前栈中的数据量。
PUSH操作会压入一个数据，所以势差为1，而且实际代价为1，所以PUSH操作的摊还代价为2。
POP/MULTIPOP操作会弹出x个对象，所以势差为-x，而实际代价为x，所以POP/MULTIPOP操作的摊还代价为0。

栈每个操作的摊还代价都是O(1)，所以n个操作的总摊还代价为O(n)。由于初始的栈为空，而结束时栈的数据量大于等于0，也就是说$\Phi(D_n) - \Phi(D_{0})$为正数，所以n个操作的总摊还代价为总实际代价的上界O(n)。

二进制计数器递增：
将势函数定义为计数器中1的个数$b_i$。
假设第i个操作将$t_i$个位复位，则实际代价为$t_i + 1$，因为还有一次置位操作。
再来考虑摊还代价，如果$b_i$为0，则将k位都复位了，因此$b_{i-1} = t_i = k$。
如果$b_i>0$，则$b_i = b_{i-1} - t_i + 1$。无论那种情况都有$b_i \le b_{i-1} - t_i + 1$，所以摊还代价为：

$\hat{c}_i = c_i + \Phi(D_i) - \Phi(D_{i-1}) \le (t_i + 1) + (-t_i + 1) = 2$

所以总摊还代价的上界为O(n)，由于计数器初始状态为0，而1的个数不会为负数，所以n个操作的总摊还代价为总实际代价的上界O(n)。

势能法的还有个优势就是，我们可以简单的得出当初始状态不是0时实际代价的上界。由公式2可得：

$\sum_{i=1}^n c_i \le \sum_{i=1}^n 2 - b_n + b_0 = 2n -b_n + b_0 \le 2n + b_0 \le 2n + k$

所以只要$k=O(n)$那么总实际代价就是O(n)。

动态表

11章介绍了散列表，我们知道散列表的大小和数据量应该是线性相关的，数据太稀疏会浪费空间，数据太稠密会浪费时间。当我们不知道会有多少数据插入时，就没办法确定散列表的大小，所以要引进动态表。
表扩张——当插入数据时发现散列表的装载因子为1（表的大小和数据量相同），对表进行扩张，比如扩张为两倍大小，首先创建一个两倍大小的表，将原来散列表中的数据插入新创建的表中，然后插入新的数据，最后不能忘了把老的表析构。

表收缩——当删除数据时发现散列包的装载因子小于某个值，对表进行收缩，与表的扩张类似，比如收缩为原来表大小的一半，首先删除这个数据并创建一般大小的新表，将原来散列表中的数据插入新创建的表中，然后把老的表析构。

表收缩时判断的装载因子的选择很重要，如果表扩张选择的两倍大小，而表收缩是选择装载因子为二分之一，这时如果在表扩张后立马删除数据，又反复的插入，删除数据，这样会导致表大小的震荡。这种情况下时间机器浪费时间。
所以我们通常可以选择该装载因子为四分之一来避免这个问题。

无论是表扩张还是表收缩，我们都可以用平摊分析来证明动态表的插入和删除操作的平均代价还是O(1)。详细情况参考书中内容，这里不做赘述。

留言與分享

16. 贪心算法

求解最优化问题的算法通常需要经过一系列的步骤，在每个步骤都面临多种选择。对于许多最优化问题，使用动态规划算法来求最优解有些杀鸡用牛刀了，可以使用更简单、更高效的算法．贪心算法就是这样的算法，它在每一步都做出当时看起来最佳的选择。也就是说，它总是做出局部最优的选择，寄希望这样的选择能导致全局最优解。

活动选择问题

我们的第一个例子是一个调度竞争共享资源的多个活动的问题，目标是选出一个最大的互相兼容的活动集合。假定有一个 n 个活动的集合$S={a_1, a_2, ... , a_n}$，这些活动使用同一个资源（例如一个阶梯教室），而这个资源在某个时刻只能供一个活动使用。每个活动$a_i$都有一个开始时间$s_i$和一个结束时间$f_i$ ，其中$0\le s_i < f_i < \infty$。如果被选中，任务$a_i$发生在半开时间区间$[s_i, f_i)$期间．如果两个活动$a_i$和$a_j$满足$[s_i, f_i)$和$[s_j, f_j)$不重叠，则称它们是兼容的．也就是说，若$s_i\ge f_j$或$s_j\ge f_i$，则$a_i$和$a_j$是兼容的．在活动选择问题中，我们希望选出一个最大兼容活动集。假定活动已按结束时间的单调递增顺序排序：

$f_1 \le f_2 \le f_3 \le ... \le f_{n-1} \le f_n$

例如下面的活动集合S：

对于这个例子，子集$S={a_3, a_9,a_11}$由相互兼容的活动组成。但它不是一个最大集，因为子集$S={a_1, a_4, a_8,a_11}$更大。实际上$S={a_1, a_4, a_8,a_11}$是一个最大兼容活动子集，另一个最大子集是
$S={a_2, a_4, a_9,a_11}$。

最优子结构

假定集合$S_{ij}$（$a_i$结束之后，$a_j$开始之前的活动）的一个最大相互兼容的活动子集为$A_{ij}$，而且$A_{ij}$包含了活动$a_k$。 由于最优解包含了活动$a_k$，我们得到两个子问题：$S_{ik}$（$a_k$开始之前的活动）的兼容活动 和 $S_{kj}$（$a_k$结束之后的活动）的兼容活动。由剪切-粘贴法可以证明$S_{ij}$的最优解为$S_{ik}$的最优解、$a_k$、 $S_{kj}$的最优解组成。即$A_{ij} = A_{ik} \cup a_k \cup A_{kj}$。
所以这时我们可以用动态规划方法求出所有k的取值情况对应的解，然后找出一个最大的作为最优解。
如果用$c[i, j]$表示$S_{ij}$的最优解的大小，递归式为：

$c[i, j] = \left \{ \begin{aligned} & 0 & if \quad S_{ij} = \varnothing, \\ & \max_{a_k \in S_{ij}}\{c[i, k] + c[k, j] + 1\} & if \quad S_{ij} \ne \varnothing, \end{aligned} \right.$

贪心选择

动态规划需要考虑每种子问题，比较之后得出哪种选择是最优解。假如我们无需求解每种子问题，直接就可以找出一个最优选择。那么将大大减少计算过程。对于活动选择问题，我们可以只考虑一个选择：贪心选择。

直观上，我们应该选择一个活动，使得选择它后剩下的资源可以尽可能多的被其他活动利用。直觉告诉我们选择S中最早结束的活动，因为他可以剩下尽可能多的资源。由于我们活动时按照结束时间排序的，也就是说贪心选择就是$a_1$。

当做出贪心选择后，只剩下一个子问题需要求解，因为不会有在$a_1$开始前结束的活动。加上前面的最优子结构性质，所以S的最优解就是$a_1 \bigcup A_{1n}$。

证明贪心选择——最早结束的活动总是最优解的一部分：

定理： 考虑任意非空子问题$S_k$，另$a_m$是$S_k$种最早结束的活动，则$a_m$在$S_k$的某个最大兼容活动子集中。

证明： $A_k$是$S_k$的某个最大兼容活动子集，$a_j$是$A_k$种最早结束的活动。如果$a_j = a_m$，则已经证明$a_m$在$S_k$的某个最大兼容活动子集中；如果$a_j \ne a_m$，那么用$a_m$替换$A_k$中的$a_j$产生$A_k^\prime$，由于$f_m \le f_j$而且$A_k$中的活动都是不相交的，所以$A_k^\prime$也是$S_k$的最大兼容活动子集，得证。

递归贪心算法

为了方便算法初始化，添加一个虚拟活动$a_0$，其结束时间$f_0 = 0$，这样子问题$S_0$就是完整的活动集S。求解原问题即可调用 RECURSIVE-ACTIVITY-SELECTOR(s, f, 0, n)。

2,3行是为了找到第一个，与$a_k$不相交的最先结束的活动$a_m$。

迭代贪心算法

贪心算法原理

前面我们讨论活动选择问题的贪心算法的过程比通常的过程繁琐一些，我们当时经过了如下几个步骤：

1. 确定问题的最优子结构．
2. 设计一个递归算法（动态规划）。
3. 证明如果我们做出一个贪心选择，则只剩下一个子问题。
4. 证明贪心选择总是安全的（步骤 3 、 4 的顺序可以调换）。
5. 设计一个递归算法实现贪心策略。
6. 将递归算法转换为迭代算法。

这个过程详细的看到贪心算法是如何以动态规划为基础的。
更一般地，我们可以简单按如下步骤设计贪心算法：

1. 将最优化问题转化为这样的形式：对其做出一次选择后，只剩下一个子问题需要求解。
2. 证明做出贪心选择后，原问题总是存在最优解，即贪心选择总是安全的。
3. 证明做出贪心选择后，剩余的子问题满足性质：其最优解与贪心选择组合即可得到原问题的最优解，这样就得到了最优子结构。

在本章剩余部分中，我们将使用这种更直接的设计方法．但我们应该知道，在每个贪心算法之下，几乎总有一个更繁琐的动态规划算法。

在证明一个问题可以用贪心算法解决过程中，证明一个问题的贪心选择性质 和 最优子结构性质 是非常重要的。

赫夫曼编码

赫夫曼编码k而已有效的压缩数据：通常可以节省20%~90%的空间，具体压缩率依赖于数据的特型。我们将待压缩数据看做字符序列。根据每个字符出现的频率，赫夫曼贪心算法构造出字符的最优二进制表示。

下面给出一个例子，待压缩文件有10万个字符数据。文件中出现了6种不同的字符，出现的次数见表格：

如果用定长编码来存储信息，那么需要用3位来表示z和6个字符。这种方法总共需要300000个二进制位来存储。
但如果用变长编码来存储信息，如表格中所示，出现频率较高字符的使用短一点的二进制编码表示。这种方法需要$(45 \times 1 + 13 \times 3 + 12 \times 3 + 16 \times 3 + 9 \times 4 + 5 \times 4) \times 1000 = 224000$位。与定长编码相比节约了25%的空间。

我们这里只考虑前缀码，即没有任何码字是其他码字的前缀。只考虑前缀码的原因是，前缀码与任何其他编码方式比都可以保证最优的压缩率，但此处不予证明。
前缀码的作用是简化解码过程，由于没有任何码字是其他码字的前缀，所以文件开头的码字是无歧义的。我们可以简单的识别出开始码字，将其转换为对应的字符，然后对文件剩余部分重复的做这种解码过程。
编码方案可以用一个二叉树表示——编码树。下图a表示了定长编码的编码树，图b表示了变长编码的编码树。叶子结点为字符，并存储字符的出现频率，字符的二进制编码是从根节点到叶子节点的路径，向左表示0，向右表示1。内部结点不存储字符，存储子节点的频率和。

给定一棵对应前缀码的树T，我们可以容易地计算出编码一个文件需要多少个二进制位。对于字母表C中的每个字符c，令属性c.freq表示c在文件中出现的频率，令$d_T(c)$表示c的叶结点在树中的深度。注意，$d_T(c)$也是字符c的码字的长度。则编码文件需要习：

$B(T) = \sum_{c \in C}c.freq\cdot d_T(c)$

个二进制位，我们将B(T)定义为T的代价。

构造赫夫曼编码

赫夫曼设计了一个贪心算法得到最优前缀码，被称为赫夫曼编码。

下面是哈夫曼算法的伪代码，输入一个字母表C，输出一个最优前缀码对应的编码树。

构造编码树的过程如下：
算法中使用一个优先队列Q，初始化时将字母表C全部插入Q。每次循环时，创造一个结点z，从优先队列取出两个最小值作为结点z的孩子，并使得z.freq为两个孩子的freq之和，然后将z结点重新插入优先队列Q中。直到优先队列中只剩下一个元素，这时这个元素为编码树的root结点，并且所有字符都已加入编码树。

下面是一个例子：

证明赫夫曼算法的正确性

贪心选择

贪心选择性质：若x,y是C中频率最低的两个字符。那么C存在一个最优前缀码，x,y的码字长度相等，且只有最后一个二进制位不同。

证明思路，假设有一个最优前缀码的编码树$T$，把x,y与编码树中深度最深的两个兄弟节点进行交换产生$T^\prime$。通过计算得出$B(T^\prime) \le B(T)$。

最优子结构

最优子结构性质：若x,y是C中频率最低的两个字符。令$C^\prime$为$C$去除x和y，然后加入字符z，且$z.freq = x.freq + y.freq$。$T^\prime$为$C^\prime$的一个最优前缀码的编码树，将$T^\prime$中的叶子节点z替换为一个以x,y为子节点的内部节点得到树$T$。则T为C的一个最优前缀码的编码树。

证明思路，反证法，假设C存在一个$T^{\prime\prime}$比T更优，由贪心选择性质可知x,y为兄弟节点（不一定吧？书中用不失一般性，不太理解），将$T^{\prime\prime}$中的x,y替换为z得到$T^{\prime\prime\prime}$，结果得到$B(T^{\prime\prime\prime}) < B(T^{\prime})$，与$T^\prime$为$C^\prime$的一个最优前缀码的编码树的假设相悖。

涉及算法

点击查看实现

• 活动选择问题
• 赫夫曼编码

留言與分享

15. 动态规划

动态规划和分治法相似，都是通过组合子问题的解来求解原问题的解。但不同的是分治法将问题划分为不相交的子问题，递归的求解子问题，再将它们的解组合起来，求出原问题的解。相反，动态规划应用于子问题相互重叠的情况，即不同的子问题具有公共的子子问题。这种情况下，用分治法会反复的求解相同的子问题，造成时间浪费。动态规划对每个子问题只求解一次，并且保存在表格中，从而避免重复计算。

动态规划通常用于求解最优化问题。

我们通常按照如下步骤设计一个动态规划算法：

1. 刻画一个最优解的结构特点
2. 递归的定义最优解的值
3. 计算最优解的值，通常采用自底而上的方法
4. 利用计算出的信息构造一个最优解

如果只需要最优解的值，而非最优解本身，可忽略第四步。如果需要做第四步，往往需要在执行第三步时维护一些额外信息。

钢条切割

钢条切割问题：给定一段长度为n英寸的钢条和一个价格表p_i(i=1,2…n)，求切割方案使得收益r_n最大。

考虑当n=4时，所有的切割方案如下，我们发现当把4英寸的钢条切割为两段2英寸长的钢条时产生的收益最大。

当区分1+3切法和3+1切法时，n英寸的钢条有2^n-1种切割方案。
我们可以手动得到r_i，及对应切割方案：

对于r_n，我们可以用更短的钢条的最优切割收益来描述它：
$r_n = max(p_n, r_1+r_{n-1}, r_2+r_{n-2},...,r_{n-1}+r_1) .$

这时我们找到了钢条切割问题的最优子结构。每次切割后我们可以把切开的两段钢条看做独立的钢条切歌问题。一段钢条的最优解为所有切割方案中最大的。
我们可以将r_n的式子简化为：$r_n=\max\limits_{1\le i\le n}(p_i+i_{n-i})$ .

自顶而下的递归实现

根据上面的式子，我们可以得到下面这种自顶而下的递归实现（p为价格数组）：

这种实现的递归式为：$T(n) = 1 + \sum\limits_{j=0}^{n-1}T(j).$
我们可以求得$T(n) = 2^i$，这个时间是非常惊人的，而这么大的原因就是因为包含了非常多的重复计算。下图是递归树，红线圈起来的都是重复计算部分。

使用动态规划方法

• 带备忘的自顶而下法
  与前面的递归方法类似，但是用了一个数组r[1…n]记录最优解，如果发现r[n]已经有记录，则无需重复计算，直接返回该值。
  
• 自底而上法
  自底而上法采用子问题的自然顺序，依次求解规模为j=1,2…,n的子问题。
  

无论是带备忘的自顶而下法还是自底而上法，其中都是用了一个额外数组来存储子问题的最优解，从而避免重复计算子问题的解。动态规划付出额外的内存空间来节省时间，是典型的时空权衡的例子。

子问题图

当思考动态规划问题时，弄清所涉及子问题及子问题之间的依赖关系很关键。
问题的子问题图准确的表达了这些信息：

途中结点表示钢条切割的子问题规模，箭头表示子问题之间的依赖关系。
有了子问题图，我们便可轻易的得知自底而上的顺序应该如何。
而且子问题图G=(V, E)还可以帮助我们确定动态规划算法的运行时间。由于每个子问题只求解一次，因此算法运行时间等于每个子问题求解时间之和。通常，一个子问题的求解时间与子问题图中对应顶点的度（射出边的数目）成正比，而子问题的数目等于子问题图的顶点数。因此，通常情况下，动态规划算法的运行时间与顶点和边数量呈线性关系。

重构解

前面的方法只给出的如何求出最优解的值，但并未返回解本身。
下面是基于自底而上动态规划算法的扩展，其中用数组s[0…n]存储n英寸钢条的最优切割方案，s[i]存储的数字表示，i英寸的钢条切割为s[i]英寸和i-s[i]英寸时最优，其中s[i]英寸部分不再切割，i-s[i]英寸部分的最优解需要递归查看。
所以可以通过PRINT-CUT-ROD-SOLUTION来打印出n英寸钢条的一个最优切割方案。

矩阵链乘法

给定一个n个矩阵的序列（矩阵链）$<A_1, A_2, ..., A_n>$，我们希望求他们的乘积$A_1A_2...A_n$，由矩阵乘法的结合律我们可知给乘积链随意的加上括号不会改变最终的结果，但是不同的括号会影响矩阵链乘法的时间。
以矩阵链$<A_1, A_2, A_3>$为例，三个矩阵的规模分别是10×100、100×5和5×50。如果按照$((A_1A_2)A_3)$的顺序计算，$A_1A_2$（规模10×5）需要做10×100×5=5000次标量乘法，再与$A_3$相乘又需要10×5×50=2500次标量乘法，总共需要7500次标量乘法。如果按照$(A_1(A_2A_3))$的顺序计算，$A_2A_3$（规模100×50）需要做100×5×50=25000次标量乘法，A₁再与之相乘又需要10×100×50=50000次标量乘法，总共需要75000次标量乘法。因此，按照$((A_1A_2)A_3)$的顺序比按照$(A_1(A_2A_3))$的顺序快了10倍。
所以引入了矩阵链乘法问题：给定n个矩阵的链$<A_1, A_2, ..., A_n>$，矩阵i的规模为$p_{i-1} \times p_i$，求一个完全括号化方案，使得计算乘积$A_1A_2...A_n$所需的标量乘法次数最少。

我们可以用下面这个递归式表示矩阵链完全括号化方案的个数，这个递归式的解为$\Theta(2^n)$，所以穷举法是不可行的。

$P(n) = \left \{ \begin{aligned} & 1 & if \quad n = 1, \\ & \sum_{k=1}^{n-1}P(k)P(n-k) & if \quad n \ge 2, \end{aligned} \right.$

动态规划法

1. 刻画一个最优解的结构特点
   用$A_{i..j}$表示$A_iA_{i+1}...A_j$。我们可以将$A_{i..j}$分解为子问题$A_{i..k}$和$A_{k+1..j}$，即$A_iA_{i+1}...A_j -> (A_iA_{i+1}...A_k)(A_{k+1}A_{k+2}...A_j)$。然后在对子问题递归求解。
2. 递归的定义最优解的值
   用$m[i,j]$表示$A_{i..j}$所需标量乘法次数的最优值。假设$A_{i..j}$的最优分割点为$k(i<=k<j)$，那么$m[i,j] = m[i,k] + m[k+1,j] + p_{i-1} p_{k} p_{j}$，然而这个分割点我们是不知道的，但是肯定在i和j之间，我们只需找到所有可能中的最小值，所以我们有如下递归式：