# 1 * X' * X is not symmetric

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.

After some discussion in the slack chat (see https://julialang.slack.com/archives/C67910KEH/p1528205759000199), it seems that the culprit code is here:

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

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

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

1 Like

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.

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

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

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

A good starting point would be

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

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

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?

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

3 Likes

I see. Thanks!