Given a square matrix A, how can I generate a basis for the generalized eigenspace corresponding to all eigenvectors $\lambda_i$ such that $\vert \lambda_i \vert > 1$? I.e., if $A$ is $n \times n$ and acts on $\mathbb{R}^n$, then $\mathbb{R}^n = E^u \oplus E^1 \oplus E^s$, corresponding to the action of eigenvalues of modulus greater than / less than / equal to 1.

If A is diagonalizable, this can easily be accomplished with the Eigensystem command:

{evals, evecs}=Eigensystem[A];
Take[evecs, Length@Select[evals, Abs[#] > 1 &]]

This approach relies on the Eigenvalue command returning the eigenvalues (and hence eigenvectors) in decreasing order of their moduli. Also note that MMA returns the eigenvalues with repetition according to their algebraic multiplicity, and that the eigenvectors are returned as "row vectors" (technically, a list of vectors).

My problem arises when A is not diagonalizable. The natural function to look at would be JordanDecomposition, which retuns the generalized eigenvectors (N.B. as columns in a matrix!) and the Jordan Canonical Form. However, the ordering is a bit...abstract. As best as I've been able to tell from taking the JordanDecomposition of random matrices, the ordering of the Jordan blocks seems to be some combination of the order in which roots are returned by the Root command for irreducible factors of the characteristic polynomial, the ordinary Greater sorting, the size of the Jordan block, and whether it corresponds to a complex eigenvalue.

The kludge that I've been able to come up with is the following:

{s, j} = JordanDecomposition[A];
bigevals=Union@Select[Eigenvalues[A], Abs[#] > 1 &];
cols= Flatten[ Position[Diagonal[j], #]& /@ bigevals];
s[[All, cols]]

This seems to return the appropriate generalized eigenvectors.

The Questions

  1. Surely there's a better way? The above kludge seems...kludgy. Perhaps something like NullSpace[MatrixPower[A-lambda IdentityMatrix[First@Dimensions@A], n]] upon detecting that an eigenvalue lambda has algebraic multiplicity n, i.e., going directly to the definition?
  2. For my application, the matrix A will only ever be real valued. Thus all complex eigenvalues occur as complex conjugates (perhaps with multiplicity!). I would prefer to have a real basis. For example, if A={{0,2},{-2,0} then the eigenvalues are $\lambda_{\pm} = \pm 2 i$. Then $\left\{ \begin{bmatrix} 1 \\ 0 \end{bmatrix} , \begin{bmatrix} 0 \\ 1 \end{bmatrix} \right\}$ is a basis for the unstable manifold of $A$. Nullspace[MatrixPower[A,2]+4 IdentityMatrix[2]] gives the correct result here. Note that this question is a tertiary concern at best; there are other ways of dealing with this in "post-production."
  • 1
    $\begingroup$ Are you dealing with exact or approximate matrices? I assume approximate, in which case the Eigenvalues and seemingly the diagonal elements of j are sorted by largest Abs and so your reordering can be done much more cleanly (and quickly) by this: s[[All, Ordering[TakeWhile[Diagonal@j, Abs[#] > 1 &]]]] $\endgroup$ – b3m2a1 Aug 7 '18 at 19:17
  • 1
    $\begingroup$ That's because it can't exactly determine the Abs of a Root object and those only appear when the matrix is exact. $\endgroup$ – b3m2a1 Aug 7 '18 at 19:53
  • 1
    $\begingroup$ Seems odd that Eigenvalues still sorts when given an exact matrix though...? $\endgroup$ – erfink Aug 7 '18 at 20:16
  • 1
    $\begingroup$ The method you show looks okay to me. An alternative avoiding Position would be to Order eigenvalues by absolute value, and use that ordering to pick the desired generalized eigenvectors. $\endgroup$ – Daniel Lichtblau Aug 7 '18 at 22:37
  • 1
    $\begingroup$ @erfink The order of the generalized eigenvectors seems to be such that OrderedQ[Transpose[j, {1, 1}]] yields True. $\endgroup$ – Hector Aug 7 '18 at 23:29

Remember that the diagonal elements of the Jordan normal form are the eigenvalues. There is no need to compute the eigenvalues and the decomposition. Instead, the following code creates a picker list from the diagonal elements.

A = 1/8 {{27, 48, 81}, {-6, 0, 0}, {1, 0, 3}};
{s, j} = JordanDecomposition[A];
Pick[Transpose[s], Map[Abs[N[#]] > 1 &, Transpose[j, {1, 1}]]]

From there, you might want to use Orthogonalize


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.