Diagonalization of symmetric matrices in julia yields strange results

Question

I would like to understand what is going on with my program in julia. The issue is the following: I have a symmetric non-negative matrix which I diagonalize using

egvals, egvecs = eig(H_mat)

By a theorem, my matrix should have a maximum eigenvalue which is associated a non-negative eigenvector. The H_mat has a further trick, its first column and row has an entry filled with zeros.

Diagonalization yields a maximum positive eigenvalue E_max, in fact it is the last eigenvalue because julia arranges the eigenvalues in order up to the biggest, but my eigenvector associated with E_max does not have all its entries zero or positive (i.e. they have negative entries)

egvecs[:,end] # Some or several components ii, egvecs[ii,end]<0

This is the matrix from I don't get proper results, for instance:

[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
 0.0 0.0 2.0 0.0 0.0 1.414213562373095 0.0 0.0 0.0 0.0 0.0 1.414213562373095 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 3.0 0.0 0.0 1.7320508075688774 0.0 0.0 0.0 0.0 0.0 1.7320508075688774 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 1.414213562373095 0.0 0.0 2.0 0.0 1.414213562373095 0.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 1.7320508075688774 0.0 0.0 3.0 0.0 2.0 0.0 0.0 0.0 1.0 0.0 1.414213562373095 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1.414213562373095 0.0 2.0 0.0 0.0 0.0 0.0 0.0 1.414213562373095 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 3.0 1.7320508075688774 0.0 0.0 0.0 0.0 1.414213562373095 1.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.7320508075688774 3.0 0.0 0.0 0.0 0.0 0.0 1.7320508075688774 0.0 0.0 0.0 0.0
0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 1.414213562373095 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 1.0 0.0 0.0 1.414213562373095 0.0 0.0 0.0
0.0 0.0 0.0 1.7320508075688774 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 3.0 0.0 1.414213562373095 0.0 0.0 2.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.414213562373095 0.0 0.0 0.0 1.0 0.0 2.0 0.0 0.0 1.414213562373095 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 1.414213562373095 0.0 1.414213562373095 0.0 0.0 0.0 1.414213562373095 0.0 3.0 1.414213562373095 0.0 1.414213562373095 1.414213562373095 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 1.7320508075688774 0.0 0.0 0.0 0.0 1.414213562373095 3.0 0.0 0.0 2.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.414213562373095 0.0 1.414213562373095 0.0 0.0 2.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 2.0 0.0 1.414213562373095 0.0 0.0 3.0 1.0 1.7320508075688774
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.414213562373095 2.0 0.0 1.0 3.0 1.7320508075688774
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.7320508075688774 1.7320508075688774 3.0]

(Should be a 20X20 matrix)

Tried a random H_mat with the properties you state and the max eigenvalue is associated with an eigenvector with all non-negative values. Can you provide your H_mat? — rickhg12hs
Very curious. Rather than H_mat, have a look at the eig results for s=svd(H_mat);eig(s[1] * diagm(s[2]) * s[3]'). Other than small precision issues, the property you state seems to hold. — rickhg12hs
Also have a look at eig(H_mat ^ 2). Property seems to hold here as well. — rickhg12hs
Hmmm, seems my knowledge of Linear Algebra is fading and/or has some holes in it. Eigenvectors are defined up to a phase, i.e. if $Av=\lambda v then e^{i\theta}v$ is also an eigenvector, and specifically so is $-v$.. So feel free to rotate/negate. Julia's negated eig-vectors are legit and fir your theory. [Kinda wished I'd have researched/relearned this before I chimed in.] 8-) — rickhg12hs

Andreas Noack Andreas Noack · Accepted Answer · 2015-04-10T14:52:31

Eigenvectors are only determined up to proportionality because the requirement is that they solve (A - lambda*I)v = 0 and if v1 solves the equation, then also v2 = -v1 does. In the symmetric case, it is standard to normalize all the vectors to one, but that still leaves the sign undetermined. Therefore the theorem you are using must say something like that the vector corresponding to the largest values can be chosen such that it has non-negative elements. Indeed, I get that the eigenvector corresponding to the largest eigenvalue of your matrix has non-positive elements. julia> eig(A) |> t -> all(t[2][:,indmax(t[1])] .<= 0) true

Diagonalization of symmetric matrices in julia yields strange results

1 Answers