# Stable method to compute $A^n$ for this defective matrix $A$?

I'm looking for a stable method to compute $A^n$, where $A$ is the following defective $12 \times 12$ matrix:

$$A = \left(\begin{array}{cccc|cccc|cccc} \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline \frac{1}{4} & \frac{1}{4} & 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & 0 & 0 & 0 & 0 & 0\\ \frac{1}{16} & \frac{3}{8} & \frac{1}{16} & 0 & \frac{1}{16} & \frac{3}{8} & \frac{1}{16} & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{1}{4} & \frac{1}{4} & 0 & 0 & \frac{1}{4} & \frac{1}{4} & 0 & 0 & 0 & 0 & 0\\ 0 & \frac{1}{16} & \frac{3}{8} & \frac{1}{16} & 0 & \frac{1}{16} & \frac{3}{8} & \frac{1}{16} & 0 & 0 & 0 & 0\\ \hline \frac{1}{16} & \frac{1}{16} & 0 & 0 & \frac{3}{8} & \frac{3}{8} & 0 & 0 & \frac{1}{16} & \frac{1}{16} & 0 & 0\\ \frac{1}{64} & \frac{3}{32} & \frac{1}{64} & 0 & \frac{3}{32} & \frac{9}{16} & \frac{3}{32} & 0 & \frac{1}{64} & \frac{3}{32} & \frac{1}{64} & 0\\ 0 & \frac{1}{16} & \frac{1}{16} & 0 & 0 & \frac{3}{8} & \frac{3}{8} & 0 & 0 & \frac{1}{16} & \frac{1}{16} & 0\\ 0 & \frac{1}{64} & \frac{3}{32} & \frac{1}{64} & 0 & \frac{3}{32} & \frac{9}{16} & \frac{3}{32} & 0 & \frac{1}{64} & \frac{3}{32} & \frac{1}{64} \end{array}\right).$$

Since $A$ is defective, it doesn't have $12$ linearly independent eigenvectors. It is therefore not possible to use eigendecomposition (i.e. write $A$ as $V \Lambda V^{-1}$ and therefore $A^n$ as $V \Lambda^n V^{-1}$).

An alternative would be to use the Jordan normal form (i.e. write $A$ as $W J W^{-1}$ and therefore $A^n$ as $W J^n W^{-1}$), but I don't have any experience with this method. Moreover, I read that the numerical approach for finding $J$ can be very unstable. Is it manageable to analytically compute $J$ for this matrix $A$?

 Notice that $A$ is a block lower triangular matrix (or is it called a lower triangular block matrix). Can we use this structure to compute $J$ in a block-wise fashion?

-
Singularity has nothing to do with eigendecomposition. Singularity means zero is an eigenvalue. $A$ is diagonalizable if and only if it has 12 linearly independent eigenvectors. These two properties are independent of each other --- knowing a matrix does or doesn't have one property does not tell you whether or not it has the other. – Gerry Myerson Jan 15 '13 at 10:04
@GerryMyerson, you're completely right! I meant to say that $A$ is defective, so it doesn't have $12$ linearly independent eigenvectors. Therefore, $V$ is singular and cannot be inverted, so the eigendecomposition cannot be used. I will update my question accordingly :) – Ailurus Jan 15 '13 at 10:11
Which matrix $V$? Since there do not exist $12$ linearly independent eigenvalues you can only find fewer, and what do you put in the remaining columns? (And even if $A$ were diagonlisable, one could not speak of the matrix $V$, as one can scale and permute eigenvectors.) – Marc van Leeuwen Jan 15 '13 at 10:35
@MarcvanLeeuwen, I never thought about that. Most of the times I just use MATLAB to compute the eigenvectors of an $n \times n$ matrix, and it always returns $n$ eigenvectors, even if the matrix is defective. – Ailurus Jan 15 '13 at 11:11

The eigenvalues $\lambda_k$ of this matrix are $$1,\ {1\over2},\ {1\over2},\ {1\over2},\ {1\over4},\ {1\over8},\ {1\over8},\ {1\over16},\ {1\over16},\ {1\over32},\ {1\over64},\ 0\ .$$ This implies that you can find a suitable Jordan basis by purely rational calculations. For this you have to determine the kernels of the maps $$A-I,\quad (A-{1\over2} I)^3,\quad A-{1\over4} I, \quad(A-{1\over8} I)^2,\quad (A-{1\over16} I)^2, \quad A-{1\over32}I, \quad A-{1\over64}I, \quad A$$ which are "independent" and whose dimensions add up to 12. For the kernels of dimension $>1$ you have to adjust the basis such that $A-\lambda_k I$ has ones in the second diagonal and all the rest zeros.

In the end you have a basis of all of ${\mathbb R}^{12}$ with respect to which $A$ assumes Jordan form.

All this can be handled by Mathematica, Maple or similar without going into floating point numerics.

-
Thanks! I just figured out these eigenvalues by exploiting the block triangular structure of the matrix (i.e. they are the union of the eigenvalues of the three $4 \times 4$ blocks on the diagonal). Is there any particular reason to put the eigenvalues in this order (or did you intend to start with $1$ and end with $0$)? – Ailurus Jan 15 '13 at 13:31
Could you perhaps elaborate a little on how to construct the matrix $W$, such that we have $A=WJW^{-1}$? Constructing $J$ is relatively simple, but I'm not so sure about $W$. Additionally, with regard to your remark about floating point numerics, would it be all right to use MATLAB's null() function to determine the kernel of a matrix? – Ailurus Jan 15 '13 at 14:28
@Ailurus: See $\S$ 57 in Halmos' Finite dimensional vector spaces. Your $A-\lambda_k I$ are nilpotent on the kernels of $(A-\lambda_k I)^r$. - I'm not fluent with contemporary Matlab. – Christian Blatter Jan 15 '13 at 14:36
Will do! Mathematica code would also be fine, I'm just more comfortable with Matlab. – Ailurus Jan 15 '13 at 15:22

I think Cayley–Hamilton theorem will work for your problem. But it only works with square matrice as in your case it is square.

-
If you would like to solve it block-by-block, then Kalman Decomposition Theorem can work for that purpose but I don't think so that you may need it because Cayley-Hamilton Theorem will work for you. – Salman Jan 15 '13 at 10:33
Thanks for suggesting the Cayley-Hamilton theorem. However, I should probably have mentioned that these $A^n$ computations are happening at runtime. Therefore, computing $A^n$ should be as efficient as possible, which is the main reason for me to propose the Jordan normal form. Once $J$ is known, $A^n$ can be computed quite efficiently for any $n$. – Ailurus Jan 15 '13 at 11:27
Ohh.. I look for it. – Salman Jan 15 '13 at 14:33