英文:
Is there a way for python to perform a matrix inversion at 500 decimal precision
问题
有一个对输出精度敏感的算法。特别是，论文要求执行一个500位小数精度的实数矩阵求逆运算。我想写一个Python脚本来检查结果。然而，在NumPy中使用的最大浮点数据类型是np.float128，而decimal包似乎没有矩阵求逆函数。
我找到了这篇帖子https://stackoverflow.com/questions/32685280/matrix-inverse-with-decimal-type-numpy，它介绍了将decimal作为NumPy函数中的对象。但我在Stack Exchange上看到了关于Python将对象传递为float(32)对象的评论，这是由于%运算符的原因。
Python是否有一种方法可以以500位小数精度执行矩阵求逆运算？
英文:
There's an algorithm sensitive to the precision of the output. Especially, the paper required a real matrix inverse to be performed at 500 decimals precision. I wanted to write a script to check the result with python. However, the largest float data type used in numpy was np.float128, and the deicmal package does not seem to have a matrix inverse function.
I found the post https://stackoverflow.com/questions/32685280/matrix-inverse-with-decimal-type-numpy which introduce the decimal as an object in the numpy's function. But I saw comments on stack exchange that python transfer the object as float(32) object because of the % operators.
Is there a way for python to perform a matrix inversion at 500 decimal precision?
答案1
得分: 1
根据https://mpmath.org/doc/current/matrices.html#advanced-methods，我导入了mpmath，并创建了一个(3,3)的矩阵：
A = mpmath.matrix([[1,3,4],[0,1,2],[3,0,1]])
并计算了它的逆矩阵：
A**-1
numpy 中的等效操作：
M = np.array([[1,3,4],[0,1,2],[3,0,1]])
np.linalg.inv(M)
以及比较的运行时间：
timeit np.linalg.inv(M)
timeit A**-1
在这个示例中，dps（显示精度）设置为30（默认值为15，类似于float64）：
mpmath.mp.dps
有时用户担心逆矩阵与数组的点乘不完全产生单位矩阵，但请注意接近零的项非常小：
Out[295] @ M
以及30位小数精度：
Out[291] * A
英文:
Following https://mpmath.org/doc/current/matrices.html#advanced-methods,
I imported mpmath, and created (3,3) matrix:
In [289]: A =mpmath.matrix([[1,3,4],[0,1,2],[3,0,1]])
In [290]: A
Out[290]:
matrix(
[['1.0', '3.0', '4.0'],
['0.0', '1.0', '2.0'],
['3.0', '0.0', '1.0']])
And calculated its inverse:
In [291]: A**-1
Out[291]:
matrix(
[['0.142857142857142857142857142857', '-0.428571428571428571428571428571', '0.285714285714285714285714285714'],
['0.857142857142857142857142857143', '-1.57142857142857142857142857143', '-0.285714285714285714285714285714'],
['-0.428571428571428571428571428571', '1.28571428571428571428571428571', '0.142857142857142857142857142857']])
The equivalent numpy operation:
In [293]: M = np.array([[1,3,4],[0,1,2],[3,0,1]])
In [294]: M
Out[294]:
array([[1, 3, 4],
[0, 1, 2],
[3, 0, 1]])
In [295]: np.linalg.inv(M)
Out[295]:
array([[ 0.14285714, -0.42857143, 0.28571429],
[ 0.85714286, -1.57142857, -0.28571429],
[-0.42857143, 1.28571429, 0.14285714]])
And the comparative times:
In [296]: timeit np.linalg.inv(M)
17.2 µs ± 31 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)
In [297]: timeit A**-1
1.19 ms ± 1.26 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
In this case I have the dps (display precision) set at 30 (the default is 15, similar to a float64
In [298]: mpmath.mp.dps
Out[298]: 30
SO posters often worry about the dot of inverse with array doesn't quite produce the identity matrix. But notice how small the near 0 terms are:
In [299]: Out[295]@M
Out[299]:
array([[ 1.00000000e+00, -5.55111512e-17, -5.55111512e-17],
[ 1.11022302e-16, 1.00000000e+00, 3.88578059e-16],
[ 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]])
And to 30 dps:
In [301]: Out[291]*A
Out[301]:
matrix(
[['1.0', '-7.22223729145213444895991728469e-35', '-4.81482486096808963263994485646e-35'],
['1.44444745829042688979198345694e-34', '1.0', '3.37037740267766274284796139952e-34'],
['-2.40741243048404481631997242823e-35', '9.62964972193617926527988971292e-35', '1.0']])
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论