旭穹の陋室

如果你正在读这本书，那么在你的CS旅程中的某个地方，你会有那么一刻，第一次开始关心代码的效率。

我在高中的时候，我意识到制作网站和做实用编程并不能让你获得大学的入场券，于是我进入了算法编程奥林匹克竞赛的世界。我是一个合格的程序员，尤其是对于一个高中生来说，但我之前从来没有想过我的代码需要多长时间才能执行完。但突然之间，问题开始变得重要了：现在竞赛中的每个问题都有严格的时间限制。我开始思考一秒钟能做多少操作。

在思考这个问题时，我还不太了解计算机架构。但我也不需要精确答案——我可以使用经验法则。我的想法是：“2-3GHz主频意味着每秒执行20到30亿条指令，在一个简单的循环中，对数组元素执行一些操作，我还需要增加循环计数器，检查循环结束条件，进行数组索引等这些辅助操作，因此，让我们为每一个有效指令增加3-5个指令的余量，最终我使用每秒5⋅1085\cdot10^85⋅108个操作作为估计值。上面这些陈述并不精确，但计算一下我的算法需要多少操作，并将其除以上面这个数字，对于我来说也是一个很好的经验法则。

当然，真正的答案要复杂得多，而且与操作类型本身强相关。对于指针追逐（pointer chasing，由指针链在一起的数据结构，如链表）这个估计值可以低到10710^7107，而对于SIMD加速后的线性代数运算，又可以高达101110^{11}1011。为了证明这些显著的差异，我们将使用不同的语言实现矩阵乘法，并深入研究计算机如何执行它们。

在最底层，计算机执行由二进制编码的指令（instructions）组成的机器码（machine code），这些指令用于控制CPU。它们是特定的、古怪的，需要付出大量的脑力才能使用好它，所以人们在发明计算机后首先做的事情之一就是发明编程语言（programming languages），它抽象了计算机操作的一些细节，从而简化了编程过程。

编程语言基本上只是一个接口。用它编写的任何程序都只是一个更好的更高级的表示，在某些时点仍然需要转换成在CPU上执行的机器码，下面有几种不同的方法来做这件事：

• 从程序员的角度来看，有两种类型的语言：在执行之前会进行预处理（编译）的编译语言（compiled language），和在运行时使用被称为解释器（interpreter）的单独程序进行执行的解释语言（interpreted language）。
• 从计算机的角度来看，也可分为两种类型：直接执行机器代码的原生语言（native language）和依赖某种运行时（runtime）执行的托管语言（managed language）。

由于在解释器中运行机器代码没有意义，因此总共有三种语言：

• 解释语言，如Python、JavaScript或Ruby。
• 带运行时的编译语言，如Java、C#或Erlang（以及在其VM上运行的语言，例如Scala、F#或Elixir）。
• 编译的本地语言，如C、Go或Rust。

执行计算机程序没有所谓“正确”的方法：每种方法都有其优点和缺点。解释器和虚拟机提供了灵活性，并实现了一些不错的高级编程功能，如动态类型、运行时代码更改和自动内存管理，但这些都会带来一些不可避免的性能问题，我们马上将讨论这些问题。

下面是一个纯Python实现的1024×10241024 × 10241024×1024的矩阵乘法的示例：

import time
import random

n = 1024

a = [[random.random()
      for row in range(n)]
      for col in range(n)]

b = [[random.random()
      for row in range(n)]
      for col in range(n)]

c = [[0
      for row in range(n)]
      for col in range(n)]

start = time.time()

for i in range(n):
    for j in range(n):
        for k in range(n):
            c[i][j] += a[i][k] * b[k][j]

duration = time.time() - start
print(duration)

此代码运行耗时630秒。超过10分钟！

让我们试着正确看待这个数字。运行这段代码的CPU的时钟频率是1.4GHz，这意味着它可以每秒执行1.4⋅1091.4\cdot10^91.4⋅109个时钟周期，整个计算总计101510^{15}1015个时钟周期，最内层循环中的每个乘法操作大约消耗880个时钟周期。（译者：这里的101510^{15}1015不知如何得到的）。

当我们试着分析Python解释器需要做哪些事情来理解上面这段代码的意思时，你便不会再感到奇怪了：

• 它解析表达式c[i][j] += a[i][k] * b[k][j]；

• 尝试找出a、b和c是什么，并在带有类型信息的特殊哈希表中查找它们的名称；

• 知道a是一个列表，获取其[]运算符，检索a[i]的指针，知道它也是一个列表。再次获取其[]操作符，获取a[i][k]的指针，然后获取元素本身；

• 查找元素的类型，确定是float类型，并获取*运算符的实现方法；

• 对b和c执行相同的操作，最后将结果+=给c[i][j]。

当然，Python等这些被广泛使用的语言的解释器已经做了很好的优化，在重复执行相同代码时，它们可以跳过其中的一些步骤。但是，由于语言设计本身的问题，一些相当重要的开销是不可避免的。如果我们摆脱了所有这些类型检查和指针追逐，在暂时不考虑乘法本身的“成本”是多少的前提下，也许我们可以使每个乘法所需的时钟周期数无限接近1。

相同的矩阵乘法运算，但在Java中实现：

import java.util.Random;

public class Matmul {
    static int n = 1024;
    static double[][] a = new double[n][n];
    static double[][] b = new double[n][n];
    static double[][] c = new double[n][n];

public static void main(String[] args) {
        Random rand = new Random();

for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                a[i][j] = rand.nextDouble();
                b[i][j] = rand.nextDouble();
                c[i][j] = 0;
            }
        }

long start = System.nanoTime();

for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                for (int k = 0; k < n; k++)
                    c[i][j] += a[i][k] * b[k][j];

double diff = (System.nanoTime() - start) * 1e-9;
        System.out.println(diff);
    }
}

这段代码需要10秒执行完，相当于每次乘法大约消耗13个CPU周期，比Python快63倍。考虑到我们需要从内存中非顺序地读取b的元素（译者：会有Cache Miss），运行时间大致与预期的一样。

Java是一种编译语言，但不是原生语言。程序首先编译成字节码（bytecode），然后由虚拟机（JVM）解释。为了实现更高的性能，代码中经常执行的部分，例如最内部的for循环，在运行时被编译成机器码，然后在几乎没有额外开销的情况下执行。这种技术称为实时编译（just-in-time compilation）。

JIT编译不是Java语言本身的一个特性，而是它的一种实现。还有一个名为PyPy的JIT编译的Python版本，它在不做任何更改的前提下，需要大约12秒的时间来执行上面的代码。

现在轮到C语言了：

#include <stdlib.h>
#include <stdio.h>
#include <time.h>

#define n 1024
double a[n][n], b[n][n], c[n][n];

int main() {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            a[i][j] = (double) rand() / RAND_MAX;
            b[i][j] = (double) rand() / RAND_MAX;
        }
    }

clock_t start = clock();

for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                c[i][j] += a[i][k] * b[k][j];

float seconds = (float) (clock() - start) / CLOCKS_PER_SEC;
    printf("%.4f\n", seconds);

return 0;
}

当你用gcc -O3编译它时，需要9秒。

这看起来并不是一个很大的提升——相比Java和PyPy，1-3秒的优势可以归因于JIT编译带来的额外时间——但我们还没有充分利用更好的C编译器环境。如果我们加上-march=native和-ffast-math编译选项时，时间会骤降到0.6秒！

这里发生的是，我们向编译器传达了我们正在运行的CPU的准确型号（-march=native），并赋予它重新排列浮点计算的自由（-ffast-math），因此编译器可以使用向量化实现加速。

如果不对源代码进行大的改动，仅通过调整PyPy和Java的JIT编译器以实现相同的性能并不是不可能，但对于直接编译为原生代码的语言来说，这当然更容易实现。

BLAS（Basic Linear Algebra Subprograms，基础线性代数程序集）

最后，让我们看一看经过专门优化后的实现的水平。我们将测试一个广泛使用的优化线性代数库OpenBLAS。使用它的最简单方法是回到Python，从numpy中调用它：

import time
import numpy as np

n = 1024

a = np.random.rand(n, n)
b = np.random.rand(n, n)

start = time.time()

c = np.dot(a, b)

duration = time.time() - start
print(duration)

现在它只需要约0.12秒：比自动向量化的C语言版本快约5倍，比我们最初的Python实现快约5250倍！

你通常不会看到如此显著的改善。现在，我们不准备告诉你这是如何实现的。OpenBLAS中稠密矩阵乘法的实现通常是为每个架构单独定制的5000行手写汇编。在后面的章节中，我们将逐一解释所有相关的技术，然后回到这个例子，使用不到40行的C语言代码来开发我们自己的BLAS级实现。

这里的关键点是，使用原生、低级语言并不一定能给你带来性能；但它确实能让你控制性能。

对“每秒N次操作”，再简单补充一点，许多程序员有一个误解，使用不同的编程语言会对这个N乘上某个系数。以这种方式思考并比较语言的性能并没有多大意义：编程语言从根本上来说只是一种工具，它通过方便的抽象来获取对性能的一些控制。不管运行环境如何，尽可能利用硬件提供的机会来优化性能表现，在很大程度上仍然是程序员的本职工作。