1 * X' * X is not symmetric

linearalgebra

#1

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.


#2

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.


#3

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


#4

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


#5

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.


#6

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


#7

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


#8

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


#9

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.


#10

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


#12

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?


#13

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


#14

I see. Thanks!