Should eigenvalues be ordered?

Question

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.

Answer 1

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.

Answer 2

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.