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 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]$:

Values of Q^N_i,0

Edit4:

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

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

1 Answer 1

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

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

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.