Dot product B^t.D.B doesn't return a symmetric array

Question

I'm trying to make a dot product of an expression and it was supposed to be symmetric.

It turns out that it just isn't.

B is a 4D array which I must transpose its last two dimensions to become B^t.

D is a 2D array. (It's an expression of the Stiffness Matrix known to the Finite Element Method programmers)

The numpy.dotproduct associated with numpy.transpose and as a second alternative numpy.einsum (the idea came from this topic: Numpy Matrix Multiplication U*B*U.T Results in Non-symmetric Matrix) have already been used and the problem persists.

By the end of the calculations the product B^tDB is obtained and when it's verified if it really is symmetric by subtracting its transpose B^tDB, there is still a residue.

The Dot product or the Einstein Summation are used only over the dimensions of interest (last ones).

The question is: How can these residues be eliminated?

Answer 1

You need to use arbitrary precision floating point math. Here's how you can combine numpy and the mpmath package to define an arbitrary precision version of matrix multiplication (ie the np.dot method):

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

As an example, if you then set up B, B^t, and D matrices as so:

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)

then B^tDB - (B^tDB)^t will always have a non-zero value if you calculate it using the standard matrix multiplication method from numpy:

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

but if you use the arbitrary precision version given above it won't have a residue:

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

Watch out though. Calculations using arbitrary precision math run much slower than those done using standard floating point numbers.