点乘 B^t.D.B 不能返回对称数组

Question

我想做一个表达式的点积，它应该是对称的。你知道吗

事实证明，事实并非如此

B是一个4D数组，我必须将它的最后两个维度转换成B^t

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

与numpy.transpose相关联的numpy.dot产品，以及作为第二种选择numpy.einsum（这个想法来自于这个主题：Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix）已经被使用，并且问题仍然存在。你知道吗

在计算结束时，得到了乘积B^tDB，当通过减去它的转置B^tDB来验证它是否真的是对称的时，仍然有一个余数。你知道吗

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

问题是：如何消除这些残留物？你知道吗

Answer 1

您需要使用任意精度浮点数学。下面是如何结合numpy和^{} package来定义任意精度的矩阵乘法（即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^tDB-（B^tDB）^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)

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