Is there a clean way to extract the subspaces invariant under a list of matrices?

Question

Let's say I have some $n \times n$ square matrices $A_1, A_2, \ldots, A_m$ with exact numbers for entries, and I want to find the subspaces of $V = \mathbb{C}^n$ invariant under these matrices. Is there a clean way to do this? The way I'm doing it seems to require writing ridiculous amounts of code from scratch to do things like calculate intersections, get all the invariant subspaces belonging to a particular generalized eigenspace, etc.

Answer 1

This answer is almost entirely about mathematics and algorithms, not Mathematica implementation. I'm not sure whether such answers are welcome on this site; I hope I haven't offended.

This is a very important problem in computational algebra, but is usually stated in a more sophisticated way. The usual way that one thinks about it is to consider the (noncommutative) algebra $R=\mathbb{C}\langle A_1, \ldots, A_m \rangle$ generated by the matrices $A_1$, ..., $A_m$. Think of $V$ as a module for $R$; your goal is to figure out whether the module is simple and, if not, find a submodule. The standard algorithm in this field is the Meataxe, introduced by Parker and improved by Ivanyos and Lux (among others). The problem of simply finding ONE nontrivial invariant subspace is hard enough that I will focus on that.

You REALLY don't want to solve the problem of describing all the invariant subspaces. Simply finding a way to display your output will be very deep as any projective variety can be described as a component of the variety of invariant subspaces for some module over some algebra; see this paper for a much more precise statement.

However, assuming that we restrict ourselves to finding ONE invariant subspace, I have a reasonable idea for a probabilistic algorithm which I would find reasonable to implement in a day.

Before starting, we recall the general structure theorem of $\mathbb{C}$-subalgebras of the $n \times n$ matrices (see sections 2 and 3 of Pete Clark's notes). There is always a composition $n=b_1+b_2+\cdots + b_k$, and a change of basis for $V$, such that every matrix in $R$ is block upper triangular with blocks of the size $b_i$ and the only relations between blocks on the diagonal is that some of them maybe imposed to be equal to others of them. The radical of $R$, denoted $\mathrm{Rad}(R)$, is the set of matrices which are all zeroes on the diagonal in this block presentation. Your goal is to detect whether or not this block decomposition is trivial.

Step 0 If your matrices have rational entries, I would compute the characteristic polynomial of $\sum q_i A_i$ for some random rationals $q_i$ and factor it. If you ever get an irreducible polynomial, then there is no invariant subspace over $\mathbb{Q}$, and the only way for there to be invariant subspaces over $\mathbb{C}$ is if all of $A_j$ commute with $\sum q_i A_i$ (in which case, the eigenspaces of $\sum q_i A_i$ are invariant.) So you might be able to stop here without going further.

Step 1 Let $C(R)$ be $\{ X : A_i X = X A_i \ 1 \leq i \leq m \}$. Compute a basis for $C(R)$; this is a linear algebra problem which I imagine you can make Mathematica do as well as I can. If $\mathrm{Id}_n$ doesn't span $C(R)$, let $X$ be an element of $C(R)$ not a multiple of of $\mathrm{Id}_n$. The eigenspaces of $X$ are nontrivial invariant subspaces of $R$.

From now on, we assume that $C(R)$ is trivial. This implies that either there are no invariant subspaces, or $R$ has nontrivial radical.

Step 2 Compute a basis for the algebra $R$. Here is some not very good Mathematica code for the purpose. I'm sure other members of the site could dramatically improve the implementation; I welcome them to edit in improvements.

(* n is the size of our matrices, the matrices are A[1], ..., A[m] *)
InSpan[L_, v_]:=MatrixRank[Append[L,v]] == MatrixRank[L];
OldBasis={}; NewBasis={IdentityMatrix[n]};
While[NewBasis != {},
  (NewerBasis = {};
   Do[(t=NewBasis[[i]].A[j]; 
      If[
        !InSpan[
            Map[Flatten, Join[OldBasis, NewBasis, NewerBasis]],
            Flatten[t]], 
        NewerBasis=Append[NewerBasis, t];]),
      {i,1,Length[NewBasis]}, {j,1,m}];
    OldBasis = Join[OldBasis, NewBasis]; 
    NewBasis=NewerBasis;)]
(* OldBasis now contains a basis for the ring R *)

Now, for any $v \in V$, you can test whether $v$ lies in an invariant subspace by computing MatrixRank[Map[(#.v)&, OldBasis]]. So our goal now is just to find good choices for $v$.

Step 2.5 If Length[OldBasis]==n^2, then $R$ is the entire ring of $n \times n$ matrices and there are no nontrivial invariant spaces; report this fact and stop.

Step 3 Pick a random element $\theta$ in $R$, by taking a random combination of the elements in OldBasis. Compute an eigenspace decomposition of $\theta$ and, for each eigenspace, see whether that eigenspace generates an $R$-submodule (using OldBasis as described above). Try this for a few random $\theta$'s, just in case.

If we have not terminated yet then (1) the partition $(b_1, \ldots, b_k)$ has more than one part (since $\dim R < n^2$ at Step 2.5) (2) all the $b_i$ are equal and all of the block diagonal matrices are equal (or, with high probability, Step 3 would have found an invariant subspace) and (3) the radical is nonzero (or $C(R)$ would have been nontrivial at Step 1. We assume from now on we are in this case. This is the "hard case" which Ivanyos and Lux describe how to attack. I've found a method which, to me, seems simpler than theirs, which I will now describe.

Step 4 We now want to solve the following auxiliary problem: Let $A$ be a matrix. Find the matrix $C$ which preserves the generalized eigenspaces of $A$, but acts on the generalized $\lambda$-eigenspace by the scalar $\lambda$. For example, if $A$ is diagonalizable, then $A=C$; if $A$ is nilpotent then $C=0$. This following solution returns $C$ while using exact arithmetic and never computes an eigenspace; I leave it to you as a puzzle to figure out why it works.

f=CharacteristicPolynomial[A,t];
g=D[f,t];
r=Simplify[f/PolynomialGCD[f,g]]; cr=CoefficientList[r,t];
s=D[r,t]; cs=CoefficientList[s,t];

newt[B_]:=
B - Sum[cr[[j]]*MatrixPower[B,j-1], {j,1, Length[cr]}].
   Inverse[Sum[cs[[j]]*MatrixPower[B,j-1], {j,1, Length[cs]}]]

FixedPoint[newt, A]

Step 5 Choose a random element $A$ of $R$ as before. Compute the element $C$ as above. It is highly probable that $A \neq C$ and the $R$-module generated by the image of $C-A$ is nontrivial. (More precisely, this will happen whenever the block diagonal part of $A$ is diagonalizable, but there is a nontrivial contribution from the radical.) Check that it is, to be safe.