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

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.

-
For those of us who are not mathematicians would you spell out your desired operation in elementary detail? Or at least include links to such detail? –  Mr.Wizard Jun 7 '12 at 6:51
A subspace $W \subseteq V$ is invariant under a matrix $A_i$ if, for each $w \in W$, $A_i w \in W$. I want to find all subspaces invariant under all the matrices simultaneously. In practice, you'd do this by looking at the (generalized) eigenvectors of each matrix (unless there's some nice matrix trick that I'm missing which makes it easier.) –  Daniel McLaury Jun 7 '12 at 7:03
@DanielMcLaury Suppose all your $A_i$s are diagonal matrices .. What is your expected output? –  belisarius Jun 7 '12 at 7:09
It depends on what's on the diagonal. If every diagonal entry is different, and they're 3x3, then the output should be { { { 1 , 0, 0 } }, { { 0, 1, 0 } }, { { 0, 0, 1 } } }. (The list could also contain two- and three- dimensional combinations of these if necessary, although that information would be redundant.) –  Daniel McLaury Jun 7 '12 at 8:40
For the case of one matrix, Jordan (for exact entries) and Schur (for inexact entries) are available. For two matrices, Mathematica supports computing the Schur decomposition of a matrix pencil, but does not support Kronecker decomposition (the analog of Jordan decomposition for matrix pencils). For more than two matrices, I suppose one could try assembling a generalized companion matrix with matrix entries, but that is a messy affair if memory serves... –  Ｊ. Ｍ. Jun 7 '12 at 12:45

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 != {},
Do[(t=NewBasis[[i]].A[j];
If[
!InSpan[
Flatten[t]],
{i,1,Length[NewBasis]}, {j,1,m}];
OldBasis = Join[OldBasis, NewBasis];
(* 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.

-
Thanks for this thorough answer. I don't think there is a problem given that you have provided a Mathematica implementation. Is the expected output of the first block of code a tensor? Using A[j] matrices of random normal numbers, with $n=3$ and $m=4$ I get a 9*3*3 tensor. –  Verbeia Aug 15 '12 at 4:40
For step 0: IrreduciblePolynomialQ[] is pretty convenient. –  Ｊ. Ｍ. Aug 15 '12 at 4:46
Very well written. Thanks! –  PlatoManiac Aug 15 '12 at 6:22
@Verbeia Your output sounds fine. You can call it a tensor, but it might be clearer to call it a list of matrices: It is a basis for the subspace of the $n \times n$ matrices which is spanned by all possible products $A_{i_1} A_{i_2} \cdots A_{i_r}$ (for all $r$, so this is an a priori infinite problem). Since your list has $9$ elements, that means the monomials in your four matrices span the entire space of $3 \times 3$ matrices, which I don't find surprising. There are therefore no invariant subspaces (I should insert a check for this basic case. Editing...) –  David Speyer Aug 15 '12 at 9:59