2

I have two matrixes with values differs only slightly

CloudGet["https://www.wolframcloud.com/objects/6cfbdce0-d7d4-4999-abd1-6a4c2db3bccf"];
CloudGet["https://www.wolframcloud.com/objects/299385d7-3f6f-4571-877e-e05313aa31f1"];

Hls1 - Hls2 // Flatten // Abs // Max
(* 5.46584*10^-10 *)

but EigenSystem gives large difference.

Eigensystem[Hls1] - Eigensystem[Hls2] // Flatten // Abs // Max
(* 1.99918 *)

While Matlab seems to give reasonable close results (MATLink package by Szabolcs is here).

Needs["MATLink`"]
OpenMATLAB[]

MSet["Hls1", Hls1]
MSet["Hls2", Hls2]

MEvaluate["[V1,D1] = eig(Hls1)"];
MEvaluate["[V2,D2] = eig(Hls2)"];

MEvaluate["max(max(abs(D1-D2)))"]
(* ">> 
ans =

   1.0851e-09

" *)

MEvaluate["max(max(abs(V1-V2)))"]
(* ">> 
ans =

   8.0894e-06

" *)

It appears that both values of the eigensystem for the two matrixes are correct, but there is a phase jump in the eigenvectors. So why is there this phase jump? Does Mathematica use two different methods on the two matrixes?

Note that we can always manually correct the phase by setting a convention, but sometimes this is tedious and time consuming, if one only requires the consistence of the phase.

Community
  • 1
xslittlegrass
  • 27,549
  • 9
  • 97
  • 186
  • You might try sorting the eigenvectors by length, then taking their absolute values (lengths) before comparison. That way, differences in signs (which are sensitive to both data variations and algorithms) will be ignored. – David G. Stork Feb 24 '15 at 17:09

1 Answers1

8

You are apparently looking for a way to reliably compare two numerical matrices by using their eigensystems. This can always be done for normal matrices by using the eigenvectors to construct their spectral decomposition. To do that, you shape the eigenvectors into the equivalent system of projectors. Then you can compare the projectors instead. The good thing about projectors is that they don't even contain the irrelevant phase factor you're worried about in the first place:

{eval1, evec1} = Eigensystem[Hls1];
{eval2, evec2} = Eigensystem[Hls2];

esProject1 = {eval1, Map[TensorProduct[Conjugate[#], #] &, evec1]};
esProject2 = {eval2, Map[TensorProduct[Conjugate[#], #] &, evec2]};

esProject1 - esProject2 // Flatten // Abs // Max

(* ==> 6.40185*10^-6 *)

The vectors returned by Eigensystem are turned into matrices by using TensorProduct. This is where the phase gets annihilated, and the result is a projector onto the subspace of that eigenvalue. To see that the spectral representation in this case indeed reproduces the original matrix, you would do something like this:

(Total[Times @@ esProject1] - Hls1) // Flatten // Abs // Max

(* ==> 3.44264*10^-7 *)

This approach works with your matrices because they are both Hermitian. I am assuming this is OK because you didn't say otherwise... Also, the normalization of the eigenvectors that I assume here is only true if the matrices are numerical.

Of course you can also use the construction I wrote in the last line directly to compare the matrices: i.e., reconstruct the matrices using Total[Times @@ esProject1] and then compare the results. That's probably the best way if you have degeneracies.

Jens
  • 97,245
  • 7
  • 213
  • 499
  • Thanks for the answer, indeed a Hermit matrix can always be decomposed into its spectrum form $H=\sum_n \lambda_n \left|n\right\rangle \left\langle n \right|$, and comparing the spectrum is a good way to compare the original matrix. But I am actually interested in calculating the eigenvectors of a parametric matrix which changes slightly in a range and found that this behavior of Eigensystem. I'm wondering what costs this seems random behavior? Does Eigensystem use LAPACK under the hood? – xslittlegrass Feb 24 '15 at 19:47
  • That's probably one of the mysteries that only a Wolfram employee can answer definitively... it's worth knowing that such discrepancies exist, but since there is no reason to not expect them, I would simply not depend on any properties of the eigenvectors that are mathematically not required to have any specific value. – Jens Feb 24 '15 at 20:32
  • Congrats on 50K :-) – Mr.Wizard Feb 26 '15 at 00:31
  • @Mr.Wizard Thanks - the view must be nice up there where you are at 130K... maybe I'll find out one day. – Jens Feb 26 '15 at 01:49
  • 1
    @xslittlegrass, yes, LAPACK is used under the hood. Look at the result of ??LinearAlgebra`LAPACK`*, for instance. – J. M.'s missing motivation Jul 09 '15 at 00:27
  • @J.M. Thanks! That's really helpful. Is there usage information I can found about those functions, or examples using them directly? Sometimes I have knowledge about my linear algebra problem, and a special transparent LAPACK subroutine would serve better than the general LinearSolve, for which I don't have control with. – xslittlegrass Jul 09 '15 at 04:15
  • @xslittlegrass, unfortunately, I haven't found out how exactly to use them within Mathematica; the extent of my knowledge is that I know now which algorithms they're using, per the LAPACK guide. – J. M.'s missing motivation Jul 09 '15 at 05:22