Diagonalization

This applet shows graphically the rotation and scaling operations involved in diagonalizing a matrix. Scroll below the applet to see more explanation.

Background

Imagine a sequence comprising 50 sites in which 30 sites (60%) are in state 0 and the remaining 20 sites (40%) are in state 1.

00101101101000010110000010100001010000001101110110

The fraction of sites in each of the two possible states can be represented by the (transposed) vector (0.6, 0.4). If this sequence were allowed to evolve for an infinite amount of time under an equal-rate (Cavender-Farris) model, the sequence composition would achieve an equilibrium in which half the sites were in state 0 and the other half in state 1. This ewquilibrium point is represented by the red dot on the simplex in the applet.

If, however, the amount of time and the substitution rate were such that the expected number of substitutions is 0.5, the sequence would evolve toward the equilibrium, but how close would it come to the red dot? This applet illustrates how to predict the distance a sequence would evolve toward the equilibrium point.

All possible vectors representing relative state frequencies must lie on the simplex represented in the applet by the dotted gray line. At one extreme (1.0 on the horizontal axis), the sequence comprises only state 0. At the other extreme (1.0 on the vertical axis), the sequence consists solely of sites in state 1. The dropdown list at the top allows you to select the length of the evolutionary path that the sequence will follow: higher numbers mean that the sequence will travel further toward the red equilibrium point because the substitution rate and/or time available is larger.

The transition probability matrix P = exp{Qt}, where t is time (measured in expected number of substitutions per site) and Q is the instantaneous rate matrix (scaled so that it represents one expected substitution per unit time. The exponentiation can be accomplished by first diagonalizing Q, which means Q is factored into an eigenvector matrix, a diagonal matrix of eigenvalues, and an inverse eigenvector matrix. The eigenvalues are then multiplied by t and exponentiated to yield the eigenvalues of the P matrix. The three 2x2 matrices shown below the graph in the applet represent the diagonalized version of the matrix P.

The applet shows that premultiplying a vector of relative state frequencies by the diagonalized P matrix transformed the starting relative frequency vector to the vector expected after the sequence has evolved for the evolutionary distance specified by the dropdown at the top. Premultipying by the inverse eigenvalue matrix rotates axes to align with the eigenvectors. Premultiplying by the diagonal matrix of eigenvalues scales the axes. Finally premultiplying by the eigenvector matrix rotates axes back to the original coordinate system. Note that all movement of the relative frequency vector is in the direction of the eigenvalues, with the eigenvectors determining the scale of the movement. Because one eigenvalue is 1, no movement occurs on one eigenvalue axis (multiplying by 1 does not change the value along that axis), so all movement is restricted to the line representing the simplex, which is parallel to the other eigenvector.

Acknowledgements

This applet makes use of the excellent d3js javascript library. I am also endebted to Grant Sanderson and his excellent 3Blue1Brown web site, where a much better explanation of eigenvectors and eigenvalues can be found in the Linear Algebra section, and to Mark Holder for showing me this amazing web site. Please see the GitHub site for details about licensing of other libraries that may have been used in the source code for this applet.

Licence

Creative Commons Attribution 4.0 International. License (CC BY 4.0). To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ or send a letter to Creative Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.