numpy：通过逆向计算维护decimal类

Question

我试图找到由一个特定类（decimal.Decimal）组成的矩阵的逆矩阵，并在整个过程中将值保持为十进制对象（在整个计算过程中保持精确性）。在

我的问题是numpy.linalg.inverse总是以浮点形式返回矩阵值。我已经找到了一个解决办法，在计算逆数后，将类型从float改为Decimal对象，但我更希望始终保持原始矩阵的类（我担心当数字转换为float时，我可能会失去准确性）？在

所以我想我有几个问题：（1）当矩阵的值被转换成float类型时，我是否会失去准确性/精确性（我正在处理一个11乘11的矩阵）；如果是这样，（2）在整个计算过程中使用numpy保持这些值为decimal.Decimal；如果不是，那么（3）对于这种类型，是否还有其他模块/方法需要考虑计算的结果？在

下面是我的代码的示例：

import numpy as np
from decimal import Decimal as D

a = np.array( [ [D('1'), D('2'), D('3')],
                [D('4'), D('5'), D('6')],
                [D('7'), D('8'), D('9')] ] )

print type(a[0,0])
# <class 'decimal.Decimal'>

inverse_a = np.linalg.inv(a)

print type(inverse_a[0,0])
# <type 'numpy.float64'>

inverse_a_Decimal_flat = [D(str(i)) for i in inverse_a.flat]  # change all values to decimal.Decimal
inverse_a_Decimal = np.reshape(inverse_a_Decimal_flat, [3, 3])  # reshape to 3x3

print type(inverse_a_Decimal[0,0]), d.shape
# <class 'decimal.Decimal'> (3,3)

Answer 1

您的a是一个对象数组：

In [255]: a
Out[255]: 
array([[Decimal('1'), Decimal('2'), Decimal('3')],
       [Decimal('4'), Decimal('5'), Decimal('6')],
       [Decimal('7'), Decimal('8'), Decimal('9')]], dtype=object)

在我最新更新的numpy中，我在尝试逆操作时遇到了一个错误。在

^{pr2}$

inv可能正在使用编译代码（_umath_linalg.inv），它将被编码为与标准的c浮点（float或double）一起工作。在

您询问了longdouble：

In [266]: a1=np.arange(1,10,dtype=np.longdouble).reshape(3,3)
In [267]: a1
Out[267]: 
array([[ 1.0,  2.0,  3.0],
       [ 4.0,  5.0,  6.0],
       [ 7.0,  8.0,  9.0]], dtype=float96)
In [268]: np.linalg.inv(a1)
...
TypeError: array type float96 is unsupported in linalg

同样，这是一个如何编译linalg的问题。在

In [269]: a1=np.arange(1,10,dtype=np.double).reshape(3,3)
In [270]: np.linalg.inv(a1)
Out[270]: 
array([[  3.15221191e+15,  -6.30442381e+15,   3.15221191e+15],
       [ -6.30442381e+15,   1.26088476e+16,  -6.30442381e+15],
       [  3.15221191e+15,  -6.30442381e+15,   3.15221191e+15]])

我试图dot这个，但没有得到预期的np.eye数组。但这是一个几乎是奇异的数组（行列式是0，或者差不多）。在

In [274]: np.linalg.det(a1)
Out[274]: -9.5171266700777698e-16

所以这是个坏例子。在

将值洗牌到数组不再是单数

In [288]: a1=np.array([[3,1,2],[4,6,5],[8,9,7]],np.float)
In [289]: a1.dot(np.linalg.inv(a1))
Out[289]: 
array([[  1.00000000e+00,   0.00000000e+00,   2.22044605e-16],
       [ -4.44089210e-16,   1.00000000e+00,   4.44089210e-16],
       [  0.00000000e+00,   0.00000000e+00,   1.00000000e+00]])

如果a1是np.int，我也会得到同样的结果，因为np.int被上推到{}。我可以用np.float32进行测试，以查看较低精度的效果。如前所述，我不能去longfloat。在

inv通常计算为I/a1线性方程解：

In [297]: np.linalg.solve(a1,np.eye(3))
Out[297]: 
array([[ 0.14285714, -0.52380952,  0.33333333],
       [-0.57142857, -0.23809524,  0.33333333],
       [ 0.57142857,  0.9047619 , -0.66666667]])
In [298]: np.linalg.inv(a1)
Out[298]: 
array([[ 0.14285714, -0.52380952,  0.33333333],
       [-0.57142857, -0.23809524,  0.33333333],
       [ 0.57142857,  0.9047619 , -0.66666667]])

Answer 2

使用^{}是否可以接受？在

如果我们用它：

>>> import sympy as sy
>>> a = sy.Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
>>> a.inv()
ValueError: Matrix det == 0; not invertible.

我们意识到，颠倒这个矩阵是没有意义的！在

再次尝试用10来交换9，得到

^{pr2}$

这是期望的无限精度结果。在