Eigenvalue Issues of rate matrixes


#1

Dear Stan devolopers and community,
for a validation and inference of the parameters of a Markov model I need to solve an eigenvalue problem of the ratematrix.
Browsing the old google Group.

https://groups.google.com/forum/#!msg/stan-users/QJe1TNioiyg/7VI2hBhRk8oJ.

I found some discussion which did not really solve the issue. Of course there are classes of matrixes which have only real valued Eigenvalues such as rate matrixes of time reversible markov processes. They are similar to symmetric matrix thus they share the same eigenvalues. Therefore they are all real numbers and they are all smaller then 1 etc…

https://books.google.de/books?id=j4ESDAAAQBAJ&pg=PA69&lpg=PA69&dq=Ratematrix+eigenvalues&source=bl&ots=e6GQz4fAh9&sig=chelbZJ4HRmmaKoHHBPksB1qPZs&hl=de&sa=X&ved=0ahUKEwjboqCNipPVAhUE2BoKHTl1A_IQ6AEIQjAE#v=onepage&q=Ratematrix%20eigenvalues&f=false

On page 69 of that link the talk about that. Finding that similar Matrix then requires to transform the ratematrix to its jordan normal form and I guess thats not yet written in the set algrebra functions in Stan.

It would be cool if Stan gets a function for that type of ratematrix. I guess many users would benefit since infering Markov models is a popular task :).

What do you think?

Thanks a lot in advance

Jan


#2

The Eigen library uses complex arithmetic in the asymmetric case even if the eigenvalues happen to be real. If you only need the first eigenvalue / eigenvector, you could do power iteration or something.


#3

Hi Jan,

It’s really really had to compute the Jordan normal form numerically, so that’s not a great route to go down. (To see this, note that the set of diagonalisable matrices is dense in the set of square matrices, so rounding error will make it hard to find non-trivial Jordan blocks even if they should be there)

Can you give the maths of what you need? The link to that book doesn’t work.

There are some tricks you can use (eg, can you write this as a generalised eigenproblem with an SPD pencil), but you need to work out the structure. For non-normal matrices, the eigenvalues aren’t necessarily a very smooth function of the matrix (it comes down to pseudospectra).

Best wishes,

Dan


#4

Could you share a bit more about what you need the eigenvalues for? I saw that your link was to a molecular evolution textbook - if you are working on what I think you are, you may be able to use the matrix exponential function in Stan and avoid calculating the eigenvalues directly.


#5

Hi aaronjg,
I am working with Patchclamb data from membran Proteins (Ionchannels like the HCN channel etc…)
The thing with the Matrix exponential is, that it is much slower to evaluate. For each time I have to compute the Matrix exponential. Thats bad! Though it should work.
Further more there is also the fact that Eigenvalues of the rate matrix have a physical meaning. They quantify the rate at which a disturbence out of the equilirium of the markov System decays in time…

To Daniel_Simpson I am trying to be a bit more precise with the mathematics I need in the next post. I just Need to read some stuff first.

Thanks alot

Jan


#6

If the matrix exponential works, then you can extract the matrix when the sampling is complete and calculate the eigenvalues in R or Python, which should work fine unless you need to put a prior on the eigenvalues within the model.