如何在python中做更高精度的矩阵指数？答案

【问题标题】：How to do higher precision matrix exponential in python?如何在python中做更高精度的矩阵指数？
【发布时间】：2015-07-17 23:00:35
【问题描述】：

是否可以在 Python 中做更高精度的矩阵指数？我的意思是获得比双浮点数更高的精度。

我有以下测试代码：

import sympy
from sympy import N
import random

n = 100
#A = sympy.Matrix([[random.random(),random.random()],
#                   [random.random(),random.random()]])
A = sympy.Matrix([[1,2],[3,4]])
dlt = 1000
e1 = A.exp()
e1 = N(e1, n)
ee2 = (A/dlt).exp()
ee2 = N(ee2, n)
e2 = sympy.eye(2)
for i in range(dlt):
    e2 = e2*ee2
print(N(max(e1-e2)))

理论上，最终结果应该为零。使用 scipy，误差约为 1e-14。

通过sympy，如果矩阵像[[1,2],[3,4]]，前面代码的输出大约是1e-98。然而，对于随机矩阵，误差在 1e-14 左右。随机矩阵是否有可能得到 1e-100 这样的结果？

速度不是问题。

【问题讨论】：

标签： python matrix

【解决方案1】：

一旦你使用N，你就进入了浮点运算的领域，因此你永远不能假设你会达到绝对零。正如here 和许多其他地方所讨论的，所有浮点运算都是这种情况。 The only reliable solution is to include a suitably chosen eps variable and a function to check.

所以不要检查result == 0，而是定义isZero = lambda val: abs(val) < eps并检查isZero(result)。

这是浮点运算中的普遍问题。原则上，使用sympy，您可以找到实零，因为它是代数库，而不是浮点数学库。但是，在您给出的示例中，不使用N（这是切换到浮点运算的原因），会使计算变得非常缓慢。

【讨论】：

函数 N 不返回常规浮点数。例如， a=Symbols('a');N((1/a).subs(a, 3),100) 将返回一个包含 100 个三的数字。
但是，它仍然是一个有限精度浮点数，所有常规浮点数的限制都适用。

【解决方案2】：

我在尝试 mpmath 时犯了一个错误。我再次尝试了 mpmath，它是解决这个问题的完美解决方案。

【讨论】：

我会提醒你不要依赖这个。它有时可能会起作用，但只要您处于浮动世界中，即使它们在数学上是等价的，您也不能可靠地在不同操作之间实现相等。
谢谢。我知道这一点。我想我以错误的方式描述了我的问题。我并不期待完全平等。但是双浮点的精度对我来说还不够。我需要更高的精度，比如 10 倍的双浮点。这正是 mpmath 设计的目的。