数值 Python

作者: Joey Bernard

在过去的几个月里,我一直在介绍用于科学计算的不同软件包。在接下来的几篇文章中,我将专注于使用 Python 来为您自己的科学问题提出您自己的算法。Python 似乎正在完全接管科学界以开发新代码,因此扎实掌握基础知识以便您可以为自己的问题构建解决方案是个好主意。

在本文中,我从 NumPY 开始。我的下一篇文章将介绍 SciPy,然后在接下来的文章中我将介绍一些更复杂的模块。

所以,让我们从几乎所有其他科学模块都构建于其上的根模块 NumPY 开始。Python 开箱即用支持实数和整数。您还可以使用列表、集合等创建复杂的数据结构。这使得编写算法来解决科学问题变得非常容易。但是,只是天真地投入,而不注意幕后发生的事情,将导致效率低下的代码。这对于所有编程语言都是如此,而不仅仅是 Python。大多数科学代码都需要从硬件中挤出最后可用的循环。关于 Python 需要记住的一件事是,它是一种动态语言,其中几乎所有函数和运算符都是多态的。这意味着 Python 实际上并不知道在硬件级别需要做什么,直到它遇到该操作为止。不幸的是,这排除了任何可以通过重新排列操作来利用它们在内存和缓存中的存储方式进行优化的可能性。

Python 的一个导致更大问题的特性是多态性。在这种情况下,Python 需要检查任何运算符或函数的操作数以查看其类型,确定此特定操作数或函数是否可以处理这些数据类型,然后使用操作数或函数的正确形式来执行实际操作。在大多数情况下,这实际上不是问题,因为现代计算机已经变得如此之快。但是在许多科学算法中,您最终会对成千上万甚至数百万个数据点应用相同的操作。一个简单的例子就是取前 100,000 个数字的平方


import time
a = range(100000)
c = []
starttime = time.clock()
for b in a:
   c.append(b*b)
endtime = time.clock()
print "Total time for loop: ", (endtime - starttime)

这个小程序使用 range 函数来创建前 100,000 个整数的列表。您需要导入 time 模块才能访问系统时钟以进行计时测量。在我的桌面上运行此程序大约需要 0.037775 秒(请记住始终进行多次测量,并取平均值)。之所以花费这么多时间,是因为在循环的每次迭代中,Python 都需要检查 b 变量的类型,然后才能尝试使用乘法运算符。如果您的算法甚至更复杂,您该怎么办?

这就是 NumPY 的用武之地。NumPY 引入的关键元素是 N 维数组对象。Python 列表的巨大灵活性,允许各种不同类型的元素,但这是有计算成本的。NumPY 数组通过引入限制来应对这种成本。数组可以是多维的,并且它们必须都是相同的数据类型。完成此操作后,您可以通过利用您已经知道元素类型这一事实来开始采取一些捷径。它为通用运算符和函数添加了额外的重载函数,以帮助优化数组的使用。

所有普通的算术运算符都以元素方式对 NumPY 数组进行运算。因此,要获得数组中元素的平方,只需 array1 * array1 即可。

NumPY 还使用用 C 或 FORTRAN 编写的外部标准优化库来处理对这些数组数据类型的许多实际操作。这由 BLAS 或 lapack 等库处理。Python 只是对有问题的每个数组进行初始检查,然后将它们作为单个对象传递给外部库。外部库完成所有繁重的工作,然后返回一个包含结果的单个对象。这消除了 Python 在使用上面的循环表示法时检查每个元素的需要。使用 NumPY,之前的示例变为


import numpy
import time
a = numpy.arange(1000000)
starttime = time.clock()
c = a * a
endtime = time.clock()
print "Total time used: ", (endtime - starttime)

运行此玩具代码的平均运行时间为 0.011167 秒,用于此平方运算。因此,时间缩短了三分之一,并且通过消除循环结构简化了代码的可读性。

到目前为止,我只处理了一维数组,但 NumPY 也同样容易地支持多维数组。如果您想定义一个二维数组或矩阵,您可以像这样设置一个小数组


a = numpy.array([[1,2,3,4], [2,3,4,5]])

基本上,您正在创建一个列表的列表,其中每个子列表都是矩阵的每一行。这只有在每个子列表的大小相同时才有效,也就是说,每行的列数相同。您可以使用属性 a.size 获取此矩阵中元素的总数。矩阵的维度存储在属性 a.shape 中。在这种情况下,大小为 8,形状为 (2, 4),即两行四列。前面示例中的数组是什么形状?如果您继续检查,您应该看到形状为 (1000000)。这些数组的其他关键属性是

  • ndim:维度数量。

  • dtype:元素的数据类型。

  • itemsize:每个元素的大小(字节)。

  • data:存储实际数据的缓冲区。

您还可以重塑数组。因此,如果您想采用前面 100,000 个整数的示例并将其转换为三维数组,您可以这样做


old_array = numpy.arange(100000)
new_array = old_array.reshape(10,100,100)

这将为您提供一个新的 3-D 数组,布局成 10x100x100 元素立方体。

现在,让我们看看一些可用于数组的其他函数。如果您还记得之前的知识,所有标准算术运算都被重载为一次对一个数组元素进行运算。但是,大多数矩阵编程语言使用乘法元素来表示矩阵乘法。当您开始使用 Python 时,请记住这一点。要获得真正的矩阵乘法,您需要使用 dot() 函数。如果您有两个矩阵 A 和 B,您可以使用 numpy.dot(A, B) 将它们相乘。

许多标准数学函数,如余弦、正弦、平方根等,都由 NumPY 提供,称为通用函数。与算术运算符一样,这些通用函数以元素方式应用于整个数组。在科学中,使用了几个常用函数。您可以使用对象函数 a.transpose() 获取矩阵的转置。如果您需要获取矩阵的逆矩阵,则可以使用模块函数 inv(),因此您将执行 numpy.inv(a)。迹也是一个模块函数,由 numpy.trace(a) 给出。

甚至有更复杂的函数可用。您可以使用 NumPY 求解方程组。如果您有一个由 a 给出的系数矩阵,以及一个表示方程右侧的数字向量 y,您可以使用 numpy.solve(a,y) 求解该系统。在许多情况下,您可能对查找给定系统的特征值和特征函数感兴趣。如果是这样,您可以使用 numpy.eig(array1) 来获取这些值。

我想在这里看到的最后一件事是一个类,它提供了更多的快捷方式,但代价是更多的限制。矩阵(二维数组)非常普遍,以至于 NumPY 提供了一个特殊的类来尽可能优化使用它们的操作。要创建一个新矩阵,例如 2x2 矩阵,您需要编写


A = numpy.matrix('1.0, 2.0; 3.0, 4.0')

现在,您只需使用 A.T 即可获得转置。同样,逆矩阵可以使用 A.I 找到。当您使用矩阵对象时,乘法运算将执行实际的矩阵乘法。因此,给定两个矩阵对象 A 和 B,您可以使用 A*B 执行矩阵乘法。solve 函数仍然可以按预期工作于使用矩阵对象定义的方程组。NumPY 网站上提供了许多技巧和窍门,非常值得一看,尤其是在您刚开始时。

这个简短的介绍应该让您开始考虑将 Python 作为“真实”数值计算的可行选择。NumPY 模块为构建复杂的科学工作流程提供了非常强大的基础。下个月,我将介绍一个可用的模块 SciPY。在那之前,请尽情体验 NumPY 提供的所有原始数值计算可能性。

加载 Disqus 评论