## 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 00 λ22 0 00 0 λ32 00 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:

 `An35992080.0 11831400.0 37439204.011831400.0 3941981.0 12317651.037439204.0 12317651.0 38949930.0P35992080.0 11831400.0 37439204.011831400.0 3941981.0 12317651.037439204.0 12317651.0 38949930.0L-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.01332Lexp0.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.98677Pexp0.99008 0.00331 0.00331 0.003310.00331 0.99008 0.00331 0.003310.00331 0.00331 0.99008 0.003310.00331 0.00331 0.00331 0.99008method 10.99008 0.00331 0.00331 0.003310.00331 0.99008 0.00331 0.003310.00331 0.00331 0.99008 0.003310.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 20.99008 0.00331 0.00331 0.003310.00331 0.99008 0.00331 0.003310.00331 0.00331 0.99008 0.003310.00331 0.00331 0.00331 0.99008method 30.99008 0.00331 0.00331 0.003310.00331 0.99008 0.00331 0.003310.00331 0.00331 0.99008 0.003310.00331 0.00331 0.00331 0.99008`

Listing:

 `import numpy as npdef 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 matrixA = np.array( [[2, 0, 4], [0, 3, 1], [4, 1, 2]]) evals, evecs = np.linalg.eig(A)U = evecsL = np.diag(evals)Ui = np.linalg.inv(U)# check that U Ui is close to IX = np.abs(np.dot(U,Ui) - np.eye(U.shape[0]))assert X.all() < 1e-6# let's check A^nn = 10def power(A,n): P = np.eye(A.shape[0]) for i in range(n): P = np.dot(P,A) return Pshow(power(A,n),'An')Lp = power(L,n)P = np.dot(np.dot(U,Lp),Ui)show(P,'P')#---------------------------------------# a JC69 matrixu = 0.00333A = 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 = evecsL = np.diag(evals)Ui = np.linalg.inv(U)# exponentiateimport mathshow(L,'L')for i in range(L.shape[0]): value = L[i,i] value = math.e**value L[i,i] = valueshow(L,'Lexp')show(np.dot(np.dot(U,L),Ui),'Pexp')# linalg has three methodsfrom scipy import linalgshow(linalg.expm(A),'method 1')show(linalg.expm2(A),'method 2') # eigenvalue decompositionshow(linalg.expm3(A),'method 3') # Taylor series`