二十世纪十大算法

本文同时作为北京大学《数据结构与算法》课程论文提交。转载请注明出处。

算法（Algorithm）一词来源于公元九世纪的波斯数学家、天文学家及地理学家花拉子米（Muhammad ibn Musa al-Khwarizmi）（「Algorithm」一词本身兴起于 1950 年左右。虽然对其起源的考证有很多不同的观点，但一种较为可靠的说法是「Algorism」，由此可以追溯到花拉子米。比较有趣的是，「代数」（Algebra）一词亦与花拉子米有关 —— 这来源于他的一本与代数没有什么联系的著作）。而世界上最早的算法，出现的甚至比这更早 —— 欧几里得的辗转相除法，诞生于公元三世纪。今日我们所说的算法，可以参见高德纳（Donald Ervin Knuth）在他的著作《计算机程序设计艺术（The Art of Computer Programming）》里的描述（可参阅《计算机程序设计艺术（中文版）第一卷第三版：基本算法》的第 3 至 4 页）：

• 输入：一个算法必须有零个或以上输入量。
• 输出：一个算法应有一个或以上输出量，输出量是算法计算的结果。
• 明确性：算法的描述必须无歧义，以保证算法的实际执行结果是精确地匹配要求或期望，通常要求实际运行结果是确定的。
• 有限性：依据图灵的定义，一个算法是能够被任何图灵完全系统模拟的一串运算，而图灵机只有有限个状态、有限个输入符号和有限个转移函数（指令）。而一些定义更规定算法必须在有限个步骤内完成任务。
• 有效性：又称可行性。能够实现，算法中描述的操作都是可以通过已经实现的基本运算执行有限次来实现。

本文将要介绍的二十世纪十大算法，最早由 Jack Dongarra 和 Francis Sullivan 在 2000 年提出（SIAM 上的一篇报道使其广泛传播）。他们选择了「对 20 世纪科学与工程的发展和实践影响最大的 10 种算法」，分别是：

• 蒙特卡洛方法（Monte Carlo method），John von Neumann, Stan Ulam 和 Nick Metropolis 在 1946 年提出
• 线性规划的单纯形法（Simplex method for linear programming），George Dantzig 在 1947 年提出
• Krylov 子空间迭代法（Krylov subspace iteration methods），Magnus Hestenes, Eduard Stiefel 和 Cornelius Lanczos 在 1950 年提出
• Householder 形式化矩阵计算（Decompositional approach to matrix computations），Alston Householder 在 1951 年提出
• 优化的 Fortran 编译器（Fortran optimizing compiler），John Backus 在 1957 年提出
• QR 分解（QR algorithm），J.G.F. Francis 在 1959-61 年提出
• 快速排序（Quicksort），Tony Hoare 在 1962 年提出
• 快速傅里叶变换（Fast Fourier transform），James Cooley 和 John Tukey 在 1965 年提出
• 整数关系探测算法（Integer relation detection algorithm），Helaman Ferguson 和 Rodney Forcade 在 1977 年提出
• 快速多极算法（Fast multipole algorithm），Leslie Greengard 和 Vladimir Rokhlin 在 1987 年提出

下面将逐一介绍这些算法。为了明确地表达思想，这里会给出一些示例代码。

蒙特卡洛方法

蒙特卡洛方法是 1940 年代中期由于科学技术的发展和电子计算机的发明，而提出的一种以概率统计理论为指导的数值计算方法。是指使用随机数（或更常见的伪随机数）来解决很多计算问题的方法。

20 世纪 40 年代，冯・诺伊曼，斯塔尼斯拉夫・乌拉姆和尼古拉斯・梅特罗波利斯在洛斯阿拉莫斯国家实验室为核武器计划工作时，发明了蒙特卡洛方法。乌拉姆的叔叔经常在摩纳哥的蒙特卡洛赌场输钱，而蒙特卡洛方法正是以概率为基础的方法，故因此得名。（这一段历史来自维基百科中的「蒙地卡羅方法」词条）

蒙特卡洛方法又称随机抽样或统计试验方法。其核心正是取样的思想，因而它常用来通过随机模拟计算概率或者期望。除此之外，它的另一个重要作用是计算积分：在积分区间上随机撒点，通过计算这些点上的函数值，再进行平均，得到积分结果。在计算高维的积分时，这将会特别有用 —— 例如平衡态统计物理中的伊辛模型（Ising model）。在组态足够多的情况下，常规的积分算法没有任何希望完成路径积分。而用蒙特卡洛方法做重点取样，却可以较少的代价完成计算。这时的问题就转化为如何寻找要计算的组态。

一种可行的方法是，按照一定规则以及通过随机演化构造出下一个组态，在这个过程中，每个组态只依赖于它之前的那个组态，这个过程被称为马尔可夫过程。

下面是一个用 Metropolis 算法（Metropolis algorithm）来实现马尔可夫过程的例子：

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

import numpy as np

x = np.ones(20)

def S(x):
  return np.sum(np.power(np.diff(x), 2)) \
    + np.sum(np.power(x, 2)) / 4

def update(x):
  print(x)
  for j in range(20):
    s = S(x)
    xi = np.random.random() * 2.8 - 1.4
    tmp = x[j]
    x[j] += xi
    delta = S(x) - s
    if delta < 0:
      continue
    eta = np.random.random()
    if np.exp(-delta) > eta:
      continue
    x[j] = tmp

for i in range(1000):
  update(x)

这是一个数值路径积分的案例，其物理背景是谐振子。其中为位移，为格点上的作用量，不难得到：

其中为谐振子势，可以取 。

先给定一组初始的组态 满足周期性边界条件。这里取为 20，初始组态为 。依次从到进行更新。对于第个元素的更新：在区间 之间产生一个随机数，然后更新 ，并计算 。如果 ，接收新的；如果 ，在之间产生一个随机数，如果 ，接收新的值，否则放弃。这里不能取得太大，否则接收率很低；也不能太小，否则每次更新组态变动很小。这样，用 Metropolis 算法就可以产生足够多的组态。

线性规划的单纯形法

所谓线性规划，简单的说，就是在给定一组线性约束条件（包括等式和不等式）下，最优化一个线性函数。

假设有 n 个变量和 m 个约束。线性规划的标准形式如下：

几何上，线性约束条件的集合相当于一个凸包或凸集，叫做可行域。而单纯形是维中的 个顶点的凸包。单纯形算法利用凸包的顶点构造一个可能的解，然后沿着凸包的边走到目标函数值更高的另一个顶点，不断执行此过程，直至到达最优解为止。（除了可能遇到的病态的情况，需要进行额外的判断）

单纯形法的一般解题步骤可归纳如下：

• 把线性规划问题的约束方程组表达成典范型方程组，找出基本可行解作为初始基可行解。
• 若基本可行解不存在，即约束条件有矛盾，则问题无解。
• 若基本可行解存在，从初始基可行解作为起点，根据最优性条件和可行性条件，引入非基变量取代某一基变量，找出目标函数值更优的另一基本可行解。
• 按步骤 3 进行迭代, 直到对应检验数满足最优性条件（这时目标函数值不能再改善），即得到问题的最优解。
• 若迭代过程中发现问题的目标函数值无界，则终止迭代。

Krylov 子空间迭代法

Krylov 子空间迭代法是一个求多维函数的极值的算法。它包含两部分：C. Lanczos 提出的构建 Krylov 子空间的方法；以及由 M. Hestenes 和 E. Stiefel 提出来的共轭梯度（Conjugate gradient）算法。

我们考虑线性方程组 的求解问题。假设 是阶实对称正定矩阵。求解线性方程组的其中一种思路是将它转化为在维向量空间求函数最小值的问题，这个函数是下面的 元二次函数

这个最小值点必定满足条件 ，不难证明，这样得到的 n 个方程就是线性方程组

。这说明方程组的解就是 元二次函数 的极小点，现在的问题是如何有效地寻找多元二次函数的最小值。

给定初始向量 ，第一步选负梯度方向为搜索方向，即 ，于是有

接下来，共轭梯度法的想法是，搜索方向不仅仅考虑，而是在过点 由向量 和 所张成的平面 内寻找函数的最小值，记为

直接计算可得到两个方程

其解写为。假设要求的极小值为 ，那么可以令

新的搜索方向则为

令 ，得到

那么就可以直接写成

另外，搜索步长是由 给出的，我们可以令 ，有

下一个迭代的解为

而新的残差定义为

可以继续迭代下去，直到机器精度的限制。

综上所述，Krylov 子空间迭代法用伪代码表示为

Householder 形式化矩阵计算

Householder 形式化矩阵中的一个重要内容是 Householder 变换。对于维的欧氏空间，我们可以对指定的矢量 进行镜面反射（归根结底来说，就是让矢量与「镜面」平行的分量不变，而与之垂直的分量反向），而这个「镜面」是一个 维的子空间，可以由它的法向量唯一确定。换句话说，对于 ，可以构造如下的矩阵：

其中 是矢量的欧氏模方。当我们将这个矩阵乘以时，即 ，就实现了镜面反射。

这个变换就是 Householder 变换，相应的称为 Householder 矢量，称为与相对应的 Householder 矩阵。由于它具有很好的性质，实际的计算通过矢量的内积即可完成，具有的复杂度；而无需计算更复杂的矩阵乘法。

进一步的，我们可以把 通过 Householder 变换，变换成为只有第 个方向有非零分量的矢量：

这里只需要把 取为 ，而矢量取为 。重复此操作，即可将矩阵约化为上三角形式。

下面是 Householder 变换的一种实现。

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

import numpy as np

def householder(matrix):
  N = len(matrix)
  Q = 1
  R = matrix.copy()
  for i in range(N - 1):
    x = np.array([[R[j, i] for j in range(i, N)]])
    sigma = np.linalg.norm(x)
    x[0, 0] -= sigma
    v = x / np.linalg.norm(x)
    H_hat = np.identity(N - i) - 2 * v * v.transpose()
    H = np.zeros((N, N))
    for m in range(N):
      for n in range(N):
        if m < i and n < i and m == n:
          H[m, n] = 1
        elif m >= i and n >= i:
          H[m, n] = H_hat[m - i, n - i]
    Q *= H
    R = H * R
    print(R)
  return [Q, R]

mat = np.matrix([[0, -15, 14], [4, 32, 2], [3, -1, 4]])
print(householder(mat))

不过在实际计算中，更常见的是通过 Hessenberg-Householder
约化 (Hessenberg-Householder reduction) 将矩阵变换为 Hessenberg 形式，再利用后面的 QR 分解进行进一步处理。过程如下：

优化的 Fortran 编译器

严格来说，这并不是一个「算法」，而是一套完整的编译系统。

Fortran 语言是为了满足数值计算的需求而发展出来的。早在 1953 年，IBM 的 John Backus 便提出了开发全新的计算机语言的想法。四年后，他所领导的团队成功开发出了第一套 Fortran 语言，这也是世界上第一个被正式采用并流传至今的高级编程语言。

Fortran 可以直接对矩阵和复数进行运算，这点被 Matlab 继承。对于矩阵中「行优先」还是「列优先」的问题，甚至可以在 numpy 中看到 Fortran 的影子。

在 1998 年出版的 IEEE Annals of the History of Computing 中，John Backus 回忆道：

(The compiler) produced code of such efficiency that its output would startle the programmers who studied it.

虽然今日 Fortran 已经不再是数值分析和科学计算的绝对选择（C++，Python，Matlab 和 R 等后起之秀都在诸多领域发挥着作用，一些语言在效率上与 Fortran 不相上下，而开发成本更低）；但它却是那个时代浓墨重彩的一笔，为推动计算机科学的发展做出了巨大的贡献。

QR 分解

QR 算法是目前各界应用的最为广泛的数值求解本征值和本征矢量的算法之一。为了求出一个矩阵的全部本征值和本征矢量，一种思路是这样的：对于 ，构造正交变换 ，使得是一个上三角矩阵，而是实正交矩阵。正交变换保证了的本征值与完全一致，而上三角矩阵的本征值就是对角元，由此得解。

QR 分解的两种典型方法是 Householder 变换和 Givens 旋转。Householder 变换是前述的 Householder 形式化矩阵计算的一部分。因此这里重点介绍 QR 分解的 Givens 旋转方法。

Givens 旋转的基本原理非常简单，我们取：

称为 Givens 变换矩阵。由此可以消去向量的某个分量，例如：

其中 ，且

通过不断地构造这样的 Givens 变换矩阵，我们可以对矩阵中任何一列进行消去操作。最终可以达到这样的效果：第一列中第一行以下的元素变为 0，第二列中第二行以下的元素变为 0…… 直到得到一个上三角矩阵，完成 QR 分解。

一种实现如下：

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

import numpy as np

def givens(matrix):
  N = len(matrix)
  Q = 1.0
  R = matrix.copy()
  for i in range(N - 1):
    for j in range(i + 1, N):
      if R[i, j] == 0:
        continue
      norm = (R[i, i] ** 2 + R[j, i] ** 2) ** 0.5
      cos = R[i, i] / norm
      sin = R[j, i] / norm
      G = np.identity(N)
      G[i, i] = G[j, j] = cos
      G[i, j] = sin
      G[j, i] = -sin
      Q *= G.transpose()
      R = G * R
      print(R)
  return [Q, R]

mat = np.matrix([[0, -15, 14], [4, 32, 2], [3, -1, 4]])
print(givens(mat))

考虑到在进行 Givens 变换时，许多元素是不变的，这个算法还有优化的余地，即不必直接进行矩阵乘法。（对于 Python 本身来讲，这样也不错了，毕竟 numpy 的效率还是有保证的）不过，它已经足够展现出 Givens 变换实现 QR 分解的过程。

快速排序

快速排序是对冒泡排序的一种改进，由 Tony Hoare 在 1962 年提出。快速排序使用分治法（Divide and conquer, D & C）策略来把一个序列分为两个子序列，然后递归地进行排序。其步骤为：

• 从数列中挑出一个元素，称为「基准」（pivot），
• 重新排序数列，所有比基准值小的元素移到基准前面，所有比基准值大的元素移到基准后面。在这个操作结束之后，该基准就处于数列的中间位置。
• 递归地把小于基准值元素的子数列和大于基准值元素的子数列排序。

递归到最后时，每个子数列的元素个数是零或一，也就是已经排序好了。

该算法的平均时间复杂度为，最坏情况下则为。下面是一种实现：

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

import numpy as np

def quickSort(alist, first, last):
  if first >= last:
    return

  pivot = alist[first]
  left = first + 1
  right = last
  while right >= left:
    while left <= right and alist[left] <= pivot:
      left += 1

    while right >= left and alist[right] >= pivot:
      right -= 1

    if right >= left:
      alist[left], alist[right] = alist[right], alist[left]

  alist[first], alist[right] = alist[right], alist[first]

  quickSort(alist, first, right - 1)
  quickSort(alist, right + 1, last) 

alist = [np.random.randint(0, 100) for i in range(20)]
quickSort(alist, 0, len(alist) - 1)
print(alist)

快速傅里叶变换

快速傅里叶变换，或者说 FFT，是快速计算序列的离散傅里叶变换（Discrete Fourier Transform, DFT）或其逆变换的方法。对于长度为的序列，将这个数据的相因子记为

它实际上是 1 在复平面的次方根。那么其离散傅里叶变换可写为

直接按这个定义求值的复杂度为。我们可以进一步将求和中的奇数和偶数的分开

这个公式又被称为 Danielson-Lanczos 引理 (Danielson-Lanscos lemma)，它将一个长度为的傅里叶变换改写为两个长度为的傅里叶变换的线性组合。它的最大作用在于，这种改写可以继续迭代地进行下去，直到求解的问题变为一个长度为 1 的傅里叶变换 —— 即不需要变换（值得注意的是，在某些情况下，需要提前向序列中填充 0 元素，使其元素个数为 2 的幂次，才能递归进行到长度为 1 的傅里叶变换。从数学上来讲，这不会对结果产生影响，并且由于 FFT 算法在时间复杂度上的优势，大部分情况下这仍然会远远优于传统的 DFT 算法）。这就是 Cooley-Tukey FFT 算法，由 James Cooley 和 John Tukey 提出。可见，这也是一个使用分治法的例子。

一种实现如下：

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

import numpy as np

def DFT(x):
  N = x.shape[0]
  n = np.arange(N)
  k = n.reshape((N, 1))
  M = np.exp(-2j * np.pi * k * n / N)
  return np.dot(M, x)

def FFT(x):
  N = x.shape[0]
  if N == 1:
    return x
  elif N % 2 > 0:
    raise ValueError("Size of x must be a power of 2")
  else:
    X_even = FFT(x[::2])
    X_odd = FFT(x[1::2])
    factor = np.exp(-2j * np.pi * np.arange(N) / N)
    return np.concatenate([
      X_even + factor[:int(N / 2)] * X_odd,
      X_even + factor[int(N / 2):] * X_odd
    ])

x = np.random.random(1024)
print(np.allclose(DFT(x), np.fft.fft(x)))
print(np.allclose(FFT(x), np.fft.fft(x)))

整数关系探测算法

整数关系探测的内容是：给定一组实数，是否存在不全为零的整数，使得： 。

这是一个古老的问题，其历史甚至可以追溯到欧几里德的时代 —— 对于 的情况，基于前面提到的欧几里得的辗转相除法的拓展，即可得到解决。关于任意的情况的研究，在 20 世纪 80 年代得到了大量的成果，这一时期诞生了包括 LLL 算法、HJLS 算法、PSOS 算法和 PSLQ 算法在内的诸多整数关系探测算法实现。现代的 LLL 算法已经可以解决 时的问题。

该算法可用于简化量子场论中的 Feynman 图的计算、对高次多项式进行因式分解和实现 Inverse Symbolic Calculator（即给出一个近似值，反推它可能是由哪些科学常数组合成的）。

快速多极算法

在物理学中，多极展开的方法应用极为广泛。无论是质量分布产生的重力场、电荷分布产生的电势或电场、电流分布产生的磁矢势和磁场还是电磁波的传播等问题，都可以在一些近似下使用多极展开。这些物理量都可以表达为单极项、偶极项、四极项、八极项等的叠加。在近代物理中，一个典型的范例是，通过原子核的外部多极矩与电子轨域的内部多极矩之间的交互作用能量，计算求得原子的原子核外多极矩。由于从原子核的外多极矩可以给出原子核内部的电荷分布，籍此可以研究原子核的形状。

下面给出一个在经典电动力学中，电势多极展开的例子。注意到在原点 的泰勒级数为

将这展开式代入电势的方程式，则可得到

我们关注展开式的前三项。将总电荷（电单极矩）、电偶极矩、电四极矩（Electric quadrupole moment）分别定义为

则电势的电单极矩、电偶极矩、电四极矩等「笛卡儿多极矩」项的总贡献为

多极展开在数值模拟领域用途很多。对于相互作用的粒子组成的物理系统，快速多极算法是高效率运算这系统的能量与作用力常使用的一种方法。具体来说，该方法分解所有粒子为几个小群，每一个小群内的粒子正常地互相作用（即通过全部势能），而小群与小群之间的互相作用则是由其多极矩计算求得。

无论是银河系中的星体的运行，或者蛋白质中的原子间的相互作用，都可以应用快速多极算法进行数值计算。

在撰写此文时，笔者参考了《算法导论》和《计算物理导论》中有关这些算法的介绍。当然，这是一个庞大的主题，试图用以上这些微不足道的文字去阐述其全部内涵是不现实的。它涉及到了诸多领域的内容，而这篇文章只是介绍了一些皮毛，只有通过理论与实践的结合，才能真正体会到一个算法的精妙，体会到这些前辈们的智慧。

算法，让生活更美好。