Models of sequence evolution (2): matrix exponentiation

In recent days I've posted a number of times on this topic, most of them under the title of Jukes-Cantor (XX):

• A script to simulate sequence evolution (post)
• A peek at the answers for P_XX(t) and P_XY(t) under Jukes-Cantor (post)
• A simple derivation from Higgs & Attwood (post)
• Use to get a relationship between observed changes and distance (post)
• A script that sets up a matrix form of P(t) and cranks it a bit showing the convergence to the stationary distribution (post)
• Going from the solution to the differential equations of JC69 (post)
• A first mention of the equation P(t) = e^Q(t) from Yang's book. (post)
• A test of matrix decomposition into eigenvalues and eigenvectors (post)

I think I've made some progress in understanding models of sequence evolution (though how well I've conveyed this, you'll be the judge). But there is still a way to go yet. I need to understand: (i) why this exponentiation thing seems to be important, (ii) use of the model for likelihood calculations (L ≈ P(data | model), (iii) extension to more realistic models like Kimura's 2-parameter model that treats transitions and transversions differently, and (iv) how to do practical calculations under the models.

Today, I want to do just a bit more with matrix exponentiation. Before, we talked about the idea that you can decompose a matrix (at least a square well-behaved matrix like ours) into something like this A = U Λ U-1 where, critically, the Λ matrix is a diagonal matrix so that for example, its square is:

&Lambda2 = λ12 0 0 0 0 λ22 0 0 0 0 λ32 0 0 0 0 λ42

This is supposed to be extensible to exponentiation (or any function), and I was having trouble wrapping my head around that. But then I got reminded that e^x can be written as:

ex = Σ xi / i! where the sum runs from i=0 to i=∞

So, if we can do A2, we can do eA.

I'm at my office computer this AM, and to do what follows I needed SciPy installed. So I built and installed Python 2.6, used an available install of the gfortran compiler and FFTW, ignored the sparse matrix stuff, built and installed Numpy and SciPy. All as described here, here and here.

The script (full listing at the end of the post), has two parts. First, we do the decomposition, following the same line as previously, checking that Aⁿ by the eigenvalue method matches what we get for repeated matrix multiplication.

In the second part, I do exponentiation using the eigenvalue method, and check it using three methods from SciPy (reference).

It seems to work fine!

Output:

An 35992080.0 11831400.0 37439204.0 11831400.0 3941981.0 12317651.0 37439204.0 12317651.0 38949930.0 P 35992080.0 11831400.0 37439204.0 11831400.0 3941981.0 12317651.0 37439204.0 12317651.0 38949930.0 L -0.01332 0.0 0.0 0.0 0.0 -0.0 0.0 0.0 0.0 0.0 -0.01332 0.0 0.0 0.0 0.0 -0.01332 Lexp 0.98677 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.98677 0.0 0.0 0.0 0.0 0.98677 Pexp 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 method 1 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 /Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/scipy/linalg/matfuncs.py:96: ComplexWarning: Casting complex values to real discards the imaginary part return dot(dot(vr,diag(exp(s))),vri).astype(t) method 2 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 method 3 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008 0.00331 0.00331 0.00331 0.00331 0.99008

Listing:

import numpy as np def show(A,name=None,x=5): if name: print name for row in A: L = [str(round(n,x)).rjust(7) for n in row] print ' '.join(L) print # a simple 3x3 symmetric matrix A = np.array( [[2, 0, 4], [0, 3, 1], [4, 1, 2]]) evals, evecs = np.linalg.eig(A) U = evecs L = np.diag(evals) Ui = np.linalg.inv(U) # check that U Ui is close to I X = np.abs(np.dot(U,Ui) - np.eye(U.shape[0])) assert X.all() < 1e-6 # let's check A^n n = 10 def power(A,n): P = np.eye(A.shape[0]) for i in range(n): P = np.dot(P,A) return P show(power(A,n),'An') Lp = power(L,n) P = np.dot(np.dot(U,Lp),Ui) show(P,'P') #--------------------------------------- # a JC69 matrix u = 0.00333 A = np.array( [[-3*u, u, u, u], [u, -3*u, u, u], [u, u, -3*u, u], [u, u, u, -3*u]]) evals, evecs = np.linalg.eig(A) U = evecs L = np.diag(evals) Ui = np.linalg.inv(U) # exponentiate import math show(L,'L') for i in range(L.shape[0]): value = L[i,i] value = math.e**value L[i,i] = value show(L,'Lexp') show(np.dot(np.dot(U,L),Ui),'Pexp') # linalg has three methods from scipy import linalg show(linalg.expm(A),'method 1') show(linalg.expm2(A),'method 2') # eigenvalue decomposition show(linalg.expm3(A),'method 3') # Taylor series