# Expressing the values of a matrix at pow N

I have a square matrix (that comes from a Markov Chain) that looks like that:

$$Q = \begin{bmatrix} 0 & 1& 0 & 0 & .. & 0 & 0\\ 0 & a & 1-a & 0 & .. & 0 & 0\\ 0 & 0 & b & 1 - b & .. & 0 & 0\\ .. & .. & .. & .. & .. & .. & .. \\ 0 & 0 & 0 & 0 & .. & 0 & 1 \end{bmatrix}$$

with $a, b, c,$ etc. real numbers between $0$ and $1$ included.

I am interested in the values of the first line of the matrix $Q^N$. Currently, I am using a library (numpy) that allows me to compute $Q^N$ and then I read the first line of this matrix. But with large matrices ($> 500 \times 500$) and large values of $N$ (~ 10000), it is a bit too slow for my usage.

By curiosity, I have plotted $Q^N_{0,j}$ for $N$ between 1 to 1000 and I found that they follow something that looks like a Poisson distribution or similar (but I don't know if it's only due to my specific input matrix $Q$ or not).

My question is, given such a matrix Q, is there a way to get the values of $Q^N_{0,j}$ without having to compute $Q^N$?

Edit: the terms on the diagonal are such as: 0 <= a <= b <= c <= ... < 1

Edit2:

It appears that if I can diagonalize $Q$, I can use $Q^n = P D^n P^{-1}$ which is faster for large values of $n$ than an exponentiation by squaring (as used by numpy). Problem is that I am not sure it is possible for any matrix $Q$. And if it's not possible, I'd accept a slight modification of $Q$, $P$ or $D$ if the result is close enough.

Edit3:

The values of the first line of $Q^N$ for $N \in [1,100]$:

Edit4:

it's $a <= b$ and not $a < b$, sorry!

-
If ${\bf R}^n$ has a basis consisting of eigenvectors for $Q$, then $Q^N=PD^NP^{-1}$, where $P$ is the matrix whose columns are the eigenvectors of $Q$, and $D$ is the diagonal matrix whose diagonal entries are the eigenvalues of $Q$. The point is, it is very easy to compute $D^N$, $D$ being diagonal. – Gerry Myerson Sep 20 '13 at 12:47
@GerryMyerson After writing my question, I wondered how the matrix power can be done so I'd tend to agree with your comment. But I'm wondering what is so long when I ask numpy to compute $Q^n$, could it be the time to compute the eigenvectors of Q and copying values? Maybe I could save time by doing it myself and save the $P$ and $D$ matrices (since I am computing $Q^N$ for thousands of $N$ values). Am I right to say that the eigenvalues of $Q$ are $(0, a, b, ..., 1)$ (because the matrix is triangular)? – Maxime Sep 20 '13 at 13:29
I don't know what "numpy" is, much less how it works, but, yes, you have the eigenvalues right. – Gerry Myerson Sep 20 '13 at 13:39
@GerryMyerson, Thanks. numpy is a Python library that allows to do linear algebra (and not only). numpy.org I assume that it needs to compute the matrices $P$ and $D$ in order to compute $Q^N$. Now that I ask him to compute $P$, $D$ and $P^{-1}$ before, it's 8 times faster to compute thousands of powers. – Maxime Sep 20 '13 at 14:31
@GerryMyerson One other question: Am I sure that my matrix Q will not be defective? For any dimension and values of a, b, c, etc. ? – Maxime Sep 20 '13 at 15:12

A general bidiagonal matrix is not always diagonalizable; for instance if $a=b=\ldots=\frac{1}{2}$ you will get a Jordan block of eigenvalues $\frac{1}{2}$.

But, since all of your entries are distinct, you're in luck. I'm going to change your notation a bit:

$$M = \left[\begin{array}{ccccc} a_1 & 1-a_1 & & &\\ & a_2 & 1-a_2 & &\\ & & \ddots & &\\ & & & a_{n-1} & 1-a_{n-1}\\ & & & & 1\end{array}\right].$$ (In your case, $a_1 = 0$.)

Then the eigenvalues of $M$ are simply $a_1, \ldots, a_{n-1},1$. The eigenvector corresponding to $1$ is obviously the constant 1 vector. The eigenvector for $a_i$ is the vector $v^i$, where $$v^i_j = \begin{cases}0, &j > i\\ \frac{\prod_{k=j}^{i-1} (a_k-1)}{\prod_{k=j}^{i-1} (a_k-a_i)}, & j \leq i\end{cases}$$ and $v^i_i = 1$ by the empty product convention. You can check this formula by looking at the $j$th row of $Mv^i$:

$$a_j \frac{\prod_{k=j}^{i-1} (a_k-1)}{\prod_{k=j}^{i-1} (a_k-a_i)} + (1-a_j+1)\frac{\prod_{k=j}^{i-1} (a_k-1)}{\prod_{k=j+1}^{i-1} (a_k-a_i)}=\left(a_j+\frac{(1-a_j)(a_j-a_i)}{a_j-1}\right)v^i_j = a_i v^i_j.$$

So now that you have the eigenvectors and eigenvalues, you can decompose $M = BDB^{-1}$ and compute $M^n = BD^nB^{-1}$.

EDIT: Here is the formula for $B^{-1}$. Write $$B = \left[\begin{array}{cccc} v^1 & v^2 & \ldots & v^n\end{array}\right]$$ where $v^i$ is as above and $v^n$ is the one vector. Then $$B^{-1} = \left[\begin{array}{cccc} w^1 & w^2 & \ldots & w^n\end{array}\right]$$ with $$w^i_j = \begin{cases}0, &j > i\\ \frac{\prod_{k=j}^{i-1} (a_k-1)}{\prod_{k=j+1}^{i} (a_k-a_j)}, & j \leq i\end{cases}$$ for $i<n$ and $$w^n_j = \begin{cases}1, &j = n\\ -\frac{\prod_{k=j+1}^{n-1} (a_k-1)}{\prod_{k=j}^{n-2} (a_k-a_{n-1})}, & j < n.\end{cases}$$

-
I am sorry, I've made a mistake, my diagonal entries are not distinct. But I guess I can slightly change them so they become different. Would it change a lot the result if I add something like 0.0000001 when necessary? – Maxime Sep 20 '13 at 20:18
The entries of $M^n$ are continuous functions of the entries of $M$. Maybe someone else can give more precise bounds on the error introduced by perturbing the equal values. – user7530 Sep 20 '13 at 20:23
Notice that the entries of $B$ and $B^{-1}$ become ill conditioned when two diagonal values approach each other, though. Be careful. – user7530 Sep 20 '13 at 20:25
Yes, I was also afraid of that... – Maxime Sep 20 '13 at 20:28