The idea behind the Jordan normal form does the trick, even though JordanDecomposition does not. (Incidentally, this suggests there may be a more reliable, stable algorithm to obtain Jordan decompositions than is implemented in Mathematica...) The resulting solution is very short, efficient, and numerically stable when applied to floating-point matrices. The following begins with a brief explanation of how it works, followed by several examples to show it actually does work.
Let's consider how any matrix can have an index exceeding $1$. Because the rank of $A^2$ is less than that of $A$, the nullspace of $A$ must be non-trivial. If the rank of $A^3$ further decreases, it must be because there is a sequence of vectors mapped by $A$, $v_2 \to v_1 \to 0$, with $v_1 \ne 0$. The index is the maximal length of such a sequence.
This provides a numerically stable, relatively quick way to compute the index. Beginning with computation of the nullspace $N_0$, find all vectors mapped by $A$ into the nullspace. In other words, find a basis for the solutions to $A x \in N_0.$ As far as I can tell, Mathematica does not have a direct method to do this, so let's reduce it to one it does have. NullSpace will provide a basis $e_1, \ldots, e_k$ of $N_0$. At the next step we seek a basis for the set of solutions
$$A x - \lambda_1 e_1 - \ldots - \lambda_k e_k = 0$$
This means the augmented vector $(x_1, x_2, \ldots, x_n; -\lambda_1, \ldots, -\lambda_k)$ lies in the nullspace of the matrix formed by augmenting $A$ with the columns $e_1, \ldots, e_k$. Among the solutions of this equation will be the original nullspace.
Iterating this procedure creates a flag of vector spaces $N_0 \subset N_1 \subset \cdots \subset N_i = N_{i+1} = \cdots$; $i+1$ is then the index. Provided we use a numerically stable procedure like LinearSolve or NullSpace, we can expect this algorithm to work even for large floating-point matrices. First I will provide an implementation and then illustrate it with several examples.
To help us study what's going on, this version returns a "generalized index" consisting of a sequence of bases of the flag: it gives far more information about the structure of $A$ than just its index, but obviously its index is easily derivable from the output: it is just its length. You should be able to recognize the iteration in NestWhileList, the augmentation of $A$ via Join, and the computation of nullspaces with NullSpace (including the initial nullspace, which is why this command appears twice).
index[a_] :=
With[{k = NullSpace[a]},
If[k == {}, {},
NestWhileList[
NullSpace[Join[a, Transpose[#[[All, 1 ;; Length@a]]], 2]] &, k,
Length[#1] != Length[#2] &, 2, Length@a, -1]]];
The test for termination is when the dimensions of the solutions stabilize (as computed by applying Length to their bases). Because termination occurs when $\dim(N_i) = \dim(N_{i+1})$, the final -1 in the NestWhileList command throws away the superfluous basis for $N_{i+1}$.
(Edit A special test has to be made for nonsingular matrices, for then Transpose fails when the nullspace is empty.)
To test this solution on large-ish matrices, let's generate some with known indexes. A good way to do this is to start with a bunch of Jordan blocks of zero eigenvalue: the index is one more than the longest contiguous string of superdiagonal ones (as is easily checked). Conjugating this by some random matrix (which is almost surely invertible) creates a non-sparse matrix for testing.
n = 30;
p = SparseArray[{Band[{1, 1}] -> 0, Band[{1, 2}] -> 1}, {n, n}];
j = Floor[n/3]; p[[j, j + 1]] = 0;
q = Rationalize[RandomReal[{0, 1}, {n, n}], 10^-2];
a = Inverse[q] .p . q;
ArrayPlot[p]
This is a plot of the Jordan normal form of $a$.
A direct calculation of its index will determine the ranks of the matrix powers. Because this matrix is designed to have an index of $20$, we can hard-code this into the check:
(ranks = MatrixRank[MatrixPower[a, #]] & /@ Range[20] ) // AbsoluteTiming
$\{7.1404084,\{28,26,24,22,20,18,16,14,12,10,9,8,7,6,5,4,3,2,1,0\}\}$
After seven seconds, we find that indeed the index is exactly $20$: $A^{19} \ne 0$ but $A^{20} = A^{21} = \cdots = 0$.
Of course, the corresponding numerical (floating point) calculation is far faster--but it gets the wrong answer:
(nRanks = MatrixRank[MatrixPower[N@a, #]] & /@ Range[Length@a] ) // AbsoluteTiming
$\{0.0100006,\{28,26,24,22,20,18,16,14,12,11,9,9,8,7,7,8,8,8,7,30,28,30,28,30,28,30,28,30,28,30\}\}$
Problems crop up around the $12^\text{th}$ power. This is evident in a plot of the two tests:
ListPlot[{ranks, nRanks}, PlotStyle -> PointSize[0.015], AxesLabel -> {"Power", "Rank"}]
Clearly the numerical calculations are producing garbage well before the correct index is reached.
I don't dare apply index to a itself: when using exact arithmetic, it's too slow! But let's see how it performs on the floating point version of a:
n = index[N@a]; // AbsoluteTiming
$\{0.0100006,\text{Null}\}$
It's as fast as the numerical brute-force rank-of-matrix-powers solution was. What about accuracy?
Length@n
$20$
It gets the right answer! But perhaps this was only luck? Let's check by restricting $A$ to the flag returned by index:
n0 = Reverse@(Last@n)[[All, 1 ;; Length[a]]];
u = Transpose[PseudoInverse[n0]] . a. Transpose[n0];
$u$ is just $A$ restricted to $N_{20}$, expressed in a different basis. Here are portraits of its powers through the $20^\text{th}$:
GraphicsGrid[{(MatrixPower[u, #]//Chop//ArrayPlot) & /@ Range@Length@n}, ImageSize -> 1000]
Up until the very end, the powers are nonzero, then finally the $20^\text{th}$ power is the zero matrix: this demonstrates the index of $A$ was at least $20$.
As another example, set $n=120$ in the previous one. My trials produce the generalized index (that is, the entire flag of $80$ subvectorspaces) in one-half to one second and consistently calculate the correct index of $80$. (I haven't tested with $n$ any larger than $120$ because it starts taking a long time just to create $a$.)
Finally, let's apply this solution to the matrix of the question:
a = {{1, -1, 0, 0, 0, 0}, {-1, 1, 0, 0, 0, 0}, {-1, -1, 1, -1, 0, 0},
{-1, -1, -1, 1, 0, 0}, {-1, -1, -1, 0, 2, -1}, {-1, -1, 0, -1, -1, 2}};
MatrixForm /@ index[a]
$$\left\{\left(
\begin{array}{cccccc}
0 & 0 & 1 & 1 & 1 & 1
\end{array}
\right),\left(
\begin{array}{ccccccc}
1 & 1 & 0 & 0 & 0 & 0 & 2 \\
0 & 0 & 1 & 1 & 1 & 1 & 0
\end{array}
\right)\right\}$$
The first element of the list is a basis of the nullspace of $A$, which is one-dimensional. The second element is a basis of the nullspace of an augmented version of $A$. We are interested only in the first six entries in each row. They form a basis for $N_1$. They include (at the bottom) the previous basis for $N_0$. As one last check, let's verify that $A$ sends the second basis into the span of the first:
a . Transpose@(Last@index[a])[[All, 1 ;; 6]] // Transpose // MatrixForm
$$\left(
\begin{array}{cccccc}
0 & 0 & -2 & -2 & -2 & -2 \\
0 & 0 & 0 & 0 & 0 & 0
\end{array}
\right)$$
Sure enough, the second basis vector is killed (it was in the nullspace) and the first is sent to a multiple ($-2$) of the second. This detailed additional information about precisely how $A$ achieves its index is potentially useful in applications.
MatrixRank[MatrixPower[A, k + 1]for incresingkand just compare the previous two entries. Sorry don't have time for a full answer now... $\endgroup$ – Ajasja Apr 24 '13 at 12:32MatrixRankin responses below should be better. $\endgroup$ – Daniel Lichtblau Apr 24 '13 at 14:00