Kabsch 算法的 Scipy 实现中是如何计算 RMSD 的？

Question

Scipy 像 this 一样计算 rmsd，为了方便起见，我将在这里解释一下（为了便于阅读，我删除了 weights 和 max(*, 0)）

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

对我来说这不看起来像 RMSD。

现在从 docs 可以推断出 rmsd return 值被定义为这个表达式的两倍的平方根：

后者确实是我认为的 RMSD。事实上，我继续编写代码（请注意，此函数希望我首先将估计的转换应用于一组点，而上面的代码片段没有）

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

我还绘制了随机生成的具有正态分布噪声的点对的样子（蓝色是 scipy，橙色是我的）

或将绘图扩展到 200 个点对：

总结一下：

文档中 rmsd 的定义与我认为被广泛接受的 rmsd 概念一致
rmsd的scipy代码实现与后者不一致。我什至不明白它应该用数学表示什么。
从蒙特卡洛模拟来看，很明显这两种实现有不同的结果。

所以这是怎么回事？

Answer 1

显然 SciPy 代码没有返回根-均值-平方距离。它对平方差求和，但在求平方根之前不除以向量数。 SciPy 计算结果与您的计算结果相差 sqrt(len(a)) 倍。您可以通过以下示例验证这一点。

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
     ...:

部分测试数据：

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]])

计算旋转：

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

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

这是您计算的 RMSD：

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

这是你的计算结果，乘以 sqrt(len(a))，所以它与 Rotation.align_vectors 返回的结果相匹配：

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

这看起来像是文档问题。如果您有时间，可以在 https://github.com/scipy/scipy/issues

中为此创建一个新问题

Kabsch 算法的 Scipy 实现中是如何计算 RMSD 的？

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

scipy