# Ways of computing $A^\infty$

As a follow-up on this question, I would like to ask which one is the better way of computing $A^\infty = \lim_{n \rightarrow \infty} A^n$.

1. Repeatedly square, computing $A^2, A^4, A^8$ and so on.
2. Do the Jordan decomposition, as described in the answer to this question.

The second way seems more "closed form", but I am concerned about the numerical stability of the Jordan decomposition.

Let me state that I only really care about computing $A^\infty$ approximately, in a numerical setting on a computer i.e. $A$ is a matrix which contains double-precision floating point numbers. I would also be satisfied with a close solution $A_{\text{close}}^ \infty$, provided $\| A_{\text{close}}^\infty - A^\infty\| \leq \varepsilon$, for some small $\varepsilon$.

However, the matrix may contain entries where the values are input by hand, i.e. there is no guarantee that the eigenvalues will be distinct.

-
Repeated squaring isn't a method of computing $A^\infty$ unless the eigenvalues consist solely of zeros and (power of 2) roots of unity, so that some finite number of squarings results in an idempotent matrix. It is only a method of approximating the limit. –  Glen O May 3 '13 at 12:17
@Glen: Given that the OP speaks of "numerical stability", the matrix $A$ is probably approximate to begin with, and approximating the limit is all you can do. –  Hurkyl May 3 '13 at 12:21
@Hurkyl: Numerical stability is a matter of computation, and does not require that the matrix is approximate. Repeated squaring also has numerical issues, arising from the effects of truncation. Meanwhile, there is no guarantee that repeated squaring will converge. –  Glen O May 3 '13 at 12:25
Furthermore, repeated squaring may converge where it should not, such as if the matrix has $\lambda=-1$ as an eigenvalue. –  Glen O May 3 '13 at 12:32
@Glen: Indeed, if $A^n$ does not converge, then repeated squaring may not work. But I suppose that if $A^n$ does converge, then repeated squaring should always work correctly. (And in "most" cases, $A^n$ does converge.) –  TMM May 3 '13 at 12:45
show 1 more comment

If $A$ is a matrix of real or complex numbers stored numerically (with error), with very high probability it will have $n$ distinct eigenvalues and hence be diagonalizable. In this case the Jordan decomposition will actually be a diagonal matrix consisting of the eigenvalues. You can't beat that.
The matrix $[\begin{smallmatrix}1&1\\0&1\end{smallmatrix}]$ is similar to $[\begin{smallmatrix}3&-1\\4&-1\end{smallmatrix}]$. Now multiply this entrywise by $\frac{1}{3}$. In exact arithmetic it's not diagonalizable, with Jordan form $\frac 13[\begin{smallmatrix}1&3\\0&1\end{smallmatrix}]$. But if you change the $-\frac 13$ to $-0.333333333333$, or anything similar, you will find it is diagonalizable. –  vadim123 May 7 '13 at 20:51
Ok, so we have the following: 1) $\left[\begin{matrix} 1 & 1 \\ 0 & 1 \end{matrix}\right]$ is not diagonalizable, if the entries are double precision floating point numbers; 2) $\frac 1 3 \left[ \begin{matrix} 3 & -1 \\ 4 & -1 \end{matrix}\right]$, where the entries again are double precision floating point numbers is diagonalizable because of the quantization error that comes from projecting 1/3 onto the set of numbers representable in the double precision floating point format. But 1) is a perfectly valid input, what makes you think that one will only encounter matrices of type 2)? –  ziutek May 8 '13 at 13:42