6
$\begingroup$

When I run Eigensystem on a symmetric matrix, the list of eigenvalues (and so, the corresponding eigenvectors) is ordered by absolute values, which is quite bizarre (it does make sense if you view the symmetric eigenvalue problem as a special case of the general eigenvalue problem, where the eigenvalues are complex, so modulus is as good a way to order them as any). Has it always been thus? Is this considered a feature?

If one wants the "obvious" ordering, the following function works fine, but it's a bit annoying.

sortify[es_] := Module[{pp = FindPermutation[es[[1]]]}, PermutationReplace[#, pp] & /@ es]

EDIT My personal view is that this is a bug. Mathematica (or any system, for that matter) should always return things in canonical order, which happens to be the one coming from the total ordering on the real line for real numbers.

$\endgroup$
  • $\begingroup$ Closely related: Find the eigenvector associated with the smallest eigenvalue, not smallest in magnitude $\endgroup$ – Jens Sep 8 '14 at 15:10
  • $\begingroup$ @Karsten 7. thanks for formatting the code properly. How do you actually do that? :) $\endgroup$ – Igor Rivin Sep 8 '14 at 15:11
  • 2
    $\begingroup$ Umm, the use of the "bug" word seems quite unmotivated. Also there is a dearth of actual code here to even understand the issue. For example, are we talking about exact or approximate matrices? Matrices with multiple eigenvalues (might matter, if input is approximate). $\endgroup$ – Daniel Lichtblau Sep 8 '14 at 17:25
  • 3
    $\begingroup$ (1) No, not joking. Without input (a 2x2 or 3x3, say) or other description, I cannot tell whether you refer to exact or approximate input. Which matters. (2) It's not a bug, design or otherwise. It is simply an outcome you do not happen to like. Were it to change it would be inconsistent with the handling of the nonsymmetric case (unless we change that as well). Such a change might comprise a reasonable suggestion. But that's not the same thing as the current behavior being a bug. $\endgroup$ – Daniel Lichtblau Sep 8 '14 at 17:53
  • 1
    $\begingroup$ (1) Currently the ordering or eigenvalues between real-symmetric and general cases is consistent. If we change the latter case, e.g. to the ordering you indicate, they would become inconsistent. (2) The term "fix" is out of place since it isn't at present broken. (3) The documentation states: "If they are numeric, eigenvalues are sorted in order of decreasing absolute value." What part of the semantics takes figuring out? $\endgroup$ – Daniel Lichtblau Sep 8 '14 at 20:03
10
$\begingroup$

The order of eigenvalues is the most convenient order for the algorithm, which find these eigenvalues. You can always order them as you want very simply

a = # + #\[Transpose] &@RandomReal[1, {10, 10}];
{ε, ψ} = Eigensystem[a];
{ε, ψ} = {ε[[#]], ψ[[#]]} &@ Ordering[ε];

Furthermore, the eigenvalues can be complex for non-Hermitian matrices. There is no native order of complex values.

$\endgroup$
  • $\begingroup$ The way you are suggesting is a somewhat more elegant variant of what I am suggesting in the question, but the lexicographic ordering of complex numbers is a natural extension of the "usual" ordering of real numbers. $\endgroup$ – Igor Rivin Sep 8 '14 at 16:23
  • 2
    $\begingroup$ @IgorRivin It is not just more elegant. It works much faster then FindPermutation and about 10% faster then Jens's solution. $\endgroup$ – ybeltukov Sep 8 '14 at 16:42
  • 1
    $\begingroup$ @ybeltukov Yes indeed, this is very nice and captures the spirit of what Igor was suggesting (+1). $\endgroup$ – Jens Sep 8 '14 at 16:44
3
$\begingroup$

As I mention in the linked answer, you can always translate the spectrum of the matrix by the Hilbert-Schmidt norm to be certain that Mathematica's ordering of the eigenvalues will coincide with the natural order on the real line. For a large matrix, this is more efficient that FindPermutations.

Edit

Here is an implementation of the shift approach:

Clear[sortedEigensystem];
sortedEigensystem[
   matrix_?MatrixQ] :=
  (Eigensystem[
       matrix - # IdentityMatrix[Dimensions[matrix]]] + {#, 0}) &@
   Norm[Flatten[matrix]];

Responding to the comment, as a test, consider the matrix

d = (# + Transpose@#) &@N@RandomInteger[{-10, 10}, {2000, 2000}];
AbsoluteTiming[eigs = sortedEigensystem[d];]

(* ==> {2.220592, Null} *)

On the other hand, trying the approach in the question for the large matrix d runs for a very long time and ultimately fails to produce an ordered list of eigenvalues and eigenvectors.

Edit 2

When doing machine arithmetic computations, it's worth keeping in mind that Eigensystem is not the only normal form that sorts eigenvalues according to absolute value. Another function that does the same thing is JordanDecomposition. But its output is not in the same form as Eigensystem. The eigenvalues appear on the diagonal of a generally upper-triangular matrix in block form. In order to sort the output consistently, one would have to sort not only these blocks but also find the corresponding similarity transformation matrix.

On the other hand, the approach based on shifting the spectrum can be directly applied to this normal form as well:

Clear[sortedJordanDecomposition];
sortedJordanDecomposition[
   matrix_?MatrixQ] := (JordanDecomposition[matrix - #] + {0, #}) &[
   Norm[Flatten[matrix]] IdentityMatrix[Dimensions[matrix]]];

{s, j} = sortedJordanDecomposition[d];

Diagonal[j] == eigs[[1]]

(* ==> True *)

The last line shows that the Jordan blocks are now automatically sorted the same way as for orderedEigensystem above.

$\endgroup$
  • $\begingroup$ Actually, I disagree with this. for an $n\times n$ matrix, computing the norm is $O(n^2),$ finding the permutation is $O(n \log n).$ But, more importantly, both these are insignificant compared to the cost of actually computing the eigensystem. $\endgroup$ – Igor Rivin Sep 8 '14 at 15:49
  • $\begingroup$ I checked the timings for the example from the linked question. See my update. $\endgroup$ – Jens Sep 8 '14 at 16:08
  • $\begingroup$ Interesting! The natural question is why the FindPermutation[] thing fails (looks like a bug to me), but in the meantime I agree that yours is the way to go (since it actually works...) $\endgroup$ – Igor Rivin Sep 8 '14 at 16:12
  • $\begingroup$ @IgorRivin Yes - it certainly shouldn't fail without warning in your approach. Maybe it's a memory issue. $\endgroup$ – Jens Sep 8 '14 at 16:16
  • $\begingroup$ In principle my approach should also work for JordanDecomposition, but I have to verify how the latter orders its Jordan blocks under different conditions. This may be an advantage of my approach, in that it doesn't need different sorting strategies for different normal forms. $\endgroup$ – Jens Sep 8 '14 at 17:28

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.