开源图书《Python完全自学教程》12.4科学计算

老齐

发布于 2022-12-09 20:23:44

1.4K0

发布于 2022-12-09 20:23:44

12.4 科学计算

科学计算是科学、工程等项目中必不可少的，MATLAB 曾风光一时，但它是收费的，并且有“被禁”的风险——坚决反对用盗版软件，“被禁”不是盗版的理由。其实，Python ——开源、免费——是做科学计算的选择之一，它不仅能做 MATLAB 所能做的一切，还能做它不能做的。所以隆重推荐，在科学计算上选用 Python 。

12.4.1 Jupyter

Jupyter 是一款基于浏览器的开源的交互开发环境，常用于科学计算、数据科学、机器学习等业务中。其官方网站是：https://jupyter.org/ ，目前有 JupyterLab 和 Jupyter Notebook ，一般认为 JupyterLab 更趋近于 IDE，功能多于 Jupyter Notebook。下面以 JupyterLab 为例，演示安装和使用方法。

% pip install jupyterlab

使用 pip 安装，安装之后，执行如下指令启动 JupyterLab ：

% jupyter-lab
... （省略部分显示信息）
    Or copy and paste one of these URLs:
        http://localhost:8888/?token=549de63a2fcd53325c5d17b8dd7a00ce4b63766aa368efb1
     or http://127.0.0.1:8888/?token=549de63a2fcd53325c5d17b8dd7a00ce4b63766aa368efb1

一般情况下，会自动打开本地计算机操作系统默认的浏览器，并显示图12-4-1所示效果。

图12-4-1 JupyterLab 默认页面

如果没有打开图12-4-1所示页面，或者打算用其他浏览器，可以复制终端所显示的地址（如：http://localhost:8888/?token=549de63a2fcd53325c5d17b8dd7a00ce4b63766aa368efb1 ）到指定浏览器地址栏，同样能打开图12-4-1所示界面（注意，地址中的 token 项不能缺少）。

JupyterLab 所提供的各项功能，与通常 IDE 类似，此处不做详细说明，有意逐项了解的读者可以自行查阅专门资料。下面仅就如何使用它编写程序和执行所编写的代码给予说明，这是最基本的应用。

在图12-4-1所示界面的 Launcher 标签下，选择 Notebook 中的第一项“Python 3”（读者的开发环境可能与图中所示不同，只要选择“Python 3 ”作为程序执行的驱动即可。然后进入图12-4-2所示页面。

图12-4-2 编辑页面

参照图12-4-3所示的功能，将当前文件改名为 scicomputing.ipynb ，并将左侧文件栏折叠收起。此外，从图示中还可看到其他关于文件和目录的基本操作，供读者参阅。

图12-4-3 基本操作

然后点击菜单 View ，并选中如图12-4-4所示的项目，即可在代码块中显示行号。

图12-4-4 显示代码块中的行号

将鼠标移动到代码块中并单击，如图12-4-5所示，开始输入一行代码，然后回车，输入第二行——注意，这里与 Python 交互模式不同，回车意味着换行，而不是执行当前行代码。输入完图12-4-5所示的所有代码。

图12-4-5 编写代码

按照图12-4-5所示，点击该按钮，可以执行当前代码块；或者按组合键 shift + Return/Enter 执行，效果如图12-4-6所示。

图12-4-6 执行代码块

执行了代码块 [1] 之后，执行结果显示在代码之后，并自动进入下一个代码块。在编辑界面顶部，有各种常见的编辑、控制按钮，把鼠标滑动到按钮上，会显示该按钮的作用。

以上是 JupyterLab 的最基本使用方法，其他复杂功能，可以在应用过程中学习。

12.4.2 第三方库

Python 生态中拥有非常丰富的支持科学计算的第三方库，常用的有 NumPy 、Pandas 、SciPy 、Matplotlib 、SymPy 等，建议读者将这些库依次安装。

% pip install numpy pandas sympy scipy matplotlib

除了能在终端执行安装指令之外，在 JupyterLab 的代码块中也可以执行终端指令，如图12-4-7所示，在代码块中输入 !pip install numpy 并执行——注意前面的 ! 符号。其效果等同于以往在终端执行安装 pip install numpy 指令。

图12-4-7 在代码块中执行安装指令

安装好之后，在代码块中输入如下代码，并执行，即可查看所安装的 NumPy 的版本。

[3]: import numpy as np
     np.__version__
[3]: '1.19.4'       # 上述代码块的输出结果

在数据科学中，引入一些常用的第三方库时，习惯于再命名一个别称（或简称），例如以 np 作为 numpy 的别称，并且这是一种习惯命名，如果非要以 ny 为别称，语法上没有问题，但不是大多数人的习惯——代码更多时候是给人看的。

安装好基础库之后，再列举几个示例（随后几个小节内容），体会 Python 在科学计算中的应用。

12.4.3 矩阵

矩阵不仅在线性代数中占据重要地位，也是科学计算的主角。如果读者还忌惮于当初用纸笔完成有关矩阵计算的痛苦，现在使用 Python 中科学计算的工具包则会体验到无比的畅快。

在 JupyterLab 代码块中输入如下代码（如无特别声明，本节的代码均在 JupyterLab 中调试）。

[4]: import numpy as np
     I_m = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
     I_m
[4]: array([[1, 0, 0],
            [0, 1, 0],
            [0, 0, 1]])

变量 I_m 引用的对象是用 np.array() 创建的二维数组，可以用二维数组表示矩阵—— I_m 表示一个

3\times3

的单位矩阵。

此外，用 NumPy 中的 np.mat() 也能生成矩阵对象：

[5]: np.mat("1 2 3; 4 5 6")
[5]: matrix([[1, 2, 3],
             [4, 5, 6]])

以这两种形式表示的矩阵，在矩阵加法和数量乘法上没有任何差别，例如：

# 二维数组表示的矩阵相加
[6]: a = np.array([[1, 4], [8, 9]])
     b = np.array([[2, 3], [7, 11]])
     a + b
[6]: array([[ 3,  7],
            [15, 20]])

# 矩阵对象相加
[7]: ma = np.mat("1 4; 8 9")
     mb = np.mat("2 3; 7 11")
     ma + mb
[7]: matrix([[ 3,  7],
             [15, 20]])

# 矩阵数量乘法
[8]: print('array: ', 7 * a)
     print("")
     print('mat: ', 7 * ma)
[8]: array:  [[ 7 28]
              [56 63]]

     mat:  [[ 7 28]
            [56 63]]

但是，矩阵乘法会有所不同。首先来看数组 a 和 b ，如果计算 a * b ，其所得并不是它们所表征的矩阵乘法结果。

[9]: a * b    # 数组相乘，不是矩阵乘法
[9]: array([[  2,  12],
            [ 56,  99]])

而矩阵对象 ma 和 mb 直接相乘 ma * mb 的结果是根据矩阵相乘规则所得。

[10]: ma * mb
[10]: matrix([[ 30,  47],
              [ 79, 123]])

如果用数组表示矩阵，且计算矩阵乘法，可以通过 np.dot() 函数或者数组的方法 dot() 实现。

[11]: np.dot(a, b)   # 或者 a.dot(b)
[11]: array([[ 30,  47],
             [ 79, 123]])

在实际的项目中，经常会遇到稀疏矩阵——详细内容请参阅拙作《机器学习数学基础》（电子工业出版社）。例如：

[12]: sm = np.array([[1, 0, 2, 0, 0, 0, 0, 0],
                     [0, 0, 0, 0, 0, 0, 0, 0],
                     [0, 0, 3, 0, 0, 0, 0, 0]])
      sm
[12]: array([[1, 0, 2, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 3, 0, 0, 0, 0, 0]])

用 sm 表示一个稀疏矩阵，由于行列数还不是很多，像上面那样写出来还能忍受。稀疏矩阵包含了很多零，虽然耗费了时间和精力，但它们事实上没有什么意义，还占用内存空间。为此可以对此类矩阵进行压缩，并直接生成压缩矩阵对象。

[13]: from scipy.sparse import csr_matrix
      row = np.array([0, 0, 2])
      col = np.array([0, 2, 2])
      data = np.array([1, 2, 3])
      csrm = csr_matrix((data, (row, col)), shape=(3, 8))
      csrm
[13]: <3x8 sparse matrix of type '<class 'numpy.longlong'>'
             with 3 stored elements in Compressed Sparse Row format>

现在得到的压缩矩阵 csrm 如果转化为数组，与前面创建的 sm 一样。

[14]: csrm.toarray()
[14]: array([[1, 0, 2, 0, 0, 0, 0, 0],
             [0, 0, 0, 0, 0, 0, 0, 0],
             [0, 0, 3, 0, 0, 0, 0, 0]], dtype=int64)

本书作者在《机器学习数学基础》(http://math.itdiffer.com)一书中，专门介绍了矩阵及其有关计算的 Python 实现方法，推荐有兴趣的读者深入学习参考。

12.4.4 解线性方程组

最一般的解线性方程组的方法是高斯消元法，在传统的数学教材中，还会列出其他巧妙的方法。这个工作如果教给 Python ，就是找到适合的函数——背后的事情教给了发明者，如果读者立志创新研究，可以自己成为发明者。

例如，求方程组的解：

\begin{cases}-x_1 + 3x_2 - 5x_3 &= -3 \\2x_1 -2x_2 + 4x_3 &= 8 \\ x_1 + 3x_2 &= 6\end{cases}

用矩阵的方式，可以将方程组表示为：

\begin{bmatrix}-1&3&-5\\2&-2&4\\1&3&0\end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}-3\\8\\6\end{bmatrix}

然后，编写如下代码：

[15]: A = np.mat("-1 3 -5; 2 -2 4; 1 3 0")    # 系数矩阵
      b = np.mat('-3 8 6').T                  # 常数项
      r = np.linalg.solve(A, b)               # 求解
      print(r)
[15]: [[ 4.5]
       [ 0.5]
       [-0. ]]

这里使用 np.linalg.solve() 函数求解线性方程组，从输出结果可知，其解为

x_1 = 4.5, x_2 = 0.5, x_3 =0

。

但是，如果遇到这样的方程组：

\begin{cases}x_1+3x_2-4x_3+2x_4&=0\\3x_1-x_2+2x_3-x_4&=0\\-2x_1+4x_2-x_3+3x_4&=0\\3x_1+9x_2-7x_3+6x_4&=0\end{cases}

还是使用前面的函数，对此方程组求解。

[17]: A = np.mat("1 3 -4 2;3 -1 2 -1;-2 4 -1 3;3 0 -7 6")
      b = np.mat("0 0 0 0").T
      r = np.linalg.solve(A, b)
      print(r)
[17]: [[ 0.]
       [ 0.]
       [-0.]
       [ 0.]]

当所有变量的解都是 0 ，原线性方程组成立，但这仅仅是一个特解。根据线性代数的知识可以判断，此方程组有无穷多个解（参阅《机器学习数学基础》2.4.2节），还能用程序计算吗？

[18]: from sympy import *
      from sympy.solvers.solveset import linsolve
      x1, x2, x3, x4 = symbols("x1 x2 x3 x4")
      linsolve([x1 + 3*x2 - 4*x3 + 2*x4, 
                3*x1 - x2 + 2*x3 - x4, 
                -2*x1 + 4*x2 - x3 + 3*x4, 
                3*x1 +9*x2 - 7*x3 + 6*x4], 
               (x1, x2, x3, x4))

输出：

这就是该线性方程组的通解，如果对应到未知量，即：

\begin{cases}x_1=\frac{x_4}{10}\\x_2=-\frac{7}{10}x_4\\x_3=0\\x_4=x_4\end{cases}

其中

x_4

是自由变量。

12.4.5 假设检验

假设检验是统计学中的重要内容，也广泛地应用在科学研究和工程项目中。下面选用《机器学习数学基础》6.5.1节的一个案例，读者从中感受 Python 在这方面的应用。

假设有人制造了一个骰子，他声称是均匀的，也就是假设分布律为：

H_0:P(X=i)=\frac{1}{6}, (i=1,2,3,4,5,6)

为了证明自己的判断，他做了

n=6\times10^{10}

次投掷试验，并将各个点数出现的次数记录下来（为了便于观察，出现次数用

10^{10}

加或减一个数表示。

点数	1	2	3	4	5	6
次数

接下来利用

\chi^2

检验法，通过如下程序对此人的结论——假设——进行检验。

[19]: from scipy.stats import chisquare
      chisquare([1e10-1e6, 1e10+1.5e6, 
                 1e10-2e6, 1e10+4e6, 1e10-3e6, 1e10+0.5e6])
[19]: Power_divergenceResult(statistic=3250.0, pvalue=0.0)

输出结果显示，检验统计量的值

\chi^2 =3250.0

。根据如下公式（对此公式详细讲解，请参阅《机器学习数学基础》6.5.1节）：

\chi^2 = \sum_{i=1}^k\frac{(np_i-f_i)^2}{np_i}

此处

𝑘=6

，则在

\alpha = 0.05

的显著水平下，

\chi^2_{0.05}(6-1)

的值为：

[20]: from scipy.stats import chi2 
      chi2.isf(0.05, (6-1))
[20]: 11.070497693516355

由输出值可知，在显著水平

\alpha = 0.05

下，

\chi^2\gt\chi^2_{0.05}(6-1)=11.07

，故试验数据不支持“骰子均匀”这个假设。

还可以根据 p 值检验此假设。

[21]: p_value = 1 - chi2.cdf(3250.0, (6-1)) 
      print(p_value)
[21]: 0.0

得到的 p 值结果说明拒绝原假设犯错误的概率是 0.0% 。

★自学建议 科学计算不仅仅是科研和工程项目的必备技能，也是后面所列举的数据分析、机器学习的基础。所以，有意在数据科学方向发展的读者务必要学习科学计算。首先，要具备相关的数学基础，其次要掌握相关计算的第三方包的引用，最后要将前两者应用到实际问题之中。也正是基于这些思考，我出版了《机器学习数学基础》（网址：http://math.itdiffer.com），在这本书中，不强调传统数学教材中的“纸笔计算”，重点是在理解有关数学原理之后，用程序工具完成计算，并以贴近真实的问题为案例。有兴趣的读者不妨参阅。”

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2022-08-17，如有侵权请联系 cloudcommunity@tencent.com 删除

python