正交 Procrustes 問題

本文的閱讀等級：高級

普洛克路斯忒斯 (Procrustes) 是希臘神話中海神波塞頓 (Poseidon) 的兒子。他在雅典到埃萊夫西納 (Eleusis) 的神聖之路 (The Sacred Way) 上開設一間黑店，向路過的旅人謊稱店內設有一張適合所有人的鐵床。旅客投宿時，普洛克路斯忒斯將身高者截斷雙足，身矮者則強行拉長，使之與床的長短相同。從來沒有一個人的身長與鐵床的長度相同而免於凌遲，因為他暗地裡準備了兩張床^[1]。後人於是以 Procrustean 表示「削足適履，殺頭便冠」，意思是將不同的長度、大小或屬性安裝到一個任意的標準。

正交 Procrustes 問題：給定兩個 $m\times n$ 階實矩陣 $A$ 和 $B$ ，求一個 $n\times n$ 階實正交矩陣 $Q$ ， $Q^T=Q^{-1}$ ，使得 $\Vert B-AQ\Vert_F$ 具有最小值^[2]，其中 $\Vert\cdot\Vert_F$ 是 Frobenius 範數。

正交 Procrustes 問題是一個最小平方矩陣近似問題，可以這麼解讀： $A$ 是旅人， $B$ 是鐵床，正交 Procrustean 變換 (包含旋轉和鏡射) $Q$ 即為施予旅人的肢體酷刑。我們的問題是求出一酷刑使旅人變形後的身長與鐵床的長度最為吻合。以下討論限定於實矩陣。

Procrustes 和他的客人(點圖看動畫) From http://www.mythweb.com/today/media/procrustes07.gif

對於一個 $m\times n$ 階矩陣 $A=[a_{ij}]$ ，Frobenius 範數定義為 (見“矩陣範數”)

$\displaystyle \Vert A\Vert_{F}=\sqrt{\mathrm{trace}(A^TA)}=\sqrt{\sum_{i=1}^m\sum_{j=1}^n a_{ij}^2}$ 。

使用跡數循環不變性和轉置不變性 (見“跡數的性質與應用”)，

$\displaystyle\begin{aligned} \Vert B-AQ\Vert_F^2&=\hbox{trace}\left((B-AQ)^T(B-AQ)\right)\\ &=\hbox{trace}(B^TB-B^TAQ-Q^TA^TB+Q^TA^TAQ)\\ &=\hbox{trace}(B^TB)-\hbox{trace}(B^TAQ)-\hbox{trace}(Q^TA^TB)+\hbox{trace}(Q^TA^TAQ)\\ &=\hbox{trace}(B^TB)-\hbox{trace}(QB^TA)-\hbox{trace}(A^TBQ^T)+\hbox{trace}(A^TAQQ^T)\\ &=\hbox{trace}(B^TB)-\hbox{trace}(QB^TA)^T-\hbox{trace}(A^TBQ^T)+\hbox{trace}(A^TA)\\ &=\Vert B\Vert_F^2-2\hbox{trace}(A^TBQ^T)+\Vert A\Vert_F^2. \end{aligned}$

因為 $A$ 和 $B$ 是給定的矩陣，最小化 $\Vert B-AQ\Vert_F$ 等價於最大化 $\hbox{trace}(A^TBQ^T)$ 。另外，設 $B=0$ ，上式表明 $\Vert AQ\Vert_F=\Vert A\Vert_F$ ，正交矩陣乘法不改變 Frobenius 範數。

正交 procrustes 問題可轉換為底下的約束最佳化問題：

$\displaystyle \max_{QQ^T=I}\hbox{trace}(A^TBQ^T)$

使用 Lagrange 乘數法 (見“Lagrange 乘數法”)。令 Lagrangian 函數為

$\displaystyle L=\hbox{trace}(A^TBQ^T)-\hbox{trace}(\Lambda(QQ^T-I))$ ，

其中 $\Lambda$ 是 $n\times n$ 階 Lagrange 乘數矩陣。利用跡數的導數公式 (見“跡數與行列式的導數”，tr-6，tr-13)，計算 $L$ 對於 $Q$ 的導數：

$\displaystyle\begin{aligned} \frac{\partial L}{\partial Q}&=\frac{\partial\hbox{trace}(A^TBQ^T)}{\partial Q}-\frac{\partial\hbox{trace}(Q^T\Lambda Q)}{\partial Q}+\frac{\partial\hbox{trace}(\Lambda)}{\partial Q}\\ &=A^TB-(\Lambda+\Lambda^T)Q. \end{aligned}$

設上式為零，可得 $A^TB=(\Lambda+\Lambda^T)Q$ ，或 $A^TBQ^T=\Lambda+\Lambda^T$ 。設 $n\times n$ 階矩陣 $A^TB$ 的奇異值分解為 $A^TB=U\Sigma V^T$ ，其中 $U$ 和 $V$ 是 $n\times n$ 階正交矩陣， $\Sigma=\hbox{diag}(\sigma_1(A^TB),\ldots,\sigma_n(A^TB))$ ，這裡 $\sigma_i(A^TB)\ge 0$ 是 $A^TB$ 的奇異值 (見“奇異值分解 (SVD)”)。因為 $\Lambda+\Lambda^T$ 是對稱矩陣，可知 $A^TBQ^T=Q(A^TB)^T$ 。代入 $A^TB$ 的奇異值分解，

$\displaystyle U\Sigma V^TQ^T=QV\Sigma^TU^T=QV\Sigma U^T$ 。

比較等號兩邊可得 $QV=U$ ，故 $Q=UV^T$ 為一個最佳解，而且

$\displaystyle\begin{aligned} \hbox{trace}(A^TBQ^T)&=\hbox{trace}(U\Sigma V^TVU^T)=\hbox{trace}(U\Sigma U^T)\\ &=\hbox{trace}(\Sigma UU^T)=\hbox{trace}\Sigma=\sum_{i=1}^n\sigma_i(A^TB).\end{aligned}$

矩陣 $B$ 的最小平方近似矩陣是 $AQ=AUV^T$ ，最小誤差為

$\displaystyle \Vert B-AQ\Vert_F=\Vert B-AUV^T\Vert_F=\sqrt{\Vert A\Vert_F^2+\Vert B\Vert_F^2-2\sum_{i=1}^n\sigma_i(A^TB)}$ 。

上式提供了 $B=AQ$ 的一個充要條件。若 $B=AQ$ ，則 $\Vert B\Vert_F=\Vert AQ\Vert_F=\Vert A\Vert_F$ ，就有 $\Vert A\Vert_F^2=\Vert B\Vert_F^2=\sum_{i=1}^n\sigma_i(A^TB)$ 。反向的推論十分明顯。

補充一個較快捷的推導方式。如果不採行 Lagrange 乘數法，直接將 $A^TB=U\Sigma V^T$ 代入目標函數，如下：

$\hbox{trace}(A^TBQ^T)=\hbox{trace}(U\Sigma V^TQ^T)=\hbox{trace}(\Sigma V^TQ^TU)=\hbox{trace}(\Sigma P)$ ，

上面令 $P=[p_{ij}]=V^TQ^TU$ 。因為 $U, V, Q$ 都是正交矩陣，不難證明 $P^T=P^{-1}$ ， $P$ 也是一個正交矩陣。明確寫出

$\displaystyle \hbox{trace}(\Sigma P)=\sum_{i=1}^n\sigma_i(A^TB)p_{ii}$ ，

但 $\sum_{i=1}^np^2_{ij}=1$ 且 $\sum_{j=1}^np^2_{ij}=1$ ，上式的最大值必發生於 $p_{ii}=1$ ， $i=1,\ldots,n$ ，即 $P=I$ 。因此， $V^TQ^TU=I$ ，解得 $Q=UV^T$ 。

舉一個例子^[3]，若

$\displaystyle A=\begin{bmatrix} 1.2&2.1\\ 2.9&4.3\\ 5.2&6.1\\ 6.8&8.1 \end{bmatrix},~~ B=\begin{bmatrix} 1&2\\ 3&4\\ 5&6\\ 7&8 \end{bmatrix}$ ，

則正交矩陣

$\displaystyle Q=\left[\!\!\begin{array}{cr} 0.9999 &-0.0126\\ 0.0126&0.9999 \end{array}\!\!\right]$

最小化 $\Vert B-AQ\Vert_F$ ，最小誤差為 $0.4661$ 。

接下來我們探討正交 Procrustes 問題的一個變形問題，稱為兩面正交 Procrustes 問題。

兩面正交 Procrustes 問題：給定兩個 $m\times n$ 階實矩陣 $A$ 和 $B$ ，求一個 $n\times n$ 階實正交矩陣 $Q$ 和一個 $m\times m$ 階實正交矩陣 $P$ ，使得 $\Vert B-PAQ\Vert_F$ 具有最小值。

仿造前面的作法，計算目標函數

$\displaystyle\begin{aligned} \Vert B-PAQ\Vert_F^2&=\hbox{trace}\left((B-PAQ)^T(B-PAQ)\right)\\ &=\hbox{trace}(B^TB-B^TPAQ-Q^TA^TP^TB+Q^TA^TP^TPAQ)\\ &=\hbox{trace}(B^TB)-\hbox{trace}(B^TPAQ)-\hbox{trace}(Q^TA^TP^TB)+\hbox{trace}(Q^TA^TAQ)\\ &=\hbox{trace}(B^TB)-\hbox{trace}(QB^TPA)-\hbox{trace}(A^TP^TBQ^T)+\hbox{trace}(A^TAQQ^T)\\ &=\hbox{trace}(B^TB)-\hbox{trace}(QB^TPA)^T-\hbox{trace}(A^TP^TBQ^T)+\hbox{trace}(A^TA)\\ &=\Vert B\Vert_F^2-2\hbox{trace}(A^TP^TBQ^T)+\Vert A\Vert_F^2, \end{aligned}$

所以最小化 $\Vert B-PAQ\Vert_F$ 等價於最大化 $\hbox{trace}(A^TP^TBQ^T)$ 。

設 $A$ 和 $B$ 的奇異值分解為 $A=U_1\Sigma_1V_1^T$ 和 $B=U_2\Sigma_2V_2^T$ ，其中 $U_1$ 和 $U_2$ 是 $m\times m$ 階正交矩陣， $V_1$ 和 $V_2$ 是 $n\times n$ 階正交矩陣， $\Sigma_1=\begin{bmatrix} S_1&0\\ 0&0 \end{bmatrix}$ 和 $\Sigma_2=\begin{bmatrix} S_2&0\\ 0&0 \end{bmatrix}$ 是 $m\times n$ 階奇異值矩陣 (零列和零行必須隨 $m$ 和 $n$ 的大小調整)， $S_1=\hbox{diag}(\sigma_1(A),\ldots,\sigma_s(A))$ ， $S_2=\hbox{diag}(\sigma_1(B),\ldots,\sigma_s(B))$ ， $s=\min\{m,n\}$ 。代入奇異值分解，使用跡數循環不變性，

$\displaystyle\begin{aligned} \hbox{trace}(A^TP^TBQ^T)&=\hbox{trace}(V_1\Sigma_1^TU_1^TP^TU_2\Sigma_2V_2^TQ^T)\\ &=\hbox{trace}(\Sigma_1^TU_1^TP^TU_2\Sigma_2V_2^TQ^TV_1)\\ &=\hbox{trace}(\Sigma_1^TX^T\Sigma_2Y^T), \end{aligned}$

上面令 $m\times m$ 階正交矩陣 $X^T=U_1^TP^TU_2$ ， $n\times n$ 階正交矩陣 $Y^T=V_2^TQ^TV_1$ 。兩式取轉置，解得 $P=U_2XU_1^T$ 和 $Q=V_1YV_2^T$ 。剩下的工作是找出正交矩陣 $X$ 和 $Y$ 使得 $\hbox{trace}(\Sigma_1^TX^T\Sigma_2Y^T)$ 有最大值。

考慮約束最佳化問題

$\displaystyle \max_{XX^T=I\atop{YY^T=I}}\hbox{trace}(\Sigma_1^TX^T\Sigma_2Y^T)$

使用 Lagrange 乘數法。寫出 Lagrangian 函數

$\displaystyle L=\hbox{trace}(\Sigma_1^TX^T\Sigma_2Y^T)-\hbox{trace}(\Lambda(XX^T-I))-\hbox{trace}(M(YY^T-I))$ ，

其中 $\Lambda$ 是 $m\times m$ 階且 $M$ 是 $n\times n$ 階 Lagrange 乘數矩陣。利用跡數的導數公式，計算 $L$ 對於 $X$ 和 $Y$ 的導數：

$\displaystyle\begin{aligned} \frac{\partial L}{\partial X}&=\frac{\partial\hbox{trace}(X^T\Sigma_2Y^T\Sigma_1^T)}{\partial X}-\frac{\partial\hbox{trace}(X^T\Lambda X-\Lambda)}{\partial X}-\frac{\partial\hbox{trace}(Y^TMY-M)}{\partial X}\\ &=\Sigma_2Y^T\Sigma_1^T-(\Lambda+\Lambda^T)X, \end{aligned}$

$\displaystyle\begin{aligned} \frac{\partial L}{\partial Y}&=\frac{\partial\hbox{trace}(\Sigma_1^TX^T\Sigma_2Y^T)}{\partial Y}-\frac{\partial\hbox{trace}(X^T\Lambda X-\Lambda)}{\partial Y}-\frac{\partial\hbox{trace}(Y^TMY-M)}{\partial Y}\\ &=\Sigma_1^TX^T\Sigma_2-(M+M^T)Y. \end{aligned}$

設上面兩式為零，可得最佳解的必要條件：

$\displaystyle\begin{aligned} \Sigma_2Y^T\Sigma_1^TX^T&=\Lambda+\Lambda^T\\ \Sigma_1^TX^T\Sigma_2Y^T&=M+M^T. \end{aligned}$

因為 $\Lambda+\Lambda^T$ 和 $M+M^T$ 是對稱矩陣，

$\displaystyle\begin{aligned} \Sigma_2Y^T\Sigma_1^TX^T&=X\Sigma_1Y\Sigma_2^T\\ \Sigma_1^TX^T\Sigma_2Y^T&=Y\Sigma_2^TX\Sigma_1 .\end{aligned}$

第二式右乘 $\Sigma_1^TX^T$ ，

$\displaystyle \Sigma_1^TX^T\Sigma_2Y^T\Sigma_1^TX^T=Y\Sigma_2^TX\Sigma_1\Sigma_1^TX^T$ 。

將第一式代入等號左邊，

$\displaystyle \Sigma_1^TX^T\Sigma_2Y^T\Sigma_1^TX^T=\Sigma_1^TX^TX\Sigma_1Y\Sigma_2^T=\Sigma_1^T\Sigma_1Y\Sigma_2^T$ 。

合併上面兩式，右乘 $X$ ，就有

$\displaystyle (Y\Sigma_2^TX)(\Sigma_1\Sigma_1^T)=(\Sigma_1^T\Sigma_1)(Y\Sigma_2^TX)$ ，

其中 $\Sigma_1\Sigma_1^T$ 是 $m\times m$ 階對角矩陣， $\Sigma_1^T\Sigma_1$ 是 $n\times n$ 階對角矩陣，兩者有相同的主對角元 (扣除零元) $\sigma_i^2(A)$ ， $1\le i\le s$ 。比較等號兩邊的 $(i,j)$ 元，可得

$\displaystyle (Y\Sigma_2^TX)_{ij}\sigma_j^2(A)=\sigma_i^2(A)(Y\Sigma_2^TX)_{ij}$ 。

對於 $i\neq j$ ，若 $\sigma_i(A)\neq\sigma_j(A)$ ，則 $(Y\Sigma_2^TX)_{ij}=0$ 。這個結果提示 $X=I_m$ 和 $Y=I_n$ 為一組最佳解 (事實上，所有的非零純量矩陣 $X=cI_m$ ， $Y=(1/c)I_n$ 皆為最佳解)，也就有 $P=U_2U_1^T$ ， $Q=V_1V_2^T$ ，目標函數可化簡成

$\displaystyle \hbox{trace}(A^TP^TBQ^T)=\hbox{trace}(\Sigma_1^T\Sigma_2)=\sum_{i=1}^s\sigma_i(A)\sigma_i(B)$ 。

根據排序不等式 (rearrangement inequality)^[4]，上式的最大值發生於 $\sigma_1(A)\ge\cdots\ge\sigma_s(A)$ 且 $\sigma_1(B)\ge\cdots\ge\sigma_s(B)$ 。因此，矩陣 $B$ 的最小平方近似矩陣為

$\displaystyle PAQ=U_2U_1^TU_1\Sigma_1V_1^TV_1V_2^T=U_2\Sigma_1V_2^T$ ，

最小誤差則為

$\displaystyle\begin{aligned} \Vert B-PAQ\Vert_F&=\Vert B-U_2\Sigma_1V_2^T\Vert_F=\sqrt{\Vert A\Vert_F^2+\Vert B\Vert_F^2-2\sum_{i=1}^s\sigma_i(A)\sigma_i(B)}\\ &=\sqrt{\sum_{i=1}^s\sigma_i^2(A)+\sum_{i=1}^s\sigma_i^2(B)-2\sum_{i=1}^s\sigma_i(A)\sigma_i(B)}\\ &=\sqrt{\sum_{i=1}^s\left(\sigma_i(A)-\sigma_i(B)\right)^2}, \end{aligned}$

上面我們使用了 Frobenius 範數的奇異值恆等式 (見“SVD 於矩陣近似的應用”)

$\displaystyle \Vert A\Vert_F=\sqrt{\sum_{i=1}^s\sigma_i^2(A)}$ 。

在兩面正交 Procrustes 問題中，從最小誤差表達式可知 $B=PAQ$ 的充分必要條件是 $A$ 和 $B$ 有相同的奇異值。

直觀上，兩面正交 Procrustes 問題 $\min_{PP^T=I,QQ^T=I}\Vert B-PAQ\Vert_F$ 的極小值應比正交 Procrustes 問題 $\min_{QQ^T=I}\Vert B-AQ\Vert_F$ 的極小值來得小，這意味

$\displaystyle \sum_{i=1}^n\sigma_i(A^TB)\le\sum_{i=1}^{\min\{m,n\}}\sigma_i(A)\sigma_i(B)$ 。

這個不等式的證明過程相當冗長，在此從略。

參考來源：
[1] 維基百科：Procrustes
[2] 正交 Procrustes 問題也可以表示為求一個 $m\times m$ 階正交矩陣 $Q$ 使最小化 $\Vert B-QA\Vert_F$ 。因為 $\Vert A\Vert_F=\hbox{trace}(A^TA)=\hbox{trace}(AA^T)=\Vert A^T\Vert_F$ ，推得 $\Vert B-QA\Vert_F=\Vert B^T-A^TQ^T\Vert$ ，其中 $Q^T$ 是一個正交矩陣。
[3] 取自 G. H. Golub and C. F. Van Loan, Matrix Computations, 2nd edition, 1989, pp 582.
[4] 維基百科：排序不等式

8 Responses to 正交 Procrustes 問題

張盛東 says:

11/21/2013 at 1:21 pm

老師，這種矩陣近似一般在什麼場合下使用呢？

- ccjou says:
  
  11/21/2013 at 2:10 pm
  
  我知道的主要應用有二個。
  (1) 幾何形態測量 Geometric Morphometric，下文圖一包含縮放、平移和旋轉。旋轉是正交變換。
  http://www.intechopen.com/books/sexual-dimorphism/sexual-dimorphism-using-geometric-morphometric-approach
  (2) 主成分因素分析 principal factor analysis，見下文
  http://support.sas.com/documentation/cdl/en/statug/63347/HTML/default/viewer.htm#statug_factor_sect029.htm
  
  第二張圖顯示PCA得到的五個變數的factor loading (見
  https://ccjou.wordpress.com/2013/04/18/%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90%E8%88%87%E5%A5%87%E7%95%B0%E5%80%BC%E5%88%86%E8%A7%A3/
  )
  從圖看不出這五個變數的分類，所以我們可以旋轉座標，常用的方法有二：Varimax Prerotation，結果見第三張圖，以及Promax Rotation (其中包含 Procrustes 旋轉)，結果見最後一張圖。從旋轉後的座標可知五個變數可以分成三類。
  
  P.S. 貼一張動畫上來。
  
徐君 says:

01/09/2017 at 6:21 pm

老师，你在文中说：比較等號兩邊可得 QV=U，故 Q=UV^T 為一最佳解。那这个解是不是唯一的吗？或者说正交 Procrustes 問題是不是一个凸问题？

徐君 says:

01/09/2017 at 6:51 pm

老师，你在文中说: 設上式為零，可得 A^TB=(\Lambda+\Lambda^T)Q，或 A^TBQ^T=\Lambda+\Lambda^T。这里用到了QQ^T=I这个条件。可是，在Lagrangian 函數里还能用到之前有的约束条件（比如QQ^T=I）吗？

徐君 says:

01/09/2017 at 7:18 pm

老师，我已经弄明白了，读者可以参考这篇论文
http://nemo.nic.uoregon.edu/wiki/images/0/07/Psychometrika_1966_Sch%C3%B6nemann_A_generalized_solution_of_the.pdf

徐君 says:

01/09/2017 at 7:36 pm

正交 Procrustes 問題里的Q可以有更直观的解决方法，请看https://jianfengwang.files.wordpress.com/2015/07/proof_of_the_orthogonal_procrustes_problem.pdf

- ccjou says:
  
  01/09/2017 at 8:11 pm
  
  謝謝，確實是有更簡潔的推導法，有空我再補上來。
  
  (見“補充一個較快捷的推導方式”該段)
  
徐君 says:

01/09/2017 at 9:40 pm

谢谢老师！老师这篇帖子里提到的技巧帮助我证明了我论文中的一个定理。

	TIGER on 特殊矩陣 (6)：正定矩陣
	tsao on 線性變換與矩陣的用語比較
	123 on 奇異值分解 (SVD)
	小鎮馬 on Lagrange 乘數法
	johnson on 特殊矩陣 (3)：么正矩陣 (酉矩陣)
	AnneLULU on Gapminder 讓資料說故事