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

Question

有没有可能在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。

根据渐近性，如果矩阵类似于[1,2,3,4]，则先前代码的输出约为1e-98。然而，对于随机矩阵，误差在1e-14左右。对于随机矩阵，有可能得到像1e-100这样的结果吗？

速度不是问题。

Answer 1

一旦您使用了N，您就进入了浮点运算的领域，因此您永远不能假设您将达到绝对零。这是所有浮点算术的情况，正如here和许多其他地方所讨论的那样。唯一可靠的解决方案是包含一个适当选择的eps变量和一个要检查的函数。

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

这是浮点运算中的一个普遍问题。原则上，使用sympy可以找到真正的零，因为它是一个代数库，而不是浮点数学库。但是，在您给出的示例中，如果不使用N (它将切换到浮点算法)，则会使计算极其缓慢。

Answer 2

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