向量运算入门

本节首先介绍向量操作的基本概念，这些概念对于理解高性能计算是极为关键的。

以Python中的循环完成向量计算

当今的科学和工程问题中，待处理的数据通常表示向量或矩阵，或者说“多维数组”。

假设我们有一个地理信息数据集，它包含大量位置的经纬度信息，我们现在要从中查找出离指定位置最佳的地点。若这些信息以tuple列表保存，我们可以用循环来完成上述功能：

def closest(position, positions):
    x0, y0 = position
    dbest, ibest = None, None

    for i, (x, y) in enumerate(positions):
        d = (x - x0) ** 2 + (y - y0) ** 2
        if dbest is None or d < dbest:
            dbest, ibest = d, i

    return ibest

这是最“原始”的Python实现，通过%timeit来测试它的性能：

import random
positions = [(random.random(), random.random()) for _ in xrange(10000000)]
%timeit closest((.5, .5), positions)

1 loops, best of 3: 7.19 s per loop

整个循环需要7秒多，来看看有无更快的方法。

使用NumPy中的数组

上面的循环中包含大量相同操作（计算与固定点的距离），NumPy提供的多维数组类型非常适合这种场景的需求。

数组的元素类型必须是相同的（data type, dtype）。在多维情形下，元素的存储顺序可有不同选择，如row-major order（C-order）、column-major order（Fortran-order）。NumPy默认采用row-major，但可以在创建数组时指定order。

在性能方面，NumPy由C实现，不需要依赖于Python的循环。此外，现代的CPU实现了向量指令，可使用更大的寄存器（register），从而包含更多的浮点数。当然，这取决于NumPy的编译选项。现在用数组和向量操作重新实现：

positions = rand(10000000, 2)
x, y = positions[:, 0], positions[:, 1]

%timeit distances = (x - 0.5) ** 2 + (y - .5) ** 2
# 10 loops, best of 3: 201 ms per loop

%timeit ibest = distances.argmin()
# 100 loops, best of 3: 10.3 ms per loop

只需要0.21秒左右。

我们有足够的理由说，在数值计算中，尽可能用数组来代替Python循环。