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.