Multiplying by 1 changes the symmetry, probably due to floating-point arithmetic issues. But I fail to understand where the failure happens.

Let A = X' * X. Naively I expect that:

A[i,j] == sum(X[k,i] * X[k,j] for k = 1 : N)
A[j,i] == sum(X[k,j] * X[k,i] for k = 1 : N)

Since multiplication is commutative (this holds in floating-point arithmetic as well), it follows that A[i,j] == A[j,i], so A is symmetric. I don’t see how multiplying by 1 changes this.

Thanks. However I am not trying to fix the issue. I want to understand where it comes from to understand the kind of floating point errors I might be committing here.

Those are different BLAS routines. If not mistaken syrk (symmetric rank-k update) takes advantage of crossproducts while gemm (General Matrix Multiply) is general.

TL;DR: order matters, if you know (from theory) that two things are equal and want to use it like that, don’t calculate it twice, or discard one of them. Hence Symmetric.

Because checking === is very fast: it only requires comparing two pointers. Checking == is expensive: in general, it requires looking at all elements of both matrices (which may be very large).