1 * X' * X is not symmetric

[e3c6]

I find the following snippet unexpected.

srand(1);
X = rand(1000, 1000) - 0.5;
(X' *X) == (X' *X)' #true
(1 * X' *X) == (1 * X' *X)' #false

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.

[e3c6]

After some discussion in the slack chat (see Slack), it seems that the culprit code is here:

https://github.com/JuliaLang/julia/blob/2c5aec90c9f9b58158eb61a516f6c2565c25eb9a/stdlib/LinearAlgebra/src/matmul.jl#L216

1 * A creates a copy of A, which the fails the identity test A === B, which results in a call to gemm_wrapper!…

But I don’t understand why that affects the result and breaks the symmetry.

[Tamas_Papp]

I think that 1 * X' * X is evaluated as (1 * X') * X, and 1 * X' cannot be recognized as coming from the same matrix X.

You should use 1 * (X' * X), or even better, Symmetric(X' * X) * c for some generic real c (I am assuming the 1 is just for the MWE).

[kristoffer.carlsson]

@e3c6

But I don’t understand why that affects the result and breaks the symmetry.

Different algorithms computing the same thing will give different results in floating point arithmetic. Don’t use == to compare floating point values.

[e3c6]

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.

[e3c6]

What exactly is the difference between syrk_wrapper! and gemm_wrapper!?

[Nosferican]

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

[e3c6]

Is there a reason for === instead of == in the snippet of code of mul!?

[Tamas_Papp]

@e3c6

I want to understand where it comes from to understand the kind of floating point errors I might be committing here.

A good starting point would be

PSA: floating-point arithmetic

Sometimes people are surprised by the results of floating-point calculations such as julia> 5/6 0.8333333333333334 # shouldn't the last digit be 3? julia> 2.6 - 0.7 - 1.9 2.220446049250313e-16 # shouldn't the answer be 0? These are not bugs in Julia. They’re consequences of the IEEE-standard 64-bit binary representation of floating-point numbers that is burned into computer hardware, which Julia and many other languages use by default. Brief explanation You can t…

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.

[StefanKarpinski]

@e3c6

Is there a reason for === instead of == in the snippet of code of mul! ?

Yes, A === B checks if A and B are the same matrix rather than checking if they are equal.

[e3c6]

Wouldn’t the same thing work with == instead of ===? Why do you need them to be the same matrix, instead of just having the same numerical values?

[StefanKarpinski]

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

[e3c6]

I see. Thanks!