3
$\begingroup$

I stumbled upon an issue with computing eigenvalues of a matrix which I discovered while debugging some code which uses reconstructing the matrix from its eigenvalues/eigenvectors.

First of all, here is the (anti-hermitian) matrix in question:

mat=Uncompress["1:eJxTTMoPShNiYGAoZgESPpnFJWg8JhCPHUg45+cW5KRWFDFAAZwxGJT8y70gF+/F7IBHic3fJV9OcxzfPxicSzclKprPnolK/MXnaVNIuNgPBufClVAxDR4u7N6WuqeT9vFuSjiBEREd9A9rWiuxIZzAiMi+9PcRfcvB18fT1pm4/6Awaey7or17Ff6ALM+4bbF7pJWD1FFC/9J00NXF5+5drlxbhTeZNhBOYMuj9dfViY+wcpA6SuhfmhKbBqnTBqNOCfZqlui+ObFGBwdDjI0qobw0JTYNUqcNRrcSbKG678HRZDpEStNBVxfTTQl12qbUUUKdWoaI7AstBIZtv3j4KSGibUodi6hTy5QTTmBENIboHwGjaZAapSnHgeHT06d/aTqaBmmuJDZYVuRjAt5kOph6+vQvTUduGhxMPf2LTzuc5g2rLhQppenITYNDrKffUXTAdQglU1JKU2LT4FCdixxM7R5aKxlM46ak1DLEpsGhOhc5mNo9g0AJ3cZNSall6FsXj5amg99H1GmbktJOpm8aHC1NB7+PiFBCRNuUlHYyfdMgdZSMlqaDv6cPaZs6EVWaDsU0OFqaDouevlvixv2gZAoAconf8g=="];

You will notice, that the elements of mat differ in magnitude by six orders of magnitude. We compute eigenvalues and eigenvectors

{eigVals, eigVecs} = Eigensystem[mat];

and inspect what they look like:

Chop[eigVals, 10^-7]//Column

Now, we make a sanity check:

Chop[eigVecs.mat.ConjugateTranspose[eigVecs], 10^-7]//MatrixForm

and observe that there are off-diagonal elements which can get close to 100 in magnitude. This should definitely not be the case. It also appears that the eigenvectors are not orthonormalized properly, although they should according to the documentation of Eigensystem where it says:

For approximate numerical matrices m, the eigenvectors are normalized.

Chop[eigVecs.eigVecs\[ConjugateTranspose], 10^-7] // MatrixForm

reveals that there are off-diagonal elements ~0.25.

I am aware that the huge difference in scales does not render this an easy task, though I need to understand what is going wrong here and find a workaround. Can anyone give me some insight/fix/explanation?

$\endgroup$
  • $\begingroup$ Did you try to take the MatrixLog of your matrix (which I obviously cannot access), and to diagonalize it? The set of normalized eigenvectors of MatrixLog[mat] should coincide with those of mat. $\endgroup$ – user46676 Feb 21 '17 at 16:41
  • $\begingroup$ @marmot No, haven't tried that but I will. Thanks for the suggestion. But why can't you access the matrix? mat should be accessible via the Uncompress line at the beginning of the question. $\endgroup$ – Lukas Feb 21 '17 at 16:45
  • $\begingroup$ Sorry, didn't know what Uncompress does. The MatrixLog tricks appears to work. $\endgroup$ – user46676 Feb 21 '17 at 16:58
  • 3
    $\begingroup$ Lukas, I suspect that you may be suffering from loss of precision in your numerical calculations. Consider for instance setting the precision of your mat higher and repeating your calculations: {eigVals, eigVecs} = Eigensystem[SetPrecision[mat, 30]]; Chop[eigVecs.ConjugateTranspose[eigVecs] == IdentityMatrix[Length[mat]] returns True. $\endgroup$ – MarcoB Feb 21 '17 at 18:00
  • $\begingroup$ @MarcoB Thanks for this suggestion. In fact, the original problem comes from a Matlab simulation and I believe there is nothing like SetPrecision there. But glad to see it is such a simple fix in Mathematica. $\endgroup$ – Lukas Feb 21 '17 at 20:34
4
$\begingroup$

A possible workaround is based on MatrixLog

mat = Uncompress[
   "1:eJxTTMoPShNiYGAoZgESPpnFJWg8JhCPHUg45+cW5KRWFDFAAZwxGJT8y70gF+/\
F7IBHic3fJV9OcxzfPxicSzclKprPnolK/\
MXnaVNIuNgPBufClVAxDR4u7N6WuqeT9vFuSjiBEREd9A9rWiuxIZzAiMi+9PcRfcvB18fT1pm4/\
6Awaey7or17Ff6ALM+4bbF7pJWD1FFC/9J00NXF5+5drlxbhTeZNhBOYMuj9dfViY+\
wcpA6SuhfmhKbBqnTBqNOCfZqlui+\
ObFGBwdDjI0qobw0JTYNUqcNRrcSbKG678HRZDpEStNBVxfTTQl12qbUUUKdWoaI7AstBIZtv3j4KSGibUod\
i6hTy5QTTmBENIboHwGjaZAapSnHgeHT06d/aTqaBmmuJDZYVuRjAt5kOph6+\
vQvTUduGhxMPf2LTzuc5g2rLhQppenITYNDrKffUXTAdQglU1JKU2LT4FCdixxM7R5aKxlM46ak1DLEpsGhO\
hc5mNo9g0AJ3cZNSall6FsXj5amg99H1GmbktJOpm8aHC1NB7+\
PiFBCRNuUlHYyfdMgdZSMlqaDv6cPaZs6EVWaDsU0OFqaDouevlvixv2gZAoAconf8g=="];

lmat = MatrixLog[mat];
ev = Orthogonalize[Eigenvectors[lmat]];
Chop[ev.mat.ConjugateTranspose[ev], 10^-7] // MatrixForm
Chop[ev.ConjugateTranspose[ev]] // MatrixForm
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.