Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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$?

[Edit] 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?

share|improve this question
    
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
add comment

2 Answers

up vote 1 down vote accepted

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.

share|improve this answer
    
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
add comment

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

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.