I have a symmetric matrix B, and I am applying a quadratic operation to the matrix with a second matrix f, of the form C = f’ B f. The result of this operation should also be a symmetric matrix. However, due to numeric precision, my implementation gives a result that is not symmetric, but deviates from symmetry very slightly. This is an issue for later processes that require that the resulting matrix is symmetric (actually, symmetric positive definite).
I wonder if there’s an elegant way to have the result stay symmetric, i.e. not having to copy-paste the lower-triangular matrix into the upper triangular part of the matrix.
Code that replicates this problem:
# Generate a SPDM
import torch
n = 400
d = 5
A = torch.randn(n, n)
B = torch.matmul(A, A.transpose(0,1)) * 0.01 # Original symmetric matrix
f = torch.randn(d, n) * 0.1 # Random filters matrix
C = torch.einsum('kd,db,mb->km', f, B, f) # Compute the quadratic form
print(f'Symmetric residue of original matrix: \n {B - B.transpose(0,1)}')
print(f'Symmetric residue of quadratic formula result: \n {C - C.transpose(0,1)}')
Symmetric residue of original matrix:
tensor([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])
Symmetric residue of quadratic formula result:
tensor([[ 0.0000e+00, 3.5763e-07, -1.1921e-07, 2.3842e-07, -4.7684e-07],
[-3.5763e-07, 0.0000e+00, -2.9802e-07, -4.7684e-07, -1.1921e-07],
[ 1.1921e-07, 2.9802e-07, 0.0000e+00, -1.7881e-07, 7.1526e-07],
[-2.3842e-07, 4.7684e-07, 1.7881e-07, 0.0000e+00, -4.7684e-07],
[ 4.7684e-07, 1.1921e-07, -7.1526e-07, 4.7684e-07, 0.0000e+00]])