Kabsch Algorithm

作者：Eric

创建时间：2026-04-16 11:06

这篇记什么

• 量化同一序列多肽（蛋白质）相似性，使用RMSD+Kabsch算法
• 进行Kabsch算法的推导及揭示其对应的数学原理
• 给出算法对应的python实现

主要内容

问题提出（刚体配准问题）

设有两组组成完全相同对应点（pdb文件中的ATOM行）：

$$p_i \in P,\qquad q_i \in Q,\qquad i=1,2,\dots,N$$

在经过配准后

$$P_{fit} = PR + T$$

希望求旋转矩阵 $R$ 和平移向量 $T$ ，使得目标函数

$$\min_{R,T}\sum_{i=1}^N \|p_iR+T-q_i\|^2$$

达到最小，其中 $R$ 满足

$$R^TR=I$$

在蛋白质结构中，一般还需令 $det(R)=1$ 来保证该蛋白质不存在镜面对称构象。

求解过程

1. 消除平移向量 $T$

   定义目标函数

   $$F(R,T)=\sum_{i=1}^N \|p_iR+T-q_i\|^2$$

   对 $T$ 求偏导，并令其为零：

   $$\frac{\partial F}{\partial T}=2\sum_{i=1}^N (p_iR+T-q_i)=0$$

   因此

   $$NT=\sum_{i=1}^N (q_i-p_iR)$$

   从而

   $$T=\frac{1}{N}\sum_{i=1}^N q_i-\frac{1}{N}\sum_{i=1}^N p_iR$$

   记两组点的质心分别为

   $$\bar p=\frac{1}{N}\sum_{i=1}^N p_i\qquad \bar q=\frac{1}{N}\sum_{i=1}^N q_i$$

   则平移项可改写为

   $$T=\bar q-\bar pR$$

   将 $T$ 代入目标函数，得

   $$F(R)=\sum_{i=1}^N \|p_iR+\bar q-\bar pR-q_i\|^2$$

   整理后可得

   $$F(R)=\sum_{i=1}^N \|(p_i-\bar p)R-(q_i-\bar q)\|^2$$

2. 优化目标函数 $F(R)$

   定义去中心化后的坐标

   $$x_i=p_i-\bar p,\qquad y_i=q_i-\bar q$$

   则目标函数可以改写为

   $$F(R)=\sum_{i=1}^N \|x_iR-y_i\|^2$$

   将所有 $x_i$ 和 $y_i$ 按行堆叠成矩阵：

   $$X=\begin{bmatrix} x_1\\ x_2\\ \vdots\\ x_N \end{bmatrix}, \qquad Y=\begin{bmatrix} y_1\\ y_2\\ \vdots\\ y_N \end{bmatrix}$$

   则有

   $$F(R)=\|XR-Y\|_F^2$$

   由 $F范数$ 与迹的关系（[[F范数与迹]]）可得

   $$F(R)=\|XR-Y\|_F^2=\operatorname{tr}\big((XR-Y)^T(XR-Y)\big)$$

   展开后

   $$F(R)=\operatorname{tr}(R^TX^TXR)-\operatorname{tr}(R^TX^TY)-\operatorname{tr}(Y^TXR)+\operatorname{tr}(Y^TY)$$

   又由于 $R^TR=I$ 与迹的循环不变性

   $$\operatorname{tr}(Y^TXR)=\operatorname{tr}(R^TX^TY)$$

   所以

   $$F(R)=\operatorname{tr}(X^TX)-2\operatorname{tr}(R^TX^TY)+\operatorname{tr}(Y^TY)$$

   其中 $\operatorname{tr}(X^TX)$ 和 $\operatorname{tr}(Y^TY)$ 都与 $R$ 无关，因此

   $$\boxed{\text{最小化} F(R) \text{等价于最大化} \operatorname{tr}(R^TX^TY)}.$$

   记 $C=X^TY$，则问题转化为

   $$\max_R \operatorname{tr}(R^TC),\qquad \text{s.t. } R^TR=I$$

3. 求解最优旋转矩阵 $R^\ast$ 和最优平移 $T^\ast$

   对矩阵 $C$ 做 $SVD$ 分解

   $$C=U\Sigma V^T$$

   其中

   $$\Sigma=\operatorname{diag}(\sigma_1,\sigma_2,\dots,\sigma_d),\qquad \sigma_i\ge 0$$

   则

   $$\operatorname{tr}(R^TC)=\operatorname{tr}(R^TU\Sigma V^T)$$

   利用迹的循环不变性

   $$\operatorname{tr}(R^TU\Sigma V^T)=\operatorname{tr}(V^TR^TU\Sigma)$$

   记 $M = V^TR^TU$ ，且由于 $R,U,V$ 都是正交矩阵，$M$ 也为正交矩阵。
   故

   $$\operatorname{tr}(R^TC)=\operatorname{tr}(M\Sigma)$$

   在三维情景下

   $$\Sigma=\operatorname{diag}(\sigma_1,\sigma_2,\sigma_3)$$

   则

   $$\operatorname{tr}(M\Sigma)=m_{11}\sigma_1+m_{22}\sigma_2+m_{33}\sigma_3$$

   由于

   $$|m_{ii}| \le 1, \qquad \sigma_i \ge 0, \qquad i=1,2,3$$

   因此

   $$m_{11}\sigma_1+m_{22}\sigma_2+m_{33}\sigma_3 \le \sigma_1 + \sigma_2 + \sigma_3$$

   当且仅当 $m_{11}=m_{22}=m_{33}=1$ 时，等号成立
   此时

   $$M=I$$

   于是有

   $$V^TR^TU=I$$

   从而得到最优旋转矩阵

   $$R^\ast=UV^T$$

   代回 $T$ 所对应的方程，得

   $$T^\ast=\bar q-\bar pR^\ast$$

4. 基于现实情况对 $R^\ast$ 修正

   上式得到的 $R^\ast$ 虽然满足正交性，但对于蛋白质来说，其可能不是真正的旋转矩阵。 由于

   $$det(R^\ast) = det(U)det(V^T)$$

   而正交矩阵的行列式只能取 $\pm 1$ ，因此可能出现

   $$det(R^\ast)=-1$$

   此时意味着在这个变换中发生了 镜面反射 ，即蛋白质经过变换后发生了手性改变。

   镜面反射指在三维空间的 $x , y, z$ 三轴中，有一条轴被翻转，此时 $R^\ast$ 代表的不只是一个单纯的旋转矩阵。

   因此，需要对 $det(R^\ast)$ 的结果进行修正。

   修正：在三个奇异方向中，找到翻转代价最小的那一条轴，再进行一次符号翻转

   由于奇异值满足

   $$\sigma_1\ge \sigma_2\ge \sigma_3\ge 0$$

   故选择对应最小奇异值 $\sigma_3$ 的方向进行进行修正。 引入

   $$D=\operatorname{diag}(1,1,-1)$$

   综合来看

   $$D = \begin{cases} \operatorname{diag}(1,1,1), & \det(VU^T)=1， &（无需修正）,\\[6pt] \operatorname{diag}(1,1,-1), & \det(VU^T)=-1， &（需要修正）. \end{cases}$$

   更一般的

   $$D=\operatorname{diag}\bigl(1,1,\det(VU^T)\bigr)$$

   于是修正后的 $R^\ast$ 为

   $$R^\ast=VDU^T$$

结果汇总

经修正后原问题（刚体配准问题

$$\min_{R,T}\sum_{i=1}^N \|p_iR+T-q_i\|^2$$

的解为

$$R^\ast=VDU^T,\qquad T^\ast=\bar q-\bar pR^\ast$$

其中

$$D=\operatorname{diag}\bigl(1,1,\det(VU^T)\bigr)$$

且

$$X=\begin{bmatrix} p_1-\bar p\\ p_2-\bar p\\ \vdots\\ p_N-\bar p \end{bmatrix}, \qquad Y=\begin{bmatrix} q_1-\bar q\\ q_2-\bar q\\ \vdots\\ q_N-\bar q \end{bmatrix}, \qquad C=X^TY=U\Sigma V^T$$

代码实现

import numpy as np
def Kabsch(P, Q):
	P = np.asarray(P, dtype = float)
	Q = np.asarray(Q, dtype = float)
	
	# Calc X & Y & C
	p_cent = P.mean(axis=0)
	q_cent = Q.mean(axis=0)
	X = P - p_cent
	Y = Q - p_cent
	C = X.T @ Y
	
	# SVD
	U, Sig, Vt = np.linalg.svd(C)
	V = Vt.T
	
	# Add correction D (3-Dim)
	D = np.diag([1, 1, np.sign(np.linalg.det(V @ U.T))])
	
	# Calc R* & T*
	R_star = V @ D @ U.T
	T_star = q_cent - p_cent @ R_star
	
	return R_star, T_star

若是想要返回配准后的RMSD（[[RMSD]]）的计算结果

def RMSD_Kabsch(P, Q):
	...
	P_fit = P @ R_star + T_star
	rmsd = np.sqrt(np.mean(np.sum((P_fit - Q) ** 2, axis=1)))
	return rmsd

相关

• 点云配准算法(kabsch) - silvergp