How is RMSD calculated in the Scipy implementation of the Kabsch algorithm?

user_4391249 0
Viewed 82

Scipy calculates the rmsd like this, and I'll paraphrase it here for convenience (for readability I drop the weights and the max(*, 0))

rmsd = np.sqrt(np.sum(b ** 2 + a ** 2) - 2 * np.sum(s))

To me this does not look like RMSD.

Now from the docs one would infer that the rmsd return value is defined as the square root of double this expression:

enter image description here

The latter is indeed what I would consider to be the RMSD. In fact I went ahead and coded it up (note that this function expects me to apply the estimated transformation to one of the sets of points first whereas the snippet above does not)

def _calc_rmsd(a: np.ndarray, b_transformed: np.ndarray) -> float:
    distances = np.linalg.norm(a - b_transformed, axis=-1)
    rmsd = np.sqrt((distances ** 2).sum() / len(distances))
    return rmsd

I also plotted out what these would look like for randomly generated point pairs with normally distributed noise (blue is scipy, orange is mine)

enter image description here

Or extending the plot out to 200 point pairs:

enter image description here

So to sum it up:

  • The definition of rmsd in the docs is in agreement with what I believe to be the widely accepted notion of rmsd
  • The scipy code implementation of rmsd disagrees with the latter. I don't even understand what it's supposed to mathematically represent.
  • From monte-carlo simulations, clearly the two implementations have different outcomes.

So what's going on?

scipy
1 Answers

Apparently the SciPy code is not returning the root-mean-squared distance. It sums the squared differences, but it does not divide by the number of vectors before taking the square root. The difference between the SciPy calculation and yours is a factor of sqrt(len(a)). You can verify this with an example such as the following.

In [157]: from scipy.spatial.transform import Rotation

In [158]: def _calc_rmsd(a: np.ndarray, b_transformed: np.ndarray) -> float:
     ...:     distances = np.linalg.norm(a - b_transformed, axis=-1)
     ...:     rmsd = np.sqrt((distances ** 2).sum() / len(distances))
     ...:     return rmsd
     ...:

Some test data:

In [159]: a = np.array([[0, 1, 1], [1, 1, 1.5], [2.0, -1.0, 4.0], [-1, 0, 5]])

In [160]: b = np.array([[0, 1, 1.5], [2, 2, 2], [1, -1, 5], [-3, 0.1, 1]])

Compute the rotation:

In [161]: R, rmsd = Rotation.align_vectors(a, b)

In [162]: rmsd
Out[162]: 3.8753534834716685

Here's your calculation of the RMSD:

In [163]: _calc_rmsd(a, R.apply(b))
Out[163]: 1.9376767417358356

And here is your calculation, multiplied by sqrt(len(a)), so it matches the result returned by Rotation.align_vectors:

In [164]: _calc_rmsd(a, R.apply(b)) * np.sqrt(len(a))
Out[164]: 3.875353483471671

This looks like a documentation issue. If you have a moment, you could create a new issue for this over in https://github.com/scipy/scipy/issues

user_12173580
Related