1

我正在尝试制作表达式的点积，它应该是对称的。

事实证明，事实并非如此。

B是一个 4D 数组，我必须将其最后两个维度转置为B ^t。

D是一个二维数组。（它是有限元方法程序员已知的刚度矩阵的表达式）

与第二个选择相关的numpy.dot产品（这个想法来自这个主题：Numpy Matrix Multiplication U*B*UT Results in Non-symmetric Matrix）已经被使用并且问题仍然存在。numpy.transposenumpy.einsum

在计算结束时，获得了乘积B ^t DB，当通过减去它的转置B ^t DB来验证它是否真的是对称的时，仍然有一个残差。

点积或爱因斯坦求和仅用于感兴趣的维度（最后一个）。

问题是：如何消除这些残留物？

4

1 回答 1

1

您需要使用任意精度的浮点数学。以下是您如何组合numpy和mpmath包来定义矩阵乘法的任意精度版本（即np.dot方法）：

from mpmath import mp, mpf
import numpy as np

# stands for "decimal places". Larger values 
# mean higher precision, but slower computation
mp.dps = 75

def tompf(arr):
    """Convert any numpy array to one of arbitrary precision mpmath.mpf floats
    """
    if arr.size and not isinstance(arr.flat[0], mpf):
        return np.array([mpf(x) for x in arr.flat]).reshape(*arr.shape)
    else:
        return arr

def dotmpf(arr0, arr1):
    """An arbitrary precision version of np.dot
    """
    return tompf(arr0).dot(tompf(arr1))

例如，如果您随后设置B、B ^t 和D矩阵，如下所示：

bshape = (8,8,8,8)
dshape = (8,8)

B = np.random.rand(*bshape)
BT = np.swapaxes(B, -2, -1)

d = np.random.rand(*dshape)
D = d.dot(d.T)

那么B ^t DB - ( B ^t DB )^t 如果您使用标准矩阵乘法方法从以下位置计算它，它将始终具有非零值numpy：

M = np.dot(np.dot(B, D), BT)
np.sum(M - M.T)

但如果您使用上面给出的任意精度版本，它不会有残留：

M = dotmpf(dotmpf(B, D), BT)
np.sum(M - M.T)

不过要小心。使用任意精度数学的计算比使用标准浮点数的计算运行得慢得多。

于 2018-04-18T02:35:30.660 回答