# How to find the index of a square matrix in Mathematica quickly?

Let $A$ be an $n\times n$ complex matrix. The smallest nonnegative integer $k$ such that $\mathrm{rank}(A^{k+1})=\mathrm{rank}(A^{k})$, is the index fo $A$ and denoted by $\mathrm{Ind}(A)$. I would like to compute $\mathrm{Ind}(A)$ quickly in Mathematica (I am using V8).

Let us as a very simple example consider

$$A=\left(\begin{array}{rrrrrr} 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 \end{array}\right)$$

then it is clear that $\mathrm{Ind}(A)=2$. For computing this in Mathematica, I used

Solve[MatrixRank[MatrixPower[A, k + 1]] == MatrixRank[MatrixPower[A, k]] && k > 0,
k, Integers]


But unfortunately, I cannot get the result. I will be grateful if anyone could design a way to compute $\mathrm{Ind}(A)$ for random matrices of the size $n\times n=200\times 200$ in Mathematica.

-
Hmmm, can't you just compute MatrixRank[MatrixPower[A, k + 1] for incresing k and just compare the previous two entries. Sorry don't have time for a full answer now... –  Ajasja Apr 24 '13 at 12:32
In theory, one could look carefully at the Jordan block that goes with the eigenvalue of zero. In practice this will not scale very well, and the numeric methods involving MatrixRank in responses below should be better. –  Daniel Lichtblau Apr 24 '13 at 14:00
@Daniel That's correct. I implemented the Jordan block method (with JordanDecomposition): it processes a $200$ by $200$ floating matrix in under 0.1 seconds but the numerical error is apparent even in tiny matrices; just a little bit of error destroys the calculation of the index (all the eigenvalues become distinct under the tiny perturbations caused by FP arithmetic and the index always seems like it's equal to $1$). But why should MatrixRank perform any better with FP matrices? My experiments show it fails as badly, even on $3$ by $3$ matrices. –  whuber Apr 24 '13 at 14:28
@whuber MatrixRank should use singular values under the hood if I am not mistaken. I think this would be fairly reliable. Can you post or send me the 3x3 where it fails? –  Daniel Lichtblau Apr 24 '13 at 15:30
@whuber (part 2) JordanDecomposition is not a good thing to try on non-exact input. I was never strongly in favor of extending it to approximate number input. To this day we find and fix bugs from that functionality. I'm not sure if we gained any advantage from having it. –  Daniel Lichtblau Apr 24 '13 at 15:33

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.

-
Overkill for the problem at hand but as the first to upvote this I will say it has an aesthetic appeal as well as independent interest. Could perhaps be adapted to a Jordan decomposition algorithm, for example, since finding those nontrivial block eigenspaces is the hard part. –  Daniel Lichtblau Apr 25 '13 at 13:36
Excellent idea. Can you check that why your technique cannot produce the result for a simple matrix in floating points arithmetic as follows: $n=100; a = RandomReal[{}, {n, n}]$. I faced with some errors such as "Transpose::nmtx: "The first two levels of the one-dimensional list {} cannot be transposed."". –  Fazlollah Soleymani Apr 25 '13 at 16:28
Thanks: I tacitly assumed $A$ is singular; if not, the index is $0$. I will add a test for this. However, the expression RandomReal[,n,n] is not even syntactically correct. The message you got indicates you created a vector rather than a matrix. An expression that will give you an actual $n$ by $n$ matrix would be RandomReal[{0,1},{n,n}]. For almost all such random matrices the index is $0$, so they don't make a very good test. Besides, you don't know in advance what the index is! Test instead with the code I provide, which produces exact matrices of known index. –  whuber Apr 25 '13 at 17:39

Try this:

Length@NestWhileList[A.# &, A, MatrixRank[#] != MatrixRank[#2] &, 2] - 1


This'll keep multiplying $A$ together and checking the matrix rank at each step.

-
f@x is a shorthand for f[x] and # and & are used in writing pure functions. The rest should be easy to find in the documentation. –  Szabolcs Apr 24 '13 at 12:37
Since MatrixRank is relatively costly (compared to Dot) this could be improved by caching the rank of of the matrix between steps. –  Mr.Wizard Apr 24 '13 at 12:42
@Mr.Wizard Yes, there's often a tradeoff between simplicity and performance. Why don't you post the faster solution? –  Szabolcs Apr 24 '13 at 12:43
@Szabolcs try testing your code with A={{1, 1, -1, 2}, {-1, 1, 0, 0}, {0, 1, -1, 2}, {1, 2, -1, 0}}. In my machine the kernel keeps on running.... –  PlatoManiac Apr 24 '13 at 12:57
I think the code is not working well, Mr. Wizard is right down here. For a random matrix, it sould give 0 while produces 5 as the $Ind(A)$. –  Fazlollah Soleymani Apr 24 '13 at 13:25

Here is the bad type of example mentioned in comments by @whuber.

SeedRandom[1111];
n = 3; p =
SparseArray[{Band[{1, 1}] -> 0, Band[{1, 2}] -> 1}, {n, n}]; q =
Rationalize[RandomReal[{0, 1}, {n, n}], 10^-3]; a = Inverse[q].p.q;
na = N[a];


MatrixRank[MatrixPower[na, #]] & /@ Range[Length@na]

(* Out[318]= {2, 1, 3} *)


To see the source of trouble we'll have a look at singular values.

svlist =
Map[SingularValueList[#, Tolerance -> 10^(-6)] &,
FoldList[Dot, na, Table[na, {n - 1}]]]

9* Out[319]= {{13.1845622926,
1.63446803072}, {21.5497455662}, {2.03322204001*10^-14,
2.76424896137*10^-15, 1.86215833311*10^-16}} *)


Notice that the tolerance setting does nothing for that last list (in this example it does nothing for the others either, but that need not be the case in general). The reason is that all of them are small, and none are a factor of 10^6 or more smaller than the largest.

What we want in this situation is to have also an absolute threshold for removing them, as Tolerance only uses a relative threshold. Chop can do this. I will illustrate with a larger example.

SeedRandom[1111];
n = 100;
p = SparseArray[{Band[{1, 1}] -> 0, Band[{1, 2}] -> 1}, {n, n}]; q =
Rationalize[RandomReal[{0, 1}, {n, n}], 10^-3]; a = Inverse[q].p.q;
na = N[a];

Timing[
svlist = Map[SingularValueList[#, Tolerance -> 10^(-6)] &,
FoldList[Dot, na, Table[na, {n - 1}]]];]

(* Out[340]= {0.180000, Null} *)

Map[Length, Chop[svlist, 10^(-6)] /. 0 :> Sequence[]]

(* Out[341]= {99, 98, 97, 96, 95, 94, 93, 92, 91, 90, 89, 88, 87, 86, \
85, 84, 83, 82, 81, 80, 79, 78, 77, 76, 75, 74, 73, 72, 71, 70, 69, \
68, 67, 66, 65, 64, 63, 62, 61, 60, 59, 58, 57, 56, 55, 54, 53, 52, \
51, 50, 49, 48, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, \
34, 33, 32, 31, 30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, \
17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0} *)


The upshot is you might want to experiment with setting Tolerance, and postprocessing using Chop, on SingularValueList[] of these powers.

-

I believe this is correct but I'm running out of time to check it or develop it further. I believe Szabolcs intented to write != where he wrote ==, but there is still the matter from starting with power zero.

Update #3:

f1 = {#, MatrixRank@#} &;
f2[m_][{_, {a2_, r2_}, n_}] := {{a2, r2}, f1[m.a2], n + 1}
f3[m_] := NestWhile[f2[m], {{-1, -1}, f1 @ MatrixPower[m, 0], -1}, #[[1, 2]] != #[[2, 2]] &]


Use:

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}};

f3[A]


Output is in the form: {{matrix1, rank1}, {matrix2, rank2}, power}
You can use Part to extract whatever you need.

-
I think Mr. Wizard answer needs some revision. I mean the answers are not the same. For instance, check the following sample examples: Clear["Global*"] SeedRandom[1234]; n = 500; A = RandomReal[{}, {n, n}]; Length@NestWhileList[A.# &, A, MatrixRank[#] == MatrixRank[#2] &, 2] // AbsoluteTiming Length@NestWhileList[{A.#[[1]], MatrixRank@#[[1]]} &, {A, -1}, #[[2]] == #2[[2]] &, 2] // AbsoluteTiming –  Fazlollah Soleymani Apr 24 '13 at 12:56
@FazlollahSoleymani Sorry, let me see if I can fix that. –  Mr.Wizard Apr 24 '13 at 13:00
@FazlollahSoleymani I don't know that I understand the question. Can you confirm that the output of Szabolcs's code is in fact the result that you want? –  Mr.Wizard Apr 24 '13 at 13:07
@FazlollahSoleymani please take for example SeedRandom[12345]; A = RandomReal[{}, {500, 500}]; then Table[MatrixRank[MatrixPower[A, i]], {i, 0, 5}] is {500, 500, 500, 500, 500, 494}` -- wouldn't the index be either zero or one? Szabolcs's code yields five. –  Mr.Wizard Apr 24 '13 at 13:11
@FazlollahSoleymani I think my updated answer now agrees with the first one, but I still don't think I understand the question. If the Rank of both A^1 and A^2 is 500, why isn't the answer one? –  Mr.Wizard Apr 24 '13 at 13:17