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