深度对比Python(Numpy，Scipy)与Matlab的数值精度

Matlab一度被认为是最专业的数值计算工具之一，相信许多同学都或多或少用过这个工具。相比而言，Python作为一种胶水式的语言，其设计之初就不是为科学计算服务的。之前也看到许多人在吐槽说用Python去复现一些计算过程时经常失败，因此（包括本人）也怀疑过是Python本身数值精度不够导致的。那么Python的精度究竟如何，本文就来一探究竟。为了方便，我们就用线性方程组的求解来对比这一事实。

1、实验设计

• 基本思路：

本文考虑的线性方程组：
$$Ax=b$$

• 关键问题一：保证两种语言处理的数据完全一致

这里第二步其实很重要，由于Matlab默认数据类型为double，而Python的默认类型为float，如果直接粘贴或者导入就有可能在Python中被截断，从而导致在一开始使用的数据就存在误差，使得实验结果不可信。

• 关键问题二：选取合适的指标来进行对比

第四步主要是计算两种程序跑完后均方根误差：
$$\Vert Ax-b\Vert$$
这种验证方式比较简单。之所以这么做是因为数值计算过程中不一定能够得到100%的精度，这个问题的主要原因是在数值计算时数字不会一直保持为真实值。尤其在除法时由于数值长度的限制会被强行截断一部分，因此用这个指标来验证方法是否正确是很有必要的。当然，根据它的定义不难看出，这个值最小精度就越高。

2、数据生成和处理

首先由Matlab随机生成一个矩阵和一个向量：

A = rand(5,5)
b = rand(5,1)

然后在Matlab中将数据粘至txt文件，再整理成Python的格式粘生成Numpy数组，同时也将数据重新粘回Matlab代码。这里原始数据是这个样子：

A的值

0.814723686393179	0.0975404049994095	0.157613081677548	0.141886338627215	0.655740699156587
0.905791937075619	0.278498218867048	0.970592781760616	0.421761282626275	0.0357116785741896
0.126986816293506	0.546881519204984	0.957166948242946	0.915735525189067	0.849129305868777
0.913375856139019	0.957506835434298	0.485375648722841	0.792207329559554	0.933993247757551
0.632359246225410	0.964888535199277	0.800280468888800	0.959492426392903	0.678735154857774

b的值

0.757700000000000
0.743100000000000
0.392200000000000
0.655500000000000
0.171200000000000

当然也可以采用round函数一类的方式进行处理。

注意 这种操作虽然看似简单并且有点笨拙，但能很好地保证二者在初始化时得到的结果是一样的。

Matlab的操作相对简单，直接将txt中的内容首尾加上[]，再粘回代码即可。对于Python要稍加处理，最终的代码是：

import numpy as np
A = np.array([0.814723686393179,0.905791937075619,0.126986816293506,0.913375856139019,0.632359246225410,0.0975404049994095,0.278498218867048,0.546881519204984,0.957506835434298,0.964888535199277,0.157613081677548,0.970592781760616,0.957166948242946,0.485375648722841,0.800280468888800,0.141886338627215,0.421761282626275,0.915735525189067,0.792207329559554,0.959492426392903,0.655740699156587,0.0357116785741896,0.849129305868777,0.933993247757551,0.678735154857774]).reshape([5,5]).T

b = np.array([0.757700000000000,0.743100000000000,0.392200000000000,0.655500000000000,0.171200000000000])

注意这里Python的代码最后要进行一下转置，否则与Matlab的数据就就不一致。

整理完成看看数据和txt中的形状就完全一致了：

A

array([[0.81472369, 0.0975404 , 0.15761308, 0.14188634, 0.6557407 ],
       [0.90579194, 0.27849822, 0.97059278, 0.42176128, 0.03571168],
       [0.12698682, 0.54688152, 0.95716695, 0.91573553, 0.84912931],
       [0.91337586, 0.95750684, 0.48537565, 0.79220733, 0.93399325],
       [0.63235925, 0.96488854, 0.80028047, 0.95949243, 0.67873515]])

b

array([0.7577, 0.7431, 0.3922, 0.6555, 0.1712])

这里b变成了一维数组但不影响计算。

3、求解方程组

• Matlab部分

用Matlab求解线性方程组十分简单，直接一个反除号就可以搞定：

>> x = A \ b

x =

   -0.8409
    3.7701
    3.8572
   -8.0050
    2.4445

注意这里x跑出来的值只是显示值，实际值的位数是远大于这个值的：

-0.840928749143009	3.77006972003583	3.85716025378275	-8.00495386938441	2.44448015270832

• Python部分

用Python求解线性方程组的方式其实也很多，最简单的办法就是直接用scipy库的linalg模块中的solve方法：

from scipy.linalg import solve

x = solve(A,b)

x

array([-0.84092875,  3.77006972,  3.85716025, -8.00495387,  2.44448015])

为了更准确地对比数值的精度，我们将x的值用np.savetxt函数保存至文件再看具体数字：

-8.409287491430091910e-01
3.770069720035832184e+00
3.857160253782746739e+00
-8.004953869384410226e+00
2.444480152708315313e+00

仔细观察发现，Python输出的数值位数要稍多一些，这是因为我们用Matlab粘贴的时候它默认显示的位数不太一样。

注：这也是为什么我们一定要在最开始进行数据统一规范的原因。Matlab在存储变量时数据是保存了完整的位数的，但即便从变量查看窗口粘出来的数值也不一定是完整的数字。

4、精度验证

• Matlab部分

error = norm(A*x-b)

error =

   1.5801e-15

% 查看15位
% 1.580118849929844e-15

• Python部分

np.sqrt(np.sum((A.dot(x)-b)**2))

1.580118849929844e-15

注意，由于两个数的量级都是1e-15，因此这里只需要看它们各自前16位即可。仔细观察不难发现，两个程序得到的结果完全一致。

这里我们把有效数字从第一位到最后一位排列直接放出来对比

# python
1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625

% matlab
1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625

肉眼观察信不过，全部放入Python字符串试试对不对：

#第一行数字(python结果)：
s1 = '1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625'

#第二行数字(matlab结果)：
s2 = '1580118849929843973242195223453751011511480618378089335607228349545039236545562744140625'

判断一下是否相等：

s1 == s2

True

搞定！

5、结论

Python(Numpy, Scipy)的数值精度与Matlab是完全一致的。

其本质的原因来自于二者在矩阵计算都是使用的LAPACK, 详情可见：
Matlab和Python(Numpy,Scipy)与Lapack的关系