eig function returning complex eigen vector values for a symmetric matrix, why?

I have a matrix of size 50 * 4, I stored it as H. To find the eigen values of H I took the gramian of H by the following code

Ht = H'
X= H*Ht
[V,D] = eig(X)

by doing this I am getting the eigen vector V as complex. But when I use [V,D] = eig(H*H') I am getting all real values. May I know why there is such a difference?

Answer 1

[Christine Tobler]

EIG checks whether the input matrix is symmetric, in which case it chooses an algorithm that always returns real matrices. It checks for exact equality, so a matrix that is non-symmetric on the level of round-off error is seen as not symmetric.

The reason that H*Ht and H*H' are not the same is more subtle: If matrix multiplication and transposition is combined, MATLAB does both operations at once. This is faster, and also means it can detect cases where both matrices are the same (X*X' or X'*X). In this case, a specific algorithm guarantees that the output matrix is exactly symmetric.

However, this only works if MATLAB recognizes that the two arrays are the same:

% Exact same matrix
>> H = randn(50, 4); A = H*H'; norm(A - A')
ans =
     0

% Direct copy of H
>> H = randn(50, 4); H2 = H; A = H*H2'; norm(A - A')
ans =
     0

% Modified copy of H - H and H2 are not recognized as the same
>> H = randn(50, 4); H2 = 1*H; A = H*H2'; norm(A - A')
ans =
   8.3299e-16

% Transpose of H: H and H2 not recognized as the same, and also
% the transpose is computed in a separate step from the matrix multiply.
>> H = randn(50, 4); H2 = H'; A = H*H2; norm(A - A')
ans =
   8.3751e-16

[Walter Roberson]

Oh, subtle!!

[Jaya Harsha Thippana]

Yes, I too noticed that the difference is zero but i got complex values for eigen vectors when i take Ht = H' here I am attaching a text file which contains my dataset and the commands I executed. Please take a look at that. Thank you.

[Christine Tobler]

The difference is that

issymmetric(H*H') 

is true, but

issymmetric(H*Ht)

is false. EIGS checks issymmetric as a first step to choose between algorithms, and (typically) returns complex values if the input is not exactly symmetric.

In general, for a nearly-symmetric matrix A, call eig((A+A')/2) instead of eig(A) to get real eigenvalues.

Is this question related to some code where you need to do the transpose and the matrix multiplication separately? If yes, note that this will not be faster: A = H*H' does not compute X' at all, it just applies the matrix multiplication and transposition in one pass. So the call Ht = H; A = H*Ht is probably slower than A = H*H', even if you are using Ht in several places.