Wednesday, March 3, 2010

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 PXX(t) and PXY(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) = eQ(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 ex 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 An 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