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
Best Answer
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: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
):