Matlab code for computing index of matrix

Question

I need help to make Matlab code for computing index of matrix.

The index of matrix $A$ is the least non negative integer $k$ such that $\operatorname{rank}(A^{k+1}) = \operatorname{rank}(A^k)$.

If it is not possible to make Matlab code for this problem. What other ways can I use to find the index of matrix.

Please help me. I would be very much grateful to you.

Thanks

Answer 1

I'm new to this topic, so any critical reading of this answer would be appreciated – but I think one can do without looping over powers.

Every square matrix $A$ has a Jordan normal (or canonical) form $J =$

where $\lambda_i$, $i = 1 \ldots n$, are the eigenvalues of $A$, and all entries shown as blank are zero. Since $J$ is related to $A$ by a similarity transformation, $J = P^{-1} A P$, and this extends to integer powers, $J^k = P^{-1} A^k P$, it holds $$ \DeclareMathOperator{\rank}{rank} \rank(A^k) = \rank(J^k). $$

$J$ has a block-diagonal form, where the size of each Jordan block corresponds to the algebraic multiplicity $a_i$ of the respective eigenvalue $\lambda_i$. The total rank of $J$ is therefore the sum of the ranks of the Jordan blocks. Blocks with $\lambda_i \neq 0$ are of full rank, with $\lambda_i = 0$ of rank $a_i - 1$: $$ \rank(J) = \sum_{\lambda_i \neq 0} a_i + \sum_{\lambda_i = 0} (a_i - 1). $$

The powers $J^k$ also have block-diagonal form, where each block of $J^k$ is the $k$th power of the corresponding block of $J$. Full-rank blocks remain of full rank under finite powers $k$, while blocks with $\lambda_i = 0$ are nilpotent matrices of degree $a_i$ whose rank decreases with increasing $k$ until it reaches 0. Therefore: $$ \rank(J^k) = \sum_{\lambda_i \neq 0} a_i + \sum_{\lambda_i = 0} \max(a_i - k, 0) $$

The rank of $J^k$ and with it $A^k$ therefore decreases with increasing $k$ until the power of each Jordan block with $\lambda_i = 0$ has reached $a_i$. Consequently, the index of $A$ is the maximum of the degrees of the nilpotent blocks of $J$, or

$$ \mathrm{index} (A) = \max_{\lambda_i = 0} ~ a_i . $$

---

Matlab implementation using jordan from the Symbolic Math Toolbox:

J = jordan(A);
d = diag(J);                            % diagonal
sd = diag(J, 1);                        % superdiagonal
lambda = d([sd ; 0] == 0);              % eigenvalues
a = diff([0 ; find([sd ; 0] == 0)]);    % arithmetic multiplicities
index = max(a(lambda == 0));

Note however that this procedure is numerically problematic. Almost all square matrices are diagonalizable, which means a small numeric error in the representation of a non-diagonalizable matrix will often make it diagonalizable (see TMW's text about the Jordan Canonical Form). Since the rank of a diagonalizable matrix does not change under finite powers $k$, its index is 1. All of this means that the index of a matrix is not a numerically stable property. It is no mistake that the function jordan, though it also works with numerical matrices, is included with Matlab's symbolic math functions.

A related problem is that the code above uses a comparison of eigenvalues with zero. This might be fixed to some degree by using a finite threshold, as in the following code (adapted from rank.m):

tol = max(size(A)) * eps(max(lambda));
index = max(a(abs(lambda) < tol));

Answer 2

Naïve proposal: (Only recommended for small, stable matrices):

k = rank(A); 
index = 0;
A_copy = A;
A = A*A_copy;
aux = rank(A);
while(aux == k)
    index = index+1;
    k = rank(A);
    A = A*A_copy;
    aux = rank(A);
end

Where the final result is index. I hurried this, so any corrections will be welcomed.