Received 15 April 2005
Rapid calculation of RMSDs using a quaternion-based characteristic polynomial
A common measure of conformational similarity in structural bioinformatics is the minimum root mean square deviation (RMSD) between the coordinates of two macromolecules. In many applications, the rotations relating the structures are not needed. Several common algorithms for calculating RMSDs require the computationally costly procedures of determining either the eigen decomposition or matrix inversion of a or matrix. Using a quaternion-based method, here a simple algorithm is developed that rapidly and stably determines RMSDs by circumventing the decomposition and inversion problems.
Orthogonal rotations are commonly used for comparing macromolecular structures, and the root mean square deviation (RMSD) is a natural metric for quantitating the similarity of two optimally rotated structures (Flower, 1999). Several least-squares algorithms have been proposed for finding the rotation that minimizes the squared distances between corresponding atoms in the structures. The most efficient of these methods requires the eigen decomposition (or inversion) of a key matrix constructed from the sums and squares of the atomic coordinates in the structures being compared. The popular method of Kabsch (1978) involves a key matrix, whereas the methods of Diamond (1988), Horn (1987) and Kearsley (1989) use quaternion representations of the rotations with key matrices. In each of these methods, the minimizing rotation is formed from the eigenvectors of the key matrices and the minimizing RMSD is found from the eigenvalues of the key matrices.
In Kabsch's method, improper rotation matrices (rotoinversions) may arise and one must verify that the determinant of the rotation matrix is positive (Kabsch, 1978). Quaternion-based methods lack this complication since quaternion rotations are always proper. Thus, with quaternion-based methods the minimum RMSD can be calculated using the eigenvalues of the key matrix alone, without knowing the eigenvectors. If only RMSDs are desired, the quaternion-based methods can save a considerable amount of computation by avoiding a complete eigen decomposition of the key matrix.
Diamond (1988) proposed the fastest method known for determining minimum RMSDs using an iterative algorithm. Diamond's method requires one matrix inversion per cycle in order to find the largest eigenvalue of a quaternion-based matrix. While matrix inversions are generally expensive, the inversion of a matrix can be calculated analytically and relatively quickly. Diamond's method, however, is unstable when the minimizing rotation is near 180°. This problem is especially grave when the coordinate differences are of magnitude similar to the coordinates themselves, such as is common when superimposing relatively small fragments of proteins or nucleic acids.
Based on Horn's quaternion superposition method (Horn, 1987), which is relatively unknown in structural biology, an extremely fast, numerically stable and accurate algorithm is presented here for determining minimum RMSDs. This algorithm does not require separate consideration of special problematic cases and it bypasses the computationally costly diagonalization and inversions commonly used to find the needed eigenvalues. Rather, a quick and simple Newton-Raphson root-finding method is used to determine the appropriate eigenvalue from the characteristic polynomial of the key matrix. This quaternion-based characteristic polynomial (QCP) algorithm requires significantly less computation than the alternate methods, including Diamond's.
A molecular structure consisting of N atoms can be represented as an N×3 matrix, where the three columns correspond to the x, y and z coordinates. The superposition objective is to find the orthogonal rotation and translation that minimizes the squared Euclidean distance between the rows of two matrices corresponding to the two structures being compared. The translational component of the problem can be removed from the outset by translating each structure so that its respective centroid is at the origin. When superimposing structure onto structure , in matrix notation the problem is to minimize the sum of squared errors E with respect to the orthogonal rotation :
where denotes the Frobenius (or Euclidean) norm of the matrix ,
where GA is the inner product of structure ,
and M is the inner product of the matrices A and B. The inner product matrix M is given by
The eigenvector corresponding to the largest eigenvalue of K is a quaternion equivalent to the optimal rotation, where
Eigenvalues are usually determined by diagonalization of the matrix. Alternatively, one can find the eigenvalues of a symmetric matrix by locating the roots of the characteristic polynomial of the matrix. With denoting the determinant of matrix , the characteristic equation for the key matrix is a quartic given by
Using these coefficients, the characteristic equation can be reduced to the following simple form:
Although the coefficient equations may appear somewhat horrendous, coding them is straightforward. With careful programming, by avoiding duplicate evaluations of several terms (such as Syz - Szy which appears four times), calculating the coefficients from the sums found in equation (4) involves at most only 66 floating point operations (FLOPs).
Horn (1987) proposed solving this characteristic equation using Ferrari's method (Abramowitz & Stegun, 1972) for quartics. While analytically exact, Ferrari's method has seen little use as it is computationally slow and notoriously numerically unstable. Additionally, Horn (1987) omits expressions for C0 and an equation given for C2 contains an important typographical error.
Fortunately in this case, the largest positive root of the characteristic polynomial given in equation (8) can be found very quickly and with excellent numerical stability using the Newton-Raphson algorithm. Successful use of the Newton-Raphson method requires a good initial guess for the desired root. From equation (3) and noting that the largest positive eigenvalue of the key matrix K is equal to with the optimal rotation, the largest positive root of the characteristic polynomial attains its maximum possible value when the two structures are in identical conformations (i.e. when RMSD = 0). Thus, the maximum possible value of the largest positive root is constrained to be less than or equal to the average of the inner products of the the two structures, (GA + GB) / 2. Using this value as an initial guess in the Newton-Raphson algorithm is essentially fool-proof, since the algorithm will converge first on the largest root of the polynomial. The Newton-Raphson QCP algorithm can be described by the following simple pseudocode:
where the first derivative of the polynomial is given by .
The Newton-Raphson QCP algorithm given above is much faster than other common methods that require determination of eigenvectors, and it is also faster than methods using eigen decomposition for only eigenvalues. The efficiency and speed of the QCP Newton-Raphson algorithm were compared to the popular method using Householder reduction to tridiagonal form followed by QL decomposition with implicit shifts (H-QL method; Golub & Van Loan, 1996). Each iteration of the QCP algorithm requires merely 11 FLOPS, of which none are square roots. In contrast, the H-QL algorithm requires a comparable number of total FLOPs but over two dozen costly square roots. Diamond's (1988) method is likewise much more efficient than the eigen decomposition methods and also requires no square roots. Each cycle of the Diamond iteration uses FLOPS (accounting for symmetry considerations), plus 15 FLOPS for set-up. Since both the QCP and Diamond's methods take an average of about five iterations to converge to a relative precision of 10-6, the QCP method uses only about one third of the operations of Diamond's method for eigenvalue determination. Thus, essentially all of the computational effort in the QCP method is expended in calculating the summations of the key matrix given in equation (4), a step that is common to all methods and dependent on the number of atoms compared.
The QCP algorithm has been benchmarked on several computing platforms (OSX, Linux and Irix) for independent superposition calculations with protein fragments ranging from 5 to 500 residues. In all cases, the same RMSD solution was found, within error, as the Kabsch and Kearsley methods. In contrast, especially with smaller fragments, Diamond's method converged on an incorrect solution (corresponding to the second largest eigenvalue) in a significant fraction of cases. In the author's experience, if the construction of the key matrices is ignored, the QCP algorithm is over four times as fast as Diamond's method and 30 to 70 times faster than the eigen decomposition methods. An ANSI C library implementing the QCP method is available at http://monkshood.colorado.edu/QCP/ or by request from the author.
I thank B. K. P. Horn for helpful discussion, an anonymous referee for alerting me to Diamond's scalar method and the American Cancer Society for funding.
Abramowitz, M. & Stegun, I. A. (1972). Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th ed., pp. 17-18. New York: Dover.
Diamond, R. (1988). Acta Cryst., A44, 211-216.
Flower, D. R. (1999). J. Mol. Graph. Model. 17, 238-244.
Golub, G. H. & Van Loan, C. F. (1996). Matrix Computations, 3rd ed. Baltimore: Johns Hopkins University Press.
Horn, B. K. P. (1987). J. Opt. Soc. Am. A Opt. Image Sci. Vis. 4, 629-642.
Kabsch, W. (1978). Acta Cryst. A34, 827-828.
Kearsley, S. K. (1989). Acta Cryst. A45, 208-210.