Точечный продукт B^tDB не возвращает симметричный массив

Question

Я пытаюсь сделать точечное произведение выражения, и оно должно было быть симметричным.

Оказывается, это просто не так.

B - это четырехмерный массив, который я должен перенести в два последних измерения, чтобы он стал B^ t.

D - это двумерный массив. (Это выражение матрицы жесткости, известное программистам метода конечных элементов)

numpy.dotпродукт, связанный с numpy.transpose и как вторая альтернатива numpy.einsum (идея пришла из этой темы: умножение матриц Numpy U*B*UT Результаты в несимметричной матрице) уже были использованы, и проблема остается.

К концу вычислений получается произведение B^ tDB, и когда проверяется, действительно ли оно симметрично, вычитая его транспонированную B^ tDB, остается остаток.

Продукт Точки или Суммирование Эйнштейна используются только по интересующим измерениям (последние).

Вопрос в том, как эти остатки можно устранить?

Answer 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^ 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)

Берегись, хотя. Расчеты с использованием математических вычислений произвольной точности выполняются намного медленнее, чем при использовании стандартных чисел с плавающей запятой.